基于STM32F4的心电监护仪

1. 处理板的选择

本系统为了更好的人机交互，采用4.3英寸触摸屏搭配开源图形库LVGL，一方面，将显示波形和数据与心电信号的收集和处理隔离开来，另一方面，系统运行界面如下图所示，便于交互和美观。整个界面主要包括菜单、返回、图表、数据栏、导联状态灯、打开心电采集按钮。

通过系统菜单LVGL的roller绘制控件，roller内部选项的事件以回调函数的形式调用。因为选择roller的选项，LVGL会回到选项值，所以我自己设计了函数指针数组注册回调函数，并通过数组调用选定的序号。代码如下：

void (*oper_fuc[4])();//函数指针数组 void Menuitem_Init(void) { 
             oper_fuc[0]=send_type_server;     oper_fuc[1]=Set_chart_div_line;     oper_fuc[2]=clear_step;  oper_fuc[3]=smooth_filter; } static void roller_event_handler(lv_obj_t * obj, lv_event_t event) { 
             static unsigned char count=0;   if(event==LV_EVENT_VALUE_CHANGED)     { 
                 count=
       
        lv_roller_get_selected
        (obj
        )
        ; 
        } 
        if
        (event
        ==
        LV_EVENT_CLICKED
        ) 
        { 
          oper_fuc
        [count
        ]
        (
        )
        ; 
        } 
        } 
       

roller_event_handler是选中roller中的事件函数，在事件函数里面来回调选项的处理函数。roller中总共写了4个选项，分别为
send_type选择发送类型（支持发送到本地显示或者串口发送给上位机）、set_div_line是否设置图表的等分线、
clear_step清除界面上的数据、
smooth_filter是否进行平滑滤波

本系统还设计导联状态指示灯，前面讨论过ADS1292R可以检测电联的脱落状态，因而这里用LVGL的led控件作为导联的状态指示，当检测到导联接入人体，led控件就会点亮。设计了红心周期性跳动，当检测到导联接入人体后，红心就会周期性跳动，当心电数据采样开始后，红心随着心率值的改变而跳动着。同时还设计了采样开始/停止按钮，可以随时暂停和开始采样心电信号。
除了以上看得见的设计之外，还创建了四个周期性的任务，任务优先级从高到低分别为：更新数据栏里的数据、更新导联状态、检查心电信号的纵轴范围、系统状态的检查。

心电信号本质上是测量人体体表的电信号，将电极通过一定的导联体系就可以记录到心电图，因而选择合适的电极是观察心电图至关重要的选择。在医学上常见的导联体系分别为标准12导联体系、Lewis导联、Fontaine导联、Cabrera导联、Nehb导联、frank导联、Mason-Likar导联等。标准12导联体系是医院所使用的，它由3个双极肢体导联、6个单极胸前导联、3个单极加压肢体导联所组成。
该系统的主要目的是实时检测心率和QRS宽度，因此选择的导联应该基于能观察心电中R波较大的原则。因而选择标准12导联中标准肢体导联I（见图左），或者Mason-Likar导联（见图右）。

人体的内阻很高，因而心电信号是一个高内阻且幅度很低的信号，如果处理不好就会造成心电信号的衰减，因此就需要从两方面解决：
（1）降低与电极的接触阻抗（2）提高采集电路的输入阻抗。

目前，市面上有三种电极，分别为湿电极、干电极和非接触式电极，这三种电极中湿电极的接触电阻最小，因而对于模拟前端的输入电阻不需要太大。湿电极主要由电极片、Ag/AgCl 涂层、导电胶等物质组成。 医学电极贴片与身体接触的是水凝胶（亲水化合物），“黑色”部分为Ag/AgCl，使用导电金属和导线与仪器连接，实物如图所示。

人体的心电信号是一种非平稳、非线性、随机性比较强的微弱生理信号,幅值约为毫伏(mV)级,频率在0.05-100Hz之间。心电信号的每一个心跳循环由一系列有规律的波形组成，它们分别是P波、QRS复合波和T波，而这些波形的起点、终点、波峰、波谷、以及间期分别记录着心脏活动状态的详细信息

心电信号各个波段的详细说明如下：

心电各个波段的功率谱如下：

心电信号的噪声分析如下：

读者想对心电信号进一步了解可以参考如下链接：http://www.mythbird.com/ecgxin-hao-te-zheng/。

