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
采纳的回答
Star Strider
2020-8-8
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.
.
4 个评论
Alex Alex
2020-8-8
Star Strider
2020-8-8
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
2020-8-8
Star Strider
2020-8-8
As always, my pleasure!
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Spectral Measurements 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
