How to find average heart beat using fast Fourier transform function

71 visualizzazioni (ultimi 30 giorni)
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
  3 Commenti

Accedi per commentare.

Risposta accettata

Meg Noah
Meg Noah il 9 Gen 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 Commenti
Meg Noah
Meg Noah il 9 Gen 2020
Modificato: Meg Noah il 9 Gen 2020
Note: there are strong components in the higher frequencies indicating that the heart might be in some sort of irregular beat.
To make the frequency from Hz to BPM just multiply the frequency axis by 60 since there are 60 seconds in a minute.
The Amplitude should be squared for your problem. Just plot Amp.^2 instead of Amp.
Eva Saleme
Eva Saleme il 16 Mag 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))']';

Accedi per commentare.

Più risposte (0)

Categorie

Scopri di più su ECG / EKG in Help Center e File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by