TEAP (Toolbox for Emotion Analysis using Physiological Signals) doc
rpeakdetect.m
Go to the documentation of this file.
1 %> @file rpeakdetect.m
2 function [hrv, R_t, R_amp, R_index, S_t, S_amp] = rpeakdetect(data, samp_freq, thresh, testmode)
3 
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.
13 %
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.
20 %
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.
24 %
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://www.gnu.org/
28 %
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.
33 %
34 % gari@ieee.org
35 
36 %%%%%%%%%%% make threshold default 0.2 -> this was 0.15 on MIT data
37 if nargin < 4
38  testmode = 0;
39 end
40 %%%%%%%%%%% make threshold default 0.2 -> this was 0.15 on MIT data
41 if nargin < 3
42  thresh = 0.2;
43 
44 end
45 %%%%%%%%%%% make sample frequency default 256 Hz
46 if nargin < 2
47  samp_freq = 256;
48  if(testmode==1)
49  fprintf('Assuming sampling frequency of %iHz\n', samp_freq);
50  end
51 end
52 
53 %%%%%%%%%%% check format of data %%%%%%%%%%
54 [a b] = size(data);
55 if(a > b)
56  len = a;
57 end
58 if(b > a)
59  len = b;
60 end
61 
62 %%%%%%%%%% if there's no time axis - make one
63 if (a || b == 1)
64 % make time axis
65  tt = [1/samp_freq:1/samp_freq:ceil(len/samp_freq)];
66  t = tt(1:len);
67  x = data;
68 end
69 %%%%%%%%%% check if data is in columns or rows
70 if (a == 2)
71  x=data(:,1);
72  t=data(:,2);
73 end
74 if (b == 2)
75  t=data(:,1);
76  x=data(:,2);
77 end
78 
79 %%%%%%%%% bandpass filter data - assume 256hz data %%%%%
80  % remove mean
81  x = x-mean(x);
82 
83  % FIR filtering stage
84  bpf=x; %Initialise
85 if( (samp_freq == 128) && (exist('filterECG128Hz') ~= 0) )
86  bpf = filterECG128Hz(x);
87 end
88 if( (samp_freq == 256) && (exist('filterECG256Hz') ~= 0) )
89  bpf = filterECG256Hz(x);
90 end
91 
92 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
93 
94 %%%%%%%%% differentiate data %%%%%%%%%%%%%%%%%%%%%%%%%%%
95  dff = diff(bpf); % now it's one datum shorter than before
96 
97 %%%%%%%%% square data %%%%%%%%%%%%%%%%%%%%%%%%%%%
98  sqr = dff .* dff; %
99  len = len - 1; % how long is the new vector now?
100 
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))];
105  end
106  % integrate
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)));
114 
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);
117 
118  %%%% then build an array of segments to look in %%%%%
119  %thresh = 0.2;
120  poss_reg = mdfint>(thresh*max_h);
121 
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
126 
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
133  end
134 
135  R_index = maxloc;
136  R_t = t(maxloc);
137  R_amp = maxval;
138  S_amp = minval; %%%% Assuming the S-wave is the lowest
139  %%%% amp in the given window
140  S_t = t(minloc);
141 
142 %%%%%%%%%% check for lead inversion %%%%%%%%%%%%%%%%%%%
143  % i.e. do minima precede maxima?
144  if (minloc(length(minloc))<maxloc(length(minloc)))
145  R_t = t(minloc);
146  R_amp = minval;
147  S_t = t(maxloc);
148  S_amp = maxval;
149  end
150 
151 %%%%%%%%%%%%
152 hrv = diff(R_t);
153 resp = R_amp-S_amp;
154 
155 %%%%%%%%%%%%%%%%%%%%
156 if (testmode ~= 0)
157  figure(1)
158  hold off
159  subplot(4, 1, 1)
160  plot(t, x);
161  hold on;
162  plot(t, bpf, 'r')
163  title('raw ECG (blue) and zero-pahse FIR filtered ECG (red)')
164  ylabel('ECG')
165  hold off;
166  subplot(4, 1, 2)
167 
168  plot(t(1:length(mdfint)), mdfint);
169  hold on;
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')
174  ylabel('Int ECG')
175  hold off;
176 
177  subplot(4, 1, 3)
178  plot(t, bpf, 'r');
179  hold on;
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')
183  ylabel('ECG+R+S')
184  hold off;
185  subplot(4, 1, 4)
186  hold off
187  plot(R_t(1:length(hrv)), hrv, 'r+')
188  hold on
189  title('RR intervals')
190  ylabel('RR (s)')
191  hold off
192  fprintf('Press any key for next block of data\n');
193  pause
194 end
filterECG128Hz
function filterECG128Hz(in data)
filterECG256Hz
function filterECG256Hz(in data)
rpeakdetect
function rpeakdetect(in data, in samp_freq, in thresh, in testmode)