时域模态参数识别方法，有效的识别时域信号 并提取模态参数

%STD法模态参数识别 clear clc close all hidden format long %fni=input('STD法模态参数识别－输入数据文件名','s'); %fid=fopen(fni,'r') %mn=fscanf(fid,'%f',1); %模态阶数 %ig=fscanf(fid,'%d',1); %ig=1 时域数据 ig=2 频域数据 %f=fscanf(fid,'%f',1); %ig=1时f为采样频率sf,ig=2时f为频率间隔df %fno=fscanf(fid,'%d',1); %读入输出数据文件名 %b=fscanf(fid,'%f',[ig inf]); %读入数据 %status=fclose(fid); mn=3;ig=2;f=0.1953125;fno='out9_4.mat'; load y b=[y(1:100);y(101:200)]; nm=2*mn; if ig==1 sf=f; n=fix(length(b)/2); h=b(1,1:2*n)'; 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)); h=real(ifft(H)); t=linspace(0,1/df,2*n); dt=t(2)-t(1); end M=length(h)-nm; for k=1:M x1(k,:)=h(k:k+nm-1)'; x2(k,:)=h(k+1:k+nm)'; end B=zeros(nm,nm); B(2:nm,1:nm-1)=eye(nm-1,nm-1); B(:,nm)=x1\x2(:,nm);%B(:,nm)=inv(x1'*x1)*x1'*x2(:,nm); [A,V]=eig(B); 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=1:nm if abs(real(U(k)))<=1&abs(imag(U(k)))<=1 V0(l)=U(k); l=l+1; end end for k=0:(2*n-1) Va(k+1,:)=[conj(U).^k]; end S1=(inv(conj(Va')*Va)*conj(Va')*h); 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: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);