File content as of revision 62:8e2fbe131b53:

/*****************************************************************************
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) ;
    }