Issue with FFT/Power Spectral Density

3 views (last 30 days)
Sebastian Daneli
Sebastian Daneli on 4 Nov 2021
Commented: Paul on 8 Nov 2021
Can someone explain what I'm doing wrong here? I wanna evaluate the PSD of ECG data (electrocardiogram), but I end up with this enormous peak in the fourier domain that I don't know what to do with. Where does this come from? What to do about it?
clear
load ecg_stress_test.mat
n=length(ecg_data);
fhat=fft(ecg_data,n);
PSD=fhat.*conj(fhat)/n;
figure
subplot(2,1,1), plot(ecg_data), xlim([0 length(ecg_data)])
subplot(2,1,2), plot(PSD), xlim([0 length(PSD)])

Accepted Answer

Star Strider
Star Strider on 5 Nov 2021
Subtract the mean of the EKG data before calculating the fft. Also, it likely needs to be rescaled, because the amplitude should be between ±1 mV.
The signal needs to be filtered. Some of the noise is broadband, however much can be filtered with a frequency-selective filte —
LD = load('Sebastian Daneli ecg_stress_test.mat')
LD = struct with fields:
ecg_data: [450001×1 double] fs: 250
EKG = LD.ecg_data;
Fs = LD.fs;
L = numel(EKG);
t = linspace(0, 1, L)/Fs;
[EKGfilt,df] = bandpass(EKG, [1 25], Fs, 'ImpulseResponse','iir');
figure
subplot(3,1,1)
plot(t, EKG)
grid
subplot(3,1,2)
idx = 200000:201000;
plot(t(idx), EKG(idx))
grid
subplot(3,1,3)
idx = 200000:201000;
plot(t(idx), EKGfilt(idx))
grid
Fn = Fs/2;
NFFT = 2^nextpow2(L);
FTEKGfilt = fft(EKGfilt,NFFT)/L;
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
[ampmax,idx] = max(abs(FTEKGfilt(Iv))*2);
fprintf('Amplitude = %.2f mV at %.2f Hz',ampmax,Fv(idx))
Amplitude = 1.20 mV at 1.92 Hz
figure
plot(Fv, abs(FTEKGfilt(Iv))*2)
grid
The pwelch function is probably the mose frequently used method to determine the power spectral density. To see the the signal spectrum and power spectrum as a function of time, use the pspectrum or spectrogram functions.
.

More Answers (1)

Paul
Paul on 4 Nov 2021
That peak is at dc because the mean of the data is very large. Try
fhat=fft(detrend(ecg_data),n);
to take out the mean before taking the fft.
  6 Comments
Paul
Paul on 8 Nov 2021
I was confused abou the use of detrend. It should have been
fhat = fhat=fft(detrend(ecg_data,'constant'),n);
to subtract the mean.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!