%  eegsim.m
%  Generates a simulated eeg signal and leaves it in the workspace
%  Does time-frequency analysis and creates two figures
%	This program is distributed as a supplement to the book
%	"Biomedical Signal Processing and Signal Modeling" by E. N. Bruce,
%	published by Wiley, 2000.  It is provided for educational use only.
%  While every effort has been made to insure its suitability to illustrate
%  principles described in the above book, no specific feature or capability 
%  is implied or guaranteed.
rand('uniform')
%generate 31-interval PRB sequence using shift register
reg=round(rand(1,5));
prbase=zeros(1,31);
for j=1:31
	reg1=reg(1,3)|reg(1,5);
	if reg(1,3)==reg(1,5)
		reg1=0;
		end
	for i=5:-1:2
		reg(1,i)=reg(1,i-1);
		end
	reg(1,1)=reg1;
	prbase(j)=reg(1,5);
	end
pr=[1,prbase,0];
%  Generate filters for 3 eeg bands
%  and 80 seconds of eeg data for each band
[awb,awa]=butter(8,[.33 .5]);
[s2b,s2a]=butter(8,[.16 .33]);
[s3b,s3a]=butter(8,.16);
in=randn(8192,1);
out=filter(awb,awa,in);
alpha=0.71*out(4097:8192);
in=randn(8192,1);
out=filter(s2b,s2a,in);
theta=out(4097:8192);
in=randn(8192,1);
out=filter(s3b,s3a,in);
delta=out(4097:8192);
%Switch between states randomly
%Set transition states at 0 sec
eeg=zeros(3*4096,1); ke=1; nst(1)=1; delst=1000;
eeg(1:delst)=alpha(1:delst);jg=1;
gan=[0,1,0;1,0,0;0,0,1];
for je=2:10
	if nst(je-1)>=3*4096, break,end;
	addon=0;
	for jf=1:3
		if pr(ke)==1,
		addon=addon+delst;
		end % if
		ke=ke+1;
	end % jf for
	nst(je)=nst(je-1)+addon;
	if nst(je)>3*4096, nst(je)=3*4096; end;
	ginx=3-rem(jg,3);
	g1=gan(ginx,1);g2=gan(ginx,2);g3=gan(ginx,3);
eeg(nst(je-1)+1:nst(je))=g1*alpha(1:(nst(je)-nst(je-1)))+g2*theta(1:(nst(je)-nst(je-1)))+g3*delta(1:(nst(je)-nst(je-1)));
	jg=jg+1;
end % je for
spec=tfar(eeg,15,length(eeg)/50,length(eeg)/50,.5,15,50);