a problem with power spectrum calculation (FFT)

Hi,

I try to calculate power spectrum of the brain EEG signal. There are 1000 samples, sampling rate is 250. I see a clear 50 Hz in this signal (attached). However, when I do FFT, I cannot see 50Hz but only power in the low frequecnies. What I am doing wrong? I attach my code and mat file.

Thak you a lot for the help!

Alex

Answer 1

[Star Strider]

The mean of ‘data’ is -16657.4418730469. This is the D-C component, so it appears at 0 Hz, and completely prevents the details of the Fourier transform of the signal from being visible. I have no idea what the amplitude units are, however the 50 Hz peak is only 70.59, so , so it is completely hidden.

The solution is to first subtract the mean of ‘data’, then do the fft:

D = load('data_1000samples.mat');
data = D.data;
Fs = 250;
Fn = Fs/2;
L = numel(data);
t = linspace(0, 1, L)/Fs;

q1 = mean(data)
datam = data - mean(data);
FTdata = fft(datam)/L;
Fv = linspace(0, 1, fix(L/2)+1)*Fn;
Iv = 1:numel(Fv);

figure
plot(Fv, abs(FTdata(Iv))*2)
grid

producing:

And now the 50 Hz peak is clearly visible.

If you want to filter it out, that is straightfforward.

.

[Alex Alex]

Cool, thanks!

Can you explain me why this DC masks the effect?

Interestingly, after positng the question I also found that if I take 20*log(powerVec) I also can see the the 50Hz (without subtracting the mean). Do you know what is the rule of thumb when to use log-scale and when not ?

[Star Strider]

My pleasure!

The D-C value of -16657 appears at 0 Hz, so the automatic scaling that plot uses stretches the y-axis limits (ylim) to accommodate it. The rest of the Fourier transform (note the small amplitude in the plot) is then essentially squashed at the upper y-axis limit, so it is not visible with a linear scale. Eliminating the D-C offset allows the rest of the fft output to be visible.

I plotted the amplitude scale here because that is easiest. Plotting the power (and seeing any details) requires the y-scale to be in dB. That version is then:

figure
plot(Fv, 20*log10(abs(FTdata(Iv))*2))
grid
axis([min(Fv)  max(Fv)    -20  +40])
xlabel('Frequency (Hz)')
ylabel('Power (dB)')

The rest of the code is unchanged.

.

[Alex Alex]

Thanks!

[Star Strider]

As always, my pleasure!