function []=DeNoise_Demo();
% *************************************************************************
% 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 frontal; global dfrontal; 
global parietal; global dparietal;
global temporal; global dtemporal;
global MT; global dMT;

good_input=0;
while ~good_input
    display('Regions: frontal-->1; parietal-->2; temporal-->3; medial-temporal-->4'); 
                      
    region=input('Choose a brain region to work with (1-4):  '); 
    if (region>0 & region<5)
        good_input=1;
        if (rem(region,1)==0)
            good_input=1;
        else
            good_input=0;
        end    
    end
end

switch region
    case 1
        useregion=frontal;
    case 2
        useregion=parietal;
    case 3
        useregion=temporal;
    case 4
        useregion=MT;
end    

YN=1; godetrend=0;
while ~ischar(YN)
    YN=input('Should we detrend the data (Y/N)?  ','s');
    if ((YN=='Y')|(YN=='y'))
        godetrend=1;
    elseif ((YN=='N')|(YN=='n'))
        godetrend=0;
    else
        YN=1;
    end
end

if godetrend
    duration=length(useregion(:,1))./Fs;
    good_input=0;
    display(['length of recorded segments(secs): ' int2str(duration) ]);
    while ~good_input
        window=input('Choose a window length (secs):  '); 
        if (window>0 & window<duration)
            good_input=1;
        end
        if isempty(window)
            good_input=1; window=1;
        end    
    end
    good_input=0;
    while ~good_input
        winstep=input('Choose a window step size (secs); < window length:  '); 
        if (winstep>0 & winstep<=window)
            good_input=1;
        end
        if isempty(winstep)
            good_input=1; winstep=window/2;
        end    
    end
    data=locdetrend(useregion,Fs,[window winstep]);
elseif godetrend==0
    data=useregion;
end    

YN=1; gormlines=0;
while ~ischar(YN)
    YN=input('Should we remove line-noise (Y/N)?  ','s');
    if ((YN=='Y')|(YN=='y'))
        gormlines=1;
    elseif ((YN=='N')|(YN=='n'))
        gormlines=0;
    else
        YN=1;
    end
end

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];

if gormlines
    good_input=0;
    while ~good_input
        tapers=input('Choose the number of tapers:  '); 
        if tapers>0
            good_input=1;
        end
        if isempty(tapers)
            good_input=1; tapers=Ktapers;
        end    
    end
    good_input=0;
    while ~good_input
        padding=input('Choose a padding factor:  '); 
        if padding>0
            good_input=1;
        end
        if isempty(padding)
            good_input=1; padding=params.pad;
        end    
    end
    good_input=0;
    while ~good_input
        pvalue=input('Choose a significance level for F-test:  '); 
        if pvalue>0
            good_input=1;
        end
        if isempty(pvalue)
            good_input=1; pvalue=0.05;
        end    
    end
    NW=(tapers+1)/2;
    params.tapers = [ NW tapers ];
    params.pad = padding;
    Cdata=rmlinesc(data,params,pvalue,'n');
elseif gormlines==0
    Cdata=data;
end    

[CU,CS,CV]=svd(Cdata,0);        %SVD

figure                          %figure of singular values of data
plot(log10(diag(CS)))
hold
plot(log10(diag(CS)),'*')
ylabel('log10'); xlabel('row in S');
title('Singular values for svd of chosen data')
set(gcf,'Units','pixels','Position',sschop);

wh=waitbar(0,'Spectral Analysis of Temporal Modes');
for ispec=1:size(CU,2)
    waitbar(ispec/size(CU,2),wh);
    [S,f,Serr]=mtspectrumc(CU(:,ispec),params);
    SCU(:,ispec)=S(1:200:end);
    SerrCU(:,ispec)=Serr(1,1:200:end)';
    clear S; clear Serr;
end
close(wh)
fsub=f(1:200:end);

figure
plot(fsub,log10(SCU(:,1)))
hold
plot(fsub,log10(SCU(:,2)),'r')
plot(fsub,log10(SCU(:,3)),'k')
plot(fsub,log10(SCU(:,4)),'m')
plot(fsub,log10(SCU(:,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')
set(gcf,'Units','pixels','Position',sschop);


display(['Type ''return'' when ready to exit']);
keyboard


return;