% *****************************************************************
% 求频谱
%
%
% 对trial先平均后做频谱
%
% *****************************************************************

%%

clear;clc;

addpath ..\..\tools
addpath ..\..\UnderlyingFunction

load ..\NN_dssall03_27_f6c6

triLen = 10;
DIMs=[triLen+1 2*triLen+1];

E1yb = yb;
clear yb
%%
E1yb = E1yb(2*fs+[1:triLen*fs],:,:,:);
f=1:fs*triLen;f=f-1;f=f/triLen;

E1spec = [2*abs(fft(E1yb))/size(E1yb,1)].^2;
E1topo = (E1spec(:,2:3:end,:,:)+E1spec(:,3:3:end,:,:))/2;

E1spec = mean(E1spec,4);
E1topo = mean(E1topo,4);

%%

E1spec = double(E1spec);
E1topo = double(E1topo);

clear pps
for dm=1:2
    DIM=DIMs(dm);
    for blo = 1:size(E1spec,3)
        for ch = 1:size(E1spec,2)
            ss = E1spec(DIM,ch,blo)/mean(E1spec(DIM+[-1 1],ch,blo));
            p  = fcdf(ss,2*16,4*16,'upper');
            pps(ch,blo,dm) = p;
        end
        disp(['E1dm' num2str(dm) 'blo' num2str(blo)]);
    end
end
for dm = 1:size(pps,3)
    for blo = 1:size(pps,2)
        [h, crit_p, pps(:,blo,dm)] = fdr_bh(pps(:,blo,dm));
    end
end
save('E1_sensor','pps')

clear pps
for dm = 1:2
    DIM=DIMs(dm);
    for blo = 1:size(E1topo,3)
        for ch = 1:size(E1topo,2)
            ss = E1topo(DIM,ch,blo)/mean(E1topo(DIM+[-1 1],ch,blo));
            p  = fcdf(ss,2*16,4*16,'upper');
            pps(ch,blo,dm) = p;
        end
        disp(['E1dm' num2str(dm) 'blo' num2str(blo)]);
    end
end
for dm = 1:size(pps,3)
    for blo = 1:size(pps,2)
        [h, crit_p, pps(:,blo,dm)]=fdr_bh(pps(:,blo,dm));
    end
end
save('E1_channel','pps')


%%

E1topom = squeeze(mean(E1topo,2));

DIMs = [triLen/2+1 triLen+1 2*triLen+1];
clear pps
for dm = 1:3
    DIM = DIMs(dm);
    for blo = 1:size(E1topo,3)
        ss  = E1topom(DIM,blo)/mean(E1topom(DIM+[-1 1],blo));
        p   = fcdf(ss,2*16,4*16,'upper');
        pps(blo,dm) = p;
    end
end
for blo = 1:size(pps,1)
    [h, crit_p, pps(blo,:)] = fdr_bh(pps(blo,:));
end
pps








