MATLAB Answers

## How to do FFT Analysis to EEG signals Using Matlab

Asked by Ouamar ferhani

### Ouamar ferhani (view profile)

on 16 Nov 2012
Accepted Answer by Wayne King

### Wayne King (view profile)

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

## Products

No products are associated with this question.

## 6 Answers

### Wayne King (view profile)

Answer by Wayne King

### Wayne King (view profile)

on 17 Nov 2012
Edited by Wayne King

### Wayne King (view profile)

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.

### Wayne King (view profile)

Answer by Wayne King

### Wayne King (view profile)

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)
```

### Ouamar ferhani (view profile)

Answer by Ouamar ferhani

### Ouamar ferhani (view profile)

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

### Wayne King (view profile)

Answer by Wayne King

### Wayne King (view profile)

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.

### Ouamar ferhani (view profile)

Answer by Ouamar ferhani

### Ouamar ferhani (view profile)

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

### Ouamar ferhani (view profile)

Answer by Ouamar ferhani

### Ouamar ferhani (view profile)

on 18 Nov 2012
Edited by Ouamar ferhani

### Ouamar ferhani (view profile)

on 18 Nov 2012

Thank you for your answer and explanations!

#### Join the 15-year community celebration.

Play games and win prizes!

Learn more

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

### MATLAB Academy

New to MATLAB?

Learn MATLAB today!