3DOF.rar_newmark

clear all; load QQ M=[4 0 0; 0 2 0; 0 0 4]; K=100*[7 -1 -1; -1 3 -1; -1 -1 7]; N=3 a=0.02; b=0.004; e=1 [c d]=eig(K,M) C=(a*M+b*K)*e I=eye(N); D=(a*I+b*d); wn=sqrt(d); drr=D/wn/2 wn=diag(sqrt(d))./(2*pi) damping=diag(drr) AA1=[1;2;1]./4 AA2=[-1;2;-1]./4 AA3=[-1;0;1]./sqrt(8) AA=[AA1 AA2 AA3] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% WW=0.2; for j=1:N for k=1:N XN(k,:,j)=X(k,:,j)+std(X(k,:,j))*WW*randn(1,length(X(k,:,j))); end end r=1; f_max=5; df=1/(RecordLength*dt); f_r=round(f_max/df); for j=1:N for k=1:N %X_(1,:,k,j)=fft(X(k,:,j));%未加雜訊 X_(1,:,k,j)=fft(XN(k,:,j));%加雜訊 end end for p=1:N F_(p,:)=fft(F(p,:)); end for i=1:f_r H11(1,i)=X_(1,i,1,1)/F_(1,i); H21(1,i)=X_(1,i,2,1)/F_(1,i); H31(1,i)=X_(1,i,3,1)/F_(1,i); H12(1,i)=X_(1,i,1,2)/F_(2,i); H22(1,i)=X_(1,i,2,2)/F_(2,i); H32(1,i)=X_(1,i,3,2)/F_(2,i); H13(1,i)=X_(1,i,1,3)/F_(3,i); H23(1,i)=X_(1,i,2,3)/F_(3,i); H33(1,i)=X_(1,i,3,3)/F_(3,i); end f=0 for w=0:df*2*pi:df*(f_r-1)*2*pi; f=f+1; H(:,:,f)=[H11(1,f) H12(1,f) ; H21(1,f) H22(1,f) ; H31(1,f) H32(1,f) ]; [u,s,v]=svd(H(:,:,f)); U1(:,f)=u(:,1); V1(:,f)=v(:,1); U2(:,f)=u(:,2); V2(:,f)=v(:,2); s1(1,f)=s(1,1); s2(1,f)=s(2,2); end; s=0:df:(f_r-1)*df; figure(1); semilogy(s,abs(H11)); xlabel('Frequency(Hz)'); ylabel('Amplitude(log)'); title('H11 (20% noise)'); figure(2); semilogy(s,s1,'k'); hold on; semilogy(s,s2,'b'); xlabel('Frequency(Hz)'); ylabel('Amplitude(log)'); hold off; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure(3); semilogy(s,s1,'r'); xlabel('Frequency(rad/s)'); ylabel('Amplitude(log)'); title('SVD Plot'); [A1,B1]=ginput(2) figure(4); semilogy(s,s2,'r'); xlabel('Frequency(rad/s)'); ylabel('Amplitude(log)'); title('SVD Plot'); [A2,B2]=ginput(1) %%%%%%%%%%%%%%%%%%%%%%%%%%% x1=round(A1(1,1)/df); x2=round(A1(2,1)/df); x3=round(A2(1,1)/df); for j=-r:1:r peak1(1,r+1+j)=s1(1,x1+j); peak2(1,r+1+j)=s1(1,x2+j); peak3(1,r+1+j)=s2(1,x3+j); end; [C1,D1]=max(peak1); [C2,D2]=max(peak2); [C3,D3]=max(peak3); pp(1,1)=D1+x1-r-1; pp(1,2)=D2+x2-r-1; pp(1,3)=D3+x3-r-1; modeshape1=U1(:,pp(1,1)); modeshape2=U1(:,pp(1,2)); modeshape3=U2(:,pp(1,3)); modalpart1=V1(:,pp(1,1)); modalpart2=V1(:,pp(1,2)); modalpart3=V2(:,pp(1,3)); for k=1:f_r EFRF1(1,k)=modeshape1'*H(:,:,k)*modalpart1; EFRF2(1,k)=modeshape2'*H(:,:,k)*modalpart2; EFRF3(1,k)=modeshape3'*H(:,:,k)*modalpart3; end; for l=-r:1:r peak_1(1,r+1+l)=EFRF1(1,x1+l); peak_2(1,r+1+l)=EFRF2(1,x2+l); peak_3(1,r+1+l)=EFRF3(1,x3+l); end; [E1,F1]=max(peak_1); [E2,F2]=max(peak_2); [E3,F3]=max(peak_3); figure(6); semilogy(s,abs(EFRF1)); title('Enhance FRF For Mode1'); xlabel('Frequency(Hz)'); ylabel('Amplitude(log)'); figure(7); semilogy(s,abs(EFRF2)); title('Enhance FRF For Mode2'); xlabel('Frequency(Hz)'); ylabel('Amplitude(log)') figure(8); semilogy(s,abs(EFRF3)); title('Enhance FRF For Mode3'); xlabel('Frequency(Hz)'); ylabel('Amplitude(log)') MAC11=(abs(modeshape1'*AA1)^2)/(abs(modeshape1'*modeshape1)*abs(AA1'*AA1)) MAC22=(abs(modeshape2'*AA2)^2)/(abs(modeshape2'*modeshape2)*abs(AA2'*AA2)) MAC33=(abs(modeshape3'*AA3)^2)/(abs(modeshape3'*modeshape3)*abs(AA3'*AA3)) save QQ1 s1 s2 pp df wn damping EFRF1 EFRF2 EFRF3 U1 U2 modeshape1 modeshape2 modeshape3