MATLAB在振动信号处理中的应用 配套源码

%导纳圆法模态参数识别 clear clc close all hidden format long global mn %fni=input('导纳圆法模态参数识别－输入数据文件名','s'); %fid=fopen(fni,'r') %mn=fscanf(fid,'%f',1); %模态阶数 %df=fscanf(fid,'%d',1); %频率间隔 %f0=fscanf(fid,'%d',[1 mn]); %模态频率初值数组 %fno=fscanf(fid,'%d',1); %读入输出数据文件名 %H0=fscanf(fid,'%f',[2 inf]); %实测频响函数实部虚部数据 %status=fclose(fid); mn=3;df=0.02;f0=[5 10 20];fno='out8_1.mat'; %tt=0:0.01:20; %y=(1-0.2*tt).*sin(2*pi*5*tt); %yy=fft(y);y=yy(1:1000); %H0=[real(y);imag(y)]; load y H0=[y(1:1000);y(1001:2000)]; 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=1: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=1: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.^2+x.*y.^2);-sum(x.^2.*y+y.^3);-sum(x.^2+y.^2)]; c=inv(A)*B; x0=-c(1)/2; y0=-c(1)/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,S(k)); end status=fclose(fid);