How to find average heart beat using fast Fourier transform function

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

[dpb]

What units is fs?

[Anon]

it is hz

[Star Strider]

Should have posted this as a follow-up to: How to work out BPM from ECG data

Answer 1

[Meg Noah]

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:

[Meg Noah]

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]

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))']';