Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

How to do FFT Analysis to EEG signals Using Matlab

Asked by Ouamar ferhani on 16 Nov 2012

I'm looking for FFT analysis to EEG or ECG signals in MATLAB.

I have tested some codes. Nevertheless I don't know how to obtain a normalized result.

I have read in the documentation that I need to normalize by sample length and sampling rate as well as the window energy, but I don't know if MATLAB functions applied a normalization before...

This is my code for ECG analysis of RR intervals:

{Hwelch=spectrum.welch('Hann',winwidth,50);
psdperiodo = psd(Hwelch,RR,'nfft',nfft/2,'Fs',freqinterpol,'SpectrumType','onesided');
PSD = psdperiodo.data /(size(RR,1)*(freqinterpol * size(psdperiodo.data,1)));
% find the indexes corresponding to the VLF, LF, and HF bands
    iVLF= find((psdperiodo.frequencies>VLFmin) & (psdperiodo.frequencies<=VLFmax));
    iLF = find((psdperiodo.frequencies>LFmin) & (psdperiodo.frequencies<=LFmax));
    iHF = find((psdperiodo.frequencies>HFmin) & (psdperiodo.frequencies<=HFmax));
% calculate raw areas (power under curve), within the freq bands (ms^2)
    aVLF=trapz(psdperiodo.frequencies(iVLF),PSD(iVLF))*10^6;
    aLF=trapz(psdperiodo.frequencies(iLF),PSD(iLF))*10^6;
    aHF=trapz(psdperiodo.frequencies(iHF),PSD(iHF))*10^6;
}

0 Comments

Ouamar ferhani

Tags

Products

No products are associated with this question.

6 Answers

Answer by Wayne King on 17 Nov 2012
Edited by Wayne King on 17 Nov 2012
Accepted answer

I do not see the need for this line:

   psdest = psdperiodo.data /(size(psdperiodo.data,1));

All the normalization for the PSD estimate is taken care of inside of spectrum.welch. That includes dividing by the window energy, which depends on the window used, e.g. Hamming, Bartlett, etc.

I think you should use the one-sided estimate not the two-sided based on the way you are using avgpower(). For example, if you want to do

 aVLF = avgpower(psdest,[VLFmin VLFmax]);

for a two-sided estimate, you really need to do

Fs = 1000;
t = 0:0.001:1-0.001;
x = cos(2*pi*100*t)+randn(size(t));
psd2 = psd(spectrum.periodogram,x,'SpectrumType','twosided','Fs',1000);
psd1 = psd(spectrum.periodogram,x,'Fs',1000);
avgpower(psd2,[98 102])+avgpower(psd2,[898 902])
avgpower(psd1,[98 102])

Or use the 'centerdc' option

    psd2 =  psd(spectrum.periodogram,x,'SpectrumType','twosided','Fs',1000,'centerDC',true);
  avgpower(psd2,[-102 -98])+avgpower(psd2,[98 102])

The way you are using it, you are missing the energy in a real-valued signal which is necessarily "mirrored" for the negative and positive frequencies.

0 Comments

Wayne King
Answer by Wayne King on 16 Nov 2012

You should format your code when you post. Since you are using spectrum.welch, you should use the avgpower() method to calculate the average power in specific frequency intervals.

For example:

    Fs = 1000;
    t = 0:0.001:1-0.001;
    x = cos(2*pi*100*t)+randn(size(t));
    hs = spectrum.welch;
    hs.SegmentLength = 200;
    psdest = psd(hs,x,'Fs',Fs);

Now to determine the average power in an interval around 100 Hz and compare that to the total power over the Nyquist range.

avgpower(psdest,[90 110])/avgpower(psdest)

0 Comments

Wayne King
Answer by Ouamar ferhani on 16 Nov 2012

Thank you for your answer. Sorry for the format... For the normalization of the result is it necessary to divide "psdest" by sample length and sampling rate in matlab ? Thank you for your help

0 Comments

Ouamar ferhani
Answer by Wayne King on 16 Nov 2012

The nonparametric PSD estimates in MATLAB like the periodogram and Welch estimator already "normalize" the result to create the PSD estimate. For example, in the periodogram with the default rectangular window, the magnitude squared DFT values are "normalized" by dividing by the length of the input and the sampling rate as you state in your post, although for a one-sided PSD estimate, the frequencies other than 0 and the Nyquist are multiplied by 2 for energy conversion. You can see this worked out explicitly here:

https://www.mathworks.com/help/signal/ug/psd-estimate-using-fft.html

In the case of the Welch estimate, it's more complicated since you have to divide also by the energy of the window.

0 Comments

Wayne King
Answer by Ouamar ferhani on 17 Nov 2012

Thank you for your answer. I have also read the post. But I have two more questions : Should I calculate psd in oneesided or two sided option in my case ? How could I know it ?

And the second one, if I understand correctly your explanations my code should be as following :

psdperiodo = psd(Hwelch,RR,'nfft',nfft/2,'Fs',freqinterpol,'SpectrumType','twosided'); 
% selected in twosided option
psdest = psdperiodo.data /(size(psdperiodo.data,1)); 
% divided by the window energy
% Values of interest in pourcent
aVLF = avgpower(psdest,[VLFmin VLFmax])/avgpower(psdest)*100;
aLF = avgpower(psdest,[LFmin LFmax])/avgpower(psdest)*100;
aHF = avgpower(psdest,[HFmin HFmax])/avgpower(psdest)*100;
% Values of interest in absolute values
aVLF = avgpower(psdest,[VLFmin VLFmax]);
aLF = avgpower(psdest,[LFmin LFmax]);
aHF = avgpower(psdest,[HFmin HFmax]);

Are you agree with my code ? Thanks a lot

0 Comments

Ouamar ferhani
Answer by Ouamar ferhani on 18 Nov 2012
Edited by Ouamar ferhani on 18 Nov 2012

Thank you for your answer and explanations!

0 Comments

Ouamar ferhani

Contact us