zhengjiao.rar_matlab 频响函数_固有频率_固有频率识别_模态阻尼_频响

clear clc close all hidden format long %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% fni=input('input file name:','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*mn; %定义正交多项式的阶次 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,1:n)+i*h(2,1:n)]; %建立对应（从负到正）的实测频响数据 [k,l]=size(w); %根据行、列数调整，这里把w改为n*1的列向量 if k<l w=w.'; end [k,l]=size(H); %根据行、列数调整，这里把H改为n*1的列向量 if k<l H=H.'; end [p,cp]=fun85(H,w,1,nm); %建立分子正交多项式【和】幂多项式转换矩阵 [q,cq]=fun85(H,w,2,nm); %建立分母正交多项式【和】幂多项式转换矩阵 P=p; %%%%%%%%%%%%%%%%最小二乘法求解分子系数C和分母向量D%%%%%%%%%%%%%%%%%%%%%% for k=1:nm 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; %%%%%%%%%%%%%%%%最小二乘法求解分子系数C和分母向量D%%%%结束结束结束%%%%%%%% A=cp*C; %通过转换矩阵求分子幂多项式的系数向量A A=A(nm+1:-1:1).'; %将A降幂排，且为列向量 B=cq*D; %通过转换矩阵求分母幂多项式的系数向量B B=B(nm+1:-1:1).'; %将B降幂排，且为列向量 [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 H=h(1,1:n)+i*h(2,1:n); %建立正频率的的频响函数向量 for l=1:n %计算生成的频响函数 H1(l)=sum(C'.*p(n+1,:))/sum(D'.*q(n+1,:)); 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'); %输出结果到文件output fprint(fid,' 频率（Hz） 阻尼比(%%) 振型系数\n'); for k=1:m fprintf(fid,'%10.4f %10.4f %10.6\n',F(k),D(k)*100,imag(S(k))); end status=fclose(fid);