系统先从硬件初始化开始，其中包括串口初始化、触摸屏初始化、外部SRAM初始化、ADS1292R初始化、LMT70初始化、LVGL心跳定时器初始化。
其次就是LVGL初始化，主要是一些主题和变量的初始化。然后创建系统的UI界面和一些定时的任务。
最后初始化心电数据缓存、 数字低通滤波器初始化、心率数据缓存初始化。
完成以上的初始化，系统便进入主循环，等待心电数据输入缓存中出现数据，随后开始滤波，将滤波之后的数据写入心电输出缓存中，然后轮询LVGL的任务和触摸屏扫描。就这样不停地循环。其中心电输入缓存中的数据是通过中断从ADS1292R的输出引脚中读取，而心电输出缓存则是原始数据经过低通处理后的数据，等待LVGL显示任务的到来并显示在触摸屏上。系统总体框图和软件框图如下所示

在前面讨论过心电信号频谱和噪声，因而要对心电信号进行滤波，为了同时实现心电信号的实时滤波和心电波形实时显示，所以有必要设计一个缓存区来解决这个难题。这里我打算用我自己设计的两个循环队列解决这个难题。

为了使得在滤波的时候，心电数据依然能够采集，设计两个循环队列，如上图所示，其中IN_Buffer和OUT_Buffer的每个矩形框表示25x4个字节的空间，这取决一次需要多少字节的数据滤波。这里一次滤波需要25个int型的数据，因而每个缓存需要25x4字节。图中的蓝色填充表示缓存区中填满了数据，每次读完数据之后都需要切换缓存区，且IN_Buffer和OUT_Buffer的读写操作相反，即IN_Buffer的读操作是OUT_Buffer的写操作，程序框图如下图所示。

图上所示的三个程序均是并行处理的，
程序1是通过外部中断的服务函数调用的，
程序2则是在UI画图程序里面通过定时器周期性的调用，
程序3则是在主程序中的滤波函数里面调用
程序1代码如下（ADS1292R采用中断方式读取数据）：

void EXTI9_5_IRQHandler(void)
{ 
        
    if(EXTI->IMR&EXTI_Line5 && ADS_DRDY==0)//数据接收中断
    { 
        
        ADS1292_Read_Data(ads1292_Cache);//数据存到9字节缓冲区
        Update_ECG_Data(ads1292_Cache);
        Cheack_lead_stata(ads1292_Cache);
		if(state_pcb.SampleStartFlag==true)
			WriteAdsInBuffer(ecg_info.ecg_data);//数据写入缓存区
    } 
	EXTI_ClearITPendingBit(EXTI_Line5);
}

程序2代码如下（LVGL的心跳在定时器中周期调用，同时程序2也在其中运行，主要从滤波后的数据缓存中取出数据进行波形显示）：

void Wave_show(void)
{ 
        
    int value=0;
    if(ReadEcgOutBuffer(&value)!=0) { 
        
        if(ecg_graph.send_type==GRAPH) { 
        
         ecg_graph.y_pose=Transf_EcgData_To_Vert(value,ecg_graph.sacle);
			  chart_add_data(ecg_graph.y_pose);
            set_data_into_heart_buff(ecg_graph.y_pose);
        } else if(ecg_graph.send_type==USART) { 
        
            //EcgSendByUart(value);
			printf("%d\r\n",(int)alg(value/200));
        }
    }
}
//定时器3中断服务程序
void TIM3_IRQHandler(void)
{ 
         	static u8 show_cnt=0;   		  			    
	if(TIM3->SR&TIM_IT_Update)//溢出中断
	{ 
        	show_cnt++;
		lv_tick_inc(1);//lvgl的1ms心跳
		if(show_cnt==3){ 
        
		show_cnt=0;
		 Wave_show();
		}
	}				   
	TIM3->SR = (uint16_t)~TIM_IT_Update;
}

程序3代码如下（在滤波函数中调用，用于承上启下，即从IN缓存中取出数据，滤波之后写入OUT缓存中）：

void arm_fir_f32_lp(void)
{ 
        
	float32_t *inputf32, *outputf32;
	if(ReadAdsInBuffer() && WriterEcgOutBuffer()){ 
        //指针定位成功
		/* 初始化输入输出缓存指针 */
		inputf32 = (float32_t *)InFifoDev.rp;
		outputf32 =(float32_t *)OutFifoDev.wp;
		
	  /* 实现FIR滤波 */
		arm_fir_f32(&S, inputf32, outputf32, BLOCK_SIZE);
		//my_memcpy(OutFifoDev.wp,InFifoDev.rp,BLOCK_SIZE*4);
		InFifoDev.state[InFifoDev.read_front]=Empty;
		InFifoDev.read_front=(InFifoDev.read_front+1)%PACK_NUM;//切换读缓存块
	
		OutFifoDev.state[OutFifoDev.writer_rear]=Full;
		OutFifoDev.writer_rear=(OutFifoDev.writer_rear+1)%PACK_NUM;//切换写缓存块
	}

}

