How can I get the power of a specific frequency band after FFT?

I need help with my matlab code. I have recorded EEG data and now I want to analyse it in matlab. I want to find the mean power of specific frequency ranges (theta: 4-8 Hz,alpha: 8-12Hz etc.). For that I performed a fft. Next I used the interp1 function to find the amplitudes at specific frequency ranges. But I am not sure if I have done it correctly...It would be great if someone can look at the code.

%p2 is the data in the time domain
% Variables

fftlength = 1024;                                % Set FFT length
stepsize = 32;                                   % Set step size
samples = fftlength:stepsize:length(p2);    % create an index array
c = 1:fftlength;

%Specify hanning window

hanning = [1:fftlength]';
hanning_in = 2* pi() * (hanning - (fftlength+1)/2)/(fftlength+1);

%rescaled x-axis to match sample length?

hanning1 = (sin(hanning_in)./hanning_in).^2;

f=[128/fftlength:128/fftlength:128];   % frequency index for the spectral
array

bp_p2=cell(1,1);

channel=1;

for kk = 1:length(samples)

k = samples(kk);

spectrum = fft(p2(k-fftlength+1:k,:) .* hanning1)/fftlength; % apply window
to filtered data

spectrum2 =(sqrt(spectrum .* conj(spectrum)));  % get magnitude uV

% theta (4-8),alpha (8-12),beta low(12-16), beta high(16-25)

theta=mean(interp1(f,spectrum2(:,channel),[4:f:8]));
alpha=mean(interp1(f,spectrum2(:,channel),[8:f:12]));
betalow=mean(interp1(f,spectrum2(:,channel),[12:f:16]));
betahigh=mean(interp1(f,spectrum2(:,channel),[16:f:25]));

%creating an array
%store all iteration loop outputs in a matrix

m(kk,:)=[theta alpha betalow betahigh];

end


To get the total power across bins, sum the power in each bin. Also you need to compensate for your window loss if you want an accurate result.

For a rectangular window, the power in each DFT bin is given as:

$$P = \frac{Y[k]Y[k]^*}{N^2} = \bigg(\frac{|Y[k]|}{N}\bigg)^2$$

and $$P_{dB} = 20Log_{10}\bigg(\frac{|Y[k]|}{N}\bigg)$$

The total power over a range of such bins would be:

$$P_{k} = \frac{\Sigma_k(|Y[k]|^2)}{N^2}$$

Windowing Loss

When we apply a window, the coherent gain and non-coherent gain are affected differently. If the N samples are coherent (such as a frequency completely contained in a bin), the output grows at rate N (direct summation). If the N samples are non-coherent (spread across bins), the output grows at the rate of $$\sqrt{N}$$ (rms summation).

If we are only concerned with the signal power in each bin, we would use the coherent gain which is given as

$$G_c = \frac{\Sigma w[n]}{N}$$

where $$w[n]$$ are the window coefficients. (For example the Blackman(30) window has a coherent gain of 0.406 or -7.83 dB):

Likewise non-coherent samples grow at $$\sqrt{N}$$. The non-coherent gain refers to the change in DFT output level of the non-coherent samples (such as white noise) due to the window function.

The normalized non-coherent gain is the rms value of the window weights given as:

$$G_{nc} = \sqrt{\frac{\Sigma w[n]^2}{N}}$$

The SNR change or processing gain of the window would be the ratio of the two, in dB this would be:

$$PG = 20log_{10}\bigg( \frac{G_c}{G_{nc}} \bigg)$$

If you are only interested in the signal power level and not signal to noise ratio, so you would simply compensate for the power computed by using $$G_c$$.

It may confusing to some how the power for a single tone can be effected differently by the window than the power in a signal that is spread evenly across multiple tones as described above. An intuitive explanation of this is that windowing the time domain function results in each bin in the FFT having an equivalent noise bandwidth (ENBW) that exceeds one bin. A rectangular window in comparison has an ENBW of exactly one bin, so when we sum the power of all the bins in the frequency domain, Parseval's theorem proves correct in that the total power from that summation is equal to the total power in the time domain signal. However if we window in time with something other than a rectangular window, the effect of the wider ENBW means that each bin in frequency is now reporting the power for that bin PLUS adjacent bins. So in this case when we sum the reported power over all bins, the total summation will over-estimate the actual noise power. In a case of a single tone, or tones that are spaced wider than the ENBW, each bin will accurately report the power for each tone when we have scaled the result by the recripocal of $$\sum W_n$$ where $$W$$ is the windowing funciton. (for the rectangular window this scaling results to $$1/N$$).