This preview shows page 1. Sign up to view the full content.
Unformatted text preview: BME 6360
Neural Engineering Homework #9 Spring 2011
Prof. Wheeler The goal of this exercise is to practice using the autoregressive (AR) modeling technique. Students
should start with an understanding of the frequency content of signals and methods to extract it. This
exercise uses sample sections of recorded EEG signals which can be analyzed for frequency content in
order to identify different kinds of biological activity. The difficulties of doing this will also be apparent
as the exercise progresses.
Notes are included suggesting what should be submitted in your write up. This is not an easy exercise
and will require multiple sessions for its understanding.
The following files are available for this exercise:
ifftplot.m, fftplot.m, ftp.m, xcorcirc.m, ar2spectrum.m
The following data files are available:
norm.mat, muscle.mat, spike.mat;
The following filter parameter files are available:
lowpassa.mat
lowpassb.mat
All the data and coefficient vectors may be loaded by executing: load hw9data. If you execute this,
ignore the references to loading individual signals.
Part I. Exploration of the frequency domain differences among three sample EEG signals.
1. Display the real EEG signals. There are three EEG files for you to examine: norm, muscle, spike.
Respectively, they come from relatively normal recordings from an awake subject, from a recording
dominated by muscle potentials (e.g. during a scalp twitch or contraction), and from a recording in
which epileptic spikes, perhaps presaging a seizure, were occurring.
a. load the signals by executing successively the functions:
load norm
load muscle
load spike
Use the function ftp to investigate the nature of each signal and its Fourier transform.
ftp(norm);
% normal display
ftp(norm,1);
% set the DC component to zero  often the largest
component is the DC or average value, which is of little interest
Repeat for the other signals. Can you distinguish one from the other? In the time domain? In the
frequency domain? Is phase information useful? For your writeup indicate the differences among the
signals and the features useful for distinguishing them. 1 b. A low pass filter has been created for you and the coefficients are saved in lowpassa.cof and
lowpassb.cof. This is digital, 4th order, low pass Butterworth filter, with a cutoff of 0.1 times the
Nyquist rate, or at 1 cycle per 20 samples. Load these coefficients and filter the signals. (The design of
digital filters is material for another course.) Filter the EEG signals and investigate the spectra.
load lowpassa
load lowpassb;
normf=filter(lowpassb,lowpassa,norm);
ftp(normf,1);
% other EEG signals:
musclef=filter(lowpassb,lowpassa,muscle);
spikef=filter(lowpassb,lowpassa,spike);
Is it easier or harder to distinguish the signals after they have been filtered?
(If you want to try other filter settings, filter creation is easy:
[b,a]=butter(N,wn);
% the prestored coefficients were generated using: [lowpassb,lowpassa]=butter(4,0.1);
For your writeup indicate what the filter process has done and whether or not it makes the task of
distinguishing signals easier.
Include a few example plots illustrating the features which help you distinguish among the three signals.
You may choose either the filtered or unfiltered signals, and either the frequency or time domain
representations.
2. Correlation and power spectra of the EEG signals. Compute the autocorrelation of the filtered and
unfiltered versions of each signal.
a. ftp(xcorcirc(norm), 1);
ftp(xcorcirc(muscle), 1);
ftp(xcorcirc(spike), 1); b. ftp(xcorcirc(normf), 1);
ftp(xcorcirc(musclef), 1);
ftp(xcorcirc(spikef), 1); c. Compare the power spectra with Matlab supplied functions spa and ffplot.
ffplot(spa(normf));
ffplot(spa(norm));
ffplot(spa(spike));
ffplot(spa(muscl));
For your writeup: how do the power spectra of (a), (b), and (c) differ? Do these plots make it any easier
(as compared to part 1) to distinguish among the three signals?
Comment: The modeling process below will attempt to utilize the differences in these power spectra.
Considerable work can go into choosing exactly the right filter to do bring the differences out. At the
end of this exercise, you'll see the choices I made to exaggerate the differences.
2 II. Autoregressive (AR) Modeling of the EEG Signal Processes.
Matlab provides functions which do all the work for an autoregressive (AR) model of the EEG signals:
ar determines the AR model parameters; compare computes the output of the AR filter on the given
data and plots the results along with the original data; and present interprets the AR model information
in terms of the coefficients ai which have been discussed in class. Note that the signals have been
modified so as to exaggerate the differences in the predictive quality of the various AR models.
One goal of this part of the exercise is to indicate the duality of the derivation of the AR model. In part
1, you develop three different models which claim to be best at fitting the time domain data for three
different types of EEG signals. In 2 you demonstrate that this claim is reasonably met, while in 3 and 4
you demonstrate that the time domain predictions of epileptic spikes using the normal and muscle model
are not as good as when using the spike model. In 5 you demonstrate the implicit claim in AR modeling:
that the each model is not only the best fit to the time domain signal but also to the frequency spectrum
of the signal.
1. Determine the AR parameters for a 5th order model derived from a single epileptic spike, and from
the corresponding sequences from a normal eeg pattern and from muscle. I have chosen the data points
with indices from 170 to 190 in spike because they include three obvious epileptic spikes. The same
relative points were taken from norm and muscle only as a matter of convenience for this exercise.
ths=ar(spike(170:190), 5);
thn=ar(norm(170:190), 5);
thm=ar(muscle(170:190), 5);
present(ths); %model parameters for the spike sequence
%model parameters for a sequence from normal EEG data
%model parameters for a sequence from muscle activity %result: the parameters a(j) for the AR model: a(0)y(k)=Σy(kj)a(j)
%
remember Matlab indices start at 1 not 0 as in the theory ths, thn, and thm will be called the models for the spike, normal and muscle EEG data.
You should get:
0 = 1.000 y(k)  1.4495 y(k1) + 0.2868 y(k2) + 0.4012 y(k3)  0.0333 y(k4)  0.0793 y(k5)
or
y(k) = 1.4495 y(k1)  0.2868 y(k2)  0.4012 y(k3) + 0.0333 y(k4) + 0.0793 y(k5)
For your writeup, give the predictive equations for the normal and muscle data. 2. Test the results. You will plot out both the actual data, the model estimated data, and the error. Note
the FIT value, which is the mean squared error, a measure of how well the model predicts the data.
subplot(2,1,1)
estimate = compare(spike(120:250),ths,1,);
% this is a good fit  it should be
error = estimatespike(120:250);
% plots the estimate over a longer range (includes the EEG spikes
3 subplot(2,1,2)
plot(error);
max(error) % overlay the error signal For your writeup, submit the plots generated and the maximum errors in the estimation. Annotate the
graphs to indicate what you understand. 3. Repeat with the "wrong" model  the model for the normal signal thn. Can you see differences in
the quality of the fit?
figure
subplot(2,1,1)
estimate2 = compare(spike(120:250),thn,1);
% fitting spike data with a normal model  not quite as good
error2 = estimate2spike(120:250);
% plots the estimate over a longer range (includes the EEG spikes)
subplot(2,1,2)
plot(error2);
% overlay the error signal
max(error2)
For your writeup, submit the plots generated and the maximum errors in the estimation. Annotate the
graphs and show how this estimate is worse than the previous.
4. Repeat for the muscle signal model thm.
5. Compare the frequency responses of the three filters implied by the ths, thn, and the thm models, as
well as the spectrum of the signal from which the ths model was derived. This is easily done with the
course function ar2spectrum.
ar2spectrum(thm);
ar2spectrum(thn);
ar2spectrum(ths);
ftp(spike(170:190)); % this plots the transform of the data from which the
% AR model was derived; this should look like ar2spectrum(ths) Submit the four spectra you have plotted. Note the similarities and differences between the two spike
spectra, and among the three AR spectra. Annotate your plots. 4 ...
View
Full
Document
This note was uploaded on 01/23/2012 for the course EEL 6502 taught by Professor Principe during the Spring '08 term at University of Florida.
 Spring '08
 PRINCIPE
 Frequency, Signal Processing

Click to edit the document details