最近一直在搞R波检测算法，对OSEA代码主要是对注释做一个翻译，增加注释，使代码更容易理解。
一、首先看QRSDE.H
1/***************************************************************************** 2FILE: qrsdet.h 3AUTHOR: Patrick S. Hamilton 4REVISED: 4/16/2002 5 ___________________________________________________________________________ 6 7qrsdet.h QRS detector parameter definitions 8Copywrite (C) 2000 Patrick S. Hamilton 9 10This file is free software; you can redistribute it and/or modify it under 11the terms of the GNU Library General Public License as published by the Free 12Software Foundation; either version 2 of the License, or (at your option) any 13later version. 14 15This software is distributed in the hope that it will be useful, but WITHOUT ANY 16WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A 17PARTICULAR PURPOSE. See the GNU Library General Public License for more 18details. 19 20You should have received a copy of the GNU Library General Public License along 21with this library; if not, write to the Free Software Foundation, Inc., 59 22Temple Place - Suite 330, Boston, MA 02111-1307, USA. 23 24You may contact the author by e-mail (pat@eplimited.com) or postal mail 25(Patrick Hamilton, E.P. Limited, 35 Medford St., Suite 204 Somerville, 26MA 02143 USA). For updates to this software, please visit our website 27(http://www.eplimited.com). 28 __________________________________________________________________________ 29 Revisions: 30 4/16: Modified to allow simplified modification of digital filters in 31 qrsfilt(). 32*****************************************************************************/ 33 34 35#define SAMPLE_RATE 200 /* Sample rate in Hz. */ 36#define MS_PER_SAMPLE ( (double) 1000/ (double) SAMPLE_RATE) 37#define MS10 ((int) (10/ MS_PER_SAMPLE + 0.5)) 38#define MS25 ((int) (25/MS_PER_SAMPLE + 0.5)) 39#define MS30 ((int) (30/MS_PER_SAMPLE + 0.5)) 40#define MS80 ((int) (80/MS_PER_SAMPLE + 0.5)) 41#define MS95 ((int) (95/MS_PER_SAMPLE + 0.5)) 42#define MS100 ((int) (100/MS_PER_SAMPLE + 0.5)) 43#define MS125 ((int) (125/MS_PER_SAMPLE + 0.5)) 44#define MS150 ((int) (150/MS_PER_SAMPLE + 0.5)) 45#define MS160 ((int) (160/MS_PER_SAMPLE + 0.5)) 46#define MS175 ((int) (175/MS_PER_SAMPLE + 0.5)) 47#define MS195 ((int) (195/MS_PER_SAMPLE + 0.5)) 48#define MS200 ((int) (200/MS_PER_SAMPLE + 0.5)) 49#define MS220 ((int) (220/MS_PER_SAMPLE + 0.5)) 50#define MS250 ((int) (250/MS_PER_SAMPLE + 0.5)) 51#define MS300 ((int) (300/MS_PER_SAMPLE + 0.5)) 52#define MS360 ((int) (360/MS_PER_SAMPLE + 0.5)) 53#define MS450 ((int) (450/MS_PER_SAMPLE + 0.5)) 54#define MS1000 SAMPLE_RATE 55#define MS1500 ((int) (1500/MS_PER_SAMPLE)) 56#define DERIV_LENGTH MS10 57#define LPBUFFER_LGTH ((int) (2*MS25)) 58#define HPBUFFER_LGTH MS125 59 60#define WINDOW_WIDTH MS80 // Moving window integration width. 61#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 过滤器的延迟加200毫秒消隐延迟 62#define DER_DELAY WINDOW_WIDTH + FILTER_DELAY + MS100
二、QRSDET2.CPP
1/***************************************************************************** 2FILE: qrsdet2.cpp 3AUTHOR: Patrick S. Hamilton 4REVISED: 7/08/2002 5 ___________________________________________________________________________ 6 7qrsdet2.cpp: A QRS detector. 8Copywrite (C) 2002 Patrick S. Hamilton 9 10This file is free software; you can redistribute it and/or modify it under 11the terms of the GNU Library General Public License as published by the Free 12Software Foundation; either version 2 of the License, or (at your option) any 13later version. 14 15This software is distributed in the hope that it will be useful, but WITHOUT ANY 16WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A 17PARTICULAR PURPOSE. See the GNU Library General Public License for more 18details. 19 20You should have received a copy of the GNU Library General Public License along 21with this library; if not, write to the Free Software Foundation, Inc., 59 22Temple Place - Suite 330, Boston, MA 02111-1307, USA. 23 24You may contact the author by e-mail (pat@eplimited.edu) or postal mail 25(Patrick Hamilton, E.P. Limited, 35 Medford St., Suite 204 Somerville, 26MA 02143 USA). For updates to this software, please visit our website 27(http://www.eplimited.com). 28 __________________________________________________________________________ 29 30This file contains functions for detecting QRS complexes in an ECG. The 31QRS detector requires filter functions in qrsfilt.cpp and parameter 32definitions in qrsdet.h. QRSDet is the only function that needs to be 33visable outside of these files. 34 35Syntax: 36 int QRSDet(int ecgSample, int init) ; 37 38Description: 39 QRSDet() implements a modified version of the QRS detection 40 algorithm described in: 41 42 Hamilton, Tompkins, W. J., "Quantitative investigation of QRS 43 detection rules using the MIT/BIH arrhythmia database", 44 IEEE Trans. Biomed. Eng., BME-33, pp. 1158-1165, 1987. 45 46 Consecutive ECG samples are passed to QRSDet. QRSDet was 47 designed for a 200 Hz sample rate. QRSDet contains a number 48 of static variables that it uses to adapt to different ECG 49 signals. These variables can be reset by passing any value 50 not equal to 0 in init. 51 52 Note: QRSDet() requires filters in QRSFilt.cpp 53 54Returns: 55 When a QRS complex is detected QRSDet returns the detection delay. 56 57****************************************************************/ 58 59#include <mem.h> /* For memmov. */ 60#include <math.h> 61#include "qrsdet.h" 62 63#define PRE_BLANK MS195 64#define MIN_PEAK_AMP 7 // Prevents detections of peaks smaller than 150 uV. 65 66// External Prototypes. 67 68int QRSFilter(int datum, int init) ; 69int deriv1( int x0, int init ) ; 70 71// Local Prototypes. 72 73int Peak( int datum, int init ) ; 74int mean(int *array, int datnum) ; 75int thresh(int qmean, int nmean) ; 76int BLSCheck(int *dBuf,int dbPtr,int *maxder) ; 77 78double TH = .3125 ; 79 80int DDBuffer[DER_DELAY], DDPtr ; /* Buffer holding derivative data. */ 81int Dly = 0 ; 82 83const int MEMMOVELEN = 7*sizeof(int); 84 85int QRSDet( int datum, int init ) 86 { 87 static int det_thresh, qpkcnt = 0 ; 88 static int qrsbuf[8], noise[8], rrbuf[8] ; 89 static int rsetBuff[8], rsetCount = 0 ; 90 static int nmean, qmean, rrmean ; 91 static int count, sbpeak = 0, sbloc, sbcount = MS1500 ; 92 static int maxder, lastmax ; 93 static int initBlank, initMax ; 94 static int preBlankCnt, tempPeak ; 95 96 int fdatum, QrsDelay = 0 ; 97 int i, newPeak, aPeak ; 98 99/* Initialize all buffers to 0 on the first call. */ 100 101 if( init ) 102 { 103 for(i = 0; i < 8; ++i) 104 { 105 noise[i] = 0 ; /* Initialize noise buffer */ 106 rrbuf[i] = MS1000 ;/* and R-to-R interval buffer. */ 107 } 108 109 qpkcnt = maxder = lastmax = count = sbpeak = 0 ; 110 initBlank = initMax = preBlankCnt = DDPtr = 0 ; 111 sbcount = MS1500 ; 112 QRSFilter(0,1) ; /* initialize filters. */ 113 Peak(0,1) ; 114 } 115 116 fdatum = QRSFilter(datum,0) ; /* Filter data. 滤波*/ 117 118 119 /* Wait until normal detector is ready before calling early detections. */ 120 /* 直到正常检测器准备好了之后开始早起检测。此处我理解的是等到心电波形稳定之后开始采集。 */aPeak = Peak(fdatum,0) ;if(aPeak < MIN_PEAK_AMP)aPeak = 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./**/保存200ms内所有的波峰，以防止后面出现更大的波峰，只所以这样，是因为我们认为前后200ms内只能有一个QRS波。*/newPeak 121 = 0 ;if(aPeak && !preBlankCnt) // If there has been no peak for 200 ms,save this one and start counting. 122{ // 如果在一个波峰之后200ms内没有波峰，则保存这个波峰并开始计数。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)//如果我们保存了一个波峰后过了200ms的检测依然没有检测到新的符合条件的波峰，newPeak 123 = tempPeak ; //则认为这个波峰为新的QRS波 }else if(aPeak) // If we were holding a peak, but{ // this ones bigger, save it andif(aPeak > tempPeak) // start counting to 200 ms again.{// 如果我们检测出一个波峰，并且比以前的最大波峰还大，则保存下来重新开始计数。tempPeak 124 = aPeak ;preBlankCnt = PRE_BLANK ; // MS200}else if(--preBlankCnt == 0)newPeak = tempPeak ;}/* Save derivative of raw signal for T-wave and baseline shift discrimination. * 保存原始信号的导数用于识别T波和基线漂移 */ 125 126DDBuffer[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 127 = 0 ;qrsbuf[qpkcnt] = initMax ;initMax = 0 ;++qpkcnt ;if(qpkcnt == 8){qmean = mean( qrsbuf, 8 ) ;nmean = 0 ;rrmean = MS1000 ;sbcount = MS1500+MS150 ;det_thresh = thresh(qmean,nmean) ;}}if( newPeak > initMax )initMax = newPeak ;}else /* Else test for a qrs. 128 */{++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 129 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 ;qmean = mean(qrsbuf,8) ;det_thresh = thresh(qmean,nmean) ;memmove(&rrbuf[1], rrbuf, MEMMOVELEN) 130 ;rrbuf[0] = count - WINDOW_WIDTH ;rrmean = mean(rrbuf,8) ;sbcount = rrmean + (rrmean >> 1) + WINDOW_WIDTH ;count = WINDOW_WIDTH ;sbpeak = 0 ;lastmax = maxder ;maxder = 0 ;QrsDelay = WINDOW_WIDTH + FILTER_DELAY ;initBlank = initMax = rsetCount = 0 ;}// If a 131 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 ;nmean = mean(noise,8) ;det_thresh = thresh(qmean,nmean) ;// Don't include early peaks (which might be 132 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 133 back update the QRS buffer and det_thresh. */if((count > sbcount) && (sbpeak > (det_thresh >> 1))){memmove(&qrsbuf[1],qrsbuf,MEMMOVELEN) ;qrsbuf[0] = sbpeak ;qmean = mean(qrsbuf,8) ;det_thresh = thresh(qmean,nmean) ;memmove(&rrbuf[1],rrbuf,MEMMOVELEN) ;rrbuf[0] 134 = sbloc ;rrmean = mean(rrbuf,8) ;sbcount = rrmean + (rrmean >> 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 135 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 136 == 8){for(i = 0; i < 8; ++i){qrsbuf[i] = rsetBuff[i] ;noise[i] = 0 ;}qmean = mean( rsetBuff, 8 ) ;nmean = 0 ;rrmean = MS1000 ;sbcount = MS1500+MS150 ;det_thresh = thresh(qmean,nmean) ;initBlank = initMax = rsetCount = 0 ;}}if( newPeak > initMax )initMax = 137 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 138 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 < 139 (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) ;}/********************************************************************mean returns the mean 140 of an array of integers. It uses a slowsort algorithm, but these arrays are small, so it hardly matters.********************************************************************/int mean(int *array, int datnum){long sum ;int i ;for(i = 0, sum = 0; i < datnum; ++i)sum 141 += array[i] ;sum /= datnum ;return(sum) ;}/**************************************************************************** thresh() calculates the detection threshold from the qrs mean and noise mean estimates.****************************************************************************/int 142 thresh(int qmean, int nmean){int thrsh, dmed ;double temp ;dmed = qmean - nmean ;/* thrsh = nmean + (dmed>>2) + (dmed>>3) + (dmed>>4); */temp = dmed ;temp *= TH ;dmed = temp ;thrsh = nmean + dmed ; /* dmed * THRESHOLD */return(thrsh) ;}/***********************************************************************BLSCheck() 143 reviews data to see if a baseline shift has occurred.This is done by looking for both positive and negative slopes ofroughly the same magnitude in a 220 ms window.***********************************************************************/int BLSCheck(int *dBuf,int 144 dbPtr,int *maxder){int max, min, maxt, mint, t, x ;max = min = 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 145 if a maximum and minimum pair are foundwhere the interval between them is less than 150 ms. */ if((max > (min>>3)) && (min > (max>>3)) &&(abs(maxt - mint) < MS150))return(0) ;elsereturn(1) ;}