Diff: qrsfilt.cpp

Revision:

62:8e2fbe131b53

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/qrsfilt.cpp	Sun Jun 28 03:06:00 2015 +0000
@@ -0,0 +1,242 @@
+/*****************************************************************************
+FILE:  qrsfilt.cpp
+AUTHOR: Patrick S. Hamilton
+REVISED:    5/13/2002
+  ___________________________________________________________________________
+
+qrsfilt.cpp filter functions to aid beat detecton in electrocardiograms.
+Copywrite (C) 2000 Patrick S. Hamilton
+
+This file is free software; you can redistribute it and/or modify it under
+the terms of the GNU Library General Public License as published by the Free
+Software Foundation; either version 2 of the License, or (at your option) any
+later version.
+
+This software is distributed in the hope that it will be useful, but WITHOUT ANY
+WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
+PARTICULAR PURPOSE.  See the GNU Library General Public License for more
+details.
+
+You should have received a copy of the GNU Library General Public License along
+with this library; if not, write to the Free Software Foundation, Inc., 59
+Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+
+You may contact the author by e-mail (pat@eplimited.edu) or postal mail
+(Patrick Hamilton, E.P. Limited, 35 Medford St., Suite 204 Somerville,
+MA 02143 USA).  For updates to this software, please visit our website
+(http://www.eplimited.com).
+  __________________________________________________________________________
+
+    This file includes QRSFilt() and associated filtering files used for QRS
+    detection.  Only QRSFilt() and deriv1() are called by the QRS detector
+    other functions can be hidden.
+
+    Revisions:
+        5/13: Filter implementations have been modified to allow simplified
+            modification for different sample rates.
+
+*******************************************************************************/
+
+#include <mbed.h>
+//#include <math.h>
+#include "qrsdet.h"
+
+// Local Prototypes.
+
+int lpfilt( int datum ,int init) ;
+int hpfilt( int datum, int init ) ;
+int deriv1( int x0, int init ) ;
+int deriv2( int x0, int init ) ;
+int mvwint(int datum, int init) ;
+
+/******************************************************************************
+* Syntax:
+*   int QRSFilter(int datum, int init) ;
+* Description:
+*   QRSFilter() takes samples of an ECG signal as input and returns a sample of
+*   a signal that is an estimate of the local energy in the QRS bandwidth.  In
+*   other words, the signal has a lump in it whenever a QRS complex, or QRS
+*   complex like artifact occurs.  The filters were originally designed for data
+*  sampled at 200 samples per second, but they work nearly as well at sample
+*   frequencies from 150 to 250 samples per second.
+*
+*   The filter buffers and static variables are reset if a value other than
+*   0 is passed to QRSFilter through init.
+*******************************************************************************/
+
+int QRSFilter(int datum,int init)
+    {
+    int fdatum ;
+
+    if(init)
+        {
+        hpfilt( 0, 1 ) ;        // Initialize filters.
+        lpfilt( 0, 1 ) ;
+        mvwint( 0, 1 ) ;
+        deriv1( 0, 1 ) ;
+        deriv2( 0, 1 ) ;
+        }
+
+    fdatum = lpfilt( datum, 0 ) ;       // Low pass filter data.
+    fdatum = hpfilt( fdatum, 0 ) ;  // High pass filter data.
+    fdatum = deriv2( fdatum, 0 ) ;  // Take the derivative.
+    fdatum = abs(fdatum) ;              // Take the absolute value.
+    fdatum = mvwint( fdatum, 0 ) ;  // Average over an 80 ms window .
+    return(fdatum) ;
+    }
+
+
+/*************************************************************************
+*  lpfilt() implements the digital filter represented by the difference
+*  equation:
+*
+*   y[n] = 2*y[n-1] - y[n-2] + x[n] - 2*x[t-24 ms] + x[t-48 ms]
+*
+*   Note that the filter delay is (LPBUFFER_LGTH/2)-1
+*
+**************************************************************************/
+
+int lpfilt( int datum ,int init)
+    {
+    static long y1 = 0, y2 = 0 ;
+    static int data[LPBUFFER_LGTH], ptr = 0 ;
+    long y0 ;
+    int output, halfPtr ;
+    if(init)
+        {
+        for(ptr = 0; ptr < LPBUFFER_LGTH; ++ptr)
+            data[ptr] = 0 ;
+        y1 = y2 = 0 ;
+        ptr = 0 ;
+        }
+    halfPtr = ptr-(LPBUFFER_LGTH/2) ;   // Use halfPtr to index
+    if(halfPtr < 0)                         // to x[n-6].
+        halfPtr += LPBUFFER_LGTH ;
+    y0 = (y1 << 1) - y2 + datum - (data[halfPtr] << 1) + data[ptr] ;
+    y2 = y1;
+    y1 = y0;
+    output = y0 / ((LPBUFFER_LGTH*LPBUFFER_LGTH)/4);
+    data[ptr] = datum ;         // Stick most recent sample into
+    if(++ptr == LPBUFFER_LGTH)  // the circular buffer and update
+        ptr = 0 ;                   // the buffer pointer.
+    return(output) ;
+    }
+
+
+/******************************************************************************
+*  hpfilt() implements the high pass filter represented by the following
+*  difference equation:
+*
+*   y[n] = y[n-1] + x[n] - x[n-128 ms]
+*   z[n] = x[n-64 ms] - y[n] ;
+*
+*  Filter delay is (HPBUFFER_LGTH-1)/2
+******************************************************************************/
+
+int hpfilt( int datum, int init )
+    {
+    static long y=0 ;
+    static int data[HPBUFFER_LGTH], ptr = 0 ;
+    int z, halfPtr ;
+
+    if(init)
+        {
+        for(ptr = 0; ptr < HPBUFFER_LGTH; ++ptr)
+            data[ptr] = 0 ;
+        ptr = 0 ;
+        y = 0 ;
+        }
+
+    y += datum - data[ptr];
+    halfPtr = ptr-(HPBUFFER_LGTH/2) ;
+    if(halfPtr < 0)
+        halfPtr += HPBUFFER_LGTH ;
+    z = data[halfPtr] - (y / HPBUFFER_LGTH);
+
+    data[ptr] = datum ;
+    if(++ptr == HPBUFFER_LGTH)
+        ptr = 0 ;
+
+    return( z );
+    }
+
+/*****************************************************************************
+*  deriv1 and deriv2 implement derivative approximations represented by
+*  the difference equation:
+*
+*   y[n] = x[n] - x[n - 10ms]
+*
+*  Filter delay is DERIV_LENGTH/2
+*****************************************************************************/
+
+int deriv1(int x, int init)
+    {
+    static int derBuff[DERIV_LENGTH], derI = 0 ;
+    int y ;
+
+    if(init != 0)
+        {
+        for(derI = 0; derI < DERIV_LENGTH; ++derI)
+            derBuff[derI] = 0 ;
+        derI = 0 ;
+        return(0) ;
+        }
+
+    y = x - derBuff[derI] ;
+    derBuff[derI] = x ;
+    if(++derI == DERIV_LENGTH)
+        derI = 0 ;
+    return(y) ;
+    }
+
+int deriv2(int x, int init)
+    {
+    static int derBuff[DERIV_LENGTH], derI = 0 ;
+    int y ;
+
+    if(init != 0)
+        {
+        for(derI = 0; derI < DERIV_LENGTH; ++derI)
+            derBuff[derI] = 0 ;
+        derI = 0 ;
+        return(0) ;
+        }
+
+    y = x - derBuff[derI] ;
+    derBuff[derI] = x ;
+    if(++derI == DERIV_LENGTH)
+        derI = 0 ;
+    return(y) ;
+    }
+
+
+
+
+/*****************************************************************************
+* mvwint() implements a moving window integrator.  Actually, mvwint() averages
+* the signal values over the last WINDOW_WIDTH samples.
+*****************************************************************************/
+
+int mvwint(int datum, int init)
+    {
+    static long sum = 0 ;
+    static int data[WINDOW_WIDTH], ptr = 0 ;
+    int output;
+    if(init)
+        {
+        for(ptr = 0; ptr < WINDOW_WIDTH ; ++ptr)
+            data[ptr] = 0 ;
+        sum = 0 ;
+        ptr = 0 ;
+        }
+    sum += datum ;
+    sum -= data[ptr] ;
+    data[ptr] = datum ;
+    if(++ptr == WINDOW_WIDTH)
+        ptr = 0 ;
+    if((sum / WINDOW_WIDTH) > 32000)
+        output = 32000 ;
+    else
+        output = sum / WINDOW_WIDTH ;
+    return(output) ;
+    }