-
Notifications
You must be signed in to change notification settings - Fork 0
/
instantHR_analysis.m
59 lines (44 loc) · 1.79 KB
/
instantHR_analysis.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
function [instantHR,beatStart]=instantHR_analysis(ecg,fs)
% You may delete any of the following lines of code. You may also add as
% many lines of code as you'd like.
fs = double(fs);
lower_cutoff = 0.5/(fs / 2);
upper_cutoff = 50/(fs / 2);
if upper_cutoff >=1
upper_cutoff = 0.99;
end
if lower_cutoff >= 1
lower_cutoff = 0.001;
end
% Normalize the cutoff frequencies to the Nyquist frequency
Wn = [lower_cutoff, upper_cutoff];
[b, a] = butter(3, Wn);% 3rd order Butterworth bandpass filter normalized to Nyquist frequency (applicable to multiple sampling rates
average_ecg = mean(ecg,2);
%filteredECG = filter(b, a, ecg);
filteredbnECG = filter(b, a, average_ecg);
diffECG = diff(filteredbnECG);
% Squaring
squaredECG = diffECG .^ 2;
window_size = round(0.150 * fs);
% Apply the moving average filter
smoothed_signal = movmean(squaredECG, window_size);
filteredECG = (smoothed_signal - min(smoothed_signal)) / (max(smoothed_signal) - min(smoothed_signal));
plot(filteredECG)
N = length(filteredECG);
mad_threshold = 1.4 * mad(filteredECG, 1);
% outliers based on the MAD threshold
outliers = abs(filteredECG - median(filteredECG)) > 3 * mad_threshold;
heightECG = filteredECG;
heightECG(outliers) = NaN;
% find the max value considering outliers
max_ecg = max(heightECG, [], 'omitnan');
height = max_ecg;
plot(filteredECG)
distance = round(fs/2);
[QRS_Peaks, peakindices] = findpeaks(filteredECG, 'MinPeakHeight', height, 'MinPeakDistance', distance);
beatStart = peakindices.';
numPeaks = length(QRS_Peaks);
intervals = diff(peakindices) / fs;
instantHR = round(60 ./ intervals);
%instantHR = floor(60./peak_to_peak_distance);
end