Convert sound of heart beat into beats per minute

45 次查看（过去 30 天）
Elyssa Watford
Elyssa Watford 2018-9-17
评论： Joseph Tricarico 2024-2-19
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.
  8 个评论
显示 6更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Star Strider
Star Strider 2019-2-21
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.
  4 个评论
显示 2更早的评论
Star Strider
Star Strider 2024-2-1
Ran in:
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
Joseph Tricarico 2024-2-19
Ran in:
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Audio and Video Data 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by