ERA.zip_ERA_ERA算法_模态识别_结构模态_自然激励技术

format long time=xlsread('C:\Users\何定桥\Desktop\课程文件\高等结构动力学\作业4\Ydata.csv',1,'a12000:a20000'); acc1_raw=xlsread('C:\Users\何定桥\Desktop\课程文件\高等结构动力学\作业4\Ydata.csv',1,'b12000:b20000'); acc2_raw=xlsread('C:\Users\何定桥\Desktop\课程文件\高等结构动力学\作业4\Ydata.csv',1,'c12000:c20000'); acc3_raw=xlsread('C:\Users\何定桥\Desktop\课程文件\高等结构动力学\作业4\Ydata.csv',1,'d12000:d20000'); acc4_raw=xlsread('C:\Users\何定桥\Desktop\课程文件\高等结构动力学\作业4\Ydata.csv',1,'e12000:e20000'); m1=mean(acc1_raw);m2=mean(acc2_raw);m3=mean(acc3_raw);m4=mean(acc4_raw); acc1_raw=acc1_raw-m1; acc2_raw=acc2_raw-m2; acc3_raw=acc3_raw-m3; acc4_raw=acc4_raw-m4; h=hanning(length(acc1_raw)); acc1_raw=acc1_raw.*h; acc2_raw=acc2_raw.*h; acc3_raw=acc3_raw.*h; acc4_raw=acc4_raw.*h; deltat=1/256; %对数据重叠平均 N=1000;%最终取样点数 movement=200;%每次向后移动的点数 times=fix((length(acc1_raw)-N)/movement); acc1=zeros(N,1); acc2=zeros(N,1); acc3=zeros(N,1); acc4=zeros(N,1); for i=1:times acc1=acc1_raw(movement*(i-1)+1:N+movement*(i-1))+acc1; acc2=acc2_raw(movement*(i-1)+1:N+movement*(i-1))+acc2; acc3=acc3_raw(movement*(i-1)+1:N+movement*(i-1))+acc3; acc4=acc4_raw(movement*(i-1)+1:N+movement*(i-1))+acc4; end acc1=acc1/times; acc2=acc2/times; acc3=acc3/times; acc4=acc4/times; n=0:N-1; t=n*deltat; [a11,b11]=xcorr(acc1,acc1); [a12,b12]=xcorr(acc1,acc2); [a13,b13]=xcorr(acc1,acc2); [a14,b14]=xcorr(acc1,acc4); [a22,b22]=xcorr(acc2,acc2); [a23,b23]=xcorr(acc2,acc3); [a24,b24]=xcorr(acc2,acc4); [a33,b33]=xcorr(acc3,acc3); [a34,b34]=xcorr(acc3,acc4); [a44,b44]=xcorr(acc4,acc4); alpha=600;beta=600; Ymatrix=zeros(N,0,0); for i=1:N Ymatrix(i,1)=a11(i); Ymatrix(i,2)=a12(i); Ymatrix(i,3)=a13(i); Ymatrix(i,4)=a14(i); end H0=Ymatrix(1,1,:); for i=2:beta H0=cat(2,H0,Ymatrix(i,1,:)); end for j=2:alpha temp=Ymatrix(j,1,:); for i=2:beta temp=cat(2,temp,Ymatrix(i,1,:)); end H0=cat(3,H0,temp); end H1=Ymatrix(2,1,:); for i=3:beta+1 H1=cat(2,H1,Ymatrix(i,1,:)); end for j=2:alpha+1 temp=Ymatrix(j,1,:); for i=3:beta+1 temp=cat(2,temp,Ymatrix(i,1,:)); end H1=cat(3,H1,temp); end H02d=zeros(alpha*4,beta*4); H12d=zeros(alpha*4,beta*4); for i=1:alpha*4 for j=1:beta*4 H02d(i,j)=H0(1,i,j); H12d(i,j)=H1(1,i,j); end end [R,Sigma,ST] = svd(H02d); R=R(1:4,1:4); Sigma=Sigma(1:4,1:4); ST=Sigma(1:4,1:4); Root_Sigma=sqrt(Sigma); P_alpha=R*Root_Sigma; Q_beta=Root_Sigma*ST; A=inv(P_alpha)*H12d(1:4,1:4)*inv(Q_beta); z=eig(A); s=log(z)/deltat; omega=sqrt(real(s).^2+imag(s).^2);