Patter Page 1 Patter Page 2 Patter Page 3

Patter page 1
Patter page 2
Patter page 3

A-8

B-8

C-10

D-10

E-2

F-4

Fig. G-1: RUN Macro

function run(i)

% run figure i

% bands: 1=delta, 2=theta, 3=alpha, 4=beta, 5=gamma

% i fig. description

% 1- 5 A-1- 5 FIR filters per band, mag, phase, imp, step

% 6-10 A-6-10 FIR filters per band, pole/zero plot

% 11-15 B-1- 5 IIR filters per band, mag, phase, imp, step

% 16-20 B-6-10 IIR filters per band, pole/zero plot

% 21-30 C-1-10 FIR filters per epoch, all bands

% 31-40 D-1-10 IIR filters per epoch, all bands

% 41-50 E-1-10 Power spectra per epoch

% 51 firpower(epoch,band) band/total

% 52 iirpower(epoch,band) band/total

% 53 fftpower(epoch,band) band/total

Fs=128;

bandtab ={'Delta' 'Theta' 'Alpha' 'Beta' 'Gamma'};

F0tab =[ 0.1 4 8 13 20 ];

F1tab =[ 4 8 12 20 64 ];

firerrtab=[ 1.08 1.08 1.06 1.03 1.02 ];

iirerrtab=[ 1.02 1.03 1.02 1.01 1.01 ];

rangetab =[100 20 20 20 20 ];

j = i;

if 0 < j & j <= 5

sys = designfir(Fs, F0tab(j), F1tab(j), firerrtab(j));

filtfig(Fs, sys, bandtab{j}, 'FIR', j, 'A');

end;

j = i - 5;

if 0 < j & j <= 5

sys = designfir(Fs, F0tab(j), F1tab(j), firerrtab(j));

pzfig(Fs, sys, bandtab{j}, 'FIR', j+5, 'A');

end;

j = i - 10;

if 0 < j & j <= 5

sys = designiir(Fs, F0tab(j), F1tab(j), iirerrtab(j));

filtfig(Fs, sys, bandtab{j}, 'IIR', j, 'B');

end;

j = i - 15;

if 0 < j & j <= 5

sys = designiir(Fs, F0tab(j), F1tab(j), iirerrtab(j));

pzfig(Fs, sys, bandtab{j}, 'IIR', j+5, 'B');

end;

j = i - 20;

if 0 < j & j <= 10

for k = 1:5

systab(k) = designfir(Fs, F0tab(k), F1tab(k), firerrtab(k));

end;

x=importdata(strcat('EEG_',num2str(j)));

epochfig(Fs, systab, rangetab, bandtab, 'FIR', x, j, 'C');

end;

j = i - 30;

if 0 < j & j <= 10

for k = 1:5

systab(k) = designiir(Fs, F0tab(k), F1tab(k), iirerrtab(k));

end;

x=importdata(strcat('EEG_',num2str(j)));

epochfig(Fs, systab, rangetab, bandtab, 'IIR', x, j, 'D');

end;

j = i - 40;

if 0 < j & j <= 10

x=importdata(strcat('EEG_',num2str(j)));

spectrumfig(Fs, x, j, 'E');

end;

if i == 51

for k = 1:5

systab(k) = designfir(Fs, F0tab(k), F1tab(k), firerrtab(k));

end;

firpower=[];

for j = 1:10

x=importdata(strcat('EEG_',num2str(j)));

firpower = vertcat(firpower, epoch(systab, x));

end;

firpower

end;

if i == 52

for k = 1:5

systab(k) = designiir(Fs, F0tab(k), F1tab(k), iirerrtab(k));

end;

iirpower=[];

for j = 1:10

x=importdata(strcat('EEG_',num2str(j)));

iirpower = vertcat(iirpower, epoch(systab, x));

end;

iirpower

end;

if i == 53

for j = 1:10

x=importdata(strcat('EEG_',num2str(j)));

X=fft(x);

Fn=Fs/2;

Ln=length(X)/2;

T=Ln/Fn;

Xm2=X(1:Ln).*conj(X(1:Ln));

for k=1:5

fftpower(j,k)=sum(Xm2(F0tab(k)*T+1:F1tab(k)*T))/sum(Xm2);

end;

end;

fftpower

end;

Fig. G-2: DESIGNFIR Macro

function sys = designfir(Fs, F0, F1, e)

% Design optimal fir bandpass filter