关于缓存切换代码如下：

static void WriteAdsInBuffer(int date)
{ 
        
	static u8 cnt=0;
	if(InFifoDev.state[InFifoDev.writer_rear]==Empty){ 
        //缓存块可写
		InFifoDev.wp=&AdsInBuffer[InFifoDev.writer_rear*(BLOCK_SIZE)];//将写指针定位写缓存块
		InFifoDev.wp[cnt++]=date;
		if(cnt==BLOCK_SIZE){ 
        
			cnt=0;
			InFifoDev.state[InFifoDev.writer_rear]=Full;
			InFifoDev.writer_rear=(InFifoDev.writer_rear+1)%PACK_NUM;//切换写缓存块
			
		}
	}
}

//定位读指针
//成功则返回1，不成功则返回0
u8 ReadAdsInBuffer(void)
{ 
        
	if(InFifoDev.state[InFifoDev.read_front]==Full){ 
        //缓存块可读
		InFifoDev.rp=&AdsInBuffer[InFifoDev.read_front*(BLOCK_SIZE)];//将读指针定位读缓存块
		return 1;
	}
	return 0;
}

//定位读指针
u8 WriterEcgOutBuffer(void)
{ 
        
	if(OutFifoDev.state[OutFifoDev.writer_rear]==Empty){ 
        //缓存块可写
		OutFifoDev.wp=&EcgOutBuffer[OutFifoDev.writer_rear*(BLOCK_SIZE)];//将读指针定位读缓存块
		return 1;
	}
	return 0;
}

//成功则返回1，不成功则返回0
u8 ReadEcgOutBuffer(int32_t *p)
{ 
        
	static u8 cnt=0;
	if(OutFifoDev.state[OutFifoDev.read_front]==Full){ 
        //缓存块可读
		OutFifoDev.rp=&EcgOutBuffer[OutFifoDev.read_front*(BLOCK_SIZE)];//将写指针定位读缓存块
		*p=OutFifoDev.rp[cnt++];
		if(cnt==BLOCK_SIZE){ 
        
			cnt=0;
			OutFifoDev.state[OutFifoDev.read_front]=Empty;
			OutFifoDev.read_front=(OutFifoDev.read_front+1)%PACK_NUM;//切换写读缓存块
		}
		return 1;
	}
	return 0;
}

1. 工频噪声滤除

滤除工频噪声的数字滤波算法主要有经典滤波器、小波变换、自适应滤波。小波变换能将心电信号进行多层分解，可以使得心电信号与工频噪声分离，但是计算量大，所占用的中间变量也比较多，对于单片机来说，处理的速度也不够快，因而对于系统的实时性这一指标很难实现。自适应滤波能够自动跟踪工频噪声的改变，但是需要增加一个输入信号作为参考，因而增加了系统的复杂性。在前面也讨论过心电信号95%的能量都是集中在0~40Hz,而工频噪声则在50Hz左右，过渡带比较宽，因而可以选择截止频率为40Hz的低通滤波器。
该低通滤波器利用MATLAB的FDATOOL生成，只需要选择低通滤波器是FIR结构，选择Blackman-Harris窗函数，滤波器的阶数定为50，选择采样频率为250Hz,截止频率为40Hz,参数如下图所示：

然后利用FDATOOL生成的冲激响应的数组，选择ARM官方的DSP库，调用arm_fir_f32函数，既可以完成一次滤波。但是在这之前，需要调用arm_fir_init_f32进行初始化。
滤波器系数如下：

