ITD模态参数识别matlab修改版

%ITDxx识别模态参数
clear
clc
close all hidden
format long
%% txt文件下输入
fni=input('ITD法模态参数识别-输入数据文件名：','s');
fid=fopen(fni,'r');
mn=fscanf(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=fscanf(fid,'%s',1);%输出数据文件名
b=fscanf(fid,'%f',[ig,inf]);%实测时域或频域数据
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 cancel the selection');%如果取消选择则显示提示
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=linspace(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
F1=abs(log(U'))./(2*pi*dt);%计算模态频率向量
D1=sqrt(1./(((imag(log(U'))./real(log(U'))).^2)+1));%计算阻尼比向量%计算振型系数向量
l=1;
for k=0:(2*n-1)
Va(k+1,:)=[conj(U).^k];
end
S1=(inv(conj(Va')*Va)*conj(Va')*h);%inv矩阵求逆
h1=real(Va*S1);%计算生成的脉冲响应函数%绘制脉冲响应函数拟合曲线图
figure(1);
plot(t,h,':',t,h1);
xlabel('时间（s）');
ylabel('幅值');
legend('实测','拟合');
grid on;
if ig>1
H1=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('实部');
legend('实测','拟合');
grid on;
%绘制频响函数虚部拟合曲线图
subplot(2,1,2);
plot(f,b(2,:),':',f,imag(H1(nn)),'r');
xlabel('频率(Hz)');
ylabel('虚部');
legend('实测','拟合');
grid on;
end
[F2,I]=sort(F1);%将自振频率从小到大排列
%剔除方程中的非模态项（非共轭根）和共轭项（重复项）m=0;
for k=1:1:nm-1
if F2(k)~=F2(k+1)
continue;
end
m=m+1;
l=I(k);
F(m)=F1(l);%自振频率
D(m)=D1(l);%阻尼比
S(m)=S1(l);%振型系数
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);。