%%%%%%%%%%-----FDD Programe----%%%%%%%%%%%%
clear all
close all
load matlab.mat 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; %通带20Hz，阻带30Hz
[n,Wn]=buttord(Wp,Ws,Rp,Rs); %计算滤波器的阶数n和滤波参数
[filt_num,filt_den] = butter(n,Wn); %滤波器设计
AY= filter(filt_num,filt_den,AY0); %滤波操作
AY(:,[4,6,8])= AY(:,[4,6,8]).*(-1); % 取相反数
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
评论0