const float32_t fir32LP[NUM_TAPS] = { 
        
  -7.484454468902e-22,-3.269336712398e-06,-1.365915864079e-05,-5.014073980636e-06,
  6.804735231975e-05,0.0001662336497003,7.965197426322e-05,-0.0003784662837741,
  -0.0008928563387901,-0.0005280588787408, 0.001284875839485, 0.003225662215767,
     0.0022425431358,-0.003157084585057,-0.009028737319977,-0.007219934929014,
   0.006057868257093,  0.02144319498633,  0.01971312591228,-0.009448071870685,
   -0.04806332586811, -0.05291973061693,  0.01224382260678,   0.1388254178822,
     0.2663085232723,   0.3199984843521,   0.2663085232723,   0.1388254178822,
    0.01224382260678, -0.05291973061693, -0.04806332586811,-0.009448071870685,
    0.01971312591228,  0.02144319498633, 0.006057868257093,-0.007219934929014,
  -0.009028737319977,-0.003157084585057,   0.0022425431358, 0.003225662215767,
   0.001284875839485,-0.0005280588787408,-0.0008928563387901,-0.0003784662837741,
  7.965197426322e-05,0.0001662336497003,6.804735231975e-05,-5.014073980636e-06,
  -1.365915864079e-05,-3.269336712399e-06,-7.484454468902e-22
};

static float32_t firStateF32[BLOCK_SIZE + NUM_TAPS - 1];
arm_fir_instance_f32 S;
void arm_fir_Init(void)
{ 
        
	arm_fir_init_f32(&S, NUM_TAPS, (float32_t *)&fir32LP[0], &firStateF32[0], BLOCK_SIZE);
}

滤波函数见程序3（往上找）

2. 基线漂移

基线漂移与工频噪声不同，它是由于呼吸和电极滑动变化所异致的，其频率一般低于1Hz左右。常见对于基线漂移滤除的数字算法有高通滤波器、中值滤波、小波变换、形态学滤波、曲线拟合等，其中高通滤波器可能会对心电信号的ST波段产生影响，毕竟基线漂移的频率也在ST波段里面。曲线拟合对较大的基线漂移处理能力较弱，处理的效果与处理数据的长度成正相关，因而不适用实时处理的系统。小波变换计算量大，也不适用实时处理的系统。相比之下，形态学滤波对心电信号的基线漂移滤除效果更好，计算量也比中值滤波小。但是形态学滤波要求数据长度足够长，因而会改变前面的缓存结构，并且在本系统中并未太严重的基线漂移，系统的任务也比较多，多方面权衡之下，选择不处理基线漂移。

3. 肌电噪声的抑制

肌电噪声主要是由于人体肌肉颤抖导致体表的电位发生变化，这种噪声通过电极贴传导至心电模拟前端，并且这种噪声持续时间较短，使得ECG信号波形产生细小的波纹，这种噪声频率分布比较广，前面已经将心电信号通过截止频率为40Hz的低通滤波器，因而需要5点平滑滤波将细小的波纹滤除，为了不影响心电信号的实时处理，因而改进版的平滑滤波器代码如下：

/* * 滑动平均值滤波。 * 每调用一次，就加入一个新数据，并得到当前的滤波值。 */
float alg(float new_val)
{ 
        
    /* 用一个减法，就做了"丢弃最旧的数据，加入最新的数据"这一操作 */
    sum += (new_val - buf[pos]);
    buf[pos] = new_val;
    pos = (pos + 1) % MAX_COUNT;
    /* 个数不足时，cnt是实际个数，个数足够时，cnt最多也只是MAX_COUNT */
    pcnt += (pcnt < MAX_COUNT);
    return sum / MAX_COUNT;
}

心率和QRS宽度检测作为本系统的算法核心，有了心率值和QRS宽度值才能进一步判断常见的心律失常。心率基本上都是检测两个R波之间的时隙来计算的，常见检测R的算法主要有阈值法、模板法和语句描述法。
而本系统的心率和QRS宽度检测算法是在一起检测的，所采用的算法是幅度阈值检测和差分检测相结合，因为观察心电信号的R波，发现R波是具有窄的脉冲，且脉冲的幅度是心电信号最高的，因而采用幅度和一阶差分共同约束找到R波，同时在找R波的同时还可以估计出QRS的宽度，算法的框图如图

心率检测和QRS宽度检测算法是采用状态机的编程思想，通过R波幅度大且从Q到R一直递增，并且R波到S波的一阶差分值很大，从而将R波定位出来，检测两个R波之前的时间，然后通过如下公式就可以计算出心率：
$$HR=(60\ast SampleRate)/count$$

而QRS宽度则是由
$$QRS=QRScnt\ast 2.2\ast 1000/(SampleRate)$$

上式中的2.2是估计值，因为QRS_cnt是在检测到R波之后才开始计数，并且未到S波谷停止计数，观察QRS波，发现Q到R与R到S近似对称，因而采用2.2这个估计值，这也是实时检测的缺陷，检测的样本不多。
心率算法和QRS宽度检测代码如下：