function []=Locdetrend_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 frontal;
global Fs;  global ss;

good_input=0;
while ~good_input
    channel=input('Choose a frontal channel to work with (1-24): '); 
    if (channel>0 & channel<25)
        good_input=1;
        if (rem(channel,1)==0)
            good_input=1;
        else
            good_input=0;
        end    
    end
end
duration=length(frontal(:,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=0.5;
    end    
end

y=frontal(:,channel);
n=window*Fs;
dn=winstep*Fs;
nt=length(y);
y_line=zeros(nt,1);
norm=y_line;
nwin=ceil((nt-n)/dn);
yfit=zeros(nwin,n);
xwt=((1:n)-n/2)/(n/2);
wt=(1-abs(xwt).^3).^3;
h1=figure;
h2=figure;
set(h1,'position',[ss(3)/50 ss(4)/2 ss(3)/2 ss(4)/2]);
set(h2,'position',[(ss(3)-ss(3)/2) 1 ss(3)/2 ss(4)]);

for j=1:nwin 
	tseg=y(dn*(j-1)+1:dn*(j-1)+n);
	y1=mean(tseg);
	y2=mean((1:n)'.*tseg)*2/(n+1);
	a=(y2-y1)*6/(n-1); b=y1-a*(n+1)/2;
	yfit(j,:)=(1:n)*a+b;
	y_line((j-1)*dn+(1:n))=y_line((j-1)*dn+(1:n))+(yfit(j,:).*wt)';
	norm((j-1)*dn+(1:n))=norm((j-1)*dn+(1:n))+wt';
    
    figure(h1)
    plot(tseg); hold 
    plot(1:n,y1*ones(1,n),'LineWidth',2,'Color',[1. 0 0]);
    plot(yfit(j,:),'LineWidth',2,'Color',[0 0 0]); hold
    title(['Window Number: ' int2str(j)]);
    legend('EcoG','mean','Yfit');
    
    figure(h2)
    subplot(3,1,1)
    plot(norm)
    title('WeightingFunction')
    subplot(3,1,2)
    plot(y_line); 
    if j==1, hold; end
    plot((j-1)*dn+1,yfit(j,1),'*'); plot((j-1)*dn+n,yfit(j,end),'*');
    plot((j-1)*dn+(1:n),yfit(j,1:end),'LineWidth',2);
    title('YFit * WeightingFunction')
    pause(2);
end
mask=find(norm>0); y_line(mask)=y_line(mask)./norm(mask);
indx=(nwin-1)*dn+n-1;
npts=length(y)-indx+1;
y_line(indx:end)=(n+1:n+npts)'*a+b;
figure(h2)
subplot(3,1,2)
plot(indx:length(y_line),y_line(indx:end),'LineWidth',2);
subplot(3,1,3)
plot(y_line,'LineWidth',2);
title('Finalfit = (YFit * WeightingFunction)/WeightingFunction')
pause(2)
figure(h1)
plot(y); hold 
plot(y_line,'LineWidth',2,'Color',[1 0 0]);
legend('EcoG','FinalFit');
title(['Moving Window: Window = ' num2str(window) ' secs; Winstep = ' num2str(winstep) 'secs'])
pause(2)
h3=figure;
set(h3,'position',[1 1 ss(3)/2 ss(4)/2]);
plot(y-y_line)
themax=max(abs(y-y_line))*2;
title('EcoG - FinalFit')
set(gca,'YLim',[-themax themax]);

return;