fdd.rar_fdd_测频率_频率_频率响应法_频率测

%%%%%%%%%% FDD Programe %%%%%%%%%%%% clear load sec1 data;% load data AY0 = data(:,[1 2 7 8 3:6]); % [B,A] = butter(4,0.99,'low'); % AY = filter(B,A,AY0); Wp=20/256;Ws=30/256;Rp=1;Rs=20; [n,Wn]=buttord(Wp,Ws,Rp,Rs); [filt_num,filt_den] = butter(n,Wn); AY= filter(filt_num,filt_den,AY0); AY(:,[4,6,8]) = AY(:,[4,6,8]).*(-1); % case10 AY = detrend(AY); fs = 512; % sampling frequency NFFT = 4096;% number of FFT channel = [1:8]; G0 = zeros((NFFT/2+1),length(channel)^2); j2 = 0; for i=1:length(channel) for j=1:length(channel) j2 = j2+1; [Pff2,f1] = cpsd(AY(:,channel(i)),AY(:,channel(j)),NFFT,NFFT/4,NFFT,fs); G0(:,j2) = Pff2; end end G = zeros(length(channel),length(channel)); for i=1:(NFFT/2+1) for i2=1:length(channel) G(i2,:) = G0(i,((i2-1)*length(channel)+1):((i2-1)*length(channel)+length(channel))); end [U,S,V] = svd(G); UU(((i-1)*length(channel)+1):((i-1)*length(channel)+length(channel)),:) = U; SS(((i-1)*length(channel)+1):((i-1)*length(channel)+length(channel)),:) = S; end for i=1:length(channel) for j=1:(NFFT/2+1) SS0(j,i) = SS(((j-1)*length(channel)+i),i); end end for i=1:8 for j=1:(NFFT/2+1) SS1(j,i) = 20*log10(SS0(j,i)); end end figure subplot(211) plot(f1(1:200)',SS1(1:200,:)) title('Singular value'),xlabel('Hz') subplot(212) plot(1:200,SS1(1:200,:)) title('Singular value') index = input('Enter the number of cutting off -- -1 = stop: ') index0 = 0; index1 = 0; ii = 0; while index>0 ii = ii+1; index1 = index1+index-index0; [l,index2] = max(SS1((index0+1):index1,1)); frequency(ii,1) = f1(index2+index0); shape1 = UU(((index2+index0)*length(channel)+1):((index2+index0)*length(channel)+length(channel)),:); shape2 = shape1; shape3((((ii-1)*size(channel,2)+1):((ii-1)*size(channel,2)+size(channel,2))),:) = shape2; index0 = index1; index = input('Enter the number of cutting off -- -1 = stop: ') end for i=1:ii shape(:,i) = shape3((((i-1)*size(channel,2)+1):((i-1)*size(channel,2)+size(channel,2))),1); end for i=1:ii for j=1:size(shape,1) shape(j,i) = sign(real(shape(j,i)))*abs(shape(j,i)); end end frequency for i=1:size(shape,2) figure plot(1:length(channel),shape(:,i),'-*') title(['Modal shape',',Frequency=',num2str(frequency(i))]) end