Convert sound of heart beat into beats per minute

Hello,

I have never before used matlab. I am a biology graduate student and I have a bunch of audio files(.mp3) of bird heart rate recorded on a standard audio recorder. What I would like to do is convert the audio into beats per minute, ideally for each minute in the ~48 hours of audio for each bird such that I can look at averages over select periods of time. Is this possible in matlab? If more information is needed to answer this question let me know.

[Star Strider]

I have no idea what an avian phonocardiogram contains. Most phonocardiogram records are of valve closure sounds, S1 being the first heart sound (mitral and tricuspic valve closure) and S2 being the second heart sound (pulmonary and aortic valve closure), both relatively high frequency and usually close to simultaneous (although S2 can be closely split). Healthy human hearts have no other sounds (except a systolic flow murmur, a normal variant, typically in athletes). S1 and S2 should be relatively close together with respect to the rest of the cardiac cycle, so I would begin there.

With the advent of cardiac ultrasonography, phonocardiography has largely been left by the wayside in clinical human medicine. I was unable to find any recent papers discussing automated phonocardiographic analysis. Search PubMed (link) to see what you can find.

I am not listing this as an Answer because it isn’t one.

[Naman Chaturvedi]

Can you share a sample of the audio file?

[Elyssa Watford]

Hello, thank you for your responses. It has taken me a while to find time to get back to working on this problem. I have a chunk of an audio file but I don't know how to share it as it is .mp3. If you can tell me how to share it I will supply some sample data. Thanks

[Star Strider]

You may have to load it, extract it, and then re-format it into one of the formats listed in the Algorithms (link) section of the audiowrite (link) funciton.

However, the audioread (link) function can read mp3 file on Windows 7 and later Windows machines, according to the filename (link) section of the documentation, so you may be able to upload it and attach it as a file here.

Note that with respect to Answers Attachments:

• Files must be 5 MB or smaller.
• Attachments must be one of the following formats: bmp, csv, fig, gif, jpg, jpeg, m, mat, mdl, mlx, slx, pdf, png, txt, xls, xlsx, zip.

Loading it, decoding it to a matrix, putting the matrix and sampling frequency (‘Fs’) in a .mat file and attaching the .mat file to your original Question or to a Comment would likely be the easiest option.

[Walter Roberson]

If it is more than 5 Mb after zip then you can post a link to Google Drive or equivalent .

[Elyssa Watford]

Thank you for your responses. I appologize I wasn't able to respond sooner. Here is an audio clip of the bird heatbeat if anyone is interested in listening to it or has suggestions for creating an automated counter to convert to beats per minute.

Thanks

[Iqra Chaudhary]

Do you have code?

[Joseph Tricarico]

Hi Elyssa,

I realise that this post is a few years old, but if you do this work in the future, this may be useful to you:

https://adelaide.figshare.com/articles/dataset/Matlab_code_for_determining_instantaneous_heart_rate_from_audio_signals_of_nesting_birds/6993485

Joe T.

Answer 1

[Star Strider]

I worked on that for a while, including attempting wavelet denoising, without success. The heartbeats are easily audible, however detecting them reliably using every signal processing techhnique I’m aware of is essentially impossible. There is too much noise even for wavelet decomposition, and it’s very difficult to isolate individual heartbeats.

Nevertheless, here’s my code (that works even though it doesn’t result in usable heartbeat information) so you can experiment with it:

filename = 'sample_HR_audio.mp3';
[HRA0, Fs0] = audioread(filename);

% HRA = resample(HRA0, 1, 10);
HRA = HRA0(1:10:end);                                           % Decimate
Fs = Fs0/10;

L = numel(HRA);
Ts = 1/Fs;                                                      % Sampling Interval
Fn = Fs/2;                                                      % Nyquist Frequency
t = linspace(0, L, L)/Fs;                                       % Time Vector

figure
plot(t, HRA)
grid
xlabel('Time (s)')
ylabel('Amplitude')
title('Recorded Avian Phonocardiogram')

% soundsc(HRA(1:fix(L/5)), Fs)

FT_HRA = fft(HRA)/L;                                            % Fourier Transform
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                             % Frequency Vector
Iv = 1:numel(Fv);                                               % Index Vector

[Pks,Locs] = findpeaks(abs(FT_HRA(Iv))*2, 'MinPeakHeight',2E-4, 'MinPeakDistance',2500);
Freqs = Fv(Locs);

