2 function [hrv, R_t, R_amp, R_index, S_t, S_amp] =
rpeakdetect(data, samp_freq, thresh, testmode)
4 % [hrv, R_t, R_amp, R_index, S_t, S_amp] =
rpeakdetect(data, samp_freq, thresh, testmode);
5 % R_t == RR points in time, R_amp == amplitude
6 % of R peak in bpf data & S_amp == amplitude of
7 % following minmum. sampling frequency (samp_freq = 256Hz
8 % by
default) only needed
if no time vector is
9 % specified (assumed to be 1st column or row).
10 % The
'triggering' threshold
'thresh' for the peaks in the
'integrated'
11 % waveform is 0.2 by
default. testmode = 0 (
default) indicates
12 % no graphics diagnostics. Otherwise, you get to scan through each segment.
14 % A batch QRS detector based upon that of Pan, Hamilton and Tompkins:
15 % J. Pan \& W. Tompkins - A real-time QRS detection algorithm
16 % IEEE Transactions on Biomedical Engineering, vol. BME-32 NO. 3. 1985.
17 % P. Hamilton \& W. Tompkins. Quantitative Investigation of QRS
18 % Detection Rules Using the MIT/BIH Arrythmia Database.
19 % IEEE Transactions on Biomedical Engineering, vol. BME-33, NO. 12.1986.
21 % Similar results reported by the authors above were achieved, without
22 % having to tune the thresholds on the MIT DB. An online version in C
23 % has also been written.
25 % Written by G. Clifford gari@ieee.org and made available under the
26 % GNU general
public license. If you have not received a copy of
this
27 % license, please download a copy from http:
29 % Please distribute (and modify) freely, commenting
30 % where you have added modifications.
31 % The author would appreciate correspondence regarding
32 % corrections, modifications, improvements etc.
36 %%%%%%%%%%% make threshold
default 0.2 ->
this was 0.15 on MIT data
40 %%%%%%%%%%% make threshold
default 0.2 ->
this was 0.15 on MIT data
45 %%%%%%%%%%% make sample frequency
default 256 Hz
49 fprintf(
'Assuming sampling frequency of %iHz\n', samp_freq);
53 %%%%%%%%%%% check format of data %%%%%%%%%%
62 %%%%%%%%%%
if there
's no time axis - make one
65 tt = [1/samp_freq:1/samp_freq:ceil(len/samp_freq)];
69 %%%%%%%%%% check if data is in columns or rows
79 %%%%%%%%% bandpass filter data - assume 256hz data %%%%%
86 bpf = filterECG128Hz(x);
89 bpf = filterECG256Hz(x);
92 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
94 %%%%%%%%% differentiate data %%%%%%%%%%%%%%%%%%%%%%%%%%%
95 dff = diff(bpf); % now it's one datum shorter than before
97 %%%%%%%%% square data %%%%%%%%%%%%%%%%%%%%%%%%%%%
99 len = len - 1; % how
long is the
new vector now?
101 %%%%%%%%% integrate data over window
'd' %%%%%%%%%%%%%%%%%%%%%%%%%
102 d = [1 1 1 1 1 1 1]; % window size - intialise
103 if(samp_freq >= 256) % adapt
for higher sampling rates
104 d = [ones(1, round(7*samp_freq/256))];
107 mdfint = medfilt1(filter(d, 1, sqr), 10);
108 % remove filter delay
for scanning back through ECG
109 delay = ceil(length(d)/2);
110 mdfint = mdfint(delay:length(mdfint));
111 %%%%%%%%% segment search area %%%%%%%%%%%%%%%%%%%%%%%
112 %%%% first find the highest bumps in the data %%%%%%
113 % max_h = max (mdfint(round(len/4):round(3*len/4)));
115 %GC: The maximum height is replaced by the 95% percentile to avoid outliers problemes
116 max_h = prctile(mdfint(round(len/4):round(3*len/4)),95);
118 %%%% then build an array of segments to look in %%%%%
120 poss_reg = mdfint>(thresh*max_h);
122 %%%%%%%%% and find peaks %%%%%%%%%%%%%%%%%%%%%%%%%%%%
123 %%%% find indices into boudaries of each segment %%%
124 left = find(diff([0 poss_reg
']) == 1); % remember to zero pad at start
125 right = find(diff([poss_reg' 0]) == -1); % remember to zero pad at end
127 %%%% loop through all possibilities
128 for(i=1:length(left))
129 [maxval(i) maxloc(i)] = max( bpf(left(i):right(i)) );
130 [minval(i) minloc(i)] = min( bpf(left(i):right(i)) );
131 maxloc(i) = maxloc(i)-1+left(i); % add offset of present location
132 minloc(i) = minloc(i)-1+left(i); % add offset of present location
138 S_amp = minval; %%%% Assuming the S-wave is the lowest
139 %%%% amp in the given window
142 %%%%%%%%%% check
for lead inversion %%%%%%%%%%%%%%%%%%%
143 % i.e.
do minima precede maxima?
144 if (minloc(length(minloc))<maxloc(length(minloc)))
163 title('raw ECG (blue) and zero-pahse FIR filtered ECG (red)')
168 plot(t(1:length(mdfint)), mdfint);
170 plot(t, max(mdfint)*bpf/(2*max(bpf)), 'r')
171 plot(t(left), mdfint(left), 'og')
172 plot(t(right), mdfint(right), 'om')
173 title('Integrated data with scan boundaries over scaled ECG')
180 plot(R_t, R_amp, '+k');
181 plot(S_t, S_amp, '+g');
182 title('ECG with R-peaks (black) and S-points (green) over ECG')
187 plot(R_t(1:length(hrv)), hrv, 'r+')
189 title('RR intervals')
192 fprintf('Press any key for next block of data\n');