How to find average heart beat using fast Fourier transform function

98 views (last 30 days)
Anon
Anon on 9 Jan 2020
Commented: Eva Saleme on 16 May 2022
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)

Community Treasure Hunt

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

Start Hunting!

Translated by