function []=DeNoise_SVD();
% *************************************************************************
% Written for Dynamical Neuroscience XIV, 36th Annual Meeting of the Society for Neuroscience, 2006,
% Atlanta, GA., by K.Purpura, Dept. Neurology & Neuroscience, 
% Weill Medical College of Cornell University, NYC.
% *************************************************************************

close all

global Fs
global ss;  global sschop
global ATL; global ATR;
global PTL; global PTR;
global DFL; global DFR;
global VFL; global VFR;
global APL; global APR;
global MTL; global MTR;
global frontal; global dfrontal; 
global parietal; global dparietal;
global temporal; global dtemporal;
global MT; global dMT;

figure                      %generate figure displaying all channels for the Right Hemisphere frontal lobe    
subplot(1,2,1); plot(DFR); title('DFR');
set(gca,'XTick',[0:1000:4000],'XTickLabel',[(0:1000:4000)/Fs]);
ylabel('mV'); xlabel('seconds');
subplot(1,2,2); plot(VFR); title('VFR');
set(gca,'XTick',[0:1000:4000],'XTickLabel',[(0:1000:4000)/Fs]);
xlabel('seconds');
set(gcf,'Units','pixels','Position',sschop);
h1=gcf;

display(['Type ''return'' when ready to continue']);
keyboard
                            %SVD of data organized into brain regions;
                            %multivariate analysis is in terms of time x
                            %channel
[PU,PS,PV]=svd(parietal,0); %left and right, 12 channels
[FU,FS,FV]=svd(frontal,0);  %left and right, dorsal and ventral, 24 channels
[MU,MS,MV]=svd(MTL,0);      %left and right, 12 channels
[TU,TS,TV]=svd(temporal,0); %left and right, anterior and posterior, 24 channels

                            %SVD of data organized into brain regions;
                            %multivariate analysis is in terms of time x
                            %channel; here the detrended data are used
[dPU,dPS,dPV]=svd(dparietal,0);
[dFU,dFS,dFV]=svd(dfrontal,0);
[dMU,dMS,dMV]=svd(dMT,0);
[dTU,dTS,dTV]=svd(dtemporal,0);

FrontalR=[DFR VFR];             %subset of frontal, right hemisphere channels only
[FU,FS,FV]=svd(FrontalR,0);     %SVD

figure                          %figure of singular values of FrontalR
plot(log10(diag(FS)))
hold
plot(log10(diag(FS)),'*')
ylabel('log10'); xlabel('row in S');
title('Singular values for svd of FrontalR (12 EcoG channels)')
set(gcf,'Units','pixels','Position',sschop);
h2=gcf;

display(['Hit UPARROW or type ''return'' when ready to continue']);
keyboard                        %pause

                                %generate the multi-taper spectra of each
                                %of the temporal modes in FU.
Ktapers=20;
NW=(Ktapers+1)/2;
params.tapers=[NW Ktapers];
params.Fs=500;
params.fpass=[0 params.Fs/2];
params.pad=4;
params.err=[2 .05];
data=FU;
wh=waitbar(0,'Spectral Analysis of Temporal Modes');
for ispec=1:size(data,2)
    waitbar(ispec/size(data,2),wh);
    [S,f,Serr]=mtspectrumc(data(:,ispec),params);
    SData(:,ispec)=S(1:200:end);
    SerrData(:,ispec)=Serr(1,1:200:end)';
    clear S; clear Serr;
end
close(wh)
fsub=f(1:200:end);

figure
plot(fsub,log10(SData(:,1)))
hold
plot(fsub,log10(SData(:,2)),'r')
plot(fsub,log10(SData(:,3)),'k')
plot(fsub,log10(SData(:,4)),'m')
plot(fsub,log10(SData(:,5)),'c')
legend('Temporal mode 1','mode 2','mode 3','mode 4','mode 5')
ylabel('log10(mV^2)')
xlabel('Hz')
title('Spectra of the first 5 temporal modes from FrontalR (12 EcoG channels)')
set(gcf,'Units','pixels','Position',sschop);

display(['Hit UPARROW or type ''return'' when ready to continue']);
keyboard
                                    %reconstruct the time series
FrontalRsvd=FU(3:end,3:end)*FS(3:end,3:end)*FV(3:end,3:end)';
figure                              %plot one channel of the reconstructed data
subplot(2,1,1)
plot(FrontalRsvd(:,1))
set(gca,'YLim',[-.15 .15])
set(gca,'Xtick',[0:500:4000],'XTickLabel',(0:500:4000)/params.Fs);
xlabel('seconds'); ylabel('mV'); title('Example of reconstructed time series');
set(gcf,'Units','pixels','Position',sschop);
                                    %spectrum of one channel of the
                                    %reconstructed data
Ktapers=8;
NW=(Ktapers+1)/2;
params.tapers = [ NW Ktapers ];
params.Fs=500;   
params.fpass = [ 0 params.Fs/2];
params.pad = 5;
params.err=[2 .05];
[S,f]=mtspectrumc(FrontalRsvd(:,1),params);
subplot(2,1,2); plot(f,log10(S));
set(gca,'XLim',[0 params.Fs/2]); set(gca,'YLim',[min(log10(S)) max(log10(S))])
ylabel('log10(mV^2)'); xlabel('Hz'); title('Spectrum of reconstructed time series');


display('Done');

return;