% Fs is sample frequency, F0 to F1 is passband

% e is allowed error, eg e=1.05 for 5% error

Fn=Fs/2; % Nyquist frequency

if F0/Fn < 0.01; F0= 0; end; % LPF for F0 near 0

if F1/Fn > 0.80; F1=Fn; end; % HPF for F1 near Fn

Ap=20*log10(e)*2; % passband dev (dB)

As=-20*log10((e-1)/(Fn/(F1-F0)-1)*2); % stopband dev (dB)

dp=(10^(Ap/20)-1)/(10^(Ap/20)+1); % passband dev

ds=10^(-As/20); % stopband dev

M=[0 1 0]; % magnitude of bands

F=[F0/e F0*e F1/e F1*e]; % band edges

dev=[ds dp ds]; % deviation in bands

if F0 == 0; M=M(2:3); F=F(3:4); dev=dev(2:3); end; % trim for LPF

if F1 == Fn; M=M(1:2); F=F(1:2); dev=dev(1:2); end; % trim for HPF

[N F0 M0 W]=firpmord(F, M, dev, Fs); % determine filter order

[b delta]=firpm(N, F0, M0, W); % compute coefficients

a=[1 zeros(1,N)]; % unity gain

sys = tf(b, a, 1/Fs); % package up for travel

Fig. G-3: DESIGNIIR Macro

function sys = designiir(Fs, F0, F1, e)

% Design elliptical, bilinear, iir bandpass filter

% Fs is sample frequency, F0 to F1 is passband

% e is allowed error, eg e=1.05 for 5% error

Fn=Fs/2; % Nyquist frequency

if F0/Fn < 0.01; F0= 0; end; % LPF for F0 near 0

if F1/Fn > 0.80; F1=Fn; end; % HPF for F1 near Fn

Ap=20*log10(e)*2; % passband dev (dB)

As=-20*log10((e-1)/(Fn/(F1-F0)-1)*2); % stopband dev (dB)

Wp=[F0*e F1/e]/Fn; % passband, normalized

Ws=[F0/e^6 F1*e^6]/Fn; % stopband, normalized

type='bandpass'; % bandpass type

if F0 == 0; Wp=Wp(2:2); Ws=Ws(2:2); type='low'; end; % trim for LPF

if F1 == Fn; Wp=Wp(1:1); Ws=Ws(1:1); type='high'; end; % trim for HPF

[N Fc]=ellipord(Wp, Ws, Ap, As); % determine filter order

[b a]=ellip(N, Ap, As, Fc, type); % compute coefficients

b=b*e; % scale up for unity gain

sys = tf(b, a, 1/Fs); % package up for travel

Fig. G-4: FILTFIG Macro

function filtfig(Fs, sys, band, type, j, n)

% Print filter figure

% Fs is sample frequency

% sys is transfer function

% band, type strings for figure title, fig. n-j

Fn=Fs/2;

b = sys.num{1};

a = sys.den{1};

subplot(4,1,1);

[H f]=freqz(b, a, 2^16, Fs);

plot(f, 20*log10(abs(H)))

N=length(b)-1;

suffix={'st' 'nd' 'rd' 'th' 'th' 'th' 'th' 'th' 'th' 'th'};

order=sprintf('%d%s', N, suffix{mod(N-1,10)+1});

title(sprintf('Fig. %s-%d: %s Order, %s Band, %s Filter', ...

n, j, order, band, type));

axis([0 Fn -60 20]);

faxis=linspace(0,Fn,2^3+1);

set(gca,'XTick',faxis)

set(gca,'XTickLabel',faxis)

maxis=[-60:20:20];

set(gca,'YTick',maxis)

set(gca,'YTickLabel',maxis)

xlabel('Frequency (Hz)')

ylabel('Magnitude Response (dB)')

subplot(4,1,2);

[H f]=freqz(b, a, 2^16, Fs);

plot(f, 180/pi*angle(H))

axis([0 Fn -180 180]);

set(gca,'XTick',faxis)

set(gca,'XTickLabel',faxis)

paxis=[-180:90:180];

set(gca,'YTick',paxis)

set(gca,'YTickLabel',paxis)

xlabel('Frequency (Hz)')

ylabel('Phase Response (°)')

subplot(4,1,3);

[y t] = impz(b, a, 4*Fs, Fs);

plot(t,y);

axis([0 4 -0.4 0.4]);

iaxis=[-0.4:0.2:0.4];

