Matlab code to estimate the power spectrum of the EEG signal

1.       EEG signals have been generated encompassing a 10 second epoch and numbered. Select and download the file corresponding to the number assigned to you. The sampling rate is 3000 Hz.

Estimate the power spectrum of the 10-s epoch by computing the periodogram. Compute several periodograms and compare the results. Change the window size for each periodogram, using 3 or 4 values of your own choosing. An example might be: 100 sample windows, 200 samples, 1000 samples, and 1000 samples.

2. Given three signals: x1 = cos(2*pi*4*t); x2 = cos(2*pi*20*t); x3 = cos(2*pi*40*t).Generate .3 seconds of each signal. First, combine them into new signals x4 = x1 + x2 + x3

 

x5 = [x1 x2 x3]  where x5 is concatenated from the first 3 vectors, and x4 is the sum of the first three signals.

 

 

Choose an appropriate sampling rate ( >> nyquist), and plot the spectrum of each signal. Compare the results.

 

Matlab code:

ECE/BIOM 537: Biomedical Signal Processing

Colorado State University                                                           

Student: Minh Anh Nguyen

Email: minhanhnguyen@q.com

close all; clear all; clc;

fs = 3000 % 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\Hw9\EEGsignal_1') % contains eeg1 and fs

N =length(eeg1); ls = size(eeg1); % find the length of the data per second

tx =[0:length(eeg1)-1]/fs;% Make time axis for EEG signal

figure; subplot (211), plot(tx,eeg1); xlabel('Time (s)'), ylabel('Amplitude (mV)'), title('Original EEG signal'); %EEG waveform

subplot(212), plot(tx,eeg1);

%fignum = fignum + 1; figure(fignum);

xlabel('Time (s)'), ylabel('Amplitude (mV)'), title('Zoom into original EEG signal at 1 to 2 seconds'), xlim([1,2]) % Used to zoom in on single ECG waveform

figure;

NFFT = 2 ^ nextpow2(N);  %% number points. Number of DFT points, specified as a positive integer. For a real-valued input signal, x, the PSD estimate, pxx has length (nfft/2 + 1) if nfft is even, and (nfft + 1)/2 if nfft is odd. For a complex-valued input signal,x, the PSD estimate always has length nfft. If nfft is specified as empty, the default nfft is used.

Y = fft(eeg1, NFFT)/N; %% fft of the EEG signals

f = (fs/2 * linspace(0,1,NFFT/2+1))'; % Vector containing frequencies in Hz

amp =( 2 * abs(Y(1: NFFT/2+1))); % Vector containing corresponding amplitudes

subplot(2,1,1), plot (f, amp);title ('plot single-sided amplitude spectrume of the EEG signal'); xlabel ('frequency (Hz)'); ylabel ('|y(f)|');grid on;

%Estimate the power spectrum of the 10-s epoch by computing the periodogram

%the plot using periodogram with no outputs.

periodogram_1 = periodogram(eeg1);

subplot(2,1,2), plot (f,periodogram_1);title('Periodogram power spectral Densisty Estimate of the EEG signal'); xlabel ('frequency (Hz)');ylabel ('Power/Frequency(dB)');grid on;

%%% plot the periodogram with the length of EEG signal

nfft= length(eeg1);

periodogram_2 =periodogram(eeg1,[],nfft);

f1 = (fs/2 * linspace(0,1,nfft/2+1))'; % Vector containing frequencies in Hz

figure;

subplot(2,2,1), plot(f1,periodogram_2);title('Periodogram power spectral Densisty Estimate of the EEG signal and length of EGG signal'); xlabel('frequency (Hz)');ylabel('Power/Frequency(dB)');grid on;

%%%%The signal is 30001 samples in length. Obtain the periodogram using the default rectangular window and DFT length. The DFT length is the next power of two greater than the signal length, or 32768 points.

