How to find average heart beat using fast Fourier transform function

309 views (last 30 days)
Hi,
I am trying to use the fft function to compute the power spectrum of an ECG.
ecg=load ('ecg.csv');
A=fft(ecg);
L=length(ecg);
X_mag=abs(A);
X_phase=angle(A);
fs=350/2;
Fbins=((0: 1/L: 1-1/L)*fs).';
E=(X_mag).^2;
figure
plot (Fbins, E)
So far I have done this to create a plot using the data. How do I change the frequency axis to BPM?
Also where and how do I find the average heart rate from this data?
Any help would be much appreciated as I am new to matlab and struggling.
Maybe I should add that to obtain the power spectrum, i need to compute the following equation:
? = |FFT(???)|^2

Accepted Answer

Meg Noah
Meg Noah on 9 Jan 2020
I'm going to give a solution that says your 3000 data points represent a timeseries of 6 seconds. The code can be modified if that assumption is incorrect. Just have Fs=3000/elapsed-time-in-seconds
X = dlmread('ECG.csv');
% assumed that the 3000 points are 6 seconds of data...?
Fs = 3000/6; % [samples/s] sampling frequency
T = 1/Fs; % [s] sampling period
N = 3000; % [samples] Length of signal
t = (0:N-1)*T; % [s] Time vector
deltaF = Fs/N; % [1/s]) frequency intervalue of discrete signal
figure('color','white','position',[70 100 600 900]);
subplot(3,1,1);
plot(1e3*t,X)
title({'Heartbeat data'})
xlabel('t (milliseconds)')
ylabel('X(t)')
% compute the fast fourier transform
Y = fft(X);
% manually shifting the FFT
Y = abs(Y/N);
Amp = [Y(ceil(end/2)+1:end)' Y(1) Y(2:ceil(end/2))']';
if (mod(N,2) == 0)
sampleIndex = -N/2:1:N/2-1; %raw index for FFT plot
else
sampleIndex = -(N-1)/2:1:(N-1)/2; %raw index for FFT plot
end
subplot(3,1,2);
plot(deltaF*sampleIndex, Amp);
hold on;
idx = find(Amp > 3.5);
idx(sampleIndex(idx) < 0) = [];
plot(deltaF*sampleIndex(idx), Amp(idx), '+');
for k = 1:length(idx)
if (idx(k) > (N-1)/2 && Amp(idx(k))>3.5)
h = text(deltaF*sampleIndex(idx(k)), Amp(idx(k))+0.15,...
['f=' num2str(deltaF*sampleIndex(idx(k))) ' Hz']);
set(h,'rotation',90)
end
end
xlabel('Frequency [Hz]');
ylabel('Amplitude');
title(['Heartbeat = ' num2str(deltaF*sampleIndex(idx(1))) ' Hz = ' ...
num2str(60.0*(deltaF*sampleIndex(idx(1)))) ' BPM']);
xlim([0 max(deltaF*sampleIndex)/4]);
subplot(3,1,3);
half_f = deltaF*(0:(N/2));
plot(fftshift([half_f -fliplr(half_f(2:end+mod(N,2)-1))]), ...
abs(fftshift(Y)/N));
xlabel('Frequency [Hz]');
ylabel('Amplitude');
title('Using fftshift');
This yields:
heartbeat.png
  2 Comments
Eva Saleme
Eva Saleme on 16 May 2022
Hiii, Iam writing the same code but their was an error in Amp= [Yfftm(ceil(end/2)+1:end)' Yfftm(1) Yfftm(2:ceil(end/2))']'; , how can it be fixed?
this is the complete error :
Error using horzcat
Dimensions of matrices being concatenated are not consistent.
Error in HeartRate (line 72)
Amp = [Yfftm(ceil(end/2)+1:end)' Yfftm(1) Yfftm(2:ceil(end/2))']';

Sign in to comment.

More Answers (0)

Categories

Find more on Fourier Analysis and Filtering in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by