频域分解法FDD，增强的频域分解

%%%%%%%%%%-----FDD Programe----%%%%%%%%%%%% 
 
clear all
close all

load sec1 data;% load  data 
AY0= data（:[1 2 7 8 3:6]）;   %列向量
% [BA] = butter（40.99‘low‘）; 
% AY    = filter（BAAY0）; 
Wp=20/256;Ws=30/256;Rp=1;Rs=20;   %通带20Hz，阻带30Hz
[nWn]=buttord（WpWsRpRs）;     %计算滤波器的阶数n和滤波参数
[filt_numfilt_den] = butter（nWn）;  %滤波器设计
AY= filter（filt_numfilt_denAY0）;   %滤波操作
AY（:[468]）= AY（:[468]）.*（-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; 
        [Pff2f1]= cpsd（AY（:channel（i））AY（:channel（j））NFFTNFFT/4NFFTfs）; 
        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 
    [USV]= 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）