Roy Sandberg / Mbed 2 deprecated BLE_ECG

I've got some basic filter code setup (but not yet tested).

Dependencies:   BLE_API Queue mbed nRF51822

Fork of BLE_HeartRate by Bluetooth Low Energy

qrsfilt.cpp@62:8e2fbe131b53, 2015-06-28 (annotated)

Committer:
roysandberg
Date:
Sun Jun 28 03:06:00 2015 +0000
Revision:
62:8e2fbe131b53
Working Beat Detection and Analysis

Who changed what in which revision?

UserRevisionLine numberNew contents of line
roysandberg 62:8e2fbe131b53 1 /*****************************************************************************
roysandberg 62:8e2fbe131b53 2 FILE: qrsfilt.cpp
roysandberg 62:8e2fbe131b53 3 AUTHOR: Patrick S. Hamilton
roysandberg 62:8e2fbe131b53 4 REVISED: 5/13/2002
roysandberg 62:8e2fbe131b53 5 ___________________________________________________________________________
roysandberg 62:8e2fbe131b53 6
roysandberg 62:8e2fbe131b53 7 qrsfilt.cpp filter functions to aid beat detecton in electrocardiograms.
roysandberg 62:8e2fbe131b53 8 Copywrite (C) 2000 Patrick S. Hamilton
roysandberg 62:8e2fbe131b53 9
roysandberg 62:8e2fbe131b53 10 This file is free software; you can redistribute it and/or modify it under
roysandberg 62:8e2fbe131b53 11 the terms of the GNU Library General Public License as published by the Free
roysandberg 62:8e2fbe131b53 12 Software Foundation; either version 2 of the License, or (at your option) any
roysandberg 62:8e2fbe131b53 13 later version.
roysandberg 62:8e2fbe131b53 14
roysandberg 62:8e2fbe131b53 15 This software is distributed in the hope that it will be useful, but WITHOUT ANY
roysandberg 62:8e2fbe131b53 16 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
roysandberg 62:8e2fbe131b53 17 PARTICULAR PURPOSE. See the GNU Library General Public License for more
roysandberg 62:8e2fbe131b53 18 details.
roysandberg 62:8e2fbe131b53 19
roysandberg 62:8e2fbe131b53 20 You should have received a copy of the GNU Library General Public License along
roysandberg 62:8e2fbe131b53 21 with this library; if not, write to the Free Software Foundation, Inc., 59
roysandberg 62:8e2fbe131b53 22 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
roysandberg 62:8e2fbe131b53 23
roysandberg 62:8e2fbe131b53 24 You may contact the author by e-mail (pat@eplimited.edu) or postal mail
roysandberg 62:8e2fbe131b53 25 (Patrick Hamilton, E.P. Limited, 35 Medford St., Suite 204 Somerville,
roysandberg 62:8e2fbe131b53 26 MA 02143 USA). For updates to this software, please visit our website
roysandberg 62:8e2fbe131b53 27 (http://www.eplimited.com).
roysandberg 62:8e2fbe131b53 28 __________________________________________________________________________
roysandberg 62:8e2fbe131b53 29
roysandberg 62:8e2fbe131b53 30 This file includes QRSFilt() and associated filtering files used for QRS
roysandberg 62:8e2fbe131b53 31 detection. Only QRSFilt() and deriv1() are called by the QRS detector
roysandberg 62:8e2fbe131b53 32 other functions can be hidden.
roysandberg 62:8e2fbe131b53 33
roysandberg 62:8e2fbe131b53 34 Revisions:
roysandberg 62:8e2fbe131b53 35 5/13: Filter implementations have been modified to allow simplified
roysandberg 62:8e2fbe131b53 36 modification for different sample rates.
roysandberg 62:8e2fbe131b53 37
roysandberg 62:8e2fbe131b53 38 *******************************************************************************/
roysandberg 62:8e2fbe131b53 39
roysandberg 62:8e2fbe131b53 40 #include <mbed.h>
roysandberg 62:8e2fbe131b53 41 //#include <math.h>
roysandberg 62:8e2fbe131b53 42 #include "qrsdet.h"
roysandberg 62:8e2fbe131b53 43
roysandberg 62:8e2fbe131b53 44 // Local Prototypes.
roysandberg 62:8e2fbe131b53 45
roysandberg 62:8e2fbe131b53 46 int lpfilt( int datum ,int init) ;
roysandberg 62:8e2fbe131b53 47 int hpfilt( int datum, int init ) ;
roysandberg 62:8e2fbe131b53 48 int deriv1( int x0, int init ) ;
roysandberg 62:8e2fbe131b53 49 int deriv2( int x0, int init ) ;
roysandberg 62:8e2fbe131b53 50 int mvwint(int datum, int init) ;
roysandberg 62:8e2fbe131b53 51
roysandberg 62:8e2fbe131b53 52 /******************************************************************************
roysandberg 62:8e2fbe131b53 53 * Syntax:
roysandberg 62:8e2fbe131b53 54 * int QRSFilter(int datum, int init) ;
roysandberg 62:8e2fbe131b53 55 * Description:
roysandberg 62:8e2fbe131b53 56 * QRSFilter() takes samples of an ECG signal as input and returns a sample of
roysandberg 62:8e2fbe131b53 57 * a signal that is an estimate of the local energy in the QRS bandwidth. In
roysandberg 62:8e2fbe131b53 58 * other words, the signal has a lump in it whenever a QRS complex, or QRS
roysandberg 62:8e2fbe131b53 59 * complex like artifact occurs. The filters were originally designed for data
roysandberg 62:8e2fbe131b53 60 * sampled at 200 samples per second, but they work nearly as well at sample
roysandberg 62:8e2fbe131b53 61 * frequencies from 150 to 250 samples per second.
roysandberg 62:8e2fbe131b53 62 *
roysandberg 62:8e2fbe131b53 63 * The filter buffers and static variables are reset if a value other than
roysandberg 62:8e2fbe131b53 64 * 0 is passed to QRSFilter through init.
roysandberg 62:8e2fbe131b53 65 *******************************************************************************/
roysandberg 62:8e2fbe131b53 66
roysandberg 62:8e2fbe131b53 67 int QRSFilter(int datum,int init)
roysandberg 62:8e2fbe131b53 68 {
roysandberg 62:8e2fbe131b53 69 int fdatum ;
roysandberg 62:8e2fbe131b53 70
roysandberg 62:8e2fbe131b53 71 if(init)
roysandberg 62:8e2fbe131b53 72 {
roysandberg 62:8e2fbe131b53 73 hpfilt( 0, 1 ) ; // Initialize filters.
roysandberg 62:8e2fbe131b53 74 lpfilt( 0, 1 ) ;
roysandberg 62:8e2fbe131b53 75 mvwint( 0, 1 ) ;
roysandberg 62:8e2fbe131b53 76 deriv1( 0, 1 ) ;
roysandberg 62:8e2fbe131b53 77 deriv2( 0, 1 ) ;
roysandberg 62:8e2fbe131b53 78 }
roysandberg 62:8e2fbe131b53 79
roysandberg 62:8e2fbe131b53 80 fdatum = lpfilt( datum, 0 ) ; // Low pass filter data.
roysandberg 62:8e2fbe131b53 81 fdatum = hpfilt( fdatum, 0 ) ; // High pass filter data.
roysandberg 62:8e2fbe131b53 82 fdatum = deriv2( fdatum, 0 ) ; // Take the derivative.
roysandberg 62:8e2fbe131b53 83 fdatum = abs(fdatum) ; // Take the absolute value.
roysandberg 62:8e2fbe131b53 84 fdatum = mvwint( fdatum, 0 ) ; // Average over an 80 ms window .
roysandberg 62:8e2fbe131b53 85 return(fdatum) ;
roysandberg 62:8e2fbe131b53 86 }
roysandberg 62:8e2fbe131b53 87
roysandberg 62:8e2fbe131b53 88
roysandberg 62:8e2fbe131b53 89 /*************************************************************************
roysandberg 62:8e2fbe131b53 90 * lpfilt() implements the digital filter represented by the difference
roysandberg 62:8e2fbe131b53 91 * equation:
roysandberg 62:8e2fbe131b53 92 *
roysandberg 62:8e2fbe131b53 93 * y[n] = 2*y[n-1] - y[n-2] + x[n] - 2*x[t-24 ms] + x[t-48 ms]
roysandberg 62:8e2fbe131b53 94 *
roysandberg 62:8e2fbe131b53 95 * Note that the filter delay is (LPBUFFER_LGTH/2)-1
roysandberg 62:8e2fbe131b53 96 *
roysandberg 62:8e2fbe131b53 97 **************************************************************************/
roysandberg 62:8e2fbe131b53 98
roysandberg 62:8e2fbe131b53 99 int lpfilt( int datum ,int init)
roysandberg 62:8e2fbe131b53 100 {
roysandberg 62:8e2fbe131b53 101 static long y1 = 0, y2 = 0 ;
roysandberg 62:8e2fbe131b53 102 static int data[LPBUFFER_LGTH], ptr = 0 ;
roysandberg 62:8e2fbe131b53 103 long y0 ;
roysandberg 62:8e2fbe131b53 104 int output, halfPtr ;
roysandberg 62:8e2fbe131b53 105 if(init)
roysandberg 62:8e2fbe131b53 106 {
roysandberg 62:8e2fbe131b53 107 for(ptr = 0; ptr < LPBUFFER_LGTH; ++ptr)
roysandberg 62:8e2fbe131b53 108 data[ptr] = 0 ;
roysandberg 62:8e2fbe131b53 109 y1 = y2 = 0 ;
roysandberg 62:8e2fbe131b53 110 ptr = 0 ;
roysandberg 62:8e2fbe131b53 111 }
roysandberg 62:8e2fbe131b53 112 halfPtr = ptr-(LPBUFFER_LGTH/2) ; // Use halfPtr to index
roysandberg 62:8e2fbe131b53 113 if(halfPtr < 0) // to x[n-6].
roysandberg 62:8e2fbe131b53 114 halfPtr += LPBUFFER_LGTH ;
roysandberg 62:8e2fbe131b53 115 y0 = (y1 << 1) - y2 + datum - (data[halfPtr] << 1) + data[ptr] ;
roysandberg 62:8e2fbe131b53 116 y2 = y1;
roysandberg 62:8e2fbe131b53 117 y1 = y0;
roysandberg 62:8e2fbe131b53 118 output = y0 / ((LPBUFFER_LGTH*LPBUFFER_LGTH)/4);
roysandberg 62:8e2fbe131b53 119 data[ptr] = datum ; // Stick most recent sample into
roysandberg 62:8e2fbe131b53 120 if(++ptr == LPBUFFER_LGTH) // the circular buffer and update
roysandberg 62:8e2fbe131b53 121 ptr = 0 ; // the buffer pointer.
roysandberg 62:8e2fbe131b53 122 return(output) ;
roysandberg 62:8e2fbe131b53 123 }
roysandberg 62:8e2fbe131b53 124
roysandberg 62:8e2fbe131b53 125
roysandberg 62:8e2fbe131b53 126 /******************************************************************************
roysandberg 62:8e2fbe131b53 127 * hpfilt() implements the high pass filter represented by the following
roysandberg 62:8e2fbe131b53 128 * difference equation:
roysandberg 62:8e2fbe131b53 129 *
roysandberg 62:8e2fbe131b53 130 * y[n] = y[n-1] + x[n] - x[n-128 ms]
roysandberg 62:8e2fbe131b53 131 * z[n] = x[n-64 ms] - y[n] ;
roysandberg 62:8e2fbe131b53 132 *
roysandberg 62:8e2fbe131b53 133 * Filter delay is (HPBUFFER_LGTH-1)/2
roysandberg 62:8e2fbe131b53 134 ******************************************************************************/
roysandberg 62:8e2fbe131b53 135
roysandberg 62:8e2fbe131b53 136 int hpfilt( int datum, int init )
roysandberg 62:8e2fbe131b53 137 {
roysandberg 62:8e2fbe131b53 138 static long y=0 ;
roysandberg 62:8e2fbe131b53 139 static int data[HPBUFFER_LGTH], ptr = 0 ;
roysandberg 62:8e2fbe131b53 140 int z, halfPtr ;
roysandberg 62:8e2fbe131b53 141
roysandberg 62:8e2fbe131b53 142 if(init)
roysandberg 62:8e2fbe131b53 143 {
roysandberg 62:8e2fbe131b53 144 for(ptr = 0; ptr < HPBUFFER_LGTH; ++ptr)
roysandberg 62:8e2fbe131b53 145 data[ptr] = 0 ;
roysandberg 62:8e2fbe131b53 146 ptr = 0 ;
roysandberg 62:8e2fbe131b53 147 y = 0 ;
roysandberg 62:8e2fbe131b53 148 }
roysandberg 62:8e2fbe131b53 149
roysandberg 62:8e2fbe131b53 150 y += datum - data[ptr];
roysandberg 62:8e2fbe131b53 151 halfPtr = ptr-(HPBUFFER_LGTH/2) ;
roysandberg 62:8e2fbe131b53 152 if(halfPtr < 0)
roysandberg 62:8e2fbe131b53 153 halfPtr += HPBUFFER_LGTH ;
roysandberg 62:8e2fbe131b53 154 z = data[halfPtr] - (y / HPBUFFER_LGTH);
roysandberg 62:8e2fbe131b53 155
roysandberg 62:8e2fbe131b53 156 data[ptr] = datum ;
roysandberg 62:8e2fbe131b53 157 if(++ptr == HPBUFFER_LGTH)
roysandberg 62:8e2fbe131b53 158 ptr = 0 ;
roysandberg 62:8e2fbe131b53 159
roysandberg 62:8e2fbe131b53 160 return( z );
roysandberg 62:8e2fbe131b53 161 }
roysandberg 62:8e2fbe131b53 162
roysandberg 62:8e2fbe131b53 163 /*****************************************************************************
roysandberg 62:8e2fbe131b53 164 * deriv1 and deriv2 implement derivative approximations represented by
roysandberg 62:8e2fbe131b53 165 * the difference equation:
roysandberg 62:8e2fbe131b53 166 *
roysandberg 62:8e2fbe131b53 167 * y[n] = x[n] - x[n - 10ms]
roysandberg 62:8e2fbe131b53 168 *
roysandberg 62:8e2fbe131b53 169 * Filter delay is DERIV_LENGTH/2
roysandberg 62:8e2fbe131b53 170 *****************************************************************************/
roysandberg 62:8e2fbe131b53 171
roysandberg 62:8e2fbe131b53 172 int deriv1(int x, int init)
roysandberg 62:8e2fbe131b53 173 {
roysandberg 62:8e2fbe131b53 174 static int derBuff[DERIV_LENGTH], derI = 0 ;
roysandberg 62:8e2fbe131b53 175 int y ;
roysandberg 62:8e2fbe131b53 176
roysandberg 62:8e2fbe131b53 177 if(init != 0)
roysandberg 62:8e2fbe131b53 178 {
roysandberg 62:8e2fbe131b53 179 for(derI = 0; derI < DERIV_LENGTH; ++derI)
roysandberg 62:8e2fbe131b53 180 derBuff[derI] = 0 ;
roysandberg 62:8e2fbe131b53 181 derI = 0 ;
roysandberg 62:8e2fbe131b53 182 return(0) ;
roysandberg 62:8e2fbe131b53 183 }
roysandberg 62:8e2fbe131b53 184
roysandberg 62:8e2fbe131b53 185 y = x - derBuff[derI] ;
roysandberg 62:8e2fbe131b53 186 derBuff[derI] = x ;
roysandberg 62:8e2fbe131b53 187 if(++derI == DERIV_LENGTH)
roysandberg 62:8e2fbe131b53 188 derI = 0 ;
roysandberg 62:8e2fbe131b53 189 return(y) ;
roysandberg 62:8e2fbe131b53 190 }
roysandberg 62:8e2fbe131b53 191
roysandberg 62:8e2fbe131b53 192 int deriv2(int x, int init)
roysandberg 62:8e2fbe131b53 193 {
roysandberg 62:8e2fbe131b53 194 static int derBuff[DERIV_LENGTH], derI = 0 ;
roysandberg 62:8e2fbe131b53 195 int y ;
roysandberg 62:8e2fbe131b53 196
roysandberg 62:8e2fbe131b53 197 if(init != 0)
roysandberg 62:8e2fbe131b53 198 {
roysandberg 62:8e2fbe131b53 199 for(derI = 0; derI < DERIV_LENGTH; ++derI)
roysandberg 62:8e2fbe131b53 200 derBuff[derI] = 0 ;
roysandberg 62:8e2fbe131b53 201 derI = 0 ;
roysandberg 62:8e2fbe131b53 202 return(0) ;
roysandberg 62:8e2fbe131b53 203 }
roysandberg 62:8e2fbe131b53 204
roysandberg 62:8e2fbe131b53 205 y = x - derBuff[derI] ;
roysandberg 62:8e2fbe131b53 206 derBuff[derI] = x ;
roysandberg 62:8e2fbe131b53 207 if(++derI == DERIV_LENGTH)
roysandberg 62:8e2fbe131b53 208 derI = 0 ;
roysandberg 62:8e2fbe131b53 209 return(y) ;
roysandberg 62:8e2fbe131b53 210 }
roysandberg 62:8e2fbe131b53 211
roysandberg 62:8e2fbe131b53 212
roysandberg 62:8e2fbe131b53 213
roysandberg 62:8e2fbe131b53 214
roysandberg 62:8e2fbe131b53 215 /*****************************************************************************
roysandberg 62:8e2fbe131b53 216 * mvwint() implements a moving window integrator. Actually, mvwint() averages
roysandberg 62:8e2fbe131b53 217 * the signal values over the last WINDOW_WIDTH samples.
roysandberg 62:8e2fbe131b53 218 *****************************************************************************/
roysandberg 62:8e2fbe131b53 219
roysandberg 62:8e2fbe131b53 220 int mvwint(int datum, int init)
roysandberg 62:8e2fbe131b53 221 {
roysandberg 62:8e2fbe131b53 222 static long sum = 0 ;
roysandberg 62:8e2fbe131b53 223 static int data[WINDOW_WIDTH], ptr = 0 ;
roysandberg 62:8e2fbe131b53 224 int output;
roysandberg 62:8e2fbe131b53 225 if(init)
roysandberg 62:8e2fbe131b53 226 {
roysandberg 62:8e2fbe131b53 227 for(ptr = 0; ptr < WINDOW_WIDTH ; ++ptr)
roysandberg 62:8e2fbe131b53 228 data[ptr] = 0 ;
roysandberg 62:8e2fbe131b53 229 sum = 0 ;
roysandberg 62:8e2fbe131b53 230 ptr = 0 ;
roysandberg 62:8e2fbe131b53 231 }
roysandberg 62:8e2fbe131b53 232 sum += datum ;
roysandberg 62:8e2fbe131b53 233 sum -= data[ptr] ;
roysandberg 62:8e2fbe131b53 234 data[ptr] = datum ;
roysandberg 62:8e2fbe131b53 235 if(++ptr == WINDOW_WIDTH)
roysandberg 62:8e2fbe131b53 236 ptr = 0 ;
roysandberg 62:8e2fbe131b53 237 if((sum / WINDOW_WIDTH) > 32000)
roysandberg 62:8e2fbe131b53 238 output = 32000 ;
roysandberg 62:8e2fbe131b53 239 else
roysandberg 62:8e2fbe131b53 240 output = sum / WINDOW_WIDTH ;
roysandberg 62:8e2fbe131b53 241 return(output) ;
roysandberg 62:8e2fbe131b53 242 }