daonayuan.rar_固有频率_导纳_模态图

clear clc close all hidden format long global mn; fni=input('导纳圆法模态参数识别-输入文件名:','s'); fid=fopen(fni,'r'); mn=fscanf(fid,'%d',1); df=fscanf(fid,'%f',1); f0=fscanf(fid,'%f',[1,mn]); fno=fscanf(fid,'%s',1); H0=fscanf(fid,'%f',[2,inf]); status=fclose(fid); f=0:df:(length (H0(1,:))-1)*df; w=2*pi*f; fds=zeros(1,3*mn); for n=1:4 for j=1:mn l=3*(j-1); fds(l+1:l+3)=zeros(1,3); H1=fun81(fds,w); H=H0-H1; k=round(f0(j)/df)+1; [m,kc]=max(abs(H(2,k-4:k+4))); kc=kc+k-5; v=0.15*abs(H(2,kc)); for k=2:30 l=kc-k; if abs(H(2,l))>v&abs(H(2,l))<=abs(H(2,l+1)) k1=l; else break; end end for k=2:30 l=kc+k; if abs(H(2,l))>v&abs(H(2,l))<=abs(H(2,l-1)) k2=l; else break; end end x=H(1,k1:k2); y=H(2,k1:k2); A=[sum(x.^2),sum(x.*y),sum(x); sum(x.*y),sum(y.^2),sum(y); sum(x),sum(y),length(x)]; B=[-sum(x.^3+x.*y.^2); -sum(x.^2.*y+y.^3); -sum(x.^2+y.^2)]; c=inv(A)*B; x0=-c(1)/2; y0=-c(2)/2; r=sqrt(x0^2+y0^2-c(3)); x1=f(kc-1); x2=f(kc); x3=f(kc+1); y1=H(2,kc-1); y2=H(2,kc); y3=H(2,kc+1); fc=((y2-y1)*(x3^2-x2^2)-(y3-y2)*(x2^2-x1^2))/(2*(2*y2-y1-y3)*(x2-x1)); %拟合实部 yc=(y1*(x2-fc)^2-y2*(x1-fc)^2)/((x2-fc)^2-(x1-fc)^2); %拟合虚部 if abs(y1)>abs(y3) xc=H(1,kc-1)+(H(1,kc)-H(1,kc-1))*(fc-x1)/(x2-x1); else xc=H(1,kc)+(H(1,kc+1)-H(1,kc))*(fc-x2)/(x3-x2); end F(j)=fc; [m,k]=max (abs(y)); a=(xc-x(k-1))^2+(yc-y(k-1))^2; b=(x0-xc)^2+(y0-yc)^2; c=(x0-x(k-1))^2+(y0-y(k-1))^2; a1=acos((b+c-a)/(2*sqrt(b)*sqrt(c)))/2; a=(x(k+1)-xc)^2+(y(k+1)-yc)^2; b=(x0-x(k+1))^2+(y0-y(k+1))^2; c=(x0-xc)^2+(y0-yc)^2; a2=acos((b+c-a)/(2*sqrt(b)*sqrt(c)))/2; D(j)=(f(l+1)-f(l-1))/(F(j)*(tan(a1)+tan(a2))); S(j)=yc*2*D(j); l=3*(j-1); fds(l+1:l+3)=[2*pi*F(j) D(j) S(j)]; end end H1=fun81(fds,w); %绘制频响函数实部拟合曲线 figure(1) subplot(2,1,1); plot(f,H0(1,:),':',f,H1(1,:)); xlabel('频率（HZ）'); ylabel('实部'); legend('实测','拟合'); grid on; %绘制频响函数虚部拟合曲线 subplot(2,1,2); plot(f,H0(2,:),':',f,H1(2,:)); xlabel('频率（HZ）'); ylabel('虚部'); legend('实测','拟合'); grid on; %绘制频响函数导纳拟合曲线 figure(2) plot(H0(1,:),H0(2,:),':',H1(1,:),H1(2,:)); axis('equal'); xlabel('实部'); ylabel('虚部'); legend('实测','拟合'); grid on; %打开文件输出模态参数数据 fid=fopen(fno,'w'); fprintf(fid,' 频率（Hz） 阻尼比（%%） 振型系数\n'); for k=1:mn fprintf(fid,'%10.4f %10.4f %10.6f\n',F(k),D(k)*100.0,S(k)); end status=fclose(fid);