sim_eeg.m

eeg segmentation codes...

function x=sim_eeg(t,fs,f,u)
% x=sim_eeg(t,fs,f,u)
% is a function using linear modeling for 
% generating an eeg wave specified by the user
% t  == segment-time (sec) of the eeg wave
% fs == sampling frequency (Hz)
% f  == [fl fh], the low and high frequencies
% (Hz) determined by the user for the required
% eeg wave.
% We recommend the folloing values for f:
% f=[0.50 3.50] ---> Delta  wave
% f=[4.00 7.00] ---> Theta  wave
% f=[7.5 12.00] ---> Alpha  wave
% f=[12.5 15.0] ---> Sigmma wave
% f=[15.5 20.0] ---> Beta   wave
% f=[20.0 40.0] ---> Gamma  wave  
% u == the model identifier
% u=[0 L] for Moving Average Model (MA(L))
% if u=[0], the default model is MA(64)
% u=[1 L] for Autoregressive Moving Average
% ARMA(L,L)) if u=[1] the default model is
% ARMA(4,4) 
%Reference:
%R. R. Gharieb and A. Cichocki, "Segmentation and tracking of
%the electro-encephalogram signal using an adaptive recursive
%bandpass filter, " Journal of Int. Federation for Medical & 
%Biological Engeneering & Computing, vol. 39, pp. 237-248, 
%March 2001
%*************************************************
% Please, for any comments, contact:
% R. R. Gharieb and A. Cichocki
% Lab. for Advanced Brain Signal Processing
% Brain Science Institute
% Riken, 2-1 Hirosawa, Wako-shi, Saitama 351-0198
% Japan
% e-mail {reda, cia}@bsp.brain.riken.go.jp
%*************************************************
if nargin<=3
   error('not enough arguments for sim_eeg function')
end
if u(1,1)>1
   error('model type must be 0 (i.e., MA) or 1 (i.e., ARMA)')   
end
if u(1,1)<0
   error('model type must be 0 (i.e., MA) or 1 (i.e., ARMA)')
end
if length(f)<2
   error('low and high frequencies must be determined')
   end
%number of samples n
n=fix(t*fs);
%Digital frequency
f=2*f./fs;
% adjusting the default model and its order
[r c]=size(u);
if r==c;switch u(1,1)
   case 0;u=[u 64];
   case 1;u=[u 4];
   end
end
[r c]=size(u);
if r>c;u=u'
end
switch u(1,1)
case 0
   a=fir1(u(1,2),f);
   b=1;
case 1
   [a b]=butter(u(1,2),f);
end
x=randn(2*n,1);
x=filter(a,b,x);
x=x(n+1:2*n);x=x./std(x);