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

Learn moreOpportunities for recent engineering grads.

Apply TodayTo resolve issues starting MATLAB on Mac OS X 10.10 (Yosemite) visit: http://www.mathworks.com/matlabcentral/answers/159016

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

*No products are associated with this question.*

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.

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)

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

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.

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

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