How to do FFT Analysis to EEG signals Using Matlab

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;
}

Answer 1

[Wayne King]

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.

Answer 2

[Wayne King]

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)

Answer 3

[Ouamar ferhani]

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

Answer 4

[Wayne King]

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.

Answer 5

[Ouamar ferhani]

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

Answer 6

[Ouamar ferhani]

Thank you for your answer and explanations!