Roy Sandberg
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
Revision 62:8e2fbe131b53, committed 2015-06-28
- Comitter:
- roysandberg
- Date:
- Sun Jun 28 03:06:00 2015 +0000
- Parent:
- 61:1de72bdab0ef
- Commit message:
- Working Beat Detection and Analysis
Changed in this revision
diff -r 1de72bdab0ef -r 8e2fbe131b53 Queue.lib --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Queue.lib Sun Jun 28 03:06:00 2015 +0000 @@ -0,0 +1,1 @@ +http://developer.mbed.org/users/wbasser/code/Queue/#a03810d46457
diff -r 1de72bdab0ef -r 8e2fbe131b53 analbeat.cpp
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/analbeat.cpp Sun Jun 28 03:06:00 2015 +0000
@@ -0,0 +1,384 @@
+/*****************************************************************************
+FILE: analbeat.cpp
+AUTHOR: Patrick S. Hamilton
+REVISED: 5/13/2002
+ ___________________________________________________________________________
+
+analbeat.cpp: Analyze Beat
+Copywrite (C) 2001 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 contains functions for determining the QRS onset, QRS offset,
+ beat onset, beat offset, polarity, and isoelectric level for a beat.
+
+ Revisions:
+ 4/16: Modified to prevent isoStart from being set to less than ISO_LENGTH1-1
+ 5/13/2002: Time related constants are tied to BEAT_SAMPLE_RATE in bdac.h.
+
+*****************************************************************************/
+#include <mbed.h>
+#include "bdac.h"
+//#include <stdio.h>
+//#include <stdlib.h>
+
+#define ISO_LENGTH1 BEAT_MS50
+#define ISO_LENGTH2 BEAT_MS80
+#define ISO_LIMIT 20
+
+// Local prototypes.
+
+int IsoCheck(int *data, int isoLength);
+
+/****************************************************************
+ IsoCheck determines whether the amplitudes of a run
+ of data fall within a sufficiently small amplitude that
+ the run can be considered isoelectric.
+*****************************************************************/
+
+int IsoCheck(int *data, int isoLength)
+ {
+ int i, max, min ;
+
+ for(i = 1, max=min = data[0]; i < isoLength; ++i)
+ {
+ if(data[i] > max)
+ max = data[i] ;
+ else if(data[i] < min)
+ min = data[i] ;
+ }
+
+ if(max - min < ISO_LIMIT)
+ return(1) ;
+ return(0) ;
+ }
+
+/**********************************************************************
+ AnalyzeBeat takes a beat buffer as input and returns (via pointers)
+ estimates of the QRS onset, QRS offset, polarity, isoelectric level
+ beat beginning (P-wave onset), and beat ending (T-wave offset).
+ Analyze Beat assumes that the beat has been sampled at 100 Hz, is
+ BEATLGTH long, and has an R-wave location of roughly FIDMARK.
+
+ Note that beatBegin is the number of samples before FIDMARK that
+ the beat begins and beatEnd is the number of samples after the
+ FIDMARK that the beat ends.
+************************************************************************/
+
+#define INF_CHK_N BEAT_MS40
+
+void AnalyzeBeat(int *beat, int *onset, int *offset, int *isoLevel,
+ int *beatBegin, int *beatEnd, int *amp)
+ {
+ int maxSlope = 0, maxSlopeI, minSlope = 0, minSlopeI ;
+ int maxV, minV ;
+ int isoStart, isoEnd ;
+ int slope, i ;
+
+ // Search back from the fiducial mark to find the isoelectric
+ // region preceeding the QRS complex.
+
+ for(i = FIDMARK-ISO_LENGTH2; (i > 0) && (IsoCheck(&beat[i],ISO_LENGTH2) == 0); --i) ;
+
+ // If the first search didn't turn up any isoelectric region, look for
+ // a shorter isoelectric region.
+
+ if(i == 0)
+ {
+ for(i = FIDMARK-ISO_LENGTH1; (i > 0) && (IsoCheck(&beat[i],ISO_LENGTH1) == 0); --i) ;
+ isoStart = i + (ISO_LENGTH1 - 1) ;
+ }
+ else isoStart = i + (ISO_LENGTH2 - 1) ;
+
+ // Search forward from the R-wave to find an isoelectric region following
+ // the QRS complex.
+
+ for(i = FIDMARK; (i < BEATLGTH) && (IsoCheck(&beat[i],ISO_LENGTH1) == 0); ++i) ;
+ isoEnd = i ;
+
+ // Find the maximum and minimum slopes on the
+ // QRS complex.
+
+ i = FIDMARK-BEAT_MS150 ;
+ maxSlope = maxSlope = beat[i] - beat[i-1] ;
+ maxSlopeI = minSlopeI = i ;
+
+ for(; i < FIDMARK+BEAT_MS150; ++i)
+ {
+ slope = beat[i] - beat[i-1] ;
+ if(slope > maxSlope)
+ {
+ maxSlope = slope ;
+ maxSlopeI = i ;
+ }
+ else if(slope < minSlope)
+ {
+ minSlope = slope ;
+ minSlopeI = i ;
+ }
+ }
+
+ // Use the smallest of max or min slope for search parameters.
+
+ if(maxSlope > -minSlope)
+ maxSlope = -minSlope ;
+ else minSlope = -maxSlope ;
+
+ if(maxSlopeI < minSlopeI)
+ {
+
+ // Search back from the maximum slope point for the QRS onset.
+
+ for(i = maxSlopeI;
+ (i > 0) && ((beat[i]-beat[i-1]) > (maxSlope >> 2)); --i) ;
+ *onset = i-1 ;
+
+ // Check to see if this was just a brief inflection.
+
+ for(; (i > *onset-INF_CHK_N) && ((beat[i]-beat[i-1]) <= (maxSlope >>2)); --i) ;
+ if(i > *onset-INF_CHK_N)
+ {
+ for(;(i > 0) && ((beat[i]-beat[i-1]) > (maxSlope >> 2)); --i) ;
+ *onset = i-1 ;
+ }
+ i = *onset+1 ;
+
+ // Check to see if a large negative slope follows an inflection.
+ // If so, extend the onset a little more.
+
+ for(;(i > *onset-INF_CHK_N) && ((beat[i-1]-beat[i]) < (maxSlope>>2)); --i);
+ if(i > *onset-INF_CHK_N)
+ {
+ for(; (i > 0) && ((beat[i-1]-beat[i]) > (maxSlope>>2)); --i) ;
+ *onset = i-1 ;
+ }
+
+ // Search forward from minimum slope point for QRS offset.
+
+ for(i = minSlopeI;
+ (i < BEATLGTH) && ((beat[i] - beat[i-1]) < (minSlope >>2)); ++i);
+ *offset = i ;
+
+ // Make sure this wasn't just an inflection.
+
+ for(; (i < *offset+INF_CHK_N) && ((beat[i]-beat[i-1]) >= (minSlope>>2)); ++i) ;
+ if(i < *offset+INF_CHK_N)
+ {
+ for(;(i < BEATLGTH) && ((beat[i]-beat[i-1]) < (minSlope >>2)); ++i) ;
+ *offset = i ;
+ }
+ i = *offset ;
+
+ // Check to see if there is a significant upslope following
+ // the end of the down slope.
+
+ for(;(i < *offset+BEAT_MS40) && ((beat[i-1]-beat[i]) > (minSlope>>2)); ++i);
+ if(i < *offset+BEAT_MS40)
+ {
+ for(; (i < BEATLGTH) && ((beat[i-1]-beat[i]) < (minSlope>>2)); ++i) ;
+ *offset = i ;
+
+ // One more search motivated by PVC shape in 123.
+
+ for(; (i < *offset+BEAT_MS60) && (beat[i]-beat[i-1] > (minSlope>>2)); ++i) ;
+ if(i < *offset + BEAT_MS60)
+ {
+ for(;(i < BEATLGTH) && (beat[i]-beat[i-1] < (minSlope>>2)); ++i) ;
+ *offset = i ;
+ }
+ }
+ }
+
+ else
+ {
+
+ // Search back from the minimum slope point for the QRS onset.
+
+ for(i = minSlopeI;
+ (i > 0) && ((beat[i]-beat[i-1]) < (minSlope >> 2)); --i) ;
+ *onset = i-1 ;
+
+ // Check to see if this was just a brief inflection.
+
+ for(; (i > *onset-INF_CHK_N) && ((beat[i]-beat[i-1]) >= (minSlope>>2)); --i) ;
+ if(i > *onset-INF_CHK_N)
+ {
+ for(; (i > 0) && ((beat[i]-beat[i-1]) < (minSlope>>2));--i) ;
+ *onset = i-1 ;
+ }
+ i = *onset+1 ;
+
+ // Check for significant positive slope after a turning point.
+
+ for(;(i > *onset-INF_CHK_N) && ((beat[i-1]-beat[i]) > (minSlope>>2)); --i);
+ if(i > *onset-INF_CHK_N)
+ {
+ for(; (i > 0) && ((beat[i-1]-beat[i]) < (minSlope>>2)); --i) ;
+ *onset = i-1 ;
+ }
+
+ // Search forward from maximum slope point for QRS offset.
+
+ for(i = maxSlopeI;
+ (i < BEATLGTH) && ((beat[i] - beat[i-1]) > (maxSlope >>2)); ++i) ;
+ *offset = i ;
+
+ // Check to see if this was just a brief inflection.
+
+ for(; (i < *offset+INF_CHK_N) && ((beat[i] - beat[i-1]) <= (maxSlope >> 2)); ++i) ;
+ if(i < *offset+INF_CHK_N)
+ {
+ for(;(i < BEATLGTH) && ((beat[i] - beat[i-1]) > (maxSlope >>2)); ++i) ;
+ *offset = i ;
+ }
+ i = *offset ;
+
+ // Check to see if there is a significant downslope following
+ // the end of the up slope.
+
+ for(;(i < *offset+BEAT_MS40) && ((beat[i-1]-beat[i]) < (maxSlope>>2)); ++i);
+ if(i < *offset+BEAT_MS40)
+ {
+ for(; (i < BEATLGTH) && ((beat[i-1]-beat[i]) > (maxSlope>>2)); ++i) ;
+ *offset = i ;
+ }
+ }
+
+ // If the estimate of the beginning of the isoelectric level was
+ // at the beginning of the beat, use the slope based QRS onset point
+ // as the iso electric point.
+
+ if((isoStart == ISO_LENGTH1-1)&& (*onset > isoStart)) // ** 4/19 **
+ isoStart = *onset ;
+
+ // Otherwise, if the isoelectric start and the slope based points
+ // are close, use the isoelectric start point.
+
+ else if(*onset-isoStart < BEAT_MS50)
+ *onset = isoStart ;
+
+ // If the isoelectric end and the slope based QRS offset are close
+ // use the isoelectic based point.
+
+ if(isoEnd - *offset < BEAT_MS50)
+ *offset = isoEnd ;
+
+ *isoLevel = beat[isoStart] ;
+
+
+ // Find the maximum and minimum values in the QRS.
+
+ for(i = *onset, maxV = minV = beat[*onset]; i < *offset; ++i)
+ if(beat[i] > maxV)
+ maxV = beat[i] ;
+ else if(beat[i] < minV)
+ minV = beat[i] ;
+
+ // If the offset is significantly below the onset and the offset is
+ // on a negative slope, add the next up slope to the width.
+
+ if((beat[*onset]-beat[*offset] > ((maxV-minV)>>2)+((maxV-minV)>>3)))
+ {
+
+ // Find the maximum slope between the finish and the end of the buffer.
+
+ for(i = maxSlopeI = *offset, maxSlope = beat[*offset] - beat[*offset-1];
+ (i < *offset+BEAT_MS100) && (i < BEATLGTH); ++i)
+ {
+ slope = beat[i]-beat[i-1] ;
+ if(slope > maxSlope)
+ {
+ maxSlope = slope ;
+ maxSlopeI = i ;
+ }
+ }
+
+ // Find the new offset.
+
+ if(maxSlope > 0)
+ {
+ for(i = maxSlopeI;
+ (i < BEATLGTH) && (beat[i]-beat[i-1] > (maxSlope>>1)); ++i) ;
+ *offset = i ;
+ }
+ }
+
+ // Determine beginning and ending of the beat.
+ // Search for an isoelectric region that precedes the R-wave.
+ // by at least 250 ms.
+
+ for(i = FIDMARK-BEAT_MS250;
+ (i >= BEAT_MS80) && (IsoCheck(&beat[i-BEAT_MS80],BEAT_MS80) == 0); --i) ;
+ *beatBegin = i ;
+
+ // If there was an isoelectric section at 250 ms before the
+ // R-wave, search forward for the isoelectric region closest
+ // to the R-wave. But leave at least 50 ms between beat begin
+ // and onset, or else normal beat onset is within PVC QRS complexes.
+ // that screws up noise estimation.
+
+ if(*beatBegin == FIDMARK-BEAT_MS250)
+ {
+ for(; (i < *onset-BEAT_MS50) &&
+ (IsoCheck(&beat[i-BEAT_MS80],BEAT_MS80) == 1); ++i) ;
+ *beatBegin = i-1 ;
+ }
+
+ // Rev 1.1
+ else if(*beatBegin == BEAT_MS80 - 1)
+ {
+ for(; (i < *onset) && (IsoCheck(&beat[i-BEAT_MS80],BEAT_MS80) == 0); ++i);
+ if(i < *onset)
+ {
+ for(; (i < *onset) && (IsoCheck(&beat[i-BEAT_MS80],BEAT_MS80) == 1); ++i) ;
+ if(i < *onset)
+ *beatBegin = i-1 ;
+ }
+ }
+
+ // Search for the end of the beat as the first isoelectric
+ // segment that follows the beat by at least 300 ms.
+
+ for(i = FIDMARK+BEAT_MS300;
+ (i < BEATLGTH) && (IsoCheck(&beat[i],BEAT_MS80) == 0); ++i) ;
+ *beatEnd = i ;
+
+ // If the signal was isoelectric at 300 ms. search backwards.
+/* if(*beatEnd == FIDMARK+30)
+ {
+ for(; (i > *offset) && (IsoCheck(&beat[i],8) != 0); --i) ;
+ *beatEnd = i ;
+ }
+*/
+ // Calculate beat amplitude.
+
+ maxV=minV=beat[*onset] ;
+ for(i = *onset; i < *offset; ++i)
+ if(beat[i] > maxV)
+ maxV = beat[i] ;
+ else if(beat[i] < minV)
+ minV = beat[i] ;
+ *amp = maxV-minV ;
+
+ }
+
+
+
diff -r 1de72bdab0ef -r 8e2fbe131b53 analbeat.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/analbeat.h Sun Jun 28 03:06:00 2015 +0000 @@ -0,0 +1,33 @@ +/***************************************************************************** +FILE: analbeat.h +AUTHOR: Patrick S. Hamilton +REVISED: 12/4/2001 + ___________________________________________________________________________ + +analbeat.h: Beat analysis prototype definition. +Copywrite (C) 2001 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). +******************************************************************************/ + +// External prototypes for analbeat.cpp + +void AnalyzeBeat(int *beat, int *onset, int *offset, + int *isoLevel, int *beatBegin, int *beatEnd, int *amp) ;
diff -r 1de72bdab0ef -r 8e2fbe131b53 bdac.cpp
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/bdac.cpp Sun Jun 28 03:06:00 2015 +0000
@@ -0,0 +1,268 @@
+/*****************************************************************************
+FILE: bdac.cpp
+AUTHOR: Patrick S. Hamilton
+REVISED: 5/13/2002
+ ___________________________________________________________________________
+
+bdac.cpp: Beat Detection And Classification
+Copywrite (C) 2001 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).
+ __________________________________________________________________________
+
+bdac.cpp contains functions for handling Beat Detection And Classification.
+The primary function calls a qrs detector. When a beat is detected it waits
+until a sufficient number of samples from the beat have occurred. When the
+beat is ready, BeatDetectAndClassify passes the beat and the timing
+information on to the functions that actually classify the beat.
+
+Functions in bdac.cpp require functions in the following files:
+ qrsfilt.cpp
+ qrsdet.cpp
+ classify.cpp
+ rythmchk.cpp
+ noisechk.cpp
+ analbeat.cpp
+ match.cpp
+ postclas.cpp
+
+ __________________________________________________________________________
+
+ Revisions:
+ 5/13/02:
+ Encapsulated down sampling from input stream to beat template in
+ the function DownSampleBeat.
+
+ Constants related to time are derived from SAMPLE_RATE in qrsdet
+ and BEAT_SAMPLE_RATE in bcac.h.
+
+*******************************************************************************/
+#include "qrsdet.h" // For base SAMPLE_RATE
+#include "bdac.h"
+
+#define ECG_BUFFER_LENGTH 1000 // Should be long enough for a beat
+ // plus extra space to accommodate
+ // the maximum detection delay.
+#define BEAT_QUE_LENGTH 10 // Length of que for beats awaiting
+ // classification. Because of
+ // detection delays, Multiple beats
+ // can occur before there is enough data
+ // to classify the first beat in the que.
+
+// Internal function prototypes.
+
+void DownSampleBeat(int *beatOut, int *beatIn) ;
+
+// External function prototypes.
+
+int QRSDet( int datum, int init ) ;
+int NoiseCheck(int datum, int delay, int RR, int beatBegin, int beatEnd) ;
+int Classify(int *newBeat,int rr, int noiseLevel, int *beatMatch, int *fidAdj, int init) ;
+int GetDominantType(void) ;
+int GetBeatEnd(int type) ;
+int GetBeatBegin(int type) ;
+int gcd(int x, int y) ;
+
+// Global Variables
+
+int ECGBuffer[ECG_BUFFER_LENGTH], ECGBufferIndex = 0 ; // Circular data buffer.
+int BeatBuffer[BEATLGTH] ;
+int BeatQue[BEAT_QUE_LENGTH], BeatQueCount = 0 ; // Buffer of detection delays.
+int RRCount = 0 ;
+int InitBeatFlag = 1 ;
+
+/******************************************************************************
+ ResetBDAC() resets static variables required for beat detection and
+ classification.
+*******************************************************************************/
+
+void ResetBDAC(void)
+ {
+ int dummy ;
+ QRSDet(0,1) ; // Reset the qrs detector
+ RRCount = 0 ;
+ Classify(BeatBuffer,0,0,&dummy,&dummy,1) ;
+ InitBeatFlag = 1 ;
+ BeatQueCount = 0 ; // Flush the beat que.
+ }
+
+/*****************************************************************************
+Syntax:
+ int BeatDetectAndClassify(int ecgSample, int *beatType, *beatMatch) ;
+
+Description:
+ BeatDetectAndClassify() implements a beat detector and classifier.
+ ECG samples are passed into BeatDetectAndClassify() one sample at a
+ time. BeatDetectAndClassify has been designed for a sample rate of
+ 200 Hz. When a beat has been detected and classified the detection
+ delay is returned and the beat classification is returned through the
+ pointer *beatType. For use in debugging, the number of the template
+ that the beat was matched to is returned in via *beatMatch.
+
+Returns
+ BeatDetectAndClassify() returns 0 if no new beat has been detected and
+ classified. If a beat has been classified, BeatDetectAndClassify returns
+ the number of samples since the approximate location of the R-wave.
+
+****************************************************************************/
+
+int BeatDetectAndClassify(int ecgSample, int *beatType, int *beatMatch)
+ {
+ int detectDelay, rr, i, j ;
+ int noiseEst = 0, beatBegin, beatEnd ;
+ int domType ;
+ int fidAdj ;
+ int tempBeat[(SAMPLE_RATE/BEAT_SAMPLE_RATE)*BEATLGTH] ;
+
+ // Store new sample in the circular buffer.
+
+ ECGBuffer[ECGBufferIndex] = ecgSample ;
+ if(++ECGBufferIndex == ECG_BUFFER_LENGTH)
+ ECGBufferIndex = 0 ;
+
+ // Increment RRInterval count.
+
+ ++RRCount ;
+
+ // Increment detection delays for any beats in the que.
+
+ for(i = 0; i < BeatQueCount; ++i)
+ ++BeatQue[i] ;
+
+ // Run the sample through the QRS detector.
+
+ detectDelay = QRSDet(ecgSample,0) ;
+ if(detectDelay != 0)
+ {
+ BeatQue[BeatQueCount] = detectDelay ;
+ ++BeatQueCount ;
+ }
+
+ // Return if no beat is ready for classification.
+
+ if((BeatQue[0] < (BEATLGTH-FIDMARK)*(SAMPLE_RATE/BEAT_SAMPLE_RATE))
+ || (BeatQueCount == 0))
+ {
+ NoiseCheck(ecgSample,0,rr, beatBegin, beatEnd) ; // Update noise check buffer
+ return 0 ;
+ }
+
+ // Otherwise classify the beat at the head of the que.
+
+ rr = RRCount - BeatQue[0] ; // Calculate the R-to-R interval
+ detectDelay = RRCount = BeatQue[0] ;
+
+ // Estimate low frequency noise in the beat.
+ // Might want to move this into classify().
+
+ domType = GetDominantType() ;
+ if(domType == -1)
+ {
+ beatBegin = MS250 ;
+ beatEnd = MS300 ;
+ }
+ else
+ {
+ beatBegin = (SAMPLE_RATE/BEAT_SAMPLE_RATE)*(FIDMARK-GetBeatBegin(domType)) ;
+ beatEnd = (SAMPLE_RATE/BEAT_SAMPLE_RATE)*(GetBeatEnd(domType)-FIDMARK) ;
+ }
+ noiseEst = NoiseCheck(ecgSample,detectDelay,rr,beatBegin,beatEnd) ;
+
+ // Copy the beat from the circular buffer to the beat buffer
+ // and reduce the sample rate by averageing pairs of data
+ // points.
+
+ j = ECGBufferIndex - detectDelay - (SAMPLE_RATE/BEAT_SAMPLE_RATE)*FIDMARK ;
+ if(j < 0) j += ECG_BUFFER_LENGTH ;
+
+ for(i = 0; i < (SAMPLE_RATE/BEAT_SAMPLE_RATE)*BEATLGTH; ++i)
+ {
+ tempBeat[i] = ECGBuffer[j] ;
+ if(++j == ECG_BUFFER_LENGTH)
+ j = 0 ;
+ }
+
+ DownSampleBeat(BeatBuffer,tempBeat) ;
+
+ // Update the QUE.
+
+ for(i = 0; i < BeatQueCount-1; ++i)
+ BeatQue[i] = BeatQue[i+1] ;
+ --BeatQueCount ;
+
+
+ // Skip the first beat.
+
+ if(InitBeatFlag)
+ {
+ InitBeatFlag = 0 ;
+ *beatType = 13 ;
+ *beatMatch = 0 ;
+ fidAdj = 0 ;
+ }
+
+ // Classify all other beats.
+
+ else
+ {
+ *beatType = Classify(BeatBuffer,rr,noiseEst,beatMatch,&fidAdj,0) ;
+ fidAdj *= SAMPLE_RATE/BEAT_SAMPLE_RATE ;
+ }
+
+ // Ignore detection if the classifier decides that this
+ // was the trailing edge of a PVC.
+
+ if(*beatType == 100)
+ {
+ RRCount += rr ;
+ return(0) ;
+ }
+
+ // Limit the fiducial mark adjustment in case of problems with
+ // beat onset and offset estimation.
+
+ if(fidAdj > MS80)
+ fidAdj = MS80 ;
+ else if(fidAdj < -MS80)
+ fidAdj = -MS80 ;
+
+ return(detectDelay-fidAdj) ;
+ }
+
+void DownSampleBeat(int *beatOut, int *beatIn)
+ {
+ int i ;
+
+ for(i = 0; i < BEATLGTH; ++i)
+ beatOut[i] = (beatIn[i<<2]+
+ beatIn[(i<<2)+1]+
+ beatIn[(i<<2)+2]+
+ beatIn[(i<<2)+3]
+ )>>2 ;
+ }
+
+
+void DownSampleBeatOld(int *beatOut, int *beatIn)
+ {
+ int i ;
+
+ for(i = 0; i < BEATLGTH; ++i)
+ beatOut[i] = (beatIn[i<<1]+beatIn[(i<<1)+1])>>1 ;
+ }
diff -r 1de72bdab0ef -r 8e2fbe131b53 bdac.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/bdac.h Sun Jun 28 03:06:00 2015 +0000 @@ -0,0 +1,57 @@ +/***************************************************************************** +FILE: bdac.h +AUTHOR: Patrick S. Hamilton +REVISED: 9/25/2001 + ___________________________________________________________________________ + +bdac.h: Beat detection and classification parameter definitions. +Copywrite (C) 2001 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). +******************************************************************************/ +// was 100 +#define BEAT_SAMPLE_RATE 125 + +#define BEAT_MS_PER_SAMPLE ( (double) 1000/ (double) BEAT_SAMPLE_RATE) + +#define BEAT_MS10 ((int) (10/BEAT_MS_PER_SAMPLE + 0.5)) +#define BEAT_MS20 ((int) (20/BEAT_MS_PER_SAMPLE + 0.5)) +#define BEAT_MS40 ((int) (40/BEAT_MS_PER_SAMPLE + 0.5)) +#define BEAT_MS50 ((int) (50/BEAT_MS_PER_SAMPLE + 0.5)) +#define BEAT_MS60 ((int) (60/BEAT_MS_PER_SAMPLE + 0.5)) +#define BEAT_MS70 ((int) (70/BEAT_MS_PER_SAMPLE + 0.5)) +#define BEAT_MS80 ((int) (80/BEAT_MS_PER_SAMPLE + 0.5)) +#define BEAT_MS90 ((int) (90/BEAT_MS_PER_SAMPLE + 0.5)) +#define BEAT_MS100 ((int) (100/BEAT_MS_PER_SAMPLE + 0.5)) +#define BEAT_MS110 ((int) (110/BEAT_MS_PER_SAMPLE + 0.5)) +#define BEAT_MS130 ((int) (130/BEAT_MS_PER_SAMPLE + 0.5)) +#define BEAT_MS140 ((int) (140/BEAT_MS_PER_SAMPLE + 0.5)) +#define BEAT_MS150 ((int) (150/BEAT_MS_PER_SAMPLE + 0.5)) +#define BEAT_MS250 ((int) (250/BEAT_MS_PER_SAMPLE + 0.5)) +#define BEAT_MS280 ((int) (280/BEAT_MS_PER_SAMPLE + 0.5)) +#define BEAT_MS300 ((int) (300/BEAT_MS_PER_SAMPLE + 0.5)) +#define BEAT_MS350 ((int) (350/BEAT_MS_PER_SAMPLE + 0.5)) +#define BEAT_MS400 ((int) (400/BEAT_MS_PER_SAMPLE + 0.5)) +#define BEAT_MS1000 BEAT_SAMPLE_RATE + +#define BEATLGTH BEAT_MS1000 +#define MAXTYPES 8 +#define FIDMARK BEAT_MS400 +
diff -r 1de72bdab0ef -r 8e2fbe131b53 classify.cpp
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/classify.cpp Sun Jun 28 03:06:00 2015 +0000
@@ -0,0 +1,792 @@
+/*****************************************************************************
+FILE: classify.cpp
+AUTHOR: Patrick S. Hamilton
+REVISED: 5/13/2001
+ ___________________________________________________________________________
+
+classify.cpp: Classify a given beat.
+Copywrite (C) 2001 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).
+ __________________________________________________________________________
+
+ Classify.cpp contains functions for classifying beats. The only
+ function that needs to be called externally from this file is Classify().
+
+ Functions in classify.cpp require functions in the following files:
+ match.cpp
+ rythmchk.cpp
+ classify.cpp
+ rythmchk.cpp
+ analbeat.cpp
+ postclas.cpp
+
+ __________________________________________________________________________
+
+ Revisions:
+ 5/13/02:
+ Width constants tied to BEAT_SAMPLE_RATE in bdac.h
+
+ Arrays added to track the classifications and RR intervals for the
+ most recent 8 beats, allowing GetRunCount to become a local function.
+ RR intervals and classifications are now passed to PostClassify.
+
+ Determination of whether the dominant rhythm is regular is now made
+ by examining the number of RR intervals classified as UNKNOWN in the
+ last DM_BUFFER_LENGTH beats (180). If more than 60 are UNKNOWN
+ the rhythm is too irregular to give any weight to whether the beat
+ was premature or not.
+
+*******************************************************************************/
+
+#include "ecgcodes.h"
+//#include <stdlib.h> // For abs()
+//#include <stdio.h>
+#include <mbed.h>
+#include "qrsdet.h" // For base sample rate.
+#include "bdac.h"
+#include "match.h"
+#include "rythmchk.h"
+#include "analbeat.h"
+#include "postclas.h"
+
+// Detection Rule Parameters.
+
+#define MATCH_LIMIT 1.3 // Threshold for template matching
+ // without amplitude sensitivity.
+#define MATCH_WITH_AMP_LIMIT 2.5 // Threshold for matching index that
+ // is amplitude sensitive.
+#define PVC_MATCH_WITH_AMP_LIMIT 0.9 // Amplitude sensitive limit for
+ //matching premature beats
+#define BL_SHIFT_LIMIT 100 // Threshold for assuming a baseline shift.
+#define NEW_TYPE_NOISE_THRESHOLD 18 // Above this noise level, do not create
+ // new beat types.
+#define NEW_TYPE_HF_NOISE_LIMIT 75 // Above this noise level, do not crate
+ // new beat types.
+
+#define MATCH_NOISE_THRESHOLD 0.7 // Match threshold below which noise
+ // indications are ignored.
+
+// TempClass classification rule parameters.
+
+#define R2_DI_THRESHOLD 1.0 // Rule 2 dominant similarity index threshold
+#define R3_WIDTH_THRESHOLD BEAT_MS90 // Rule 3 width threshold.
+#define R7_DI_THRESHOLD 1.2 // Rule 7 dominant similarity index threshold
+#define R8_DI_THRESHOLD 1.5 // Rule 8 dominant similarity index threshold
+#define R9_DI_THRESHOLD 2.0 // Rule 9 dominant similarity index threshold
+#define R10_BC_LIM 3 // Rule 10 beat count limit.
+#define R10_DI_THRESHOLD 2.5 // Rule 10 dominant similarity index threshold
+#define R11_MIN_WIDTH BEAT_MS110 // Rule 11 minimum width threshold.
+#define R11_WIDTH_BREAK BEAT_MS140 // Rule 11 width break.
+#define R11_WIDTH_DIFF1 BEAT_MS40 // Rule 11 width difference threshold 1
+#define R11_WIDTH_DIFF2 BEAT_MS60 // Rule 11 width difference threshold 2
+#define R11_HF_THRESHOLD 45 // Rule 11 high frequency noise threshold.
+#define R11_MA_THRESHOLD 14 // Rule 11 motion artifact threshold.
+#define R11_BC_LIM 1 // Rule 11 beat count limit.
+#define R15_DI_THRESHOLD 3.5 // Rule 15 dominant similarity index threshold
+#define R15_WIDTH_THRESHOLD BEAT_MS100 // Rule 15 width threshold.
+#define R16_WIDTH_THRESHOLD BEAT_MS100 // Rule 16 width threshold.
+#define R17_WIDTH_DELTA BEAT_MS20 // Rule 17 difference threshold.
+#define R18_DI_THRESHOLD 1.5 // Rule 18 dominant similarity index threshold.
+#define R19_HF_THRESHOLD 75 // Rule 19 high frequency noise threshold.
+
+// Dominant monitor constants.
+
+#define DM_BUFFER_LENGTH 180
+#define IRREG_RR_LIMIT 60
+
+// Local prototypes.
+
+int HFNoiseCheck(int *beat) ;
+int TempClass(int rhythmClass, int morphType, int beatWidth, int domWidth,
+ int domType, int hfNoise, int noiseLevel, int blShift, double domIndex) ;
+int DomMonitor(int morphType, int rhythmClass, int beatWidth, int rr, int reset) ;
+int GetDomRhythm(void) ;
+int GetRunCount(void) ;
+
+// Local Global variables
+
+int DomType ;
+int RecentRRs[8], RecentTypes[8] ;
+
+/***************************************************************************
+* Classify() takes a beat buffer, the previous rr interval, and the present
+* noise level estimate and returns a beat classification of NORMAL, PVC, or
+* UNKNOWN. The UNKNOWN classification is only returned. The beat template
+* type that the beat has been matched to is returned through the pointer
+* *beatMatch for debugging display. Passing anything other than 0 in init
+* resets the static variables used by Classify.
+****************************************************************************/
+
+int Classify(int *newBeat,int rr, int noiseLevel, int *beatMatch, int *fidAdj,
+ int init)
+ {
+ int rhythmClass, beatClass, i, beatWidth, blShift ;
+ static int morphType, runCount = 0 ;
+ double matchIndex, domIndex, mi2 ;
+ int shiftAdj ;
+ int domType, domWidth, onset, offset, amp ;
+ int beatBegin, beatEnd, tempClass ;
+ int hfNoise, isoLevel ;
+ static int lastIsoLevel=0, lastRhythmClass = UNKNOWN, lastBeatWasNew = 0 ;
+
+ // If initializing...
+
+ if(init)
+ {
+ ResetRhythmChk() ;
+ ResetMatch() ;
+ ResetPostClassify() ;
+ runCount = 0 ;
+ DomMonitor(0, 0, 0, 0, 1) ;
+ return(0) ;
+ }
+
+ hfNoise = HFNoiseCheck(newBeat) ; // Check for muscle noise.
+ rhythmClass = RhythmChk(rr) ; // Check the rhythm.
+
+ // Estimate beat features.
+
+ AnalyzeBeat(newBeat, &onset, &offset, &isoLevel,
+ &beatBegin, &beatEnd, &) ;
+
+ blShift = abs(lastIsoLevel-isoLevel) ;
+ lastIsoLevel = isoLevel ;
+
+ // Make isoelectric level 0.
+
+ for(i = 0; i < BEATLGTH; ++i)
+ newBeat[i] -= isoLevel ;
+
+ // If there was a significant baseline shift since
+ // the last beat and the last beat was a new type,
+ // delete the new type because it might have resulted
+ // from a baseline shift.
+
+ if( (blShift > BL_SHIFT_LIMIT)
+ && (lastBeatWasNew == 1)
+ && (lastRhythmClass == NORMAL)
+ && (rhythmClass == NORMAL) )
+ ClearLastNewType() ;
+
+ lastBeatWasNew = 0 ;
+
+ // Find the template that best matches this beat.
+
+ BestMorphMatch(newBeat,&morphType,&matchIndex,&mi2,&shiftAdj) ;
+
+ // Disregard noise if the match is good. (New)
+
+ if(matchIndex < MATCH_NOISE_THRESHOLD)
+ hfNoise = noiseLevel = blShift = 0 ;
+
+ // Apply a stricter match limit to premature beats.
+
+ if((matchIndex < MATCH_LIMIT) && (rhythmClass == PVC) &&
+ MinimumBeatVariation(morphType) && (mi2 > PVC_MATCH_WITH_AMP_LIMIT))
+ {
+ morphType = NewBeatType(newBeat) ;
+ lastBeatWasNew = 1 ;
+ }
+
+ // Match if within standard match limits.
+
+ else if((matchIndex < MATCH_LIMIT) && (mi2 <= MATCH_WITH_AMP_LIMIT))
+ UpdateBeatType(morphType,newBeat,mi2,shiftAdj) ;
+
+ // If the beat isn't noisy but doesn't match, start a new beat.
+
+ else if((blShift < BL_SHIFT_LIMIT) && (noiseLevel < NEW_TYPE_NOISE_THRESHOLD)
+ && (hfNoise < NEW_TYPE_HF_NOISE_LIMIT))
+ {
+ morphType = NewBeatType(newBeat) ;
+ lastBeatWasNew = 1 ;
+ }
+
+ // Even if it is a noisy, start new beat if it was an irregular beat.
+
+ else if((lastRhythmClass != NORMAL) || (rhythmClass != NORMAL))
+ {
+ morphType = NewBeatType(newBeat) ;
+ lastBeatWasNew = 1 ;
+ }
+
+ // If its noisy and regular, don't waste space starting a new beat.
+
+ else morphType = MAXTYPES ;
+
+ // Update recent rr and type arrays.
+
+ for(i = 7; i > 0; --i)
+ {
+ RecentRRs[i] = RecentRRs[i-1] ;
+ RecentTypes[i] = RecentTypes[i-1] ;
+ }
+ RecentRRs[0] = rr ;
+ RecentTypes[0] = morphType ;
+
+ lastRhythmClass = rhythmClass ;
+ lastIsoLevel = isoLevel ;
+
+ // Fetch beat features needed for classification.
+ // Get features from average beat if it matched.
+
+ if(morphType != MAXTYPES)
+ {
+ beatClass = GetBeatClass(morphType) ;
+ beatWidth = GetBeatWidth(morphType) ;
+ *fidAdj = GetBeatCenter(morphType)-FIDMARK ;
+
+ // If the width seems large and there have only been a few
+ // beats of this type, use the actual beat for width
+ // estimate.
+
+ if((beatWidth > offset-onset) && (GetBeatTypeCount(morphType) <= 4))
+ {
+ beatWidth = offset-onset ;
+ *fidAdj = ((offset+onset)/2)-FIDMARK ;
+ }
+ }
+
+ // If this beat didn't match get beat features directly
+ // from this beat.
+
+ else
+ {
+ beatWidth = offset-onset ;
+ beatClass = UNKNOWN ;
+ *fidAdj = ((offset+onset)/2)-FIDMARK ;
+ }
+
+ // Fetch dominant type beat features.
+
+ DomType = domType = DomMonitor(morphType, rhythmClass, beatWidth, rr, 0) ;
+ domWidth = GetBeatWidth(domType) ;
+
+ // Compare the beat type, or actual beat to the dominant beat.
+
+ if((morphType != domType) && (morphType != 8))
+ domIndex = DomCompare(morphType,domType) ;
+ else if(morphType == 8)
+ domIndex = DomCompare2(newBeat,domType) ;
+ else domIndex = matchIndex ;
+
+ // Update post classificaton of the previous beat.
+
+ PostClassify(RecentTypes, domType, RecentRRs, beatWidth, domIndex, rhythmClass) ;
+
+ // Classify regardless of how the morphology
+ // was previously classified.
+
+ tempClass = TempClass(rhythmClass, morphType, beatWidth, domWidth,
+ domType, hfNoise, noiseLevel, blShift, domIndex) ;
+
+ // If this morphology has not been classified yet, attempt to classify
+ // it.
+
+ if((beatClass == UNKNOWN) && (morphType < MAXTYPES))
+ {
+
+ // Classify as normal if there are 6 in a row
+ // or at least two in a row that meet rhythm
+ // rules for normal.
+
+ runCount = GetRunCount() ;
+
+ // Classify a morphology as NORMAL if it is not too wide, and there
+ // are three in a row. The width criterion prevents ventricular beats
+ // from being classified as normal during VTACH (MIT/BIH 205).
+
+ if((runCount >= 3) && (domType != -1) && (beatWidth < domWidth+BEAT_MS20))
+ SetBeatClass(morphType,NORMAL) ;
+
+ // If there is no dominant type established yet, classify any type
+ // with six in a row as NORMAL.
+
+ else if((runCount >= 6) && (domType == -1))
+ SetBeatClass(morphType,NORMAL) ;
+
+ // During bigeminy, classify the premature beats as ventricular if
+ // they are not too narrow.
+
+ else if(IsBigeminy() == 1)
+ {
+ if((rhythmClass == PVC) && (beatWidth > BEAT_MS100))
+ SetBeatClass(morphType,PVC) ;
+ else if(rhythmClass == NORMAL)
+ SetBeatClass(morphType,NORMAL) ;
+ }
+ }
+
+ // Save morphology type of this beat for next classification.
+
+ *beatMatch = morphType ;
+
+ beatClass = GetBeatClass(morphType) ;
+
+ // If the morphology has been previously classified.
+ // use that classification.
+ // return(rhythmClass) ;
+
+ if(beatClass != UNKNOWN)
+ return(beatClass) ;
+
+ if(CheckPostClass(morphType) == PVC)
+ return(PVC) ;
+
+ // Otherwise use the temporary classification.
+
+ return(tempClass) ;
+ }
+
+/**************************************************************************
+* HFNoiseCheck() gauges the high frequency (muscle noise) present in the
+* beat template. The high frequency noise level is estimated by highpass
+* filtering the beat (y[n] = x[n] - 2*x[n-1] + x[n-2]), averaging the
+* highpass filtered signal over five samples, and finding the maximum of
+* this averaged highpass filtered signal. The high frequency noise metric
+* is then taken to be the ratio of the maximum averaged highpassed signal
+* to the QRS amplitude.
+**************************************************************************/
+
+#define AVELENGTH BEAT_MS50
+
+int HFNoiseCheck(int *beat)
+ {
+ int maxNoiseAve = 0, i ;
+ int sum=0, aveBuff[AVELENGTH], avePtr = 0 ;
+ int qrsMax = 0, qrsMin = 0 ;
+
+ // Determine the QRS amplitude.
+
+ for(i = FIDMARK-BEAT_MS70; i < FIDMARK+BEAT_MS80; ++i)
+ if(beat[i] > qrsMax)
+ qrsMax = beat[i] ;
+ else if(beat[i] < qrsMin)
+ qrsMin = beat[i] ;
+
+ for(i = 0; i < AVELENGTH; ++i)
+ aveBuff[i] = 0 ;
+
+ for(i = FIDMARK-BEAT_MS280; i < FIDMARK+BEAT_MS280; ++i)
+ {
+ sum -= aveBuff[avePtr] ;
+ aveBuff[avePtr] = abs(beat[i] - (beat[i-BEAT_MS10]<<1) + beat[i-2*BEAT_MS10]) ;
+ sum += aveBuff[avePtr] ;
+ if(++avePtr == AVELENGTH)
+ avePtr = 0 ;
+ if((i < (FIDMARK - BEAT_MS50)) || (i > (FIDMARK + BEAT_MS110)))
+ if(sum > maxNoiseAve)
+ maxNoiseAve = sum ;
+ }
+ if((qrsMax - qrsMin)>=4)
+ return((maxNoiseAve * (50/AVELENGTH))/((qrsMax-qrsMin)>>2)) ;
+ else return(0) ;
+ }
+
+/************************************************************************
+* TempClass() classifies beats based on their beat features, relative
+* to the features of the dominant beat and the present noise level.
+*************************************************************************/
+
+int TempClass(int rhythmClass, int morphType,
+ int beatWidth, int domWidth, int domType,
+ int hfNoise, int noiseLevel, int blShift, double domIndex)
+ {
+
+ // Rule 1: If no dominant type has been detected classify all
+ // beats as UNKNOWN.
+
+ if(domType < 0)
+ return(UNKNOWN) ;
+
+ // Rule 2: If the dominant rhythm is normal, the dominant
+ // beat type doesn't vary much, this beat is premature
+ // and looks sufficiently different than the dominant beat
+ // classify as PVC.
+
+ if(MinimumBeatVariation(domType) && (rhythmClass == PVC)
+ && (domIndex > R2_DI_THRESHOLD) && (GetDomRhythm() == 1))
+ return(PVC) ;
+
+ // Rule 3: If the beat is sufficiently narrow, classify as normal.
+
+ if(beatWidth < R3_WIDTH_THRESHOLD)
+ return(NORMAL) ;
+
+ // Rule 5: If the beat cannot be matched to any previously
+ // detected morphology and it is not premature, consider it normal
+ // (probably noisy).
+
+ if((morphType == MAXTYPES) && (rhythmClass != PVC)) // == UNKNOWN
+ return(NORMAL) ;
+
+ // Rule 6: If the maximum number of beat types have been stored,
+ // this beat is not regular or premature and only one
+ // beat of this morphology has been seen, call it normal (probably
+ // noisy).
+
+ if((GetTypesCount() == MAXTYPES) && (GetBeatTypeCount(morphType)==1)
+ && (rhythmClass == UNKNOWN))
+ return(NORMAL) ;
+
+ // Rule 7: If this beat looks like the dominant beat and the
+ // rhythm is regular, call it normal.
+
+ if((domIndex < R7_DI_THRESHOLD) && (rhythmClass == NORMAL))
+ return(NORMAL) ;
+
+ // Rule 8: If post classification rhythm is normal for this
+ // type and its shape is close to the dominant shape, classify
+ // as normal.
+
+ if((domIndex < R8_DI_THRESHOLD) && (CheckPCRhythm(morphType) == NORMAL))
+ return(NORMAL) ;
+
+ // Rule 9: If the beat is not premature, it looks similar to the dominant
+ // beat type, and the dominant beat type is variable (noisy), classify as
+ // normal.
+
+ if((domIndex < R9_DI_THRESHOLD) && (rhythmClass != PVC) && WideBeatVariation(domType))
+ return(NORMAL) ;
+
+ // Rule 10: If this beat is significantly different from the dominant beat
+ // there have previously been matching beats, the post rhythm classification
+ // of this type is PVC, and the dominant rhythm is regular, classify as PVC.
+
+ if((domIndex > R10_DI_THRESHOLD)
+ && (GetBeatTypeCount(morphType) >= R10_BC_LIM) &&
+ (CheckPCRhythm(morphType) == PVC) && (GetDomRhythm() == 1))
+ return(PVC) ;
+
+ // Rule 11: if the beat is wide, wider than the dominant beat, doesn't
+ // appear to be noisy, and matches a previous type, classify it as
+ // a PVC.
+
+ if( (beatWidth >= R11_MIN_WIDTH) &&
+ (((beatWidth - domWidth >= R11_WIDTH_DIFF1) && (domWidth < R11_WIDTH_BREAK)) ||
+ (beatWidth - domWidth >= R11_WIDTH_DIFF2)) &&
+ (hfNoise < R11_HF_THRESHOLD) && (noiseLevel < R11_MA_THRESHOLD) && (blShift < BL_SHIFT_LIMIT) &&
+ (morphType < MAXTYPES) && (GetBeatTypeCount(morphType) > R11_BC_LIM)) // Rev 1.1
+
+ return(PVC) ;
+
+ // Rule 12: If the dominant rhythm is regular and this beat is premature
+ // then classify as PVC.
+
+ if((rhythmClass == PVC) && (GetDomRhythm() == 1))
+ return(PVC) ;
+
+ // Rule 14: If the beat is regular and the dominant rhythm is regular
+ // call the beat normal.
+
+ if((rhythmClass == NORMAL) && (GetDomRhythm() == 1))
+ return(NORMAL) ;
+
+ // By this point, we know that rhythm will not help us, so we
+ // have to classify based on width and similarity to the dominant
+ // beat type.
+
+ // Rule 15: If the beat is wider than normal, wide on an
+ // absolute scale, and significantly different from the
+ // dominant beat, call it a PVC.
+
+ if((beatWidth > domWidth) && (domIndex > R15_DI_THRESHOLD) &&
+ (beatWidth >= R15_WIDTH_THRESHOLD))
+ return(PVC) ;
+
+ // Rule 16: If the beat is sufficiently narrow, call it normal.
+
+ if(beatWidth < R16_WIDTH_THRESHOLD)
+ return(NORMAL) ;
+
+ // Rule 17: If the beat isn't much wider than the dominant beat
+ // call it normal.
+
+ if(beatWidth < domWidth + R17_WIDTH_DELTA)
+ return(NORMAL) ;
+
+ // If the beat is noisy but reasonably close to dominant,
+ // call it normal.
+
+ // Rule 18: If the beat is similar to the dominant beat, call it normal.
+
+ if(domIndex < R18_DI_THRESHOLD)
+ return(NORMAL) ;
+
+ // If it's noisy don't trust the width.
+
+ // Rule 19: If the beat is noisy, we can't trust our width estimate
+ // and we have no useful rhythm information, so guess normal.
+
+ if(hfNoise > R19_HF_THRESHOLD)
+ return(NORMAL) ;
+
+ // Rule 20: By this point, we have no rhythm information, the beat
+ // isn't particularly narrow, the beat isn't particulary similar to
+ // the dominant beat, so guess a PVC.
+
+ return(PVC) ;
+
+ }
+
+
+/****************************************************************************
+* DomMonitor, monitors which beat morphology is considered to be dominant.
+* The dominant morphology is the beat morphology that has been most frequently
+* classified as normal over the course of the last 120 beats. The dominant
+* beat rhythm is classified as regular if at least 3/4 of the dominant beats
+* have been classified as regular.
+*******************************************************************************/
+
+#define DM_BUFFER_LENGTH 180
+
+int NewDom, DomRhythm ;
+int DMBeatTypes[DM_BUFFER_LENGTH], DMBeatClasses[DM_BUFFER_LENGTH] ;
+int DMBeatRhythms[DM_BUFFER_LENGTH] ;
+int DMNormCounts[8], DMBeatCounts[8], DMIrregCount = 0 ;
+
+int DomMonitor(int morphType, int rhythmClass, int beatWidth, int rr, int reset)
+ {
+ static int brIndex = 0 ;
+ int i, oldType, runCount, dom, max ;
+
+ // Fetch the type of the beat before the last beat.
+
+ i = brIndex - 2 ;
+ if(i < 0)
+ i += DM_BUFFER_LENGTH ;
+ oldType = DMBeatTypes[i] ;
+
+ // If reset flag is set, reset beat type counts and
+ // beat information buffers.
+
+ if(reset != 0)
+ {
+ for(i = 0; i < DM_BUFFER_LENGTH; ++i)
+ {
+ DMBeatTypes[i] = -1 ;
+ DMBeatClasses[i] = 0 ;
+ }
+
+ for(i = 0; i < 8; ++i)
+ {
+ DMNormCounts[i] = 0 ;
+ DMBeatCounts[i] = 0 ;
+ }
+ DMIrregCount = 0 ;
+ return(0) ;
+ }
+
+ // Once we have wrapped around, subtract old beat types from
+ // the beat counts.
+
+ if((DMBeatTypes[brIndex] != -1) && (DMBeatTypes[brIndex] != MAXTYPES))
+ {
+ --DMBeatCounts[DMBeatTypes[brIndex]] ;
+ DMNormCounts[DMBeatTypes[brIndex]] -= DMBeatClasses[brIndex] ;
+ if(DMBeatRhythms[brIndex] == UNKNOWN)
+ --DMIrregCount ;
+ }
+
+ // If this is a morphology that has been detected before, decide
+ // (for the purposes of selecting the dominant normal beattype)
+ // whether it is normal or not and update the approporiate counts.
+
+ if(morphType != 8)
+ {
+
+ // Update the buffers of previous beats and increment the
+ // count for this beat type.
+
+ DMBeatTypes[brIndex] = morphType ;
+ ++DMBeatCounts[morphType] ;
+ DMBeatRhythms[brIndex] = rhythmClass ;
+
+ // If the rhythm appears regular, update the regular rhythm
+ // count.
+
+ if(rhythmClass == UNKNOWN)
+ ++DMIrregCount ;
+
+ // Check to see how many beats of this type have occurred in
+ // a row (stop counting at six).
+
+ i = brIndex - 1 ;
+ if(i < 0) i += DM_BUFFER_LENGTH ;
+ for(runCount = 0; (DMBeatTypes[i] == morphType) && (runCount < 6); ++runCount)
+ if(--i < 0) i += DM_BUFFER_LENGTH ;
+
+ // If the rhythm is regular, the beat width is less than 130 ms, and
+ // there have been at least two in a row, consider the beat to be
+ // normal.
+
+ if((rhythmClass == NORMAL) && (beatWidth < BEAT_MS130) && (runCount >= 1))
+ {
+ DMBeatClasses[brIndex] = 1 ;
+ ++DMNormCounts[morphType] ;
+ }
+
+ // If the last beat was within the normal P-R interval for this beat,
+ // and the one before that was this beat type, assume the last beat
+ // was noise and this beat is normal.
+
+ else if(rr < ((FIDMARK-GetBeatBegin(morphType))*SAMPLE_RATE/BEAT_SAMPLE_RATE)
+ && (oldType == morphType))
+ {
+ DMBeatClasses[brIndex] = 1 ;
+ ++DMNormCounts[morphType] ;
+ }
+
+ // Otherwise assume that this is not a normal beat.
+
+ else DMBeatClasses[brIndex] = 0 ;
+ }
+
+ // If the beat does not match any of the beat types, store
+ // an indication that the beat does not match.
+
+ else
+ {
+ DMBeatClasses[brIndex] = 0 ;
+ DMBeatTypes[brIndex] = -1 ;
+ }
+
+ // Increment the index to the beginning of the circular buffers.
+
+ if(++brIndex == DM_BUFFER_LENGTH)
+ brIndex = 0 ;
+
+ // Determine which beat type has the most beats that seem
+ // normal.
+
+ dom = 0 ;
+ for(i = 1; i < 8; ++i)
+ if(DMNormCounts[i] > DMNormCounts[dom])
+ dom = i ;
+
+ max = 0 ;
+ for(i = 1; i < 8; ++i)
+ if(DMBeatCounts[i] > DMBeatCounts[max])
+ max = i ;
+
+ // If there are no normal looking beats, fall back on which beat
+ // has occurred most frequently since classification began.
+
+ if((DMNormCounts[dom] == 0) || (DMBeatCounts[max]/DMBeatCounts[dom] >= 2)) // == 0
+ dom = GetDominantType() ;
+
+ // If at least half of the most frequently occuring normal
+ // type do not seem normal, fall back on choosing the most frequently
+ // occurring type since classification began.
+
+ else if(DMBeatCounts[dom]/DMNormCounts[dom] >= 2)
+ dom = GetDominantType() ;
+
+ // If there is any beat type that has been classfied as normal,
+ // but at least 10 don't seem normal, reclassify it to UNKNOWN.
+
+ for(i = 0; i < 8; ++i)
+ if((DMBeatCounts[i] > 10) && (DMNormCounts[i] == 0) && (i != dom)
+ && (GetBeatClass(i) == NORMAL))
+ SetBeatClass(i,UNKNOWN) ;
+
+ // Save the dominant type in a global variable so that it is
+ // accessable for debugging.
+
+ NewDom = dom ;
+ return(dom) ;
+ }
+
+int GetNewDominantType(void)
+ {
+ return(NewDom) ;
+ }
+
+int GetDomRhythm(void)
+ {
+ if(DMIrregCount > IRREG_RR_LIMIT)
+ return(0) ;
+ else return(1) ;
+ }
+
+
+void AdjustDomData(int oldType, int newType)
+ {
+ int i ;
+
+ for(i = 0; i < DM_BUFFER_LENGTH; ++i)
+ {
+ if(DMBeatTypes[i] == oldType)
+ DMBeatTypes[i] = newType ;
+ }
+
+ if(newType != MAXTYPES)
+ {
+ DMNormCounts[newType] = DMNormCounts[oldType] ;
+ DMBeatCounts[newType] = DMBeatCounts[oldType] ;
+ }
+
+ DMNormCounts[oldType] = DMBeatCounts[oldType] = 0 ;
+
+ }
+
+void CombineDomData(int oldType, int newType)
+ {
+ int i ;
+
+ for(i = 0; i < DM_BUFFER_LENGTH; ++i)
+ {
+ if(DMBeatTypes[i] == oldType)
+ DMBeatTypes[i] = newType ;
+ }
+
+ if(newType != MAXTYPES)
+ {
+ DMNormCounts[newType] += DMNormCounts[oldType] ;
+ DMBeatCounts[newType] += DMBeatCounts[oldType] ;
+ }
+
+ DMNormCounts[oldType] = DMBeatCounts[oldType] = 0 ;
+
+ }
+
+/***********************************************************************
+ GetRunCount() checks how many of the present beat type have occurred
+ in a row.
+***********************************************************************/
+
+int GetRunCount()
+ {
+ int i ;
+ for(i = 1; (i < 8) && (RecentTypes[0] == RecentTypes[i]); ++i) ;
+ return(i) ;
+ }
+
+
+
+
+
+
+
+
+
+
+
diff -r 1de72bdab0ef -r 8e2fbe131b53 ecgcodes.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/ecgcodes.h Sun Jun 28 03:06:00 2015 +0000 @@ -0,0 +1,52 @@ +/* file: ecgcodes.h T. Baker and G. Moody June 1981 + Last revised: 19 March 1992 dblib 7.0 +ECG annotation codes + +Copyright (C) Massachusetts Institute of Technology 1992. All rights reserved. +*/ + +#ifndef db_ECGCODES_H /* avoid multiple definitions */ +#define db_ECGCODES_H + +#define NOTQRS 0 /* not-QRS (not a getann/putann code) */ +#define NORMAL 1 /* normal beat */ +#define LBBB 2 /* left bundle branch block beat */ +#define RBBB 3 /* right bundle branch block beat */ +#define ABERR 4 /* aberrated atrial premature beat */ +#define PVC 5 /* premature ventricular contraction */ +#define FUSION 6 /* fusion of ventricular and normal beat */ +#define NPC 7 /* nodal (junctional) premature beat */ +#define APC 8 /* atrial premature contraction */ +#define SVPB 9 /* premature or ectopic supraventricular beat */ +#define VESC 10 /* ventricular escape beat */ +#define NESC 11 /* nodal (junctional) escape beat */ +#define PACE 12 /* paced beat */ +#define UNKNOWN 13 /* unclassifiable beat */ +#define NOISE 14 /* signal quality change */ +#define ARFCT 16 /* isolated QRS-like artifact */ +#define STCH 18 /* ST change */ +#define TCH 19 /* T-wave change */ +#define SYSTOLE 20 /* systole */ +#define DIASTOLE 21 /* diastole */ +#define NOTE 22 /* comment annotation */ +#define MEASURE 23 /* measurement annotation */ +#define BBB 25 /* left or right bundle branch block */ +#define PACESP 26 /* non-conducted pacer spike */ +#define RHYTHM 28 /* rhythm change */ +#define LEARN 30 /* learning */ +#define FLWAV 31 /* ventricular flutter wave */ +#define VFON 32 /* start of ventricular flutter/fibrillation */ +#define VFOFF 33 /* end of ventricular flutter/fibrillation */ +#define AESC 34 /* atrial escape beat */ +#define SVESC 35 /* supraventricular escape beat */ +#define NAPC 37 /* non-conducted P-wave (blocked APB) */ +#define PFUS 38 /* fusion of paced and normal beat */ +#define PQ 39 /* PQ junction (beginning of QRS) */ +#define JPT 40 /* J point (end of QRS) */ +#define RONT 41 /* R-on-T premature ventricular contraction */ + +/* ... annotation codes between RONT+1 and ACMAX inclusive are user-defined */ + +#define ACMAX 49 /* value of largest valid annot code (must be < 50) */ + +#endif
diff -r 1de72bdab0ef -r 8e2fbe131b53 filter.cpp
--- a/filter.cpp Sun May 17 00:27:16 2015 +0000
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,70 +0,0 @@
-// Bandpass the 1000Hz sampling rate to only allow signals that could be heart rate frequencies
-// https://www.keil.com/pack/doc/CMSIS/DSP/html/arm_fir_example_f32_8c-example.html
-
-#define TEST_LENGTH_SAMPLES 320
-#define BLOCK_SIZE 32
-#define NUM_TAPS 51
-#define PER_MINUTE 60
-
-#define MIN_HEARTRATE (30.0/PER_MINUTE)
-#define MAX_HEARTRATE (120.0/PER_MINUTE)
-
-int blockSize = BLOCK_SIZE;
-int numBlocks = TEST_LENGTH_SAMPLES/BLOCK_SIZE;
-
-static float testOutput[TEST_LENGTH_SAMPLES];
-static float firStateF32[BLOCK_SIZE + NUM_TAPS - 1];
-
-// TODO: USE CORRECT FREQUENCY FOR HEART RATE.
-// 51-tap FIR filter @1000Hz sampling rate, 30Hz-120Hz pass, 50dB attenuation: http://www.arc.id.au/FilterDesign.html
-float firCoeffs32[NUM_TAPS] = {0.000707,
-0.000317,
--0.000103,
-0.000151,
-0.001621,
-0.004056,
-0.006124,
-0.005903,
-0.002063,
--0.004881,
--0.012196,
--0.016049,
--0.013847,
--0.006569,
-0.000551,
--0.000000,
--0.013828,
--0.040240,
--0.070297,
--0.089104,
--0.082037,
--0.042270,
-0.024193,
-0.098977,
-0.157748,
-0.180000,
-0.157748,
-0.098977,
-0.024193,
--0.042270,
--0.082037,
--0.089104,
--0.070297,
--0.040240,
--0.013828,
--0.000000,
-0.000551,
--0.006569,
--0.013847,
--0.016049,
--0.012196,
--0.004881,
-0.002063,
-0.005903,
-0.006124,
-0.004056,
-0.001621,
-0.000151,
--0.000103,
-0.000317,
-0.000707};
\ No newline at end of file
diff -r 1de72bdab0ef -r 8e2fbe131b53 main.cpp
--- a/main.cpp Sun May 17 00:27:16 2015 +0000
+++ b/main.cpp Sun Jun 28 03:06:00 2015 +0000
@@ -15,6 +15,9 @@
*/
#include "mbed.h"
+#include <qrsdet.h>
+#include "queue.h"
+
#include "BLEDevice.h"
#include "HeartRateService.h"
#include "DeviceInformationService.h"
@@ -25,50 +28,90 @@
* interval.*/
#define UPDATE_PARAMS_FOR_LONGER_CONNECTION_INTERVAL 0
+// SAMPLE_RATE is defined in qrsdet.h
+#define MS_PER_SAMPLE (1000/SAMPLE_RATE)
+
extern void setup_sampler(void (*function)(uint32_t));
-
+extern int BeatDetectAndClassify(int ecgSample, int *beatType, int *beatMatch); // defined in bdac.cpp
+extern void ResetBDAC(void);
+
BLEDevice ble;
DigitalOut led1(LED1);
+PwmOut LRA(P0_15);
+InterruptIn button1(P0_17); // button 1
+InterruptIn button2(P0_18); // button 2
+InterruptIn button3(P0_19); // button 3
+InterruptIn button4(P0_20); // button 4
-volatile uint8_t hrmCounter = 100; // init HRM to 100bps
+
+//Serial pc (P0_9, P0_11); // TX, RX
+
+volatile uint8_t hrmCounter = 128;
+volatile uint16_t ecgSample = 0;
+volatile int sampleCounter=0;
+
+#define ITEMS_IN_QUEUE 50
+Queue ecgQueue ( sizeof(int), ITEMS_IN_QUEUE ); // queue of hrmCounter values
+bool itemAddedToQueue = true;
+
const static char DEVICE_NAME[] = "HRM1";
static const uint16_t uuid16_list[] = {GattService::UUID_HEART_RATE_SERVICE,
GattService::UUID_DEVICE_INFORMATION_SERVICE
};
-static volatile bool triggerSensorPolling = false;
void disconnectionCallback(Gap::Handle_t handle, Gap::DisconnectionReason_t reason)
{
ble.startAdvertising(); // restart advertising
}
-//void periodicCallback(void)
-//{
-// led1 = !led1; /* Do blinky on LED1 while we're waiting for BLE events */
-//
-// /* Note that the periodicCallback() executes in interrupt context, so it is safer to do
-// * heavy-weight sensor polling from the main thread. */
-// triggerSensorPolling = true;
-//}
+
+void power0() {
+ LRA.write(0.0);
+}
+void power1() {
+ LRA.write(0.55);
+}
+void power2() {
+ LRA.write(0.85);
+}
+void power3() {
+ LRA.write(1.0);
+}
void sampler_callback(uint32_t value)
{
- hrmCounter = (uint8_t) (value);
-// led1 = !led1;
- triggerSensorPolling = true;
+ //hrmCounter = (uint8_t) ((value+127) % 255);
+ ecgSample = (uint16_t) value;
+ itemAddedToQueue = ecgQueue.PutIrq((void*)&ecgSample);
}
+
int main(void)
{
+ uint16_t value;
+ int beatType;
+ int beatMatch;
+ int samplesSinceLastR;
+ int lastSamplesSinceLastR=0;
+ int sampleOffset=0;
+ char* beatName;
+
led1 = 1;
+
// Ticker ticker;
//ticker.attach(periodicCallback, 0.1); // blink LED every second
- printf("Executing code...\n");
-
- setup_sampler(&sampler_callback);
+ button1.rise(&power0);
+ button2.rise(&power1);
+ button3.rise(&power2);
+ button4.rise(&power3);
+ LRA.period_us(50); // 50 uS -> 20 Khz (needs to be between 10Khz and 250Khz)
+ LRA.write(0.0); // Off. < 50% = 0ff
+
+ printf("Initializing BLE...\r\n");
+
ble.init();
ble.onDisconnection(disconnectionCallback);
@@ -80,6 +123,7 @@
DeviceInformationService deviceInfo(ble, "ARM", "Model1", "SN1", "hw-rev1", "fw-rev1", "soft-rev1");
/* Setup advertising. */
+ printf("Setting up advertising...\r\n");
ble.accumulateAdvertisingPayload(GapAdvertisingData::BREDR_NOT_SUPPORTED | GapAdvertisingData::LE_GENERAL_DISCOVERABLE);
ble.accumulateAdvertisingPayload(GapAdvertisingData::COMPLETE_LIST_16BIT_SERVICE_IDS, (uint8_t *)uuid16_list, sizeof(uuid16_list));
ble.accumulateAdvertisingPayload(GapAdvertisingData::GENERIC_HEART_RATE_SENSOR);
@@ -88,27 +132,49 @@
ble.setAdvertisingInterval(1000);
ble.startAdvertising();
+ setup_sampler(&sampler_callback);
+
+ ResetBDAC(); // reset the beat detector and classifier
+
+ printf("Starting...\r\n");
+
// infinite loop
while (1) {
// check for trigger from periodicCallback()
- if (triggerSensorPolling && ble.getGapState().connected) {
- triggerSensorPolling = false;
-
- led1 = !led1;
-
- printf("v=%d\n", hrmCounter);
-// // Do blocking calls or whatever is necessary for sensor polling.
-// // In our case, we simply update the HRM measurement.
-// hrmCounter++;
-//
-// // 100 <= HRM bps <=175
-// if (hrmCounter == 175) {
-// hrmCounter = 100;
-// }
-
- // update bps
- hrService.updateHeartRate(hrmCounter);
+ if (ble.getGapState().connected) {
+ while (ecgQueue.GetNumberOfItems()>0) {
+ ecgQueue.Get(&value);
+ samplesSinceLastR = BeatDetectAndClassify(value, &beatType, &beatMatch);
+ if (samplesSinceLastR != 0) {
+ printf("[C] samplesSinceLastR=%d, type=%d\r\n", value, beatType);
+ }
+ //hrService.updateHeartRate(value);
+ led1 = !led1;
+ }
} else {
+ if (!itemAddedToQueue) {
+ printf("Queue overflow.\n\r");
+ itemAddedToQueue = true;
+ }
+ while (ecgQueue.GetNumberOfItems()>0) {
+ sampleCounter++;
+ ecgQueue.Get(&value);
+ samplesSinceLastR = BeatDetectAndClassify(value, &beatType, &beatMatch);
+ if (samplesSinceLastR != 0) {
+ sampleOffset = lastSamplesSinceLastR - samplesSinceLastR;
+
+ if (beatType == 1) {
+ beatName = "Normal";
+ } else if (beatType == 5) {
+ beatName = "Ectopic";
+ } else {
+ beatName = "Unknown";
+ }
+ printf("[NC] interval_ms=%d, bpm=%d, samplesSinceLastR=%d (%s)\r\n", ((sampleCounter-sampleOffset)*MS_PER_SAMPLE), (60000/((sampleCounter-sampleOffset)*MS_PER_SAMPLE)), value, beatName);
+ sampleCounter=0;
+ lastSamplesSinceLastR = samplesSinceLastR;
+ }
+ }
ble.waitForEvent(); // low power wait for event
}
}
diff -r 1de72bdab0ef -r 8e2fbe131b53 match.cpp
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/match.cpp Sun Jun 28 03:06:00 2015 +0000
@@ -0,0 +1,856 @@
+/*****************************************************************************
+FILE: match.cpp
+AUTHOR: Patrick S. Hamilton
+REVISED: 5/13/2002
+ ___________________________________________________________________________
+
+match.cpp: Match beats to previous beats.
+Copywrite (C) 2001 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).
+ __________________________________________________________________________
+
+Match.cpp contains functions for managing template matching of beats and
+managing of feature data associated with each beat type. These
+functions are called functions in classify.cpp. Beats are matched to
+previoiusly detected beats types based on how well they match point by point
+in a MATCH_LENGTH region centered on FIDMARK (R-wave location). The following
+is a list of functions that are available for calling by classify.
+
+ ResetMatch -- Resets global variables used in template matching.
+ CompareBeats -- Measures the difference between two beats with
+ beats scaled to produce the best match.
+ CompareBeats2 -- Measures the difference between two beats without
+ beat scaling.
+ NewBeatType -- Start a new beat type with the present beat.
+ BestMorphMatch -- Finds the beat template that best matches a new beat.
+ UpdateBeatType -- Updates existing beat template and associated features
+ based on a new beat.
+ GetDominantType -- Returns the NORMAL beat type that has occorred most often.
+ ClearLastNewType -- Removes the last new beat type from the possible beat
+ types.
+ DomCompare -- Compares the template for a given beat type to the template
+ of the dominant normal beat type.
+ DomCompare2 -- Compares a given beat template to the templat of the
+ dominant normal beat type.
+
+ PostClassify -- Classifies beats based on preceding and following beats
+ and R-to-R intervals.
+
+ ResetPostClassify -- Resets variables used for post classification.
+
+ CheckPostClass -- Check type classification based on last eight post
+ classifications.
+
+ CheckPCClass -- Check post beat rhythm classification for the last eight
+ beats.
+
+A number of simple functions allow access to beat features while maintaining
+some level of encapsulation:
+
+ GetTypesCount -- Returns number of beat types that have been detected.
+ GetBeatTypeCount -- Returns the number of beats of a given type
+ that have been detected.
+ GetBeatWidth -- Returns the width estimate for a given beat type.
+ SetBeatClass -- Associates a beat classification with a beat type.
+ GetBeatBegin -- Returns the beginning point for a given beat type.
+ GetBeatEnd -- Returns the ending point for a given beat type.
+
+******************************************************************************/
+//#include <stdlib.h>
+//#include <stdio.h>
+#include <mbed.h>
+#include "ecgcodes.h"
+
+#include "bdac.h"
+#define MATCH_LENGTH BEAT_MS300 // Number of points used for beat matching.
+#define MATCH_LIMIT 1.2 // Match limit used testing whether two
+ // beat types might be combined.
+#define COMBINE_LIMIT 0.8 // Limit used for deciding whether two types
+ // can be combined.
+
+#define MATCH_START (FIDMARK-(MATCH_LENGTH/2)) // Starting point for beat matching
+#define MATCH_END (FIDMARK+(MATCH_LENGTH/2)) // End point for beat matching.
+#define MAXPREV 8 // Number of preceeding beats used as beat features.
+#define MAX_SHIFT BEAT_MS40
+
+// Local prototypes.
+
+int NoiseCheck(int *beat) ;
+double CompareBeats(int *beat1, int *beat2, int *shiftAdj) ;
+double CompareBeats2(int *beat1, int *beat2, int *shiftAdj) ;
+void UpdateBeat(int *aveBeat, int *newBeat, int shift) ;
+void BeatCopy(int srcBeat, int destBeat) ;
+int MinimumBeatVariation(int type) ;
+
+// External prototypes.
+
+void AnalyzeBeat(int *beat, int *onset, int *offset, int *isoLevel,
+ int *beatBegin, int *beatEnd, int *amp) ;
+void AdjustDomData(int oldType, int newType) ;
+void CombineDomData(int oldType, int newType) ;
+
+// Global variables.
+
+int BeatTemplates[MAXTYPES][BEATLGTH] ;
+int BeatCounts[MAXTYPES] ;
+int BeatWidths[MAXTYPES] ;
+int BeatClassifications[MAXTYPES] ;
+int BeatBegins[MAXTYPES] ;
+int BeatEnds[MAXTYPES] ;
+int BeatsSinceLastMatch[MAXTYPES] ;
+int BeatAmps[MAXTYPES] ;
+int BeatCenters[MAXTYPES] ;
+double MIs[MAXTYPES][8] ;
+
+// Need access to these in postclas.cpp when beat types are combined
+// and moved.
+
+extern int PostClass[MAXTYPES][8] ;
+extern int PCRhythm[MAXTYPES][8] ;
+
+int TypeCount = 0 ;
+
+/***************************************************************************
+ResetMatch() resets static variables involved with template matching.
+****************************************************************************/
+
+void ResetMatch(void)
+ {
+ int i, j ;
+ TypeCount = 0 ;
+ for(i = 0; i < MAXTYPES; ++i)
+ {
+ BeatCounts[i] = 0 ;
+ BeatClassifications[i] = UNKNOWN ;
+ for(j = 0; j < 8; ++j)
+ {
+ MIs[i][j] = 0 ;
+ }
+ }
+ }
+
+/**************************************************************************
+ CompareBeats() takes two beat buffers and compares how well they match
+ point-by-point. Beat2 is shifted and scaled to produce the closest
+ possible match. The metric returned is the sum of the absolute
+ differences between beats divided by the amplitude of the beats. The
+ shift used for the match is returned via the pointer *shiftAdj.
+***************************************************************************/
+
+#define MATCH_START (FIDMARK-(MATCH_LENGTH/2))
+#define MATCH_END (FIDMARK+(MATCH_LENGTH/2))
+
+double CompareBeats(int *beat1, int *beat2, int *shiftAdj)
+ {
+ int i, max, min, magSum, shift ;
+ long beatDiff, meanDiff, minDiff, minShift ;
+ double metric, scaleFactor, tempD ;
+
+ // Calculate the magnitude of each beat.
+
+ max = min = beat1[MATCH_START] ;
+ for(i = MATCH_START+1; i < MATCH_END; ++i)
+ if(beat1[i] > max)
+ max = beat1[i] ;
+ else if(beat1[i] < min)
+ min = beat1[i] ;
+
+ magSum = max - min ;
+
+ i = MATCH_START ;
+ max = min = beat2[i] ;
+ for(i = MATCH_START+1; i < MATCH_END; ++i)
+ if(beat2[i] > max)
+ max = beat2[i] ;
+ else if(beat2[i] < min)
+ min = beat2[i] ;
+
+ // magSum += max - min ;
+ scaleFactor = magSum ;
+ scaleFactor /= max-min ;
+ magSum *= 2 ;
+
+ // Calculate the sum of the point-by-point
+ // absolute differences for five possible shifts.
+
+ for(shift = -MAX_SHIFT; shift <= MAX_SHIFT; ++shift)
+ {
+ for(i = FIDMARK-(MATCH_LENGTH>>1), meanDiff = 0;
+ i < FIDMARK + (MATCH_LENGTH>>1); ++i)
+ {
+ tempD = beat2[i+shift] ;
+ tempD *= scaleFactor ;
+ meanDiff += beat1[i]- tempD ; // beat2[i+shift] ;
+ }
+ meanDiff /= MATCH_LENGTH ;
+
+ for(i = FIDMARK-(MATCH_LENGTH>>1), beatDiff = 0;
+ i < FIDMARK + (MATCH_LENGTH>>1); ++i)
+ {
+ tempD = beat2[i+shift] ;
+ tempD *= scaleFactor ;
+ beatDiff += abs(beat1[i] - meanDiff- tempD) ; // beat2[i+shift] ) ;
+ }
+
+
+ if(shift == -MAX_SHIFT)
+ {
+ minDiff = beatDiff ;
+ minShift = -MAX_SHIFT ;
+ }
+ else if(beatDiff < minDiff)
+ {
+ minDiff = beatDiff ;
+ minShift = shift ;
+ }
+ }
+
+ metric = minDiff ;
+ *shiftAdj = minShift ;
+ metric /= magSum ;
+
+ // Metric scales inversely with match length.
+ // algorithm was originally tuned with a match
+ // length of 30.
+
+ metric *= 30 ;
+ metric /= MATCH_LENGTH ;
+ return(metric) ;
+ }
+
+/***************************************************************************
+ CompareBeats2 is nearly the same as CompareBeats above, but beat2 is
+ not scaled before calculating the match metric. The match metric is
+ then the sum of the absolute differences divided by the average amplitude
+ of the two beats.
+****************************************************************************/
+
+double CompareBeats2(int *beat1, int *beat2, int *shiftAdj)
+ {
+ int i, max, min, shift ;
+ int mag1, mag2 ;
+ long beatDiff, meanDiff, minDiff, minShift ;
+ double metric ;
+
+ // Calculate the magnitude of each beat.
+
+ max = min = beat1[MATCH_START] ;
+ for(i = MATCH_START+1; i < MATCH_END; ++i)
+ if(beat1[i] > max)
+ max = beat1[i] ;
+ else if(beat1[i] < min)
+ min = beat1[i] ;
+
+ mag1 = max - min ;
+
+ i = MATCH_START ;
+ max = min = beat2[i] ;
+ for(i = MATCH_START+1; i < MATCH_END; ++i)
+ if(beat2[i] > max)
+ max = beat2[i] ;
+ else if(beat2[i] < min)
+ min = beat2[i] ;
+
+ mag2 = max-min ;
+
+ // Calculate the sum of the point-by-point
+ // absolute differences for five possible shifts.
+
+ for(shift = -MAX_SHIFT; shift <= MAX_SHIFT; ++shift)
+ {
+ for(i = FIDMARK-(MATCH_LENGTH>>1), meanDiff = 0;
+ i < FIDMARK + (MATCH_LENGTH>>1); ++i)
+ meanDiff += beat1[i]- beat2[i+shift] ;
+ meanDiff /= MATCH_LENGTH ;
+
+ for(i = FIDMARK-(MATCH_LENGTH>>1), beatDiff = 0;
+ i < FIDMARK + (MATCH_LENGTH>>1); ++i)
+ beatDiff += abs(beat1[i] - meanDiff- beat2[i+shift]) ; ;
+
+ if(shift == -MAX_SHIFT)
+ {
+ minDiff = beatDiff ;
+ minShift = -MAX_SHIFT ;
+ }
+ else if(beatDiff < minDiff)
+ {
+ minDiff = beatDiff ;
+ minShift = shift ;
+ }
+ }
+
+ metric = minDiff ;
+ *shiftAdj = minShift ;
+ metric /= (mag1+mag2) ;
+
+ // Metric scales inversely with match length.
+ // algorithm was originally tuned with a match
+ // length of 30.
+
+ metric *= 30 ;
+ metric /= MATCH_LENGTH ;
+
+ return(metric) ;
+ }
+
+/************************************************************************
+UpdateBeat() averages a new beat into an average beat template by adding
+1/8th of the new beat to 7/8ths of the average beat.
+*************************************************************************/
+
+void UpdateBeat(int *aveBeat, int *newBeat, int shift)
+ {
+ int i ;
+ long tempLong ;
+
+ for(i = 0; i < BEATLGTH; ++i)
+ {
+ if((i+shift >= 0) && (i+shift < BEATLGTH))
+ {
+ tempLong = aveBeat[i] ;
+ tempLong *= 7 ;
+ tempLong += newBeat[i+shift] ;
+ tempLong >>= 3 ;
+ aveBeat[i] = tempLong ;
+ }
+ }
+ }
+
+/*******************************************************
+ GetTypesCount returns the number of types that have
+ been detected.
+*******************************************************/
+
+int GetTypesCount(void)
+ {
+ return(TypeCount) ;
+ }
+
+/********************************************************
+ GetBeatTypeCount returns the number of beats of a
+ a particular type have been detected.
+********************************************************/
+
+int GetBeatTypeCount(int type)
+ {
+ return(BeatCounts[type]) ;
+ }
+
+/*******************************************************
+ GetBeatWidth returns the QRS width estimate for
+ a given type of beat.
+*******************************************************/
+int GetBeatWidth(int type)
+ {
+ return(BeatWidths[type]) ;
+ }
+
+/*******************************************************
+ GetBeatCenter returns the point between the onset and
+ offset of a beat.
+********************************************************/
+
+int GetBeatCenter(int type)
+ {
+ return(BeatCenters[type]) ;
+ }
+
+/*******************************************************
+ GetBeatClass returns the present classification for
+ a given beat type (NORMAL, PVC, or UNKNOWN).
+********************************************************/
+
+int GetBeatClass(int type)
+ {
+ if(type == MAXTYPES)
+ return(UNKNOWN) ;
+ return(BeatClassifications[type]) ;
+ }
+
+/******************************************************
+ SetBeatClass sets up a beat classifation for a
+ given type.
+******************************************************/
+
+void SetBeatClass(int type, int beatClass)
+ {
+ BeatClassifications[type] = beatClass ;
+ }
+
+/******************************************************************************
+ NewBeatType starts a new beat type by storing the new beat and its
+ features as the next available beat type.
+******************************************************************************/
+
+int NewBeatType(int *newBeat )
+ {
+ int i, onset, offset, isoLevel, beatBegin, beatEnd ;
+ int mcType, amp ;
+
+ // Update count of beats since each template was matched.
+
+ for(i = 0; i < TypeCount; ++i)
+ ++BeatsSinceLastMatch[i] ;
+
+ if(TypeCount < MAXTYPES)
+ {
+ for(i = 0; i < BEATLGTH; ++i)
+ BeatTemplates[TypeCount][i] = newBeat[i] ;
+
+ BeatCounts[TypeCount] = 1 ;
+ BeatClassifications[TypeCount] = UNKNOWN ;
+ AnalyzeBeat(&BeatTemplates[TypeCount][0],&onset,&offset, &isoLevel,
+ &beatBegin, &beatEnd, &) ;
+ BeatWidths[TypeCount] = offset-onset ;
+ BeatCenters[TypeCount] = (offset+onset)/2 ;
+ BeatBegins[TypeCount] = beatBegin ;
+ BeatEnds[TypeCount] = beatEnd ;
+ BeatAmps[TypeCount] = amp ;
+
+ BeatsSinceLastMatch[TypeCount] = 0 ;
+
+ ++TypeCount ;
+ return(TypeCount-1) ;
+ }
+
+ // If we have used all the template space, replace the beat
+ // that has occurred the fewest number of times.
+
+ else
+ {
+ // Find the template with the fewest occurances,
+ // that hasn't been matched in at least 500 beats.
+
+ mcType = -1 ;
+
+ if(mcType == -1)
+ {
+ mcType = 0 ;
+ for(i = 1; i < MAXTYPES; ++i)
+ if(BeatCounts[i] < BeatCounts[mcType])
+ mcType = i ;
+ else if(BeatCounts[i] == BeatCounts[mcType])
+ {
+ if(BeatsSinceLastMatch[i] > BeatsSinceLastMatch[mcType])
+ mcType = i ;
+ }
+ }
+
+ // Adjust dominant beat monitor data.
+
+ AdjustDomData(mcType,MAXTYPES) ;
+
+ // Substitute this beat.
+
+ for(i = 0; i < BEATLGTH; ++i)
+ BeatTemplates[mcType][i] = newBeat[i] ;
+
+ BeatCounts[mcType] = 1 ;
+ BeatClassifications[mcType] = UNKNOWN ;
+ AnalyzeBeat(&BeatTemplates[mcType][0],&onset,&offset, &isoLevel,
+ &beatBegin, &beatEnd, &) ;
+ BeatWidths[mcType] = offset-onset ;
+ BeatCenters[mcType] = (offset+onset)/2 ;
+ BeatBegins[mcType] = beatBegin ;
+ BeatEnds[mcType] = beatEnd ;
+ BeatsSinceLastMatch[mcType] = 0 ;
+ BeatAmps[mcType] = amp ;
+ return(mcType) ;
+ }
+ }
+
+/***************************************************************************
+ BestMorphMatch tests a new beat against all available beat types and
+ returns (via pointers) the existing type that best matches, the match
+ metric for that type, and the shift used for that match.
+***************************************************************************/
+
+void BestMorphMatch(int *newBeat,int *matchType,double *matchIndex, double *mi2,
+ int *shiftAdj)
+ {
+ int type, i, bestMatch, nextBest, minShift, shift, temp ;
+ int bestShift2, nextShift2 ;
+ double bestDiff2, nextDiff2;
+ double beatDiff, minDiff, nextDiff=10000 ;
+
+ if(TypeCount == 0)
+ {
+ *matchType = 0 ;
+ *matchIndex = 1000 ; // Make sure there is no match so a new beat is
+ *shiftAdj = 0 ; // created.
+ return ;
+ }
+
+ // Compare the new beat to all type beat
+ // types that have been saved.
+
+ for(type = 0; type < TypeCount; ++type)
+ {
+ beatDiff = CompareBeats(&BeatTemplates[type][0],newBeat,&shift) ;
+ if(type == 0)
+ {
+ bestMatch = 0 ;
+ minDiff = beatDiff ;
+ minShift = shift ;
+ }
+ else if(beatDiff < minDiff)
+ {
+ nextBest = bestMatch ;
+ nextDiff = minDiff ;
+ bestMatch = type ;
+ minDiff = beatDiff ;
+ minShift = shift ;
+ }
+ else if((TypeCount > 1) && (type == 1))
+ {
+ nextBest = type ;
+ nextDiff = beatDiff ;
+ }
+ else if(beatDiff < nextDiff)
+ {
+ nextBest = type ;
+ nextDiff = beatDiff ;
+ }
+ }
+
+ // If this beat was close to two different
+ // templates, see if the templates which template
+ // is the best match when no scaling is used.
+ // Then check whether the two close types can be combined.
+
+ if((minDiff < MATCH_LIMIT) && (nextDiff < MATCH_LIMIT) && (TypeCount > 1))
+ {
+ // Compare without scaling.
+
+ bestDiff2 = CompareBeats2(&BeatTemplates[bestMatch][0],newBeat,&bestShift2) ;
+ nextDiff2 = CompareBeats2(&BeatTemplates[nextBest][0],newBeat,&nextShift2) ;
+ if(nextDiff2 < bestDiff2)
+ {
+ temp = bestMatch ;
+ bestMatch = nextBest ;
+ nextBest = temp ;
+ temp = minDiff ;
+ minDiff = nextDiff ;
+ nextDiff = temp ;
+ minShift = nextShift2 ;
+ *mi2 = bestDiff2 ;
+ }
+ else *mi2 = nextDiff2 ;
+
+ beatDiff = CompareBeats(&BeatTemplates[bestMatch][0],&BeatTemplates[nextBest][0],&shift) ;
+
+ if((beatDiff < COMBINE_LIMIT) &&
+ ((*mi2 < 1.0) || (!MinimumBeatVariation(nextBest))))
+ {
+
+ // Combine beats into bestMatch
+
+ if(bestMatch < nextBest)
+ {
+ for(i = 0; i < BEATLGTH; ++i)
+ {
+ if((i+shift > 0) && (i + shift < BEATLGTH))
+ {
+ BeatTemplates[bestMatch][i] += BeatTemplates[nextBest][i+shift] ;
+ BeatTemplates[bestMatch][i] >>= 1 ;
+ }
+ }
+
+ if((BeatClassifications[bestMatch] == NORMAL) || (BeatClassifications[nextBest] == NORMAL))
+ BeatClassifications[bestMatch] = NORMAL ;
+ else if((BeatClassifications[bestMatch] == PVC) || (BeatClassifications[nextBest] == PVC))
+ BeatClassifications[bestMatch] = PVC ;
+
+ BeatCounts[bestMatch] += BeatCounts[nextBest] ;
+
+ CombineDomData(nextBest,bestMatch) ;
+
+ // Shift other templates over.
+
+ for(type = nextBest; type < TypeCount-1; ++type)
+ BeatCopy(type+1,type) ;
+
+ }
+
+ // Otherwise combine beats it nextBest.
+
+ else
+ {
+ for(i = 0; i < BEATLGTH; ++i)
+ {
+ BeatTemplates[nextBest][i] += BeatTemplates[bestMatch][i] ;
+ BeatTemplates[nextBest][i] >>= 1 ;
+ }
+
+ if((BeatClassifications[bestMatch] == NORMAL) || (BeatClassifications[nextBest] == NORMAL))
+ BeatClassifications[nextBest] = NORMAL ;
+ else if((BeatClassifications[bestMatch] == PVC) || (BeatClassifications[nextBest] == PVC))
+ BeatClassifications[nextBest] = PVC ;
+
+ BeatCounts[nextBest] += BeatCounts[bestMatch] ;
+
+ CombineDomData(bestMatch,nextBest) ;
+
+ // Shift other templates over.
+
+ for(type = bestMatch; type < TypeCount-1; ++type)
+ BeatCopy(type+1,type) ;
+
+
+ bestMatch = nextBest ;
+ }
+ --TypeCount ;
+ BeatClassifications[TypeCount] = UNKNOWN ;
+ }
+ }
+ *mi2 = CompareBeats2(&BeatTemplates[bestMatch][0],newBeat,&bestShift2) ;
+ *matchType = bestMatch ;
+ *matchIndex = minDiff ;
+ *shiftAdj = minShift ;
+ }
+
+/***************************************************************************
+ UpdateBeatType updates the beat template and features of a given beat type
+ using a new beat.
+***************************************************************************/
+
+void UpdateBeatType(int matchType,int *newBeat, double mi2,
+ int shiftAdj)
+ {
+ int i,onset,offset, isoLevel, beatBegin, beatEnd ;
+ int amp ;
+
+ // Update beats since templates were matched.
+
+ for(i = 0; i < TypeCount; ++i)
+ {
+ if(i != matchType)
+ ++BeatsSinceLastMatch[i] ;
+ else BeatsSinceLastMatch[i] = 0 ;
+ }
+
+ // If this is only the second beat, average it with the existing
+ // template.
+
+ if(BeatCounts[matchType] == 1)
+ for(i = 0; i < BEATLGTH; ++i)
+ {
+ if((i+shiftAdj >= 0) && (i+shiftAdj < BEATLGTH))
+ BeatTemplates[matchType][i] = (BeatTemplates[matchType][i] + newBeat[i+shiftAdj])>>1 ;
+ }
+
+ // Otherwise do a normal update.
+
+ else
+ UpdateBeat(&BeatTemplates[matchType][0], newBeat, shiftAdj) ;
+
+ // Determine beat features for the new average beat.
+
+ AnalyzeBeat(&BeatTemplates[matchType][0],&onset,&offset,&isoLevel,
+ &beatBegin, &beatEnd, &) ;
+
+ BeatWidths[matchType] = offset-onset ;
+ BeatCenters[matchType] = (offset+onset)/2 ;
+ BeatBegins[matchType] = beatBegin ;
+ BeatEnds[matchType] = beatEnd ;
+ BeatAmps[matchType] = amp ;
+
+ ++BeatCounts[matchType] ;
+
+ for(i = MAXPREV-1; i > 0; --i)
+ MIs[matchType][i] = MIs[matchType][i-1] ;
+ MIs[matchType][0] = mi2 ;
+
+ }
+
+
+/****************************************************************************
+ GetDominantType returns the NORMAL beat type that has occurred most
+ frequently.
+****************************************************************************/
+
+int GetDominantType(void)
+ {
+ int maxCount = 0, maxType = -1 ;
+ int type, totalCount ;
+
+ for(type = 0; type < MAXTYPES; ++type)
+ {
+ if((BeatClassifications[type] == NORMAL) && (BeatCounts[type] > maxCount))
+ {
+ maxType = type ;
+ maxCount = BeatCounts[type] ;
+ }
+ }
+
+ // If no normals are found and at least 300 beats have occurred, just use
+ // the most frequently occurring beat.
+
+ if(maxType == -1)
+ {
+ for(type = 0, totalCount = 0; type < TypeCount; ++type)
+ totalCount += BeatCounts[type] ;
+ if(totalCount > 300)
+ for(type = 0; type < TypeCount; ++type)
+ if(BeatCounts[type] > maxCount)
+ {
+ maxType = type ;
+ maxCount = BeatCounts[type] ;
+ }
+ }
+
+ return(maxType) ;
+ }
+
+
+/***********************************************************************
+ ClearLastNewType removes the last new type that was initiated
+************************************************************************/
+
+void ClearLastNewType(void)
+ {
+ if(TypeCount != 0)
+ --TypeCount ;
+ }
+
+/****************************************************************
+ GetBeatBegin returns the offset from the R-wave for the
+ beginning of the beat (P-wave onset if a P-wave is found).
+*****************************************************************/
+
+int GetBeatBegin(int type)
+ {
+ return(BeatBegins[type]) ;
+ }
+
+/****************************************************************
+ GetBeatEnd returns the offset from the R-wave for the end of
+ a beat (T-wave offset).
+*****************************************************************/
+
+int GetBeatEnd(int type)
+ {
+ return(BeatEnds[type]) ;
+ }
+
+int GetBeatAmp(int type)
+ {
+ return(BeatAmps[type]) ;
+ }
+
+
+/************************************************************************
+ DomCompare2 and DomCompare return similarity indexes between a given
+ beat and the dominant normal type or a given type and the dominant
+ normal type.
+************************************************************************/
+
+double DomCompare2(int *newBeat, int domType)
+ {
+ int shift ;
+ return(CompareBeats2(&BeatTemplates[domType][0],newBeat,&shift)) ;
+ }
+
+double DomCompare(int newType, int domType)
+ {
+ int shift ;
+ return(CompareBeats2(&BeatTemplates[domType][0],&BeatTemplates[newType][0],
+ &shift)) ;
+ }
+
+/*************************************************************************
+BeatCopy copies beat data from a source beat to a destination beat.
+*************************************************************************/
+
+void BeatCopy(int srcBeat, int destBeat)
+ {
+ int i ;
+
+ // Copy template.
+
+ for(i = 0; i < BEATLGTH; ++i)
+ BeatTemplates[destBeat][i] = BeatTemplates[srcBeat][i] ;
+
+ // Move feature information.
+
+ BeatCounts[destBeat] = BeatCounts[srcBeat] ;
+ BeatWidths[destBeat] = BeatWidths[srcBeat] ;
+ BeatCenters[destBeat] = BeatCenters[srcBeat] ;
+ for(i = 0; i < MAXPREV; ++i)
+ {
+ PostClass[destBeat][i] = PostClass[srcBeat][i] ;
+ PCRhythm[destBeat][i] = PCRhythm[srcBeat][i] ;
+ }
+
+ BeatClassifications[destBeat] = BeatClassifications[srcBeat] ;
+ BeatBegins[destBeat] = BeatBegins[srcBeat] ;
+ BeatEnds[destBeat] = BeatBegins[srcBeat] ;
+ BeatsSinceLastMatch[destBeat] = BeatsSinceLastMatch[srcBeat];
+ BeatAmps[destBeat] = BeatAmps[srcBeat] ;
+
+ // Adjust data in dominant beat monitor.
+
+ AdjustDomData(srcBeat,destBeat) ;
+ }
+
+/********************************************************************
+ Minimum beat variation returns a 1 if the previous eight beats
+ have all had similarity indexes less than 0.5.
+*********************************************************************/
+
+int MinimumBeatVariation(int type)
+ {
+ int i ;
+ for(i = 0; i < MAXTYPES; ++i)
+ if(MIs[type][i] > 0.5)
+ i = MAXTYPES+2 ;
+ if(i == MAXTYPES)
+ return(1) ;
+ else return(0) ;
+ }
+
+/**********************************************************************
+ WideBeatVariation returns true if the average similarity index
+ for a given beat type to its template is greater than WIDE_VAR_LIMIT.
+***********************************************************************/
+
+#define WIDE_VAR_LIMIT 0.50
+
+int WideBeatVariation(int type)
+ {
+ int i, n ;
+ double aveMI ;
+
+ n = BeatCounts[type] ;
+ if(n > 8)
+ n = 8 ;
+
+ for(i = 0, aveMI = 0; i <n; ++i)
+ aveMI += MIs[type][i] ;
+
+ aveMI /= n ;
+ if(aveMI > WIDE_VAR_LIMIT)
+ return(1) ;
+ else return(0) ;
+ }
+
+
+
diff -r 1de72bdab0ef -r 8e2fbe131b53 match.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/match.h Sun Jun 28 03:06:00 2015 +0000 @@ -0,0 +1,51 @@ +/***************************************************************************** +FILE: match.h +AUTHOR: Patrick S. Hamilton +REVISED: 12/4/2001 + ___________________________________________________________________________ + +match.h: Beat matching prototype definitions. +Copywrite (C) 2001 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). +******************************************************************************/ + +int NewBeatType(int *beat) ; +void BestMorphMatch(int *newBeat,int *matchType,double *matchIndex, double *mi2, int *shiftAdj) ; +void UpdateBeatType(int matchType,int *newBeat, double mi2, int shiftAdj) ; +int GetTypesCount(void) ; +int GetBeatTypeCount(int type) ; +int IsTypeIsolated(int type) ; +void SetBeatClass(int type, int beatClass) ; +int GetBeatClass(int type) ; +int GetDominantType(void) ; +int GetBeatWidth(int type) ; +int GetPolarity(int type) ; +int GetRhythmIndex(int type) ; +void ResetMatch(void) ; +void ClearLastNewType(void) ; +int GetBeatBegin(int type) ; +int GetBeatEnd(int type) ; +int GetBeatAmp(int type) ; +int MinimumBeatVariation(int type) ; +int GetBeatCenter(int type) ; +int WideBeatVariation(int type) ; +double DomCompare2(int *newBeat, int domType) ; +double DomCompare(int newType, int domType) ;
diff -r 1de72bdab0ef -r 8e2fbe131b53 noisechk.cpp
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/noisechk.cpp Sun Jun 28 03:06:00 2015 +0000
@@ -0,0 +1,131 @@
+/*****************************************************************************
+FILE: noisechk.cpp
+AUTHOR: Patrick S. Hamilton
+REVISED: 5/13/2002
+ ___________________________________________________________________________
+
+noisechk.cpp: Noise Check
+Copywrite (C) 2001 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 contains functions for evaluating the noise content of a beat.
+
+*****************************************************************************/
+
+//#include <stdlib.h>
+#include <mbed.h>
+#include "qrsdet.h"
+
+#define NB_LENGTH MS1500
+#define NS_LENGTH MS50
+
+int NoiseBuffer[NB_LENGTH], NBPtr = 0 ;
+int NoiseEstimate ;
+
+/************************************************************************
+ GetNoiseEstimate() allows external access the present noise estimate.
+ this function is only used for debugging.
+*************************************************************************/
+
+int GetNoiseEstimate()
+ {
+ return(NoiseEstimate) ;
+ }
+
+/***********************************************************************
+ NoiseCheck() must be called for every sample of data. The data is
+ stored in a circular buffer to facilitate noise analysis. When a
+ beat is detected NoiseCheck() is passed the sample delay since the
+ R-wave of the beat occurred (delay), the RR interval between this
+ beat and the next most recent beat, the estimated offset from the
+ R-wave to the beginning of the beat (beatBegin), and the estimated
+ offset from the R-wave to the end of the beat.
+
+ NoiseCheck() estimates the noise in the beat by the maximum and
+ minimum signal values in either a window from the end of the
+ previous beat to the beginning of the present beat, or a 250 ms
+ window preceding the present beat, which ever is shorter.
+
+ NoiseCheck() returns ratio of the signal variation in the window
+ between beats to the length of the window between the beats. If
+ the heart rate is too high and the beat durations are too long,
+ NoiseCheck() returns 0.
+
+***********************************************************************/
+
+int NoiseCheck(int datum, int delay, int RR, int beatBegin, int beatEnd)
+ {
+ int ptr, i;
+ int ncStart, ncEnd, ncMax, ncMin ;
+ double noiseIndex ;
+
+ NoiseBuffer[NBPtr] = datum ;
+ if(++NBPtr == NB_LENGTH)
+ NBPtr = 0 ;
+
+ // Check for noise in region that is 300 ms following
+ // last R-wave and 250 ms preceding present R-wave.
+
+ ncStart = delay+RR-beatEnd ; // Calculate offset to end of previous beat.
+ ncEnd = delay+beatBegin ; // Calculate offset to beginning of this beat.
+ if(ncStart > ncEnd + MS250)
+ ncStart = ncEnd + MS250 ;
+
+
+ // Estimate noise if delay indicates a beat has been detected,
+ // the delay is not to long for the data buffer, and there is
+ // some space between the end of the last beat and the beginning
+ // of this beat.
+
+ if((delay != 0) && (ncStart < NB_LENGTH) && (ncStart > ncEnd))
+ {
+
+ ptr = NBPtr - ncStart ; // Find index to end of last beat in
+ if(ptr < 0) // the circular buffer.
+ ptr += NB_LENGTH ;
+
+ // Find the maximum and minimum values in the
+ // isoelectric region between beats.
+
+ ncMax = ncMin = NoiseBuffer[ptr] ;
+ for(i = 0; i < ncStart-ncEnd; ++i)
+ {
+ if(NoiseBuffer[ptr] > ncMax)
+ ncMax = NoiseBuffer[ptr] ;
+ else if(NoiseBuffer[ptr] < ncMin)
+ ncMin = NoiseBuffer[ptr] ;
+ if(++ptr == NB_LENGTH)
+ ptr = 0 ;
+ }
+
+ // The noise index is the ratio of the signal variation
+ // over the isoelectric window length, scaled by 10.
+
+ noiseIndex = (ncMax-ncMin) ;
+ noiseIndex /= (ncStart-ncEnd) ;
+ NoiseEstimate = noiseIndex * 10 ;
+ }
+ else
+ NoiseEstimate = 0 ;
+ return(NoiseEstimate) ;
+ }
+
diff -r 1de72bdab0ef -r 8e2fbe131b53 postclas.cpp
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/postclas.cpp Sun Jun 28 03:06:00 2015 +0000
@@ -0,0 +1,234 @@
+/*****************************************************************************
+FILE: postclas.cpp
+AUTHOR: Patrick S. Hamilton
+REVISED: 5/13/2002
+ ___________________________________________________________________________
+
+postclas.cpp: Post classifier
+Copywrite (C) 2002 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 contains functions for classifying beats based after the
+following beat is detected.
+
+ ResetPostClassify() -- Resets static variables used by
+ PostClassify()
+ PostClassify() -- classifies each beat based on six preceding
+ beats and the following beat.
+ CheckPostClass() -- classifys beat type based on the last
+ eight post classifications of that beat.
+ CheckPCRhythm() -- returns the classification of the RR interval
+ for this type of beat based its previous eight RR intervals.
+
+****************************************************************/
+
+#include <mbed.h>
+#include "bdac.h"
+#include "ecgcodes.h"
+
+// External Prototypes.
+
+double DomCompare(int newType, int domType) ;
+int GetBeatTypeCount(int type) ;
+
+// Records of post classifications.
+
+int PostClass[MAXTYPES][8], PCInitCount = 0 ;
+int PCRhythm[MAXTYPES][8] ;
+
+/**********************************************************************
+ Resets post classifications for beats.
+**********************************************************************/
+
+void ResetPostClassify()
+ {
+ int i, j ;
+ for(i = 0; i < MAXTYPES; ++i)
+ for(j = 0; j < 8; ++j)
+ {
+ PostClass[i][j] = 0 ;
+ PCRhythm[i][j] = 0 ;
+ }
+ PCInitCount = 0 ;
+ }
+
+/***********************************************************************
+ Classify the previous beat type and rhythm type based on this beat
+ and the preceding beat. This classifier is more sensitive
+ to detecting premature beats followed by compensitory pauses.
+************************************************************************/
+
+void PostClassify(int *recentTypes, int domType, int *recentRRs, int width, double mi2,
+ int rhythmClass)
+ {
+ static int lastRC, lastWidth ;
+ static double lastMI2 ;
+ int i, regCount, pvcCount, normRR ;
+ double mi3 ;
+
+ // If the preceeding and following beats are the same type,
+ // they are generally regular, and reasonably close in shape
+ // to the dominant type, consider them to be dominant.
+
+ if((recentTypes[0] == recentTypes[2]) && (recentTypes[0] != domType)
+ && (recentTypes[0] != recentTypes[1]))
+ {
+ mi3 = DomCompare(recentTypes[0],domType) ;
+ for(i = regCount = 0; i < 8; ++i)
+ if(PCRhythm[recentTypes[0]][i] == NORMAL)
+ ++regCount ;
+ if((mi3 < 2.0) && (regCount > 6))
+ domType = recentTypes[0] ;
+ }
+
+ // Don't do anything until four beats have gone by.
+
+ if(PCInitCount < 3)
+ {
+ ++PCInitCount ;
+ lastWidth = width ;
+ lastMI2 = 0 ;
+ lastRC = 0 ;
+ return ;
+ }
+
+ if(recentTypes[1] < MAXTYPES)
+ {
+
+ // Find first NN interval.
+ for(i = 2; (i < 7) && (recentTypes[i] != recentTypes[i+1]); ++i) ;
+ if(i == 7) normRR = 0 ;
+ else normRR = recentRRs[i] ;
+
+ // Shift the previous beat classifications to make room for the
+ // new classification.
+ for(i = pvcCount = 0; i < 8; ++i)
+ if(PostClass[recentTypes[1]][i] == PVC)
+ ++pvcCount ;
+
+ for(i = 7; i > 0; --i)
+ {
+ PostClass[recentTypes[1]][i] = PostClass[recentTypes[1]][i-1] ;
+ PCRhythm[recentTypes[1]][i] = PCRhythm[recentTypes[1]][i-1] ;
+ }
+
+ // If the beat is premature followed by a compensitory pause and the
+ // previous and following beats are normal, post classify as
+ // a PVC.
+
+ if(((normRR-(normRR>>3)) >= recentRRs[1]) && ((recentRRs[0]-(recentRRs[0]>>3)) >= normRR)// && (lastMI2 > 3)
+ && (recentTypes[0] == domType) && (recentTypes[2] == domType)
+ && (recentTypes[1] != domType))
+ PostClass[recentTypes[1]][0] = PVC ;
+
+ // If previous two were classified as PVCs, and this is at least slightly
+ // premature, classify as a PVC.
+
+ else if(((normRR-(normRR>>4)) > recentRRs[1]) && ((normRR+(normRR>>4)) < recentRRs[0]) &&
+ (((PostClass[recentTypes[1]][1] == PVC) && (PostClass[recentTypes[1]][2] == PVC)) ||
+ (pvcCount >= 6) ) &&
+ (recentTypes[0] == domType) && (recentTypes[2] == domType) && (recentTypes[1] != domType))
+ PostClass[recentTypes[1]][0] = PVC ;
+
+ // If the previous and following beats are the dominant beat type,
+ // and this beat is significantly different from the dominant,
+ // call it a PVC.
+
+ else if((recentTypes[0] == domType) && (recentTypes[2] == domType) && (lastMI2 > 2.5))
+ PostClass[recentTypes[1]][0] = PVC ;
+
+ // Otherwise post classify this beat as UNKNOWN.
+
+ else PostClass[recentTypes[1]][0] = UNKNOWN ;
+
+ // If the beat is premature followed by a compensitory pause, post
+ // classify the rhythm as PVC.
+
+ if(((normRR-(normRR>>3)) > recentRRs[1]) && ((recentRRs[0]-(recentRRs[0]>>3)) > normRR))
+ PCRhythm[recentTypes[1]][0] = PVC ;
+
+ // Otherwise, post classify the rhythm as the same as the
+ // regular rhythm classification.
+
+ else PCRhythm[recentTypes[1]][0] = lastRC ;
+ }
+
+ lastWidth = width ;
+ lastMI2 = mi2 ;
+ lastRC = rhythmClass ;
+ }
+
+
+/*************************************************************************
+ CheckPostClass checks to see if three of the last four or six of the
+ last eight of a given beat type have been post classified as PVC.
+*************************************************************************/
+
+int CheckPostClass(int type)
+ {
+ int i, pvcs4 = 0, pvcs8 ;
+
+ if(type == MAXTYPES)
+ return(UNKNOWN) ;
+
+ for(i = 0; i < 4; ++i)
+ if(PostClass[type][i] == PVC)
+ ++pvcs4 ;
+ for(pvcs8=pvcs4; i < 8; ++i)
+ if(PostClass[type][i] == PVC)
+ ++pvcs8 ;
+
+ if((pvcs4 >= 3) || (pvcs8 >= 6))
+ return(PVC) ;
+ else return(UNKNOWN) ;
+ }
+
+/****************************************************************************
+ Check classification of previous beats' rhythms based on post beat
+ classification. If 7 of 8 previous beats were classified as NORMAL
+ (regular) classify the beat type as NORMAL (regular).
+ Call it a PVC if 2 of the last 8 were regular.
+****************************************************************************/
+
+int CheckPCRhythm(int type)
+ {
+ int i, normCount, n ;
+
+
+ if(type == MAXTYPES)
+ return(UNKNOWN) ;
+
+ if(GetBeatTypeCount(type) < 9)
+ n = GetBeatTypeCount(type)-1 ;
+ else n = 8 ;
+
+ for(i = normCount = 0; i < n; ++i)
+ if(PCRhythm[type][i] == NORMAL)
+ ++normCount;
+ if(normCount >= 7)
+ return(NORMAL) ;
+ if(((normCount == 0) && (n < 4)) ||
+ ((normCount <= 1) && (n >= 4) && (n < 7)) ||
+ ((normCount <= 2) && (n >= 7)))
+ return(PVC) ;
+ return(UNKNOWN) ;
+ }
\ No newline at end of file
diff -r 1de72bdab0ef -r 8e2fbe131b53 postclas.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/postclas.h Sun Jun 28 03:06:00 2015 +0000 @@ -0,0 +1,5 @@ +void ResetPostClassify() ; +void PostClassify(int *recentTypes, int domType, int *recentRRs, int width, double mi2, + int rhythmClass) ; +int CheckPostClass(int type) ; +int CheckPCRhythm(int type) ; \ No newline at end of file
diff -r 1de72bdab0ef -r 8e2fbe131b53 qrsdet.cpp
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/qrsdet.cpp Sun Jun 28 03:06:00 2015 +0000
@@ -0,0 +1,452 @@
+/*****************************************************************************
+FILE: qrsdet.cpp
+AUTHOR: Patrick S. Hamilton
+REVISED: 12/04/2000
+ ___________________________________________________________________________
+
+qrsdet.cpp: A QRS detector.
+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 contains functions for detecting QRS complexes in an ECG. The
+QRS detector requires filter functions in qrsfilt.cpp and parameter
+definitions in qrsdet.h. QRSDet is the only function that needs to be
+visable outside of these files.
+
+Syntax:
+ int QRSDet(int ecgSample, int init) ;
+
+Description:
+ QRSDet() implements a modified version of the QRS detection
+ algorithm described in:
+
+ Hamilton, Tompkins, W. J., "Quantitative investigation of QRS
+ detection rules using the MIT/BIH arrhythmia database",
+ IEEE Trans. Biomed. Eng., BME-33, pp. 1158-1165, 1987.
+
+ Consecutive ECG samples are passed to QRSDet. QRSDet was
+ designed for a 200 Hz sample rate. QRSDet contains a number
+ of static variables that it uses to adapt to different ECG
+ signals. These variables can be reset by passing any value
+ not equal to 0 in init.
+
+ Note: QRSDet() requires filters in QRSFilt.cpp
+
+Returns:
+ When a QRS complex is detected QRSDet returns the detection delay.
+
+****************************************************************/
+
+#include <mbed.h>
+//#include <mem.h> /* For memmov. */
+//#include <math.h>
+#include "qrsdet.h"
+#define PRE_BLANK MS200
+
+// External Prototypes.
+
+int QRSFilter(int datum, int init) ;
+int deriv1( int x0, int init ) ;
+
+// Local Prototypes.
+
+int Peak( int datum, int init ) ;
+int median(int *array, int datnum) ;
+int thresh(int qmedian, int nmedian) ;
+int BLSCheck(int *dBuf,int dbPtr,int *maxder) ;
+
+int earlyThresh(int qmedian, int nmedian) ;
+
+
+double TH = 0.475 ;
+
+int DDBuffer[DER_DELAY], DDPtr ; /* Buffer holding derivative data. */
+int Dly = 0 ;
+
+const int MEMMOVELEN = 7*sizeof(int);
+
+int QRSDet( int datum, int init )
+ {
+ static int det_thresh, qpkcnt = 0 ;
+ static int qrsbuf[8], noise[8], rrbuf[8] ;
+ static int rsetBuff[8], rsetCount = 0 ;
+ static int nmedian, qmedian, rrmedian ;
+ static int count, sbpeak = 0, sbloc, sbcount = MS1500 ;
+ static int maxder, lastmax ;
+ static int initBlank, initMax ;
+ static int preBlankCnt, tempPeak ;
+
+ int fdatum, QrsDelay = 0 ;
+ int i, newPeak, aPeak ;
+
+/* Initialize all buffers to 0 on the first call. */
+
+ if( init )
+ {
+ for(i = 0; i < 8; ++i)
+ {
+ noise[i] = 0 ; /* Initialize noise buffer */
+ rrbuf[i] = MS1000 ;/* and R-to-R interval buffer. */
+ }
+
+ qpkcnt = maxder = lastmax = count = sbpeak = 0 ;
+ initBlank = initMax = preBlankCnt = DDPtr = 0 ;
+ sbcount = MS1500 ;
+ QRSFilter(0,1) ; /* initialize filters. */
+ Peak(0,1) ;
+ }
+
+ fdatum = QRSFilter(datum,0) ; /* Filter data. */
+
+
+ /* Wait until normal detector is ready before calling early detections. */
+
+ aPeak = Peak(fdatum,0) ;
+
+ // Hold any peak that is detected for 200 ms
+ // in case a bigger one comes along. There
+ // can only be one QRS complex in any 200 ms window.
+
+ newPeak = 0 ;
+ if(aPeak && !preBlankCnt) // If there has been no peak for 200 ms
+ { // save this one and start counting.
+ tempPeak = aPeak ;
+ preBlankCnt = PRE_BLANK ; // MS200
+ }
+
+ else if(!aPeak && preBlankCnt) // If we have held onto a peak for
+ { // 200 ms pass it on for evaluation.
+ if(--preBlankCnt == 0)
+ newPeak = tempPeak ;
+ }
+
+ else if(aPeak) // If we were holding a peak, but
+ { // this ones bigger, save it and
+ if(aPeak > tempPeak) // start counting to 200 ms again.
+ {
+ tempPeak = aPeak ;
+ preBlankCnt = PRE_BLANK ; // MS200
+ }
+ else if(--preBlankCnt == 0)
+ newPeak = tempPeak ;
+ }
+
+/* newPeak = 0 ;
+ if((aPeak != 0) && (preBlankCnt == 0))
+ newPeak = aPeak ;
+ else if(preBlankCnt != 0) --preBlankCnt ; */
+
+
+
+ /* Save derivative of raw signal for T-wave and baseline
+ shift discrimination. */
+
+ DDBuffer[DDPtr] = deriv1( datum, 0 ) ;
+ if(++DDPtr == DER_DELAY)
+ DDPtr = 0 ;
+
+ /* Initialize the qrs peak buffer with the first eight */
+ /* local maximum peaks detected. */
+
+ if( qpkcnt < 8 )
+ {
+ ++count ;
+ if(newPeak > 0) count = WINDOW_WIDTH ;
+ if(++initBlank == MS1000)
+ {
+ initBlank = 0 ;
+ qrsbuf[qpkcnt] = initMax ;
+ initMax = 0 ;
+ ++qpkcnt ;
+ if(qpkcnt == 8)
+ {
+ qmedian = median( qrsbuf, 8 ) ;
+ nmedian = 0 ;
+ rrmedian = MS1000 ;
+ sbcount = MS1500+MS150 ;
+ det_thresh = thresh(qmedian,nmedian) ;
+ }
+ }
+ if( newPeak > initMax )
+ initMax = newPeak ;
+ }
+
+ else /* Else test for a qrs. */
+ {
+ ++count ;
+ if(newPeak > 0)
+ {
+
+
+ /* Check for maximum derivative and matching minima and maxima
+ for T-wave and baseline shift rejection. Only consider this
+ peak if it doesn't seem to be a base line shift. */
+
+ if(!BLSCheck(DDBuffer, DDPtr, &maxder))
+ {
+
+
+ // Classify the beat as a QRS complex
+ // if the peak is larger than the detection threshold.
+
+ if(newPeak > det_thresh)
+ {
+ memmove(&qrsbuf[1], qrsbuf, MEMMOVELEN) ;
+ qrsbuf[0] = newPeak ;
+ qmedian = median(qrsbuf,8) ;
+ det_thresh = thresh(qmedian,nmedian) ;
+ memmove(&rrbuf[1], rrbuf, MEMMOVELEN) ;
+ rrbuf[0] = count - WINDOW_WIDTH ;
+ rrmedian = median(rrbuf,8) ;
+ sbcount = rrmedian + (rrmedian >> 1) + WINDOW_WIDTH ;
+ count = WINDOW_WIDTH ;
+
+ sbpeak = 0 ;
+
+ lastmax = maxder ;
+ maxder = 0 ;
+ QrsDelay = WINDOW_WIDTH + FILTER_DELAY ;
+ initBlank = initMax = rsetCount = 0 ;
+
+ // preBlankCnt = PRE_BLANK ;
+ }
+
+ // If a peak isn't a QRS update noise buffer and estimate.
+ // Store the peak for possible search back.
+
+
+ else
+ {
+ memmove(&noise[1],noise,MEMMOVELEN) ;
+ noise[0] = newPeak ;
+ nmedian = median(noise,8) ;
+ det_thresh = thresh(qmedian,nmedian) ;
+
+ // Don't include early peaks (which might be T-waves)
+ // in the search back process. A T-wave can mask
+ // a small following QRS.
+
+ if((newPeak > sbpeak) && ((count-WINDOW_WIDTH) >= MS360))
+ {
+ sbpeak = newPeak ;
+ sbloc = count - WINDOW_WIDTH ;
+ }
+ }
+ }
+ }
+
+ /* Test for search back condition. If a QRS is found in */
+ /* search back update the QRS buffer and det_thresh. */
+
+ if((count > sbcount) && (sbpeak > (det_thresh >> 1)))
+ {
+ memmove(&qrsbuf[1],qrsbuf,MEMMOVELEN) ;
+ qrsbuf[0] = sbpeak ;
+ qmedian = median(qrsbuf,8) ;
+ det_thresh = thresh(qmedian,nmedian) ;
+ memmove(&rrbuf[1],rrbuf,MEMMOVELEN) ;
+ rrbuf[0] = sbloc ;
+ rrmedian = median(rrbuf,8) ;
+ sbcount = rrmedian + (rrmedian >> 1) + WINDOW_WIDTH ;
+ QrsDelay = count = count - sbloc ;
+ QrsDelay += FILTER_DELAY ;
+ sbpeak = 0 ;
+ lastmax = maxder ;
+ maxder = 0 ;
+ initBlank = initMax = rsetCount = 0 ;
+ }
+ }
+
+ // In the background estimate threshold to replace adaptive threshold
+ // if eight seconds elapses without a QRS detection.
+
+ if( qpkcnt == 8 )
+ {
+ if(++initBlank == MS1000)
+ {
+ initBlank = 0 ;
+ rsetBuff[rsetCount] = initMax ;
+ initMax = 0 ;
+ ++rsetCount ;
+
+ // Reset threshold if it has been 8 seconds without
+ // a detection.
+
+ if(rsetCount == 8)
+ {
+ for(i = 0; i < 8; ++i)
+ {
+ qrsbuf[i] = rsetBuff[i] ;
+ noise[i] = 0 ;
+ }
+ qmedian = median( rsetBuff, 8 ) ;
+ nmedian = 0 ;
+ rrmedian = MS1000 ;
+ sbcount = MS1500+MS150 ;
+ det_thresh = thresh(qmedian,nmedian) ;
+ initBlank = initMax = rsetCount = 0 ;
+ sbpeak = 0 ;
+ }
+ }
+ if( newPeak > initMax )
+ initMax = newPeak ;
+ }
+
+ return(QrsDelay) ;
+ }
+
+/**************************************************************
+* peak() takes a datum as input and returns a peak height
+* when the signal returns to half its peak height, or
+**************************************************************/
+
+int Peak( int datum, int init )
+ {
+ static int max = 0, timeSinceMax = 0, lastDatum ;
+ int pk = 0 ;
+
+ if(init)
+ max = timeSinceMax = 0 ;
+
+ if(timeSinceMax > 0)
+ ++timeSinceMax ;
+
+ if((datum > lastDatum) && (datum > max))
+ {
+ max = datum ;
+ if(max > 2)
+ timeSinceMax = 1 ;
+ }
+
+ else if(datum < (max >> 1))
+ {
+ pk = max ;
+ max = 0 ;
+ timeSinceMax = 0 ;
+ Dly = 0 ;
+ }
+
+ else if(timeSinceMax > MS95)
+ {
+ pk = max ;
+ max = 0 ;
+ timeSinceMax = 0 ;
+ Dly = 3 ;
+ }
+ lastDatum = datum ;
+ return(pk) ;
+ }
+
+/********************************************************************
+median returns the median of an array of integers. It uses a slow
+sort algorithm, but these arrays are small, so it hardly matters.
+********************************************************************/
+
+int median(int *array, int datnum)
+ {
+ int i, j, k, temp, sort[20] ;
+ for(i = 0; i < datnum; ++i)
+ sort[i] = array[i] ;
+ for(i = 0; i < datnum; ++i)
+ {
+ temp = sort[i] ;
+ for(j = 0; (temp < sort[j]) && (j < i) ; ++j) ;
+ for(k = i - 1 ; k >= j ; --k)
+ sort[k+1] = sort[k] ;
+ sort[j] = temp ;
+ }
+ return(sort[datnum>>1]) ;
+ }
+/*
+int median(int *array, int datnum)
+ {
+ long sum ;
+ int i ;
+
+ for(i = 0, sum = 0; i < datnum; ++i)
+ sum += array[i] ;
+ sum /= datnum ;
+ return(sum) ;
+ } */
+
+/****************************************************************************
+ thresh() calculates the detection threshold from the qrs median and noise
+ median estimates.
+****************************************************************************/
+
+int thresh(int qmedian, int nmedian)
+ {
+ int thrsh, dmed ;
+ double temp ;
+ dmed = qmedian - nmedian ;
+/* thrsh = nmedian + (dmed>>2) + (dmed>>3) + (dmed>>4); */
+ temp = dmed ;
+ temp *= TH ;
+ dmed = temp ;
+ thrsh = nmedian + dmed ; /* dmed * THRESHOLD */
+ return(thrsh) ;
+ }
+
+/***********************************************************************
+ BLSCheck() reviews data to see if a baseline shift has occurred.
+ This is done by looking for both positive and negative slopes of
+ roughly the same magnitude in a 220 ms window.
+***********************************************************************/
+
+int BLSCheck(int *dBuf,int dbPtr,int *maxder)
+ {
+ int max, min, maxt, mint, t, x ;
+ max = min = 0 ;
+
+ return(0) ;
+
+ for(t = 0; t < MS220; ++t)
+ {
+ x = dBuf[dbPtr] ;
+ if(x > max)
+ {
+ maxt = t ;
+ max = x ;
+ }
+ else if(x < min)
+ {
+ mint = t ;
+ min = x;
+ }
+ if(++dbPtr == DER_DELAY)
+ dbPtr = 0 ;
+ }
+
+ *maxder = max ;
+ min = -min ;
+
+ /* Possible beat if a maximum and minimum pair are found
+ where the interval between them is less than 150 ms. */
+
+ if((max > (min>>3)) && (min > (max>>3)) &&
+ (abs(maxt - mint) < MS150))
+ return(0) ;
+
+ else
+ return(1) ;
+ }
+
diff -r 1de72bdab0ef -r 8e2fbe131b53 qrsdet.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/qrsdet.h Sun Jun 28 03:06:00 2015 +0000 @@ -0,0 +1,62 @@ +/***************************************************************************** +FILE: qrsdet.h +AUTHOR: Patrick S. Hamilton +REVISED: 4/16/2002 + ___________________________________________________________________________ + +qrsdet.h QRS detector parameter definitions +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.com) 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). + __________________________________________________________________________ + Revisions: + 4/16: Modified to allow simplified modification of digital filters in + qrsfilt(). +*****************************************************************************/ + + +#define SAMPLE_RATE 500 /* Sample rate in Hz. - was 200 */ +#define MS_PER_SAMPLE ( (double) 1000/ (double) SAMPLE_RATE) +#define MS10 ((int) (10/ MS_PER_SAMPLE + 0.5)) +#define MS25 ((int) (25/MS_PER_SAMPLE + 0.5)) +#define MS30 ((int) (30/MS_PER_SAMPLE + 0.5)) +#define MS80 ((int) (80/MS_PER_SAMPLE + 0.5)) +#define MS95 ((int) (95/MS_PER_SAMPLE + 0.5)) +#define MS100 ((int) (100/MS_PER_SAMPLE + 0.5)) +#define MS125 ((int) (125/MS_PER_SAMPLE + 0.5)) +#define MS150 ((int) (150/MS_PER_SAMPLE + 0.5)) +#define MS160 ((int) (160/MS_PER_SAMPLE + 0.5)) +#define MS175 ((int) (175/MS_PER_SAMPLE + 0.5)) +#define MS195 ((int) (195/MS_PER_SAMPLE + 0.5)) +#define MS200 ((int) (200/MS_PER_SAMPLE + 0.5)) +#define MS220 ((int) (220/MS_PER_SAMPLE + 0.5)) +#define MS250 ((int) (250/MS_PER_SAMPLE + 0.5)) +#define MS300 ((int) (300/MS_PER_SAMPLE + 0.5)) +#define MS360 ((int) (360/MS_PER_SAMPLE + 0.5)) +#define MS450 ((int) (450/MS_PER_SAMPLE + 0.5)) +#define MS1000 SAMPLE_RATE +#define MS1500 ((int) (1500/MS_PER_SAMPLE)) +#define DERIV_LENGTH MS10 +#define LPBUFFER_LGTH ((int) (2*MS25)) +#define HPBUFFER_LGTH MS125 + +#define WINDOW_WIDTH MS80 // Moving window integration width. +#define FILTER_DELAY (int) (((double) DERIV_LENGTH/2) + ((double) LPBUFFER_LGTH/2 - 1) + (((double) HPBUFFER_LGTH-1)/2) + PRE_BLANK) // filter delays plus 200 ms blanking delay +#define DER_DELAY WINDOW_WIDTH + FILTER_DELAY + MS100
diff -r 1de72bdab0ef -r 8e2fbe131b53 qrsfilt.cpp
--- /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) ;
+ }
diff -r 1de72bdab0ef -r 8e2fbe131b53 rythmchk.cpp
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/rythmchk.cpp Sun Jun 28 03:06:00 2015 +0000
@@ -0,0 +1,492 @@
+/*****************************************************************************
+FILE: rythmchk.cpp
+AUTHOR: Patrick S. Hamilton
+REVISED: 5/13/2002
+ 10/9/2001: 1.1 Call premature for 12.5% difference with very regular
+ rhythms. If short after VV go to QQ.
+ 5/13/2002: Check for NNVNNNV pattern when last interval was QQ.
+ ___________________________________________________________________________
+
+rythmchk.cpp: Rhythm Check
+Copywrite (C) 2001 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).
+ __________________________________________________________________________
+
+ Rythmchk.cpp contains functions for classifying RR intervals as either
+ NORMAL, PVC, or UNKNOWN. Intervals classified as NORMAL are presumed to
+ end with normal beats, intervals classified as PVC are presumed to end
+ with a premature contraction, and intervals classified as unknown do
+ not fit any pattern that the rhythm classifier recognizes.
+
+ NORMAL intervals can be part of a normal (regular) rhythm, normal beats
+ following a premature beats, or normal beats following runs of ventricular
+ beats. PVC intervals can be short intervals following a regular rhythm,
+ short intervals that are part of a run of short intervals following a
+ regular rhythm, or short intervals that are part of a bigeminal rhythm.
+
+**************************************************************************/
+
+#include <mbed.h>
+#include "qrsdet.h" // For time intervals.
+#include "ecgcodes.h" // Defines codes of NORMAL, PVC, and UNKNOWN.
+//#include <stdlib.h> // For abs()
+
+// Define RR interval types.
+
+#define QQ 0 // Unknown-Unknown interval.
+#define NN 1 // Normal-Normal interval.
+#define NV 2 // Normal-PVC interval.
+#define VN 3 // PVC-Normal interval.
+#define VV 4 // PVC-PVC interval.
+
+#define RBB_LENGTH 8
+#define LEARNING 0
+#define READY 1
+
+#define BRADY_LIMIT MS1500
+
+// Local prototypes.
+int RRMatch(int rr0,int rr1) ;
+int RRShort(int rr0,int rr1) ;
+int RRShort2(int *rrIntervals, int *rrTypes) ;
+int RRMatch2(int rr0,int rr1) ;
+
+// Global variables.
+int RRBuffer[RBB_LENGTH], RRTypes[RBB_LENGTH], BeatCount = 0;
+int ClassifyState = LEARNING ;
+
+int BigeminyFlag ;
+
+/***************************************************************************
+ ResetRhythmChk() resets static variables used for rhythm classification.
+****************************************************************************/
+
+void ResetRhythmChk(void)
+ {
+ BeatCount = 0 ;
+ ClassifyState = LEARNING ;
+ }
+
+/*****************************************************************************
+ RhythmChk() takes an R-to-R interval as input and, based on previous R-to-R
+ intervals, classifys the interval as NORMAL, PVC, or UNKNOWN.
+******************************************************************************/
+
+int RhythmChk(int rr)
+ {
+ int i, regular = 1 ;
+ int NNEst, NVEst ;
+
+ BigeminyFlag = 0 ;
+
+ // Wait for at least 4 beats before classifying anything.
+
+ if(BeatCount < 4)
+ {
+ if(++BeatCount == 4)
+ ClassifyState = READY ;
+ }
+
+ // Stick the new RR interval into the RR interval Buffer.
+
+ for(i = RBB_LENGTH-1; i > 0; --i)
+ {
+ RRBuffer[i] = RRBuffer[i-1] ;
+ RRTypes[i] = RRTypes[i-1] ;
+ }
+
+ RRBuffer[0] = rr ;
+
+ if(ClassifyState == LEARNING)
+ {
+ RRTypes[0] = QQ ;
+ return(UNKNOWN) ;
+ }
+
+ // If we couldn't tell what the last interval was...
+
+ if(RRTypes[1] == QQ)
+ {
+ for(i = 0, regular = 1; i < 3; ++i)
+ if(RRMatch(RRBuffer[i],RRBuffer[i+1]) == 0)
+ regular = 0 ;
+
+ // If this, and the last three intervals matched, classify
+ // it as Normal-Normal.
+
+ if(regular == 1)
+ {
+ RRTypes[0] = NN ;
+ return(NORMAL) ;
+ }
+
+ // Check for bigeminy.
+ // Call bigeminy if every other RR matches and
+ // consecutive beats do not match.
+
+ for(i = 0, regular = 1; i < 6; ++i)
+ if(RRMatch(RRBuffer[i],RRBuffer[i+2]) == 0)
+ regular = 0 ;
+ for(i = 0; i < 6; ++i)
+ if(RRMatch(RRBuffer[i],RRBuffer[i+1]) != 0)
+ regular = 0 ;
+
+ if(regular == 1)
+ {
+ BigeminyFlag = 1 ;
+ if(RRBuffer[0] < RRBuffer[1])
+ {
+ RRTypes[0] = NV ;
+ RRTypes[1] = VN ;
+ return(PVC) ;
+ }
+ else
+ {
+ RRTypes[0] = VN ;
+ RRTypes[1] = NV ;
+ return(NORMAL) ;
+ }
+ }
+
+ // Check for NNVNNNV pattern.
+
+ if(RRShort(RRBuffer[0],RRBuffer[1]) && RRMatch(RRBuffer[1],RRBuffer[2])
+ && RRMatch(RRBuffer[2]*2,RRBuffer[3]+RRBuffer[4]) &&
+ RRMatch(RRBuffer[4],RRBuffer[0]) && RRMatch(RRBuffer[5],RRBuffer[2]))
+ {
+ RRTypes[0] = NV ;
+ RRTypes[1] = NN ;
+ return(PVC) ;
+ }
+
+ // If the interval is not part of a
+ // bigeminal or regular pattern, give up.
+
+ else
+ {
+ RRTypes[0] = QQ ;
+ return(UNKNOWN) ;
+ }
+ }
+
+ // If the previous two beats were normal...
+
+ else if(RRTypes[1] == NN)
+ {
+
+ if(RRShort2(RRBuffer,RRTypes))
+ {
+ if(RRBuffer[1] < BRADY_LIMIT)
+ {
+ RRTypes[0] = NV ;
+ return(PVC) ;
+ }
+ else RRTypes[0] = QQ ;
+ return(UNKNOWN) ;
+ }
+
+
+ // If this interval matches the previous interval, then it
+ // is regular.
+
+ else if(RRMatch(RRBuffer[0],RRBuffer[1]))
+ {
+ RRTypes[0] = NN ;
+ return(NORMAL) ;
+ }
+
+ // If this interval is short..
+
+ else if(RRShort(RRBuffer[0],RRBuffer[1]))
+ {
+
+ // But matches the one before last and the one before
+ // last was NN, this is a normal interval.
+
+ if(RRMatch(RRBuffer[0],RRBuffer[2]) && (RRTypes[2] == NN))
+ {
+ RRTypes[0] = NN ;
+ return(NORMAL) ;
+ }
+
+ // If the rhythm wasn't bradycardia, call it a PVC.
+
+ else if(RRBuffer[1] < BRADY_LIMIT)
+ {
+ RRTypes[0] = NV ;
+ return(PVC) ;
+ }
+
+ // If the regular rhythm was bradycardia, don't assume that
+ // it was a PVC.
+
+ else
+ {
+ RRTypes[0] = QQ ;
+ return(UNKNOWN) ;
+ }
+ }
+
+ // If the interval isn't normal or short, then classify
+ // it as normal but don't assume normal for future
+ // rhythm classification.
+
+ else
+ {
+ RRTypes[0] = QQ ;
+ return(NORMAL) ;
+ }
+ }
+
+ // If the previous beat was a PVC...
+
+ else if(RRTypes[1] == NV)
+ {
+
+ if(RRShort2(&RRBuffer[1],&RRTypes[1]))
+ {
+ /* if(RRMatch2(RRBuffer[0],RRBuffer[1]))
+ {
+ RRTypes[0] = VV ;
+ return(PVC) ;
+ } */
+
+ if(RRMatch(RRBuffer[0],RRBuffer[1]))
+ {
+ RRTypes[0] = NN ;
+ RRTypes[1] = NN ;
+ return(NORMAL) ;
+ }
+ else if(RRBuffer[0] > RRBuffer[1])
+ {
+ RRTypes[0] = VN ;
+ return(NORMAL) ;
+ }
+ else
+ {
+ RRTypes[0] = QQ ;
+ return(UNKNOWN) ;
+ }
+
+
+ }
+
+ // If this interval matches the previous premature
+ // interval assume a ventricular couplet.
+
+ else if(RRMatch(RRBuffer[0],RRBuffer[1]))
+ {
+ RRTypes[0] = VV ;
+ return(PVC) ;
+ }
+
+ // If this interval is larger than the previous
+ // interval, assume that it is NORMAL.
+
+ else if(RRBuffer[0] > RRBuffer[1])
+ {
+ RRTypes[0] = VN ;
+ return(NORMAL) ;
+ }
+
+ // Otherwise don't make any assumputions about
+ // what this interval represents.
+
+ else
+ {
+ RRTypes[0] = QQ ;
+ return(UNKNOWN) ;
+ }
+ }
+
+ // If the previous beat followed a PVC or couplet etc...
+
+ else if(RRTypes[1] == VN)
+ {
+
+ // Find the last NN interval.
+
+ for(i = 2; (RRTypes[i] != NN) && (i < RBB_LENGTH); ++i) ;
+
+ // If there was an NN interval in the interval buffer...
+ if(i != RBB_LENGTH)
+ {
+ NNEst = RRBuffer[i] ;
+
+ // and it matches, classify this interval as NORMAL.
+
+ if(RRMatch(RRBuffer[0],NNEst))
+ {
+ RRTypes[0] = NN ;
+ return(NORMAL) ;
+ }
+ }
+
+ else NNEst = 0 ;
+ for(i = 2; (RRTypes[i] != NV) && (i < RBB_LENGTH); ++i) ;
+ if(i != RBB_LENGTH)
+ NVEst = RRBuffer[i] ;
+ else NVEst = 0 ;
+ if((NNEst == 0) && (NVEst != 0))
+ NNEst = (RRBuffer[1]+NVEst) >> 1 ;
+
+ // NNEst is either the last NN interval or the average
+ // of the most recent NV and VN intervals.
+
+ // If the interval is closer to NN than NV, try
+ // matching to NN.
+
+ if((NVEst != 0) &&
+ (abs(NNEst - RRBuffer[0]) < abs(NVEst - RRBuffer[0])) &&
+ RRMatch(NNEst,RRBuffer[0]))
+ {
+ RRTypes[0] = NN ;
+ return(NORMAL) ;
+ }
+
+ // If this interval is closer to NV than NN, try
+ // matching to NV.
+
+ else if((NVEst != 0) &&
+ (abs(NNEst - RRBuffer[0]) > abs(NVEst - RRBuffer[0])) &&
+ RRMatch(NVEst,RRBuffer[0]))
+ {
+ RRTypes[0] = NV ;
+ return(PVC) ;
+ }
+
+ // If equally close, or we don't have an NN or NV in the buffer,
+ // who knows what it is.
+
+ else
+ {
+ RRTypes[0] = QQ ;
+ return(UNKNOWN) ;
+ }
+ }
+
+ // Otherwise the previous interval must have been a VV
+
+ else
+ {
+
+ // Does this match previous VV.
+
+ if(RRMatch(RRBuffer[0],RRBuffer[1]))
+ {
+ RRTypes[0] = VV ;
+ return(PVC) ;
+ }
+
+ // If this doesn't match a previous VV interval, assume
+ // any new interval is recovery to Normal beat.
+
+ else
+ {
+ if(RRShort(RRBuffer[0],RRBuffer[1]))
+ {
+ RRTypes[0] = QQ ;
+ return(UNKNOWN) ;
+ }
+ else
+ {
+ RRTypes[0] = VN ;
+ return(NORMAL) ;
+ }
+ }
+ }
+ }
+
+
+/***********************************************************************
+ RRMatch() test whether two intervals are within 12.5% of their mean.
+************************************************************************/
+
+int RRMatch(int rr0,int rr1)
+ {
+ if(abs(rr0-rr1) < ((rr0+rr1)>>3))
+ return(1) ;
+ else return(0) ;
+ }
+
+/************************************************************************
+ RRShort() tests whether an interval is less than 75% of the previous
+ interval.
+*************************************************************************/
+
+int RRShort(int rr0, int rr1)
+ {
+ if(rr0 < rr1-(rr1>>2))
+ return(1) ;
+ else return(0) ;
+ }
+
+/*************************************************************************
+ IsBigeminy() allows external access to the bigeminy flag to check whether
+ a bigeminal rhythm is in progress.
+**************************************************************************/
+
+int IsBigeminy(void)
+ {
+ return(BigeminyFlag) ;
+ }
+
+/**************************************************************************
+ Check for short interval in very regular rhythm.
+**************************************************************************/
+
+int RRShort2(int *rrIntervals, int *rrTypes)
+ {
+ int rrMean = 0, i, nnCount ;
+
+ for(i = 1, nnCount = 0; (i < 7) && (nnCount < 4); ++i)
+ if(rrTypes[i] == NN)
+ {
+ ++nnCount ;
+ rrMean += rrIntervals[i] ;
+ }
+
+ // Return if there aren't at least 4 normal intervals.
+
+ if(nnCount != 4)
+ return(0) ;
+ rrMean >>= 2 ;
+
+
+ for(i = 1, nnCount = 0; (i < 7) && (nnCount < 4); ++i)
+ if(rrTypes[i] == NN)
+ {
+ if(abs(rrMean-rrIntervals[i]) > (rrMean>>4))
+ i = 10 ;
+ }
+
+ if((i < 9) && (rrIntervals[0] < (rrMean - (rrMean>>3))))
+ return(1) ;
+ else
+ return(0) ;
+ }
+
+int RRMatch2(int rr0,int rr1)
+ {
+ if(abs(rr0-rr1) < ((rr0+rr1)>>4))
+ return(1) ;
+ else return(0) ;
+ }
diff -r 1de72bdab0ef -r 8e2fbe131b53 rythmchk.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rythmchk.h Sun Jun 28 03:06:00 2015 +0000 @@ -0,0 +1,34 @@ +/***************************************************************************** +FILE: rythmchk.h +AUTHOR: Patrick S. Hamilton +REVISED: 9/25/2001 + ___________________________________________________________________________ + +rythmchk.h: Prototype definitions for rythmchk.cpp +Copywrite (C) 2001 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). +******************************************************************************/ + +// External prototypes for rythmchk.cpp + +void ResetRhythmChk(void) ; +int RhythmChk(int rr) ; +int IsBigeminy(void) ; \ No newline at end of file
diff -r 1de72bdab0ef -r 8e2fbe131b53 sampler.cpp
--- a/sampler.cpp Sun May 17 00:27:16 2015 +0000
+++ b/sampler.cpp Sun Jun 28 03:06:00 2015 +0000
@@ -1,19 +1,22 @@
#include "mbed.h"
-
+#include <qrsdet.h>
// Use a ticker (interrupt routine) to sample the ADC at a fast rate (is 200Khz the max?), and generate an average value
// Use another ticker to grab the average value at 1Khz, which is the maximum useful sampling rate for HRV data collection
// These need to divide evenly into each other
-#define SAMPLING_RATE_HZ 4000
-#define READ_OUT_RATE_HZ 50
-#define SAMPLES_PER_READOUT (SAMPLING_RATE_HZ/READ_OUT_RATE_HZ)
+// 3 means sample at 8x read-out rate - SAMPLE_RATE is defined in qrsdet.h
+#define SHIFT_FOR_READ_OUT 1
+#define ADC_RATE (SAMPLE_RATE<<SHIFT_FOR_READ_OUT)
+#define READ_OUT_RATE_HZ SAMPLE_RATE
+
+#define SAMPLES_PER_READOUT (1<<SHIFT_FOR_READ_OUT)
AnalogIn ECG(P0_1);
Ticker sampling_rate;
int NumReadings=0;
-uint32_t ReadingTotal=0;
+int32_t ReadingTotal=0;
void (*SamplingFunction)(uint32_t);
@@ -23,7 +26,7 @@
ReadingTotal += ECG.read_u16();
if (NumReadings == SAMPLES_PER_READOUT)
{
- SamplingFunction(ReadingTotal/NumReadings);
+ SamplingFunction(ReadingTotal>>SHIFT_FOR_READ_OUT);
ReadingTotal = 0;
NumReadings = 0;
}
@@ -32,5 +35,5 @@
void setup_sampler(void (*function)(uint32_t))
{
SamplingFunction = function;
- sampling_rate.attach(&ADC_read, 1.0/SAMPLING_RATE_HZ);
+ sampling_rate.attach(&ADC_read, 1.0/ADC_RATE);
}