bdac.cpp

00001 /*****************************************************************************
00002 FILE:  bdac.cpp
00003 AUTHOR: Patrick S. Hamilton
00004 REVISED:    5/13/2002
00005   ___________________________________________________________________________
00006 
00007 bdac.cpp: Beat Detection And Classification
00008 Copywrite (C) 2001 Patrick S. Hamilton
00009 
00010 This file is free software; you can redistribute it and/or modify it under
00011 the terms of the GNU Library General Public License as published by the Free
00012 Software Foundation; either version 2 of the License, or (at your option) any
00013 later version.
00014 
00015 This software is distributed in the hope that it will be useful, but WITHOUT ANY
00016 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
00017 PARTICULAR PURPOSE.  See the GNU Library General Public License for more
00018 details.
00019 
00020 You should have received a copy of the GNU Library General Public License along
00021 with this library; if not, write to the Free Software Foundation, Inc., 59
00022 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00023 
00024 You may contact the author by e-mail (pat@eplimited.edu) or postal mail
00025 (Patrick Hamilton, E.P. Limited, 35 Medford St., Suite 204 Somerville,
00026 MA 02143 USA).  For updates to this software, please visit our website
00027 (http://www.eplimited.com).
00028   __________________________________________________________________________
00029 
00030 bdac.cpp contains functions for handling Beat Detection And Classification.
00031 The primary function calls a qrs detector.  When a beat is detected it waits
00032 until a sufficient number of samples from the beat have occurred.  When the
00033 beat is ready, BeatDetectAndClassify passes the beat and the timing
00034 information on to the functions that actually classify the beat.
00035 
00036 Functions in bdac.cpp require functions in the following files:
00037         qrsfilt.cpp
00038         qrsdet.cpp
00039         classify.cpp
00040         rythmchk.cpp
00041         noisechk.cpp
00042         analbeat.cpp
00043         match.cpp
00044         postclas.cpp
00045 
00046  __________________________________________________________________________
00047 
00048     Revisions:
00049         5/13/02:
00050             Encapsulated down sampling from input stream to beat template in
00051             the function DownSampleBeat.
00052 
00053             Constants related to time are derived from SAMPLE_RATE in qrsdet
00054          and BEAT_SAMPLE_RATE in bcac.h.
00055 
00056 *******************************************************************************/
00057 #include "qrsdet.h" // For base SAMPLE_RATE
00058 #include "bdac.h"
00059 
00060 #define ECG_BUFFER_LENGTH   1000    // Should be long enough for a beat
00061                                             // plus extra space to accommodate
00062                                             // the maximum detection delay.
00063 #define BEAT_QUE_LENGTH 10          // Length of que for beats awaiting
00064                                             // classification.  Because of
00065                                             // detection delays, Multiple beats
00066                                             // can occur before there is enough data
00067                                             // to classify the first beat in the que.
00068 
00069 // Internal function prototypes.
00070 
00071 void DownSampleBeat(int *beatOut, int *beatIn) ;
00072 
00073 // External function prototypes.
00074 
00075 int QRSDet( int datum, int init ) ;
00076 int NoiseCheck(int datum, int delay, int RR, int beatBegin, int beatEnd) ;
00077 int Classify(int *newBeat,int rr, int noiseLevel, int *beatMatch, int *fidAdj, int init) ;
00078 int GetDominantType(void) ;
00079 int GetBeatEnd(int type) ;
00080 int GetBeatBegin(int type) ;
00081 int gcd(int x, int y) ;
00082 
00083 // Global Variables
00084 
00085 int ECGBuffer[ECG_BUFFER_LENGTH], ECGBufferIndex = 0 ;  // Circular data buffer.
00086 int BeatBuffer[BEATLGTH] ;
00087 int BeatQue[BEAT_QUE_LENGTH], BeatQueCount = 0 ;  // Buffer of detection delays.
00088 int RRCount = 0 ;
00089 int InitBeatFlag = 1 ;
00090 
00091 /******************************************************************************
00092     ResetBDAC() resets static variables required for beat detection and
00093     classification.
00094 *******************************************************************************/
00095 
00096 void ResetBDAC(void)
00097     {
00098     int dummy ;
00099     QRSDet(0,1) ;   // Reset the qrs detector
00100     RRCount = 0 ;
00101     Classify(BeatBuffer,0,0,&dummy,&dummy,1) ;
00102     InitBeatFlag = 1 ;
00103    BeatQueCount = 0 ;   // Flush the beat que.
00104     }
00105 
00106 /*****************************************************************************
00107 Syntax:
00108     int BeatDetectAndClassify(int ecgSample, int *beatType, *beatMatch) ;
00109 
00110 Description:
00111     BeatDetectAndClassify() implements a beat detector and classifier.
00112     ECG samples are passed into BeatDetectAndClassify() one sample at a
00113     time.  BeatDetectAndClassify has been designed for a sample rate of
00114     200 Hz.  When a beat has been detected and classified the detection
00115     delay is returned and the beat classification is returned through the
00116     pointer *beatType.  For use in debugging, the number of the template
00117    that the beat was matched to is returned in via *beatMatch.
00118 
00119 Returns
00120     BeatDetectAndClassify() returns 0 if no new beat has been detected and
00121     classified.  If a beat has been classified, BeatDetectAndClassify returns
00122     the number of samples since the approximate location of the R-wave.
00123 
00124 ****************************************************************************/
00125 
00126 int BeatDetectAndClassify(int ecgSample, int *beatType, int *beatMatch)
00127     {
00128     int detectDelay, rr, i, j ;
00129     int noiseEst = 0, beatBegin, beatEnd ;
00130     int domType ;
00131     int fidAdj ;
00132     int tempBeat[(SAMPLE_RATE/BEAT_SAMPLE_RATE)*BEATLGTH] ;
00133 
00134     // Store new sample in the circular buffer.
00135 
00136     ECGBuffer[ECGBufferIndex] = ecgSample ;
00137     if(++ECGBufferIndex == ECG_BUFFER_LENGTH)
00138         ECGBufferIndex = 0 ;
00139 
00140     // Increment RRInterval count.
00141 
00142     ++RRCount ;
00143 
00144     // Increment detection delays for any beats in the que.
00145 
00146     for(i = 0; i < BeatQueCount; ++i)
00147         ++BeatQue[i] ;
00148 
00149     // Run the sample through the QRS detector.
00150 
00151     detectDelay = QRSDet(ecgSample,0) ;
00152     if(detectDelay != 0)
00153         {
00154         BeatQue[BeatQueCount] = detectDelay ;
00155         ++BeatQueCount ;
00156         }
00157 
00158     // Return if no beat is ready for classification.
00159 
00160     if((BeatQue[0] < (BEATLGTH-FIDMARK)*(SAMPLE_RATE/BEAT_SAMPLE_RATE))
00161         || (BeatQueCount == 0))
00162         {
00163         NoiseCheck(ecgSample,0,rr, beatBegin, beatEnd) ;    // Update noise check buffer
00164         return 0 ;
00165         }
00166 
00167     // Otherwise classify the beat at the head of the que.
00168 
00169     rr = RRCount - BeatQue[0] ; // Calculate the R-to-R interval
00170     detectDelay = RRCount = BeatQue[0] ;
00171 
00172     // Estimate low frequency noise in the beat.
00173     // Might want to move this into classify().
00174 
00175     domType = GetDominantType() ;
00176     if(domType == -1)
00177         {
00178         beatBegin = MS250 ;
00179         beatEnd = MS300 ;
00180         }
00181     else
00182         {
00183         beatBegin = (SAMPLE_RATE/BEAT_SAMPLE_RATE)*(FIDMARK-GetBeatBegin(domType)) ;
00184         beatEnd = (SAMPLE_RATE/BEAT_SAMPLE_RATE)*(GetBeatEnd(domType)-FIDMARK) ;
00185         }
00186     noiseEst = NoiseCheck(ecgSample,detectDelay,rr,beatBegin,beatEnd) ;
00187 
00188     // Copy the beat from the circular buffer to the beat buffer
00189     // and reduce the sample rate by averageing pairs of data
00190     // points.
00191 
00192     j = ECGBufferIndex - detectDelay - (SAMPLE_RATE/BEAT_SAMPLE_RATE)*FIDMARK ;
00193     if(j < 0) j += ECG_BUFFER_LENGTH ;
00194 
00195     for(i = 0; i < (SAMPLE_RATE/BEAT_SAMPLE_RATE)*BEATLGTH; ++i)
00196         {
00197         tempBeat[i] = ECGBuffer[j] ;
00198         if(++j == ECG_BUFFER_LENGTH)
00199             j = 0 ;
00200         }
00201 
00202     DownSampleBeat(BeatBuffer,tempBeat) ;
00203 
00204     // Update the QUE.
00205 
00206     for(i = 0; i < BeatQueCount-1; ++i)
00207         BeatQue[i] = BeatQue[i+1] ;
00208     --BeatQueCount ;
00209 
00210 
00211     // Skip the first beat.
00212 
00213     if(InitBeatFlag)
00214         {
00215         InitBeatFlag = 0 ;
00216         *beatType = 13 ;
00217         *beatMatch = 0 ;
00218         fidAdj = 0 ;
00219         }
00220 
00221     // Classify all other beats.
00222 
00223     else
00224         {
00225         *beatType = Classify(BeatBuffer,rr,noiseEst,beatMatch,&fidAdj,0) ;
00226         fidAdj *= SAMPLE_RATE/BEAT_SAMPLE_RATE ;
00227       }
00228 
00229     // Ignore detection if the classifier decides that this
00230     // was the trailing edge of a PVC.
00231 
00232     if(*beatType == 100)
00233         {
00234         RRCount += rr ;
00235         return(0) ;
00236         }
00237 
00238     // Limit the fiducial mark adjustment in case of problems with
00239     // beat onset and offset estimation.
00240 
00241     if(fidAdj > MS80)
00242         fidAdj = MS80 ;
00243     else if(fidAdj < -MS80)
00244         fidAdj = -MS80 ;
00245 
00246     return(detectDelay-fidAdj) ;
00247     }
00248 
00249 void DownSampleBeat(int *beatOut, int *beatIn)
00250     {
00251     int i ;
00252 
00253     for(i = 0; i < BEATLGTH; ++i)
00254         beatOut[i] = (beatIn[i<<2]+
00255                       beatIn[(i<<2)+1]+
00256                       beatIn[(i<<2)+2]+
00257                       beatIn[(i<<2)+3]
00258                       )>>2 ;
00259     }
00260 
00261 
00262 void DownSampleBeatOld(int *beatOut, int *beatIn)
00263     {
00264     int i ;
00265 
00266     for(i = 0; i < BEATLGTH; ++i)
00267         beatOut[i] = (beatIn[i<<1]+beatIn[(i<<1)+1])>>1 ;
00268     }