function wc=eeg_seg(x,fs)
% function wc=eeg_seg(x,fs)
% makes segmentation of eeg signal based 
% on the ceneter frequency of each segment.
% The function uses an adaptive bandpass filter;
% the adaptive filter is implemented as 4th order
% Butterworth; the adaptive filter needs only one
% coefficient to be updated. The center frequency
% of the eeg segment is proportional to this adaptive
% coefficient
% x  --> the eeg sequence
% fs --> the sampling rate for the eeg signal
% we recommend fs to be 100 up 160 Hz 
% w --> is the adaptive coefficient that yields 
% the segmentation paprameter
% Reference:
% R. R. Gharieb and A. Cichocki, "Segmentation
% and tracking of the electo-encephalogram signal
% using and adaptive recursive bandpass filter,"
% Medical & Biological Engineering & Computing, IFMBE-2001,
% vol. 39 pp. 237-248.
% for any comments, please, contact
% R. R. Gharieb and A.Cichocki
% Lab for Advanced Brain Signal Processing
% Brain Science Institute, Riken, Wako-shi
% Saitama 351-0198, JAPAN
% {reda, cia}@bsp.brain.riken.go.jp
[lr lc]=size(x);
if lc>lr;x=x';
   end
   l=length(x);
b=0.1;
k=cot(pi*b);
a0=1/(k*k+sqrt(2)*k+1);
a4=a0;
a2=-2*a0;
a1=-4*a0;
a3=a1;
a2d=4*a0;
b1=-2*k*a0*(2*k+sqrt(2));
b2=4*k*k*a0;
b2d=2*(k*k-1)*a0;
b3=2*k*(-2*k+sqrt(2))*a0;
b4=a0*(k*k-sqrt(2)*k+1);
fc=0.34;
b=0.1;
ny=0.95;
am=0.95;
[lpf1 lpf2]=butter(8,[0.18 2*fc-3/fs]);
r=1.0e2;
y(l)=0;
w=0.0;
wc(l)=0;
ala=y;
j=1:l;
v=cos(2*pi*fc*j)';
x=x./std(x);
x=x.*v;
x=filter(lpf1,lpf2,x);
for n=5:l;
    ma=[a0 a2 a4];xm=[x(n) x(n-2) x(n-4)]';
    ar=[-b1*w -(b2*w*w+b2d) -b3*w -b4];
    yr=[y(n-1) y(n-2) y(n-3) y(n-4)]';
    y(n)=ma*xm+ar*yr;
    ma=-[b1 2*b2*w b3];xm=[y(n-1) y(n-2) y(n-3)]'; 
    ar=-[b1*w (b2d+b2*w*w) b3*w b4];
    yr=[ala(n-1) ala(n-2) ala(n-3) ala(n-4)]';
    ala(n)=ma*xm+ar*yr;
    r=ny*r+ala(n)*ala(n);
    w=w+am*ala(n)*y(n)/r;
    if abs(w)>=1
      w=wc(n-1);
    end
    wc(n)=w;
end