set(gca,'YTick',iaxis)

set(gca,'YTickLabel',iaxis)

xlabel('Time (s)')

ylabel('Impulse Response')

subplot(4,1,4);

[y t] = stepz(b, a, 4*Fs, Fs);

plot(t,y);

axis([0 4 -.5 1.5]);

saxis=[-.5:0.5:1.5];

set(gca,'YTick',saxis)

set(gca,'YTickLabel',saxis)

xlabel('Time (s)')

ylabel('Step Response')

Fig. G-5: PZFIG Macro

function pzfig(Fs, sys, band, type, j, n)

% Print pole-zero figure

% Fs is sample frequency

% sys is the transfer function

% band, type strings for figure title, fig. n-j

Fn=Fs/2;

b = sys.num{1};

a = sys.den{1};

[z p k] = zpkdata(sys);

if dot(p{1},p{1})==0; p{1}=[]; end;

zplane(z{1},p{1});

suffix={'st' 'nd' 'rd' 'th' 'th' 'th' 'th' 'th' 'th' 'th'};

N=length(b)-1;

order=sprintf('%d%s', N, suffix{mod(N-1,10)+1});

title(sprintf('Fig. %s-%d: %s Order, %s Band, %s Filter', ...

n, j, order, band, type));

Fig. G-6: EPOCHFIG Macro

function epochfig(Fs, systab, rangetab, bandtab, type, x, j, n)

% Print epoch figure

% Fs is sample frequency

% systab is transfer functions per band

% rangetab is range for y (amplitude) scale

% x is time series

% bandtab, type strings for figure title, fig. n-j

Ts=1/Fs;

Ls=length(x)*Ts;

power = epoch(systab, x);

subplot(6,1,1);

[y t] = lsim(tf(1, 1, Ts),x);

plot(t, y);

suffix={'st' 'nd' 'rd' 'th' 'th' 'th' 'th' 'th' 'th' 'th'};

number=sprintf('%d%s', j, suffix{mod(j-1,10)+1});

title(sprintf('Fig. %s-%d: %s Filters, %s Epoch\nOriginal Signal', ...

n, j, type, number));

axis([0 Ls -rangetab(1) rangetab(1)]);

ylabel('Amplitude');

for k = 1:5

subplot(6,1,k+1);

[y t] = lsim(systab(k),x);

plot(t, y);

title(sprintf('%s Band, %.2f%% Power', ...

bandtab{k}, power(k)*100));

axis([0 Ls -rangetab(k) rangetab(k)]);

ylabel('Amplitude');

if k == 5

xlabel('Time (s)');

end;

end;

Fig. G-7: SPECTRUMFIG Macro

function spectrumfig(Fs, x, j, n)

% Print frequency specturm of epoch figure

% Fs is sample frequency

% x is time series

% fig. n-j

X=fft(x);

Fn=Fs/2;

Ls=length(X);

T=Ls/Fs;

Ln=Ls/2;

Xm2=X(1:Ln).*conj(X(1:Ln));

Xm2=Xm2/sum(Xm2)*Ln;

%scatter(0:Ln-1,10*log10(Xm2),'.','SizeData',1);

scatter(0:Ln-1,10*log10(Xm2),'.','MarkerEdgeColor','yellow');

set(gca,'Box','on');

suffix={'st' 'nd' 'rd' 'th' 'th' 'th' 'th' 'th' 'th' 'th'};

number=sprintf('%d%s', j, suffix{mod(j-1,10)+1});

title(sprintf('Fig. %s-%d: Power Spectrum, %s Epoch', ...

n, j, number));

axis([0 Ln -60 20]);

xlabel('Frequency (Hz)')

ylabel('Magnitude (dB)')

set(gca,'XTick',linspace(0,Ln,1+2^3))

set(gca,'XTickLabel',linspace(0,Fn,1+2^3))

hold on

Bs=32;

Bn=Bs/2;

Bm2=sum(reshape(Xm2,Ln/Bn,Bn),1)*Bn/Ln;

Rm2=reshape(repmat(Bm2,Ln/Bn,1),Ln,1);

plot(0:Ln-1,10*log10(Rm2));

hold off

Fig. G-8: EPOCH Macro

function power = epoch(systab, x)

% measure power for bands in systab of time series x

for k = 1:5;

[y t] = lsim(systab(k),x);

power(k)=dot(y,y)/dot(x,x);

end;

G-6