MATLAB: How to do FFT Analysis to EEG signals Using Matlab

fft ecg

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

Best 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.
Related Question