Issue with FFT/Power Spectral Density

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

Answer 1

[Star Strider]

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.

.

Answer 2

[Paul]

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.

[Sebastian Daneli]

That improves it, but i still get a peak of size 2.6*10^7.

[Paul]

What do you see for

sum(detrend(ecg_data))

and what is the index of PSD wher you see the spike?

[Sebastian Daneli]

For that code i get 4.8015e-06, and the index is 1.

[Paul]

Hmm. My understanding of what should happen is:

n=length(ecg_data);
fhat=fft(detrend(ecg_data),n);

At this point we should have

fhat(1) == sum(detrend(ecg_data)) % to within rounding errors

What do you see for fhat(1) and sum(detrend(ecg_data)) ?

PSD=fhat.*conj(fhat)/n;

So you should see

PSD(1) == fhat(1)^2/n ;

Can you run this exact code and show the results:

clear 
load ecg_stress_test.mat
whos ecgdata
n=length(ecg_data);
fhat=fft(detrend(ecg_data),n);
[fhat(1) sum(detrend(ecg_data))]
PSD=fhat.*conj(fhat)/n;
[PSD(1) fhat(1).^2/n]

[Sebastian Daneli]

I will take a look at it.

[Paul]

I was confused abou the use of detrend. It should have been

fhat = fhat=fft(detrend(ecg_data,'constant'),n);

to subtract the mean.