OP.rar_振动模态分析_振型识别_模态参数_阻尼_阻尼比

%正交多项式法模态参数识别 %========================================================== clear all clc close all hidden format long %=========================================================== fni=input('正交多项式法模态参数识别-输入数据文件名:','s'); fid=fopen(fni,'r'); mn=fscanf(fid,'%d',1); %模态阶数 df=fscanf(fid,'%f',1); %频率间隔 fno=fscanf(fid,'%s',1); %输出数据文件名 h=fscanf(fid,'%f',[2,inf]); %实测频响函数实部虚部数据 status=fclose(fid); %========================================================== %定义正交多项式的阶次 nm=2*nm; %取频响函数的长度 n=length(h(1,:)); %建立离散频率向量 f=0:df:(n-1)*df; %建立从负到正的归一化频率向量 w=[-f(n:-1:1)f(1:n)]/max(f); %建立对应频率向量的实测频响函数向量 H=[h(1,n:-1:1)-i*h(2,n:-1:1)h(1,n:-1:1)+i*h(2,1:n)]; %将频率向量设置为列向量 [k,l]=size(w); if k<l w=w.'; end %将频响函数向量设置为列向量 [k,l]=size(H); if k<l H=H.'; end %建立分子正交多项式及幂多项式的转换矩阵 [p,cp]=fun85(H,w,2,nm); %最小二乘法求解正交多项式的分子系数C和分母向量D P=p; for k=1:m Q(:,k)=(H.*q(:,k)); end T=-real(P'*Q); g=real(P'*(H.*q(:,nm+1))); D=-inv(eye(nm)-T'*T)*T'*g; C=g-T*D; D(nm+1)=1; %通过转换矩阵计算分子幂多项式系数向量A A=cp*C; %将A设置成按降幂排序的列向量 A=A(nm+1:-1:1).'; %通过转换矩阵计算分母幂多项式系数向量A B=cq*D; %将B转换成按幂排序的列向量 B=B(nm+1:-1:1).'; %计算有理分式多项式的极点和留数 [R,U,K]=residue(A,B); %计算模态频率向量 F1=abs(U)*max(f); %计算阻尼比向量 D1=-real(U)./(abs(U)); %计算振型系数向量 for k=1:nm if k==1 u(1:nm-1)=U(2:nm); else u(1:k-1)=U(1:k-1); u(k:nm-1)=U(k+1:nm); end y=poly(u); S1(k)=polyval(A,U(k))/polyval(y,U(k)); end %================================================== %绘制频响函数实部拟合曲线图 subplot(2,1,1); plot(f,real(H),':',f,real(H1),'r'); xlabel('频率(Hz)'); ylabel('实部'); legend('实测'，'拟合'); grid on; %绘制频响函数虚部拟合曲线图 subplot(2,1,2); plot(f,imag(H),':',f,imag(H1),'r'); xlabel('频率(Hz)'); ylabel('虚部'); legend('实测'，'拟合'); grid on; %=================================================== %将模态频率从小到大排序 [F2,I]=sort(F1); %剔除方程解中的非模态项（非共轭根）和共轭项（重复项） m=0; for k=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);