Matlab code to study the EEG signal

ECE/BIOM 537: Biomedical Signal Processing

Colorado State University                                                           

Student: Minh Anh Nguyen

Email: minhanhnguyen@q.com

1.      This homework will demonstrate EEG signal processing techniques and interpretation. A 10 s signal, with sampling rate of 512 samples per second, has been provided. The signal was monitored and obtained using the C4 and P4 electrodes, and is a differential voltage signal (Image (Links to an external site.) ).

 

The signal in the attached file, 'EEGsig', will be used for parts a through d, while 'EEGsig_wander' will be used in part e.  Use two of the methods discussed in class to estimate the power spectral density of the 10 s epoch, and compare them with some discussion.  

 

a)  Estimate the dominant frequency region in the EEG signal using the PSD estimates. Given your knowledge of the frequency bandwidths in EEG signals and the location of the electrodes, what type of wave might this signal represent? What might it indicate?

The dominant frequency is 8.3Hz, which indicates that the EEG signal is for Theta wave (8-13Hz).

The maximum power occurs at 8.3 Hz

The power estimate is 0.67

b)  Design a low-pass filter to be used on the EEG signal. Choose a cut-off frequency that retains the energy in the dominant frequency region found in part b. Discuss your choice of cut-off frequency and how well your filter worked.

The cut-off frequency is 175Hz; because the dominant frequency is 8-10 Hz and the plot becomes flat after 175Hz.  After 175Hz signal has a lot of noise.

Low-pass filters will pass low frequencies without change, but attenuate (i.e. reduce) frequencies above the cutoff frequency

c)      Plot the low-pass filtered signal using the filter you designed in part c, then comment on the new signal.

e)      The original signal had low-frequency artifacts (baseline wander, etc) that were removed for the signal in parts a-d. Using EEGsig_wander, find a way to remove the baseline wander to approximate EEGsig. Discuss the method used and how well it worked. You can quantify the difference using a metric like mean squared error, if you prefer

. 

 

Matlab code:

ECE/BIOM 537: Biomedical Signal Processing

                                                               Colorado State University                                                           

Student: Minh Anh Nguyen

Email: minhanhnguyen@q.com

 l   close all; clear all; clc;
fs = 512 % fs — Sampling frequency, positive scalar. Sampling frequency, specified as a positive scalar. The sampling frequency is the number of samples per unit time. If the unit of time is seconds, the sampling frequency has units of hertz.
T = 1/fs;% sampling rate or frequency
load('J:\BIOM_Signal_processing\HW10\hmwk_EEGs') % contains eeg1 and fs
N =length(EEGsig); ls = size(EEGsig); % find the length of the data per second
tx =[0:length(EEGsig)-1]/fs;% Make time axis for EEG signal
fx = fs*(0:N/2-1)/N;   %Prepare freq data for plot
figure; subplot (211), plot(tx,EEGsig); xlabel('Time (s)'), ylabel('Amplitude (uV)'), title('Original EEG signal'); %EEG waveform
subplot(212), plot(tx,EEGsig);
xlabel('Time (s)'), ylabel('Amplitude (uV)'), title('Zoom into original EEG signal at 1 to 2 seconds'), xlim([1,2]) % Used to zoom in on single ECG waveformfigure;
%The mean of the PSDs of xl
mean_EEGsig = mean(EEGsig);
max_value=max(EEGsig);
mean_value=mean(EEGsig);
threshold=(max_value-mean_value)/2;

%Estimate the power spectrum of the 10-s epoch by computing the periodogram
%% this method is slide the window through the entire data at every 1/2 second, calculate the frequency, average it.
[p,f] = pwelch(EEGsig,hamming(fs),.5*fs, 2*fs,fs); %%
figure; subplot(421), plot (f,10*log10(p),'r'); xlabel('freq (hz)');
ylabel('PSD Amplitude'); title('Power SPectral Density via Welchs method and hamming window');grid on; xlabel('freq (hz)');ylabel('PSD Amplitude (dB)');
subplot(422), plot (f,10*log10(p),'g'); xlabel('freq (hz)');ylabel('PSD Amplitude (dB)'); title('Power SPectral Density via Welchs method zoom in at 60 hz'); xlim([0,60]);grid on;

x= EEGsig;
[pxx,f] = periodogram(EEGsig,hamming(length(x)),length(x),fs,'power');
[pwrest,idx] = max(pxx);
fprintf('The maximum power occurs at %3.1f Hz\n',f(idx));
fprintf('The power estimate is %2.2f\n',pwrest);
subplot(423), plot(f,10*log10(pxx));title('Power SPectral Density via periodogram method and hamming window');grid on; xlabel('freq (hz)');ylabel('PSD Amplitude (dB)');
subplot(424), plot(f,10*log10(pxx));title('Power SPectral Density via periodogram method and hamming window zoom in at 60Hz');grid on; xlabel('freq (hz)');ylabel('PSD Amplitude (dB)');xlim([0,60]);
%% low pass filter
lpfLength=127; % Order/Number of Filter coefficients
fc = 30; %% cutoff frequency
Wn=(2*fc)/fs; h1=fir1(lpfLength,Wn);
figure; plot(h1);xlabel('Time in Seconds'); ylabel('Magnitude'); title('Low-pass filter');
fi = filtfilt(h1,1,EEGsig);
figure; plot(fi); title ('filtfilt');
%Compute the Fourier transform
Tr1 = conv(EEGsig,h1);
figure; plot(Tr1);

Click on the link below to download the EEG raw data:
  hmwk_EEG.mat file