%Because the signal is real-valued and has even length, the periodogram is one-sided and there are 512/2+1 points.

[pxx,w] = periodogram(eeg1);

subplot(2,2,2), plot(w,10*log10(pxx)); title ('Periodogram using the default rectangular window and DFT length'); xlabel ('frequency(Hz)'); ylabel ('Power/Frequency(dB)');grid on;

%%Modified Periodogram with Hamming Window. Obtain the modified periodogram of an input EEG signal with no noise.  The signal is 30001 samples in length. Obtain the modified periodogram using a Hamming window and default DFT length. The DFT length is the next power of two greater than the signal length, or 32786 points.

%Because the signal is real-valued and has even length, the periodogram is one-sided and there are 32768/2+1 points.

hamming_1=periodogram(eeg1,hamming(length(eeg1)));

subplot(2,2,3), plot (f,hamming_1);title ('Periodogram using the hamming window and DFT length'); xlabel ('frequency(Hz)'); ylabel ('PSD estimate/Power/Frequency(dB) ');grid on;

%DC-Centered Periodogram

%Obtain the periodogram of EEG signal. The data are sampled at 3000 kHz. Use the 'centered' option to obtain the DC-centered periodogram and plot the result.

[periodogram_center, fcenter]= periodogram(eeg1,[],length(eeg1),fs,'centered');

subplot(2,2,4), plot(fcenter,periodogram_center); title ('Periodogram using DC-Centered and DFT length'); xlabel ('frequency(Hz)'); ylabel ('Power/Frequency(dB)');grid on;

%Change the window size for each periodogram, using 3 or 4 values of your

%own choosing. An example might be: 100 sample windows, 200 samples, 1000 samples, and 2000 samples.

%use example from this link: http://www.mathworks.com/help/signal/ref/periodogram.html

%pxx = periodogram(x,window) returns the modified periodogram PSD estimate using the window, window. window is a vector the same length as x.

%% use from this link: http://www.mathworks.com/matlabcentral/answers/73055-hai-i-am-new-to-matlab-how-to-divide-a-eeg-signal-into-frames-in-matlab

%eeg_100= eeg1(1:100,1); %% create the data sample from 0 to 100 sample from column 1.

eeg_100= eeg1(1:101); %% create the data sample from 0 to 100 sample from column 1.

nfft_100 = length(eeg_100);

[pxx_100,w1] = periodogram(eeg_100);

figure;

subplot(3,2,1), plot(w1,10*log10(pxx_100)); title ('Periodogram using the default rectangular window for 100 sample and DFT length'); xlabel ('frequency(Hz)'); ylabel ('Power/Frequency(dB)');grid on

eeg_200= eeg1(1:201); %% create the data sample from 0 to 200 sample from column 1.

nfft_200 = length(eeg_200);

[pxx_200,w2] = periodogram(eeg_200);

subplot(3,2,2), plot(w2,10*log10(pxx_200)); title ('Periodogram using the default rectangular window for 200 sample and DFT length'); xlabel ('frequency(Hz)'); ylabel ('Power/Frequency(dB)');grid on;

eeg_1000= eeg1(1:1001); %% create the data sample from 0 to 1000 sample from column 1.

nfft_1000 = length(eeg_1000);

[pxx_1000,w3] = periodogram(eeg_1000);

subplot(3,2,3), plot(w3,10*log10(pxx_1000)); title ('Periodogram using the default rectangular window for 1000 sample and DFT length'); xlabel ('frequency(Hz)'); ylabel ('Power/Frequency(dB)');grid on;

eeg_2000= eeg1(1:2001); %% create the data sample from 0 to 1000 sample from column 1.

nfft_2000 = length(eeg_2000);

[pxx_2000,w4] = periodogram(eeg_2000);

subplot(3,2,4), plot(w4,10*log10(pxx_2000)); title ('Periodogram using the default rectangular window for 2000 sample and DFT length'); xlabel ('frequency(Hz)'); ylabel ('Power/Frequency(dB)');grid on;

