ITD模态参数识别matlab修改版

%ITDxx识别模态参数
clear
clc
close all hidden
format long
%% txt 文件下输入
fni=input('ITD 法模态参数识别-输入数据文件名：','s'); fid=fopen(fni,'r'); mn二fsca nf(fid,'%d',1);%模态阶数
%定义输入实测数据类型
%ig=1 时域数据如冲击响应、自由振动、互相关函数、随机减量法处理结果
%ig=2频域数据如频响函数实部和虚部数据
ig=fscanf(fid,'%f',1);
%ig=1时，f为采样频率sf, ig=2时，f为频率间隔df
f=fscanf(fid,'%f',1);
fno=fsca nf(fid,'%s',1);%输出数据文件名
b=fsca nf(fid,'%f',[ig,i nf]);% 实测时域或频域数据
status=fclose(fid);
%%
clc;
clear all;
format long
[FileName,PathName] = uigetfile('*.mat', 'Select the Mat-files of time
signal'); % 窗口读文件,并获取包含路径的文件名
if isequal(FileName,0)
disp('User can cel the selectio n');%如果取消选择则显示提示
return;
else
FULLFILE=fullfile(PathName,FileName);
Signal_str= sprintf('User selected signal file:%s',FULLFILE);
disp(Signal_str);
Struct=load(FULLFILE);
end
c=fieldnames(Struct);%得到一个元胞数组，包含Struct中各个域名(倘若有多个的话)
b=getfield(Struct,c{1}); %获取c{1}对应的域中的内容
b=b(3601:9600);
%%
%ig=1 时域数据如冲击响应、自由振动、互相关函数、随机减量法处理结果
%ig=2 频域数据如频响函数实部和虚部数据
ig=input('数据类型ig=');
f=input('采样频率f=');%指定采样频率
mn=input('计算模态阶数mn二');％指定计算模态阶数
%建立特征方程矩阵的阶数(为模态阶数的 2 倍)
nm=2*mn;
%组织识别计算多用的时域数据及参数
if ig==1%实测时域数据
sf=f;%采样频率
n=fix(length(b)/2);%向0靠拢取整，取时域数据的长度
h=b(1,1:2* n)';%将输入时域数据赋值给列向量h
dt=1/sf;%时间间隔
t=0:dt:(2*n-1)*dt;% 建立离散时间向量
else %实测频域数据
df=f;%取频率间隔
n=length(b(1,:));
f=0:df:(n-1)*df;% 建立离散频率向量
H=b(1,:)'+b(2,:)'*i;%建立对应正负频率的实测频响函数向量
H(n+1)=real(H(n));
H( n+2:2* n)二conj(H( n:-1:2));%conj 求负数的共轭值
h=real(ifft(H));%频响函数经IFFT并取实部变换成脉冲响应函数t=li nspace(0,1/df,2* n);%建立离散时间向量
dt=t(2)-t(1);%计算时间间隔
end
%计算自由振动响应矩阵
L=length(h);
M=L/2;
for k=1:nm
x1(k,:)=h(k:L-(nm-k+1))';
x2(k,:)=h(k+1:L-(nm-k))';
end
%用最小二乘法求解特征方程矩阵B=x1\x2;%B=x2*x1'*inv(x1*x1');
[A,V]=eig(B);%计算特征值及特征向量(特征值V,特征向量A)
%变换特征值对角阵为一向量
for k=1:nm
U(k)=V(k,k);
end
F仁abs(log(U'))./(2*pi*dt);%计算模态频率向量
D仁sqrt(1./(((imag(log(U'))./real(log(U')))八2)+1));% 计算阻尼比向量%计算振型系数向量
l=1;
for k=0:(2*n-1)
Va(k+1,:)=[conj(U).A k];
end
S仁(i nv(conj(Va')*Va)*conj(Va')*h);% inv 矩阵求逆
hl二real(Va*S1);%计算生成的脉冲响应函数%绘制脉冲响应函数拟合曲线图
figure(1);
plot(t,h,':',t,h1);
xlabel('时间(s)');
ylabel('幅值');
lege nd('实测','拟合');
grid on;
if ig>1
H仁fft(Va*S1);%计算生成的频响函数%绘制频响函数实部拟合曲线图figure(2); nn=1:n;
subplot(2,1,1);
plot(f,real(H(nn)),':',f,real(H1(nn)),'r'); xlabel('频率(Hz)');
ylabel('实部');
lege nd('实测','拟合');
grid on;
%绘制频响函数虚部拟合曲线图
subplot(2,1,2);
plot(f,b(2,:),':',f,imag(H1(nn)),'r');
xlabel('频率(Hz)');
ylabel('虚部');
lege nd('实测','拟合');
grid on;
end
[F2,l]=sort(F1);%将自振频率从小到大排列
%剔除方程中的非模态项(非共轭根)和共轭项(重复项) m=0;
for k=1:1:nm-1
if F2(k)~=F2(k+1) continue;
end
m=m+1;
l=l(k);
F(m)=F1(l);%自振频率
D(m)=D1(l);%阻尼比
S(m)=S1(l);%g 型系数
end
%打开文件输出识别的模态参数数据
fno='out.txt';
fid=fopen(fno,'w');
fprintf(fid,'频率(Hz)阻尼比(％%)振型系数\n');
for k=1:m
fprintf(fid,'%10.4f %10.4f%10.6f\n',F(k),D(k)*100.0,imag(S(k))); end
status=fclose(fid);。