3971002ITD.rar_模态 参数识别_模态参数识别的方法_模态识别

%ITD法识别模态参数--王济 clear clc close all hidden format long 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); %建立特征方程矩阵的阶数（为模态阶数的2倍） nm=2*mn; %组织识别计算多用的时域数据及参数 if ig==1 %实测时域数据 sf=f; %采样频率 n=fix(length(b)/2); %向0靠拢取整，取时域数据1/2的长度 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 %打开文件输出识别的模态参数数据 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);