ITD模态参数识别matlab修改版

%ITD法识别模态参数clearclccloseallhiddenformatlong%%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为频率间隔dff=fscanf(fid,'%f',1);fno=fscanf(fid,'%s',1);%输出数据文件名b=fscanf(fid,'%f',[ig,inf]);%实测时域或频域数据status=fclose(fid);%%clc;clearall;formatlong[FileName,PathName]=uigetfile('*.mat','SelecttheMat-filesoftimesignal');%窗口读文件,并获取包含路径的文件名ifisequal(FileName,0)disp('Usercanceltheselection');%如果取消选择则显示提示return;elseFULLFILE=fullfile(PathName,FileName);Signal_str=sprintf('Userselectedsignalfile:%s',FULLFILE);disp(Signal_str);Struct=load(FULLFILE);endc=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;%组织识别计算多用的时域数据及参数ifig==1%实测时域数据sf=f;%采样频率n=fix(length(b)/2);%向0靠拢取整，取时域数据1/2的长度h=b(1,1:2*n)';%将输入时域数据赋值给列向量hdt=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;fork=1:nmx1(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）%变换特征值对角阵为一向量fork=1:nmU(k)=V(k,k);endF1=abs(log(U'))./(2*pi*dt);%计算模态频率向量D1=sqrt(1./(((imag(log(U'))./real(log(U'))).^2)+1));%计算阻尼比向量%计算振型系数向量l=1;fork=0:(2*n-1)Va(k+1,:)=[conj(U).^k];endS1=(inv(conj(Va')*Va)*conj(Va')*h);%inv矩阵求逆h1=real(Va*S1);%计算生成的脉冲响应函数%绘制脉冲响应函数拟合曲线图figure(1);plot(t,h,':',t,h1);xlabel('时间（s）');ylabel('幅值');legend('实测','拟合');gridon;ifig1H1=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('实测','拟合');gridon;%绘制频响函数虚部拟合曲线图subplot(2,1,2);plot(f,b(2,:),':',f,imag(H1(nn)),'r');xlabel('频率(Hz)');ylabel('虚部');legend('实测','拟合');gridon;end[F2,I]=sort(F1);%将自振频率从小到大排列%剔除方程中的非模态项（非共轭根）和共轭项（重复项）m=0;fork=1:1:nm-1ifF2(k)~=F2(k+1)continue;endm=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');fork=1:mfprintf(fid,'%10.4f%10.4f%10.6f\n',F(k),D(k)*100.0,imag(S(k)));endstatus=fclose(fid);