figure
plot(Fv, abs(FT_HRA(Iv))*2)
hold on
plot(Fv(Locs), Pks, '^r')
hold off
grid
xlim([0, 500])
xlabel('Frequency (Hz)')
ylabel('Amplitude')
title('Fourier Transform of Avian Phonocardiogram')
text(200, min(Pks), sprintf('%.1f Hz    %7.1E\n', [Freqs(:) Pks]'), 'VerticalAlignment','bottom')

Wp = [20 40]/Fn;                                                % Design Filter & Filter Signal
Ws = [15 45]/Fn;
Rp =  1;
Rs = 75;
[n,Wp] = ellipord(Wp, Ws, Rp, Rs);
[z,p,k] = ellip(n,Rp,Rs,Wp);
[sos,g] = zp2sos(z,p,k);
HRA_Filtered = filtfilt(sos, g, HRA);
[envh,envl] = envelope(HRA_Filtered, 5000, 'peak');             % Calculate Envelope

figure
plot(t, HRA_Filtered)
hold on
plot(t, envh,    t, envl)
hold off
grid
xlabel('Time (s)')
ylabel('Amplitude')
title('Bandpass-Filtered Avian Phonocardiogram')

env_mean = mean([envh,-envl],2);

[PksB,LocB] = findpeaks(env_mean, 'MinPeakDistance',500);
HR_mean = numel(PksB)/max(t)*60;

figure
plot(t, env_mean)
hold on
plot(t(LocB), PksB, '^r')
hold off
grid
xlabel('Time (s)')
ylabel('Amplitude')
title('Peaks and Mean Envelope of Bandpass-Filtered Avian Phonocardiogram')
text(mean(t), mean(PksB)/4, sprintf('Mean Heart Rate = %.1f / min', HR_mean), 'HorizontalAlignment','center')

The decimation step is intended to speed up the computations, since I had significant problems with the length of the signal even on my Ryzen 7 1700 machine with 64GB memory. It is not absolutely necessary.

[Elyssa Watford]

Thank you Star Strider for that code. I'll see if I can play around with it some. I would say this clip is certainly one of the less noisy ones so if there is trouble with it, the others will be more troublesome. I appreciate you working on this.

[Star Strider]

My pleasure.

[Star Strider]

I decided to take another look at that.

Try this —

Uz = unzip('sample_HR_audio.zip');
[HRA0, Fs0] = audioread(Uz{1});

L = numel(HRA0);
t = linspace(0, L-1, L).'/Fs0;                          % Create Time Vector

figure
plot(t, HRA0)
grid
xlabel('Time (s)')
ylabel('Amplitude')

[p,f,t] = pspectrum(HRA0,Fs0,'spectrogram');

figure
waterfall(f,t,p')
xlabel('Frequency (Hz)')
ylabel('Time (seconds)')
wtf = gca;
wtf.XDir = 'reverse';
colormap(turbo)
view([30 45])
xlim([0 2.5E+2])

figure
plot(t, p(1,:))
grid
xlabel('Time (s)')
ylabel('Magnitude')

p1_filt = sgolayfilt(p(1,:), 3, 51);                        % Multiband Empirical FIR Filter

[pks,locs] = findpeaks(p1_filt, 'MinPeakProminence',1E-3);

Beats = t(locs)                                             % Beat Times (s)

Beats = 4×1

13.6460 100.9650 149.1606 222.8715

IBI = diff(Beats)                                           % Inter-beat Interval

IBI = 3×1

87.3190 48.1956 73.7109

HR = 60/mean(IBI)                                           % Beats/Minute

HR = 0.8603

figure
plot(t, p1_filt)
hold on
plot(t(locs), pks, 'r^')
hold off
grid
xline((60:60:240), '--k', compose('%2d minutes',1:4))
xlabel('Time (s)')
ylabel('Magnitude')

.

[Joseph Tricarico]

Thank you for circling back to this after five years! This was a great step forward. As your work shows, there are some pretty loud spikes in the clip that make it tough to identify heart beats across the clip. I worked with one of the more steady sections (30–60s), and was able to get some of the false positives knocked out by filtering on a percentile of the magnitude. I also implemented minimum peak distance as the inverse of the maximum expected heart rate, which I assumed to be 200 BPM. I unexpectedly got a more accurate result without the Savitsky-Golay filter applied (but have only tried it on a couple of clips).

I plan to mess with some dynamic range compression and circle back to this in the future.

%% Set constants
START_S = 30;
END_S = START_S + 30;
BP_ORDER = 20;
BP_LO = 50;
BP_HI = 200;
% DOWNSAMPLE_DIVISOR = 50;
MIN_PEAK_PROMINENCE = 5E-5; 
MIN_PEAK_DISTANCE = 200^-1*60; % i.e., inverse of beats per second
SG_ORDER = 3; % polynomial order to use for Savitzky-Golay FIR filter
SG_FRAMELEN = 51; % for S-G filter; must be odd number

%% Read audio file
Uz = unzip('sample_HR_audio.zip');
[signal, fs] = audioread(Uz{1});
fprintf('Sample rate is %d Hz.\n', fs);

Sample rate is 44100 Hz.

fprintf('Band-pass FIR filter is %d-%d Hz.\n', BP_LO, BP_HI);

Band-pass FIR filter is 50-200 Hz.

fprintf('Savitzky-Golay filter frame length is %d (%.6f s).\n', SG_FRAMELEN, SG_FRAMELEN / fs);

Savitzky-Golay filter frame length is 51 (0.001156 s).

%% Clip
if exist('START_S', 'var')
    start = fs*START_S + 1;
else
    start = 1;
end
if exist('END_S', 'var')
    end_ = fs * END_S;
else
    end_ = numel(signal);
end
signal = signal(start:end_)*10; % add gain

sample_ct = numel(signal);
T = linspace(0, sample_ct-1, sample_ct).' / fs; % Create Time Vector

%% Graph unfiltered signal amplitude over time
figure('Name', 'Amplitude vs. time (unfiltered)');
plot(T, signal);
grid;
xlabel('Time (s)');
ylabel('Amplitude');

%% Apply bandpass filter (BPF)
bandpassFilter = designfilt('bandpassfir', ...
                            'FilterOrder', BP_ORDER, ...
                            'CutoffFrequency1', BP_LO, ...
                            'CutoffFrequency2', BP_HI, ...
                            'SampleRate', fs);
% signal = filtfilt(bandpassFilter, signal);
figure('Name', 'Amplitude vs. time (BPF)')

plot(T, signal);
grid;
xlabel('Time (s)');
ylabel('Amplitude');

%% TODO: Dynamic range control

%% Display power spectrum waterfall
[p,f,t] = pspectrum(signal, fs, 'spectrogram');
figure('Name', 'Power Spectrum Waterfall (BPF)')
waterfall(f,t,p')
xlabel('Frequency (Hz)')
ylabel('Time (seconds)')
wtf = gca;
wtf.XDir = 'reverse';
colormap(turbo)
view([30 45]) % something about the angle?
xlim([0 250]) % freq range to graph

%% Graph magnitude over time
titletext = 'Magnitude over time';
figure('Name', titletext);
hold on;

p1 = normalize(p(1,:));

plot(t, p1);

xlabel('Time (s)');
ylabel('Normalised magnitude');

pQ1 = prctile(p1, 25);
pQ2 = prctile(p1, 50);
pQ3 = prctile(p1, 75);
customPercentilePoint = 65;
customPercentile = prctile(p1, customPercentilePoint);
pMean = mean(p1);
% yline(pQ1, '--g', sprintf('Q1=%f', pQ1));
% yline(pQ2, '--g', sprintf('Median=%f', pQ2));
% yline(pQ3, '--g', sprintf('Q3=%f', pQ3));
% Yes I know this will say weird things like '73th percentile'
yline(customPercentile, '--r', sprintf('%dth percentile=%f', customPercentilePoint, customPercentile));
yline(pMean, '--b', sprintf('Mean=%f', pMean));

[pks, locs] = findpeaks(p1, ...
    'MinPeakProminence', MIN_PEAK_PROMINENCE, ...
    'MinPeakDistance', MIN_PEAK_DISTANCE, ...
    'MinPeakHeight', customPercentile);
plot(t(locs), pks, 'r^');

hbCt = length(t(locs));
fprintf('\nCounted %d heart beats from magnitude:\n', length(t(locs)));

Counted 55 heart beats from magnitude:

bpm = hbCt / (max(T)/60);
fprintf('Mean: %.1f BPM.\n', bpm);

Mean: 110.0 BPM.

%% Apply Savitsky-Golay filter (smooths noise); find peaks; make new fig.
p1_sg = sgolayfilt(p(1,:), SG_ORDER, SG_FRAMELEN); % Multiband Empirical FIR Filter

figure('Name', 'S-G filtered magnitude over time')
hold on; % going to plot ref lines & peaks
plot(t, p1_sg(1,:));

p1_sg_pctl25 = prctile(p1_sg(1,:), 25);
yline(p1_sg_pctl25, '--r', '25th percentile')

[pks, locs] = findpeaks(p1_sg, ...
    'MinPeakProminence', MIN_PEAK_PROMINENCE, ...
    'MinPeakDistance', MIN_PEAK_DISTANCE);
plot(t(locs), pks, 'r^');

hbCt = length(t(locs));
fprintf('\nCounted %d heart beats from S-G filtered magnitude:\n', hbCt);

Counted 52 heart beats from S-G filtered magnitude:

bpm = hbCt / (max(T)/60);
fprintf('Mean: %.1f BPM.\n', bpm);

Mean: 104.0 BPM.