function []=Remove_Lines();
% *************************************************************************
% 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 frontal; global dfrontal;
global Fs;  global ss; global sschop;

DFR=frontal(:,13);
dDFR=dfrontal(:,13);

ssrml=round([ss(3)*.3691 ss(4)*.0247 ss(3)*.6299 ss(4)*.8854]);
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];

figure
subplot(2,1,1); plot(dDFR); 
set(gca,'XTick',[0:1000:4000],'XTickLabel',[(0:1000:4000)/Fs],'YLim',[-.25 .25]);
ylabel('mV'); xlabel('seconds'); title('dDFR');

[SDFR,f,Serr]=mtspectrumc(dDFR,params);
subplot(2,1,2); plot(f,log10(SDFR));
set(gca,'XLim',[0 params.Fs/2]);
xlabel('Hz'); ylabel('log10(mV^2)'); title('SDFR');
set(gcf,'Units','pixels','Position',sschop);

                          %generate figure displaying power spectra for one channel from 
                          %each site in the Right Hemisphere;
                          %mtspectrumc(data,params)
                          %after detrending the data with locdetrend
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];
    
display(['Type ''return'' when ready to continue']);
keyboard                  %pause
close all
                          %line removal using the F-test to find all highly
                          %significant harmonic features;
                          %these may be noise or they may be physiology; using dDFR channel from
                          %the right hemisphere; the bonferroni
                          %correction is used
dnlfrontal=rmlinesc(dDFR,params,.05,'y');
set(gcf,'Units','pixels','Position',ssrml);

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

                          %line removal using the F-test to find all highly
                          %significant harmonic features; If the wrong
                          %significance level is chosen, it's possible to
                          %throw out signal, not just noise
dnlfrontal=rmlinesc(dDFR,params,.5,'y');
set(gcf,'Units','pixels','Position',ssrml);

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

f0=60;
                          %line removal investigating a narrow region in
                          %the spectrum surrounding 60 Hz
dnlfrontal=rmlinesc(dDFR,params,.05,'y',f0);
set(gcf,'Units','pixels','Position',ssrml);

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

f0=60.025;
                          %line removal investigating a narrow region in
                          %the spectrum surrounding 60.025 Hz
dnlfrontal=rmlinesc(dDFR,params,.05,'y',f0);
set(gcf,'Units','pixels','Position',ssrml);

f0=180.075;
                          %line removal investigating a narrow region in
                          %the spectrum surrounding 180.075 Hz. Use the
                          %output of the previous rmlinesc run so that both
                          %60.025 Hz and 180.075 Hz are eliminated in
                          %the output. 180.075 is the third harmonic of 60.025
                          %Hz.
dnlfrontal2=rmlinesc(dnlfrontal,params,.05,'y',f0);
set(gcf,'Units','pixels','Position',ssrml);

f0=199.88;
                          %line removal investigating a narrow region in
                          %the spectrum surrounding 199.88 Hz. Use the
                          %output of the previous two rmlinesc runs so that
                          %60.025, 180.075 and 199.88 Hz are all eliminated in
                          %the output. 
dnlfrontal3=rmlinesc(dnlfrontal2,params,.05,'y',f0);
set(gcf,'Units','pixels','Position',ssrml);

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

figure

[nolineSDFR,f,Serr]=mtspectrumc(dnlfrontal3,params);
subplot(2,1,1); plot(f,log10(SDFR));
set(gca,'XLim',[0 params.Fs/2]);
xlabel('Hz'); ylabel('log10(mV^2)'); title('SDFR');

subplot(2,1,2); plot(f,log10(nolineSDFR));
set(gca,'XLim',[0 params.Fs/2]);
xlabel('Hz'); ylabel('log10(mV^2)'); title('rmlinesc(dDFR)');
set(gcf,'Units','pixels','Position',sschop);

                          % generate figure displaying power spectra for one channel from 
                          % Right Hemisphere (DFR);
                          % compares spectrum with the results produced by
                          % successive application of rmlinesc
                          % after detrending the data with locdetrend

display(['Hit UPARROW or type ''return'' when ready to continue']);
keyboard
close all
                          %line removal investigating a region in
                          %the spectrum surrounding 60.025 Hz, but with
                          %params.pad = 1, instead of params.pad = 5 as was
                          %used in the previous examples
params.pad=1;
f0=60.025;
dnlfrontal=rmlinesc(dDFR,params,.05,'y',f0);
set(gcf,'Units','pixels','Position',ssrml);
    

return;