想问下有没有做心电算法的呢？
问下那个QRS波检测与采样率有什么关系啊？？！程序是在ourdev网友给下载的，但是运行结果对于360hz的数据检测出的QRS波个数是差不多对的，没分钟六七十个，但是对于采样率为200hz的数据检测出的QRS波个数就15个，这个差距也太大了吧？？1那我想是不是与程序的时间间隔有关呢？？？程序如下：*****************************************************************************
FILE:picqrs.c
AUTHOR: Patrick S. Hamilton
___________________________________________________________________________
picqrs.c: A PIC based QRS detector.（PIC——图像文件存储格式）
Copywrite (C) 2002 Patrick S. Hamilton
这是运行QRS波检测的程序，使用3.6864MHZ的可编程中断控制器(PIC）。
心电图样本吻合通过RS232连接的PIC（19.2K波特率），QRS探测器通过RS232
(异步通信端口设置程序)连接传回侦测。这方便PIC基础探测器的检测，它的数据来自
MIT/BIH（良性颅内高压）心率失常数据库。一个PC程序可以从MIT/BIH数据库读取数据，
然后送给PC。从PIC接收的侦测与有注解的心跳位置相对比。实施一个有效的硬件心跳探
测器需要用来自ECG放大器的A到D输入代替连续输入数据。
我们用计算机控制系统(CCS)的脉冲编码调制(PCM)C编译器3.093版本执行测试这个程序
***********************************************************************/
#include <16F877.h>
#device PIC16F877 *=16
#fuses HS,NOWDT,NOPROTECT,NOBROWNOUT,NOLVP,PUT
#use delay(clock=3686400)
#use rs232(baud=19200, xmit=PIN_C6, rcv=PIN_C7)
// Prototypes(技术原型).
byte bgetc(void) ;
int SyncRx(int8 in, int16 *out) ;
void SendInt(int16 x) ;
signed int16 hpfilt( signed int16 datum, int init ) ;
signed int16 lpfilt( signed int16 datum ,int init) ;
signed int16 deriv1( signed int16 x0, int init ) ;
int16 mvwint(signed int16 datum, int init) ;
signed int16 PICQRSDet(signed int16, int init) ;
signed int16 Peak( signed int16 datum, int init ) ;
void UpdateQ(int16 newQ) ;
void UpdateRR(int16 newRR) ;
void UpdateN(int16 newN) ;
// 连续输入缓冲器
byte buffer;
byte next_in = 0;
byte next_out = 0;
#define bkbhit (next_in!=next_out)
// 时间间隔常数
#define MS80 16
#define MS95 19
#define MS150 30
#define MS200 40
#define MS360 72
#define MS450 90
#define MS1000200
#define MS1500300
#define WINDOW_WIDTH MS80
#define FILTER_DELAY 21 + MS200
/***********************************************
主函数
***********************************************/
main()
{
char c ;
signed int16 x ;
output_low(PIN_D1) ; // Turn of the buzzer（蜂音器）
// 打开输入缓冲器的间断点
enable_interrupts(global);
enable_interrupts(int_rda);
// 初始化滤波器
PICQRSDet(0,1) ;
// 主循环
while(TRUE)
{
if(bkbhit)
{
c = bgetc() ;
if(SyncRx(c,&x) != 0)
{
output_high(PIN_D1) ;
x = PICQRSDet(x,0) ;
output_low(PIN_D1) ;
SendInt(x) ;
}
}
}
}
/***************************************************************************
syncrx处理通过RS232连接的通信同步输入收到的数据，它要求数据以21字节传输，
同步字节其次是1016位数值，先传输高位，再是低位。
当收到数据以后，syncrx返回1，数据在*out中返回
****************************************************************************/
#define SYNC_CHAR 0x55
#define RESET_CHAR 0xAA
#define FRAME_LGTH 10
int SyncRx(int8 in, int16 *out)
{
static int SyncState = 0, frameCount = 0 ;
static int16 Datum1 ;
int16 temp ;
// 等一个同步字(sync character)
if(SyncState == 0)
{
// 如果检测到一个同步字，开始接受一个frame
if(in == SYNC_CHAR)
{
SyncState = 1 ;
frameCount = 0 ;
}
// 如果检测到重置字，重置QRS探测器，开始接受一个frame
else if(in == RESET_CHAR)
{
SyncState = 1 ;
frameCount = 0 ;
PICQRSDet(0,1) ;
}
return(0) ;
}
// 保存最低有效位
else if(SyncState == 1)
{
Datum1 = in ;
Datum1 <<= 8 ;
SyncState = 2 ;
return(0) ;
}
// 将最低有效位和最高有效位放在一起，返回新值
else
{
temp = in ;
Datum1 |= temp ;
if(++frameCount == FRAME_LGTH)
SyncState = 0 ;
else
SyncState = 1 ;
*out = Datum1 ;
return(1) ;
}
return(0) ;
}
/**********************************************************************
* SendInt 传送一个16位整数，先是高位传送，再是低位
***********************************************************************/
void SendInt(int16 x)
{
int8 c ;
c = (x >> 8) & 0xFF ;
putc(c) ;
c = x & 0xFF ;
putc(c) ;
}
// QRS探测器的全值
int16 Q0 = 0, Q1 = 0, Q2 = 0, Q3 = 0, Q4 = 0, Q5 = 0, Q6 = 0, Q7 = 0 ;
int16 N0 = 0, N1 = 0, N2 = 0, N3 = 0, N4 = 0, N5 = 0, N6 = 0, N7 = 0 ;
int16 RR0=0, RR1=0, RR2=0, RR3=0, RR4=0, RR5=0, RR6=0, RR7=0 ;
int16 QSum = 0, NSum = 0, RRSum = 0 ;
int16 det_thresh, sbcount ;
int16 tempQSum, tempNSum, tempRRSum ;
int16 QN0=0, QN1=0 ;
int Reg0=0 ;
/******************************************************************************
PICQRSDet 用16位ECG样本(5 uV/LSB)作为输入，当检测到QRS波时，返回一个延迟。
经过init的一个非零值，重置QRS探测器
******************************************************************************/
signed int16 PICQRSDet(signed int16 x, int init)
{
static int16 tempPeak, initMax ;
static int8 preBlankCnt=0, qpkcnt=0, initBlank=0 ;
static int16 count, sbpeak, sbloc ;
int16 QrsDelay = 0 ;
int16 temp0, temp1 ;
if(init)
{
hpfilt(0,1) ;
lpfilt(0,1) ;
deriv1(0,1) ;
mvwint(0,1) ;
Peak(0,1) ;
qpkcnt = count = sbpeak = 0 ;
QSum = NSum = 0 ;
RRSum = MS1000<<3 ;
RR0=RR1=RR2=RR3=RR4=RR5=RR6=RR7=MS1000 ;
Q0=Q1=Q2=Q3=Q4=Q5=Q6=Q7=0 ;
N0=N1=N2=N3=N4=N5=N6=N7=0 ;
NSum = 0 ;
return(0) ;
}
x = lpfilt(x,0) ;
x = hpfilt(x,0) ;
x = deriv1(x,0) ;
if(x < 0) x = -x ;
x = mvwint(x,0) ;
x = Peak(x,0) ;
// 保持任一峰值200ms，除非更大的峰值来临。在任意200ms期间，只有唯一一个QRS波群
if(!x && !preBlankCnt)
x = 0 ;
else if(!x && preBlankCnt) // 如果我们已经抓住了一个200ms的峰值，继续，
{
if(--preBlankCnt == 0)
x = tempPeak ;
else x = 0 ;
}
else if(x && !preBlankCnt) // 如果200ms中无峰值，保存，开始计算
{
tempPeak = x ;
preBlankCnt = MS200 ;
x = 0 ;
}
else if(x) // 如果我们有一个峰值，但有更大的，保存，再开始计到200ms
{
if(x > tempPeak)
{
tempPeak = x ;
preBlankCnt = MS200 ;
x = 0 ;
}
else if(--preBlankCnt == 0)
x = tempPeak ;
else x = 0 ;
}
// 初始化检测到的首先8个最大峰值的QRS峰缓冲器
if( qpkcnt < 8 )
{
++count ;
if(x > 0) count = WINDOW_WIDTH ;
if(++initBlank == MS1000)
{
initBlank = 0 ;
UpdateQ(initMax) ;
initMax = 0 ;
++qpkcnt ;
if(qpkcnt == 8)
{
RRSum = MS1000<<3 ;
RR0=RR1=RR2=RR3=RR4=RR5=RR6=RR7=MS1000 ;
sbcount = MS1500+MS150 ;
}
}
if( x > initMax )
initMax = x ;
}
else
{
++count ;
//检查峰值是否在探测阈值之上
if(x > det_thresh)
{
UpdateQ(x) ;
// 更新RR间隔估计，寻找后边的限制
UpdateRR(count-WINDOW_WIDTH) ;
count=WINDOW_WIDTH ;
sbpeak = 0 ;
QrsDelay = WINDOW_WIDTH+FILTER_DELAY ;
}
// 如果峰值在探测阈值之下
else if(x != 0)
{
UpdateN(x) ;
QN1=QN0 ;
QN0=count ;
if((x > sbpeak) && ((count-WINDOW_WIDTH) >= MS360))
{
sbpeak = x ;
sbloc = count-WINDOW_WIDTH ;
}
}
// 测试寻找的后边的条件。如果QRS在后边的的条件中发现，更新QRS的buffer和det_thresh
if((count > sbcount) && (sbpeak > (det_thresh >> 1)))
{
UpdateQ(sbpeak) ;
// 更新RR间隔估计，寻找后边的限制
UpdateRR(sbloc) ;
QrsDelay = count = count - sbloc;
QrsDelay += FILTER_DELAY ;
sbpeak = 0 ;
}
}
return(QRSDelay) ;
}
/**************************************************************************
*UpdateQ 产生一个新的QRS峰值，更新QRS平均估计值和探测阈值
**************************************************************************/
void UpdateQ(int16 newQ)
{
QSum -= Q7 ;
Q7=Q6; Q6=Q5; Q5=Q4; Q4=Q3; Q3=Q2; Q2=Q1; Q1=Q0;
Q0=newQ ;
QSum += Q0 ;
det_thresh = QSum-NSum ;
det_thresh = NSum + (det_thresh>>1) - (det_thresh>>3) ;
det_thresh >>= 3 ;
}
/**************************************************************************
*UpdateN 产生一个新的噪声峰值，更新噪声平均值和探测阈值
**************************************************************************/
void UpdateN(int16 newN)
{
NSum -= N7 ;
N7=N6; N6=N5; N5=N4; N4=N3; N3=N2; N2=N1; N1=N0; N0=newN ;
NSum += N0 ;
det_thresh = QSum-NSum ;
det_thresh = NSum + (det_thresh>>1) - (det_thresh>>3) ;
det_thresh >>= 3 ;
}
/**************************************************************************
*UpdateRR 产生一个新的RR值，更新RR平均值
**************************************************************************/
void UpdateRR(int16 newRR)
{
RRSum -= RR7 ;
RR7=RR6; RR6=RR5; RR5=RR4; RR4=RR3; RR3=RR2; RR2=RR1; RR1=RR0 ;
RR0=newRR ;
RRSum += RR0 ;
sbcount=RRSum+(RRSum>>1) ;
sbcount >>= 3 ;
sbcount += WINDOW_WIDTH ;
}
/*************************************************************************
*lpfilt() 执行数字滤波器，由差分方程表示
* y = 2*y - y + x - 2*x + x
*
*记滤波器的延迟是5样本
*
**************************************************************************/
signed int16 lpfilt( signed int16 datum ,int init)
{
static signed int16 y1 = 0, y2 = 0 ;
static signed int16 d0,d1,d2,d3,d4,d5,d6,d7,d8,d9 ;
signed int16 y0 ;
signed int16 output ;
if(init)
{
d0=d1=d2=d3=d4=d5=d6=d7=d8=d9=0 ;
y1 = y2 = 0 ;
}
y0 = (y1 << 1) - y2 + datum - (d4<<1) + d9 ;
y2 = y1;
y1 = y0;
if(y0 >= 0) output = y0 >> 5;
else output = (y0 >> 5) | 0xF800 ;
d9=d8 ;
d8=d7 ;
d7=d6 ;
d6=d5 ;
d5=d4 ;
d4=d3 ;
d3=d2 ;
d2=d1 ;
d1=d0 ;
d0=datum ;
return(output) ;
}
/******************************************************************************
*hpfilt() 执行高通滤波器，由差分方程：
* y = y + x - x
* z = x - y ;
*滤波器延迟是15.5样本
******************************************************************************/
#define HPBUFFER_LGTH 32
signed int16 hpfilt( signed int16 datum, int init )
{
static signed int16 y=0 ;
static signed int16 data ;
static int ptr = 0 ;
signed int16 z ;
int halfPtr ;
if(init)
{
for(ptr = 0; ptr < HPBUFFER_LGTH; ++ptr)
data = 0 ;
ptr = 0 ;
y = 0 ;
return(0) ;
}
y += datum - data;
halfPtr = ptr-(HPBUFFER_LGTH/2) ;
halfPtr &= 0x1F ;
z = data ; // 补偿计算机转换故障（CCS shift bug）
if(y >= 0) z -= (y>>5) ;
else z -= (y>>5)|0xF800 ;
data = datum ;
ptr = (ptr+1) & 0x1F ;
return( z );
}
/*****************************************************************************
*deriv1 and deriv2 产生衍生的近似值，由差分方程
* y = 2*x + x - x - 2*x
*
*滤波器有2样本的延迟
*****************************************************************************/
signed int16 deriv1( signed int16 x0, int init )
{
static signed int16 x1, x2, x3, x4 ;
signed int16 output;
if(init)
x1 = x2 = x3 = x4 = 0 ;
output = x1-x3 ;
if(output < 0) output = (output>>1) | 0x8000 ; // 补偿 shift bug
else output >>= 1 ;
output += (x0-x4) ;
if(output < 0) output = (output>>1) | 0x8000 ;
else output >>= 1 ;
x4 = x3 ;
x3 = x2 ;
x2 = x1 ;
x1 = x0 ;
return(output);
}
/*****************************************************************************
* mvwint() 执行一个移动窗口积分，平均上一次16个信号值
******************************************************************************/
int16 mvwint(signed int16 datum, int init)
{
static unsigned int16 sum = 0 ;
static unsigned int d0,d1,d2,d3,d4,d5,d6,d7,d8,d9,d10,d11,d12,d13,d14,d15 ;
if(init)
{
d0=d1=d2=d3=d4=d5=d6=d7=d8=d9=d10=d11=d12=d13=d14=d15=0 ;
sum = 0 ;
}
sum -= d15 ;
d15=d14 ;
d14=d13 ;
d13=d12 ;
d12=d11 ;
d11=d10 ;
d10=d9 ;
d9=d8 ;
d8=d7 ;
d7=d6 ;
d6=d5 ;
d5=d4 ;
d4=d3 ;
d3=d2 ;
d2=d1 ;
d1=d0 ;
if(datum >= 0x0400) d0 = 0x03ff ;
else d0 = (datum>>2) ;
sum += d0 ;
return(sum>>2) ;
}
/**************************************************************
* peak() 用一个datum作为输入，返回一个峰高，当信号返回到半个峰高或者
从峰高检测到后已经经过95ms
**************************************************************/
signed int16 Peak( signed int16 datum, int init )
{
static signed int16 max = 0, lastDatum ;
static int timeSinceMax = 0 ;
signed int16 pk = 0 ;
if(init)
{
max = 0 ;
timeSinceMax = 0 ;
return(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 ;
}
else if(timeSinceMax > MS95)
{
pk = max ;
max = 0 ;
timeSinceMax = 0 ;
}
lastDatum = datum ;
return(pk) ;
}
/***********************************************************************
* 连续输入中断服务程序，输入被保存在FIFO （fist in fist out)
************************************************************************/
#int_rda
void serial_isr()
{
buffer=getc() ;
next_in=(next_in+1) & 0x3F;
}
/************************************************************************
* bgetc() (buffered getc()) 在连续输入的FIFO中返回下一个值
*************************************************************************/
byte bgetc()
{
byte c;
c=buffer;
next_out=(next_out+1) & 0x3F;
return(c);
} 楼主的R波怎么识别的？
页:
[1]