%% this method is slide the window through the entire data at every 1/2 second, calculate the frequency, average it.

[p,f] = pwelch(eeg1,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');grid on;

subplot (422), plot (f,10*log10(p),'g'); xlabel('freq (hz)');ylabel('PSD Amplitude'); title('Power SPectral Density via Welchs method zoom in at 40 hz'); xlim([0,40]);grid on;

%text(0.5,-1/3,'{\itNote the odd symmetry.}')

%%%%%%%

%%% problem 2

%2. Given three signals:

%x1 = cos(2*pi*4*t)

%x2 = cos(2*pi*20*t)

%x3 = cos(2*pi*40*t)

%Generate .3 seconds of each signal. First, combine them into new signals

%x4 = x1 + x2 + x3     x4 is the sum of the first three signals.

%x5 = [x1 x2 x3] where x5 is concatenated from the first 3 vectors

clc;clear all;close all;

t = 0:0.01:1;

x1 = cos(2*pi*4*t);  % for x1: frequency is 4hz,

x2 = cos(2*pi*20*t);%x2: frequency = 20Hz,

x3 = cos(2*pi*40*t);% is 40 Hz.

x4 = x1 + x2 + x3; 

f1 = 4;

f2 = 20;

f3= 40;

T1 = 1/f1;

T2 = 1/f2;

T3 = 1/f3;

figure;

subplot(2,2,1), plot(t,x1); xlabel('Time (s)'), ylabel('Amplitude'), title('signal: x1 = cos(2*pi*4*t) ');

subplot(2,2,2), plot(t,x2); xlabel('Time (s)'), ylabel('Amplitude'), title('signal: x2 = cos(2*pi*20*t) ');

subplot(2,2,3), plot(t,x3); xlabel('Time (s)'), ylabel('Amplitude'), title('signal: x3 = cos(2*pi*40*t) ');

subplot(2,2,4), plot(t,x4); xlabel('Time (s)'), ylabel('Amplitude'), title('signal: x4 = x1 + x2 + x3 ');

% Concatinate the three signals:

   x5 = [x1 x2 x3] ;

   figure;

   plot(x5); xlabel('Time (s)'), ylabel('Amplitude'), title('Concatinate the three signals: x1, x2, and x3 ');

%What is the Nyquist rate for this signal?

   fN1=f1/2;

   fN2 = f2/2;

   FN3= f3/2; %Fmax

 fs = 50; %The Nyquist rate is FN= 2Fmax

 %plot the spectrum of each signal. Compare the results.

NFFT_x1 = 2 ^ nextpow2(length(x1));  %% number points. Number of DFT points, specified as a positive integer. For a real-valued input signal, x, the PSD estimate, pxx has length (nfft/2 + 1) if nfft is even, and (nfft + 1)/2 if nfft is odd. For a complex-valued input signal,x, the PSD estimate always has length nfft. If nfft is specified as empty, the default nfft is used.

Y_x1 = fft(x1, NFFT_x1)/length(x1); %% fft of the EEG signals

f_x1 = (fs/2 * linspace(0,1,NFFT_x1/2+1))'; % Vector containing frequencies in Hz

amp_x1 =( 2 * abs(Y_x1(1: NFFT_x1/2+1))); % Vector containing corresponding amplitudes

NFFT_x2 = 2 ^ nextpow2(length(x2));  %% number points. Number of DFT points, specified as a positive integer. For a real-valued input signal, x, the PSD estimate, pxx has length (nfft/2 + 1) if nfft is even, and (nfft + 1)/2 if nfft is odd. For a complex-valued input signal,x, the PSD estimate always has length nfft. If nfft is specified as empty, the default nfft is used.

Y_x2 = fft(x2, NFFT_x1)/length(x2); %% fft of the EEG signals

f_x2 = (fs/2 * linspace(0,1,NFFT_x2/2+1))'; % Vector containing frequencies in Hz

amp_x2 =( 2 * abs(Y_x2(1: NFFT_x2/2+1))); % Vector containing corresponding amplitudes  

  

NFFT_x3 = 2 ^ nextpow2(length(x3));  %% number points. Number of DFT points, specified as a positive integer. For a real-valued input signal, x, the PSD estimate, pxx has length (nfft/2 + 1) if nfft is even, and (nfft + 1)/2 if nfft is odd. For a complex-valued input signal,x, the PSD estimate always has length nfft. If nfft is specified as empty, the default nfft is used.

Y_x3 = fft(x3, NFFT_x3)/length(x3); %% fft of the EEG signals

f_x3 = (fs/2 * linspace(0,1,NFFT_x3/2+1))'; % Vector containing frequencies in Hz

amp_x3 =( 2 * abs(Y_x3(1: NFFT_x3/2+1))); % Vector containing corresponding amplitudes

  

NFFT_x4 = 2 ^ nextpow2(length(x4));  %% number points. Number of DFT points, specified as a positive integer. For a real-valued input signal, x, the PSD estimate, pxx has length (nfft/2 + 1) if nfft is even, and (nfft + 1)/2 if nfft is odd. For a complex-valued input signal,x, the PSD estimate always has length nfft. If nfft is specified as empty, the default nfft is used.

Y_x4 = fft(x4, NFFT_x4)/length(x4); %% fft of the EEG signals

f_x4 = (fs/2 * linspace(0,1,NFFT_x4/2+1))'; % Vector containing frequencies in Hz

amp_x4 =( 2 * abs(Y_x3(1: NFFT_x4/2+1))); % Vector containing corresponding amplitudes 

NFFT_x5 = 2 ^ nextpow2(length(x5));  %% number points. Number of DFT points, specified as a positive integer. For a real-valued input signal, x, the PSD estimate, pxx has length (nfft/2 + 1) if nfft is even, and (nfft + 1)/2 if nfft is odd. For a complex-valued input signal,x, the PSD estimate always has length nfft. If nfft is specified as empty, the default nfft is used.

Y_x5 = fft(x5, NFFT_x5)/length(x5); %% fft of the EEG signals

f_x5 = (fs/2 * linspace(0,1,NFFT_x5/2+1))'; % Vector containing frequencies in Hz

amp_x5 =( 2 * abs(Y_x5(1: NFFT_x5/2+1))); % Vector containing corresponding amplitudes

   figure;

subplot(3,2,1), plot (f_x1, amp_x1); title ('plot single-sided amplitude spectrume of signal x1 = cos(2*pi*4*t)'); xlabel ('frequency (Hz)'); ylabel ('|y(f)|');grid on;

subplot(3,2,2), plot (f_x2, amp_x2); title ('plot single-sided amplitude spectrume of signal x2 = cos(2*pi*20*t)'); xlabel ('frequency (Hz)'); ylabel ('|y(f)|');grid on;

subplot(3,2,3), plot (f_x3, amp_x3); title ('plot single-sided amplitude spectrume of signal x3 = cos(2*pi*40*t)'); xlabel ('frequency (Hz)'); ylabel ('|y(f)|');grid on;

subplot(2,2,4),plot (f_x4, amp_x4); title ('plot single-sided amplitude spectrume of signal x4 = x1+x2+x3 '); xlabel ('frequency (Hz)'); ylabel ('|y(f)|');grid on;

figure;

subplot (211); plot(x5); xlabel('Time (s)'), ylabel('Amplitude'), title('Concatinate the three signals: x1, x2, and x3 ');

subplot (212); plot (f_x5, amp_x5); title ('plot single-sided amplitude spectrume of concatinate the three signals: x1, x2, and x3 '); xlabel ('frequency (Hz)'); ylabel ('|y(f)|');grid on;