Homework _9 Spring 11

Homework _9 Spring 11 - BME 6360 Neural Engineering...

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

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(k-j)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(k-1) + 0.2868 y(k-2) + 0.4012 y(k-3) - 0.0333 y(k-4) - 0.0793 y(k-5) or y(k) = 1.4495 y(k-1) - 0.2868 y(k-2) - 0.4012 y(k-3) + 0.0333 y(k-4) + 0.0793 y(k-5) 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 = estimate-spike(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 = estimate2-spike(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.

Ask a homework question - tutors are online