Calculate the FFT of an ECG signal

Hi all,

currently i'am trying to transform an ecg signal into frequency domain.

I use an ecg signal from MIT-BIH Arrhythmia Database (physionet). To read the ecg file i use the .dat file reader from Said Bourezg (https://de.mathworks.com/matlabcentral/fileexchange/49822-open_ecg-ecg-dat-file-reader).

After reading the signal i add a noisy sine wave to corrupt the signal. Finally I want to calculate the Fourier transform of the noisy signal.

The ECG signal is sampled at 150 Hz and has a duration of 10 seconds (3500 samples/sec).

I would like to know if my frequency plot is correct. It looks a bit weird. It has a relatively high peak at 1 Hz and no other (or small ) magnitudes. And of course you can see the 50 Hz noise. I was actually expecting more frequencies.

clear all
clc
close all

[filename, pathname] = uigetfile('*.dat', 'Open file .dat');
if isequal(filename, 0) || isequal(pathname, 0)   
    disp('File input canceled.');  
   ECG_Data = [];  
else
fid=fopen(filename,'r');
end;

time=10;  % duration of ecg signal
fs = 150; %sample frequency
 
f=fread(fid,2*fs*time,'ubit12');
len = length(f(1:2:length(f))); %length of the signal
time_step = 1/fs; 
max_time = len/fs; % duration of your signal in seconds depending on sampling rate
t = time_step:time_step:max_time;  % this is our time vector.

Orig_Sig=f(1:2:length(f))/len; %compute the clean ecg signal
subplot(4,1,1);
plot(t,Orig_Sig)


noise_coeff = 0.1;    %randomize the amplitude of the noisy ecg signal
noise_signal = 5*sin(2*pi*50*t);      %add 50Hz noise from power supply
noise_signal = noise_signal';           %transpose vector
dirty_signal = Orig_Sig + noise_coeff*noise_signal;  %mix the original ecg signal with the noisy 50 Hz sine

subplot(4,1,2);
plot(t,noise_signal);
subplot(4,1,3);
plot(t,dirty_signal);


NFFT = 2 ^ nextpow2(length(dirty_signal));  %compute FFT length depends on the signal length
Y = fft(dirty_signal,NFFT);  %compute the fft of the noisy signal
Y = Y(1:NFFT/2);  %we only need a one sided fft plot
Y_abs = 1/NFFT*abs(Y); %calculate the magnitude and normalize the spectrum
f_fft = (0:NFFT/2-1)*fs/NFFT; %scale the frequency axe and calculate the coresponding frequencys

subplot(4,1,4);
plot(f_fft,Y_abs);

I also added the used ecg signal (.dat file) and a frequency plot.

If anyone finds a mistake I would be very grateful.

Best regards.

Answer 1

[Daniel M]

I haven't ran your code but the plot looks pretty normal. The peak you see around 0 Hz is a DC component added to your signal. You can subtract the mean of your data using detrend(), and see if that helps. Otherwise, try looking at the plot using semilogy. The resolution would be better if you had more data than just 10 seconds.

[Stephan Lallinger]

Thank you so much.

I feel so stupid i forget complete about the DC.

Now it looks way better.

Just added:

Orig_Sig_DCRemoved = Orig_Sig - mean(Orig_Sig); %remove DC from ecg

Thank you !

[Nur Shahirah]

Can i know in which line you put the this code? i try to figure it out, yet i still have problems with the frequency graph. Thanks in advance.

9/5/2020 : I still can't configure this. Help?

[Arslan Khan]

😐