OSEA中QRS波检测算法代码分析

Wesley13
• 阅读 122

最近一直在搞R波检测算法,对OSEA代码主要是对注释做一个翻译,增加注释,使代码更容易理解。

一、首先看QRSDE.H

/*****************************************************************************
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    200    /* Sample rate in Hz. */
#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 过滤器的延迟加200毫秒消隐延迟
#define DER_DELAY    WINDOW_WIDTH + FILTER_DELAY + MS100

二、QRSDET2.CPP

/*****************************************************************************
FILE:  qrsdet2.cpp
AUTHOR:    Patrick S. Hamilton
REVISED:    7/08/2002
  ___________________________________________________________________________

qrsdet2.cpp: A QRS detector.
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 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 <mem.h>        /* For memmov. */
#include <math.h>
#include "qrsdet.h"

#define PRE_BLANK    MS195
#define MIN_PEAK_AMP    7 // Prevents detections of peaks smaller than 150 uV.

// External Prototypes.

int QRSFilter(int datum, int init) ;
int deriv1( int x0, int init ) ;

// Local Prototypes.

int Peak( int datum, int init ) ;
int mean(int *array, int datnum) ;
int thresh(int qmean, int nmean) ;
int BLSCheck(int *dBuf,int dbPtr,int *maxder) ;

double TH = .3125 ;

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 nmean, qmean, rrmean ;
    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) ;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
 = 0 ;if(aPeak && !preBlankCnt) // If there has been no peak for 200 ms,save this one and start counting.
{ // 如果在一个波峰之后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
 = 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
 = aPeak ;preBlankCnt = PRE_BLANK ; // MS200}else if(--preBlankCnt == 0)newPeak = tempPeak ;}/* Save derivative of raw signal for T-wave and baseline shift discrimination.     * 保存原始信号的导数用于识别T波和基线漂移 */

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){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.
 */{++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 ;qmean = mean(qrsbuf,8) ;det_thresh = thresh(qmean,nmean) ;memmove(&rrbuf[1], rrbuf, MEMMOVELEN)
 ;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
 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
 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 ;qmean = mean(qrsbuf,8) ;det_thresh = thresh(qmean,nmean) ;memmove(&rrbuf[1],rrbuf,MEMMOVELEN) ;rrbuf[0]
 = 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
 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 ;}qmean = mean( rsetBuff, 8 ) ;nmean = 0 ;rrmean = MS1000 ;sbcount = MS1500+MS150 ;det_thresh = thresh(qmean,nmean) ;initBlank = initMax = rsetCount = 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) ;}/********************************************************************mean returns the mean
 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
 += array[i] ;sum /= datnum ;return(sum) ;}/**************************************************************************** thresh() calculates the detection threshold from the qrs mean and noise mean estimates.****************************************************************************/int
 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()
 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
 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
 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) ;}
点赞
收藏
评论区
推荐文章
光头强的博客 光头强的博客
2个月前
Java面向对象试题
1、 请创建一个Animal动物类,要求有方法eat()方法,方法输出一条语句“吃东西”。 创建一个接口A,接口里有一个抽象方法fly()。创建一个Bird类继承Animal类并实现 接口A里的方法输出一条有语句“鸟儿飞翔”,重写eat()方法输出一条语句“鸟儿 吃虫”。在Test类中向上转型创建b对象,调用eat方法。然后向下转型调用eat()方
Jacquelyn38 Jacquelyn38
1年前
2020年前端实用代码段,为你的工作保驾护航
有空的时候,自己总结了几个代码段,在开发中也经常使用,谢谢。 1、使用解构获取json数据let jsonData   id: 1, status: "OK", data: ['a', 'b'] ; let  id, status, data: number   jsonData; console.log(id, status, number )
刚刚好 刚刚好
2个月前
css问题
1、 在IOS中图片不显示(给图片加了圆角或者img没有父级) <div<img src""/</div div {width: 20px; height: 20px; borderradius: 20px; overflow: h
blmius blmius
1年前
MySQL:[Err] 1292 - Incorrect datetime value: ‘0000-00-00 00:00:00‘ for column ‘CREATE_TIME‘ at row 1
文章目录 问题 用navicat导入数据时,报错: 原因这是因为当前的MySQL不支持datetime为0的情况。 解决修改sql\mode: sql\mode:SQL Mode定义了MySQL应支持的SQL语法、数据校验等,这样可以更容易地在不同的环境中使用MySQL。 全局s
小森森 小森森
2个月前
校园表白墙微信小程序V1.0 SayLove -基于微信云开发-一键快速搭建,开箱即用
后续会继续更新,敬请期待2.0全新版本 欢迎添加左边的微信一起探讨!项目地址:](https://www.aliyun.com/activity/daily/bestoffer?userCodesskuuw5n) \2. Bug修复更新日历 2. 情侣脸功能大家不要使用了,现在阿里云的接口已经要收费了(土豪请随意), \ \ 和 注意
晴空闲云 晴空闲云
2个月前
css中box-sizing解放盒子实际宽高计算
我们知道传统的盒子模型,如果增加内边距padding和边框border,那么会撑大整个盒子,造成盒子的宽度不好计算,在实务中特别不方便。boxsizing可以设置盒模型的方式,可以很好的设置固定宽高的盒模型。 盒子宽高计算假如我们设置如下盒子:宽度和高度均为200px,那么这会这个盒子实际的宽高就都是200px。但是当我们设置这个盒子的边框和内间距的时候,那
艾木酱 艾木酱
1个月前
快速入门|使用MemFire Cloud构建React Native应用程序
> MemFire Cloud是一款提供云数据库,用户可以创建云数据库,并对数据库进行管理,还可以对数据库进行备份操作。它还提供后端即服务,用户可以在1分钟内新建一个应用,使用自动生成的API和SDK,访问云数据库、对象存储、用户认证与授权等功能,可专
Stella981 Stella981
1年前
Angular material mat
Icon Icon Name mat-icon code _add\_comment_ add comment icon <mat-icon> add\_comment</mat-icon> _attach\_file_ attach file icon <mat-icon> attach\_file</mat-icon> _attach\
helloworld_28799839 helloworld_28799839
2个月前
常用知识整理
# Javascript ## 判断对象是否为空 ```js Object.keys(myObject).length === 0 ``` ## 经常使用的三元运算 > 我们经常遇到处理表格列状态字段如 `status` 的时候可以用到 ``` vue
helloworld_34035044 helloworld_34035044
4个月前
皕杰报表之UUID
​在我们用皕杰报表工具设计填报报表时,如何在新增行里自动增加id呢?能新增整数排序id吗?目前可以在新增行里自动增加id,但只能用uuid函数增加UUID编码,不能新增整数排序id。 uuid函数说明:获取一个UUID,可以在填报表中用来创建数据ID语法:uuid() 或 uuid(sep)参数说明:sep 布尔值,生成的uuid中是否包含分隔符'',缺省为