Issue with FFT/Power Spectral Density
조회 수: 2 (최근 30일)
이전 댓글 표시
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)])
댓글 수: 0
채택된 답변
Star Strider
2021년 11월 5일
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')
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))
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.
.
댓글 수: 0
추가 답변 (1개)
Paul
2021년 11월 4일
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
Paul
2021년 11월 8일
I was confused abou the use of detrend. It should have been
fhat = fhat=fft(detrend(ecg_data,'constant'),n);
to subtract the mean.
참고 항목
카테고리
Help Center 및 File Exchange에서 Parametric Spectral Estimation에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!