How to do FFT Analysis to EEG signals Using Matlab
11 views (last 30 days)
Show older comments
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
Accepted Answer
Wayne King
on 17 Nov 2012
Edited: Wayne King
on 17 Nov 2012
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
More Answers (5)
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
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:
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
See Also
Categories
Find more on Bartlett in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!