1．引言心脏是血液循环的动力器官。

心肌细胞的任何活动，都伴随着电的变化，这是一种生物电。

把特制的、有放大装置的电流计连接到体表，就可将每一心动周期内所发生的电位变化描记成连续的曲线，即心电图（简称ECG）。

由于各种病理原因引起的心脏疾病，几乎都和心脏的生物电活动相关，因此，心电图反映出心血管病人的许多病变信息，所以，它是心血管疾病诊断中十分重要的一种方法。

早期的ECG分析完全由医生用人工的方法完成。

这一过程不仅费时费力，且可靠性不高。

计算机辅助的ECG分析与诊断系统的研究始于五十年代末，在计算机辅助的ECG分析与诊断系统中，心电图中常存在由于各种干扰而造成的心电图的改变，这种改变称为心电图伪差。

伪差给心电图诊断带来一定的困难。

所以，从带有伪差的实际心电图中正确检测出我们需要的信息是很多科研工作者愿意研究的课题。

随着生活水平的提高，人们对健康的重视程度也愈来愈强。

心血管疾病是现代人患病率最高的疾病之一。

心电图能够反映出心血管病患者不少的病变信息，所以，对心电图的研究具有很重要的意义。

在心电图中，每一个周期波形代表一个心动周期，它是由以下各个波和时间段构成的（图1-1）：图1-1QRS波群：反映心室肌除极和最早复极过程的电位和时间的变化，但以心室肌除极化为主。

P波：反映心房肌除极过程的电位与时间的变化。

P—R间期：代表激动从窦房结通过心房、心室交界区到心室开始除极的时间。

S—T 间期：从QRS波群终点到T波起点间的线段。

它反映心室肌早期复极化过程的电位及时间变化。

T波：反映心室肌晚期复极化过程的电位与时间的变化。

Q—T 间期：从QRS波群起点到T波终点间的时间。

代表心室肌除极化与复极化的时间。

当心脏有病变时，将使相应的心电波形有所改变。

例如，QRS波群电压增高主要原因是心室肥大，S—T波段抬高有可能是心肌梗死，T波倒置有可能是心肌缺血等。

本设计中应用的标准心电信号ECG_X1是由UW DigiScope软件产生的，并以文本文件的形式存于Matlab的Work文件夹中。

而各种噪声是用Matlab编写程序添加的。

程序如下：clear %清除内存中的变量和函数clc %清屏幕fs=150; %设置取样频率N = 512; %设置取样点数load ECG_X1.txt %调出由UW DigiScope软件生成的标准心电图数据x=( ECG_X1/256)'; %归一化f=fs/N*(0:N/2-1); %设置频谱分辨率k=0:N-1; %设置离散频率变量z1=0.2*sin(2*pi*50*k/fs); %设置50HZ工频噪声z2=0.2*sin(2*pi*49.5*k/fs); %设置频率偏移50HZ工频噪声x1=x+z1; %合成含工频噪声心电信号x2=x+z2;z3=0.1*randn(1,N); %设置随机噪声50HZ工频噪声y1=x+z3; %合成含随机噪声心电信号save ECG_xinhao x x1 x2 y1 fs N z3 %将正常信号和含噪信号的参数存入ECG_xinhao文件中以上程序中，ecgy.txt文件是由UW DigiScope产生的标准心电图数据，x1是含50Hz工频噪声的心电信号，y1是含随机噪声的心电信号，save ECG_xinhaox x1 x2 y1 fs N z3是把工作表中的一些重要参数存入ECG_xinhao文件中，以便其它程序能调用这些参数。

2．对ECG信号噪声干扰的估计和消除从人体体表电极采集的常规心电图信号(ECG)是mV级的信号，其主要能量集中在0.05Hz 到100Hz之间，是信噪比很低的微弱信号，很容易受到周围环境的影响。

主要包括：电力系统的工频干扰，主要由50Hz及其谐波构成；由于测量电极同人体的接触不良引入的电极接触噪声包括电极同人体接触的时断时续及人体轻微运动引起的电极接触电阻的变化，这类干扰往往表现为信号基线的跳跃及其回复至基线的指数型的衰减过程，或者基线的抖动；由于人体呼吸引起的基线飘移，这类干扰的频率约为0.15Hz到0.3Hz；由于人体肌电信号（肌肉纤维的动作电位）引入的干扰，这是一种快速和随机的电变化，它的频带影响范围可以从2Hz到10Hz；此外还包括电极的极化电压，以及心电图放大电路内部噪声等。

相对于微弱的心电信号来说，环境干扰可能大几个数量级，因此为了增强心电信号的有效成分，抑制噪声，要对采集的心电信号进行数字滤波处理，以抑制噪声干扰，提高波形检测及参数分析的准确性。

2.1几种消除50Hz工频干扰的滤波方法2.1.1有限冲激响应滤波器有限冲激响应滤波器(FIR)的单位冲激响应为有限项,这意味着输出数据的计算中不含反馈项,其输出仅取决于当前的输入和过去的输入。

这种性质对数字滤波设计和应用具有重要的意义。

FIR滤波器还具有线性相位，因为在生物医学信号处理的许多应用中，通过滤波来保存信号的某些特性是重要的，例如QRS波的高度和间期。

由于线性相位滤波器的相位响应是一纯延时，因此相位失真最小，使得FIR滤波器很容易被设计为具有线性相位特性的滤波器。

FIR滤波器还有较强的稳定性。

因为滤波器未采用反馈回路，故除了那些位于z=0的极点外再没有别的极点。

不可能有极点位于单位圆之外，这意味着它是稳定的。

只要滤波器的输入是有限的，其输出也必将是有限的。

以上所有特性都使得FIR滤波器的设计变得简单。

有许多的直接设计方法能够满足任意频率响应和相位响应指标，例如窗设计法或频率取样法。

有许多软件包能够自动进行FIR设计过程，而且往往能求得在某方面最优的滤波器。

1.三点平均滤波器理论：由于心电信号的主要频率都集中在低频段，所以我们可以设计一个低通滤波器，以去除工频干扰。

我们先设计一个非常简单的低通滤波器，其输出和输入的关系为：y(n)=1/3[x(n)+x(n-1)+x(n-2)]因为当输入信号x(n)=δ(n)时，系统的输出，即h(n)=T[δ(n)]，所以y(n)就是h(n)，即h(n)=1/3[δ(n)+ δ(n-1)+ δ(n-2)]其z变换的形式为：H(z)=1/3(1+z-1+z-2)程序如下：% 三点平均滤波器clearclcload ECG_ xinhao x x1 fs N %调用心电信号参数t=(0:length(x)-1)/fs; %设置取样长度X=fft(x,N); %以512点取样点数对正常信号进行快速傅立叶变换f=fs/N*(0:N/2-1); %设置频谱分辨率figure(1) %建立图形窗口subplot(211);plot(t,x);grid;title('标准心电图时域图');subplot(212);plot(f,abs(X(1:N/2)));grid;title('标准心电图频域图');X1=fft(x1,N); %以512点取样点数对含噪信号进行快速傅立叶变换figure(2)subplot(211);plot(t,x1);grid;title('含噪心电图时域图');subplot(212);plot(f,abs(X1(1:N/2)));grid;title('含噪心电图频域图');h=1/3*[1 1 1]; %滤波器的单位脉冲响应[H,w]=freqz(h,1,512); %求h的频响y=conv(x1,h); %求含噪信号经过滤波器的输出t1=(0:length(y)-1)/fs; %设置输出信号的取样长度figure(3)subplot(211);plot(w*fs/(2*pi),abs(H));grid;title('滤波器的幅频特性');subplot(212);plot(t1,y);grid;title('滤波后心电信号时域图');在此程序中，load ECG_A x x1 fs N是调出工作表中的ECG_A.mat文件。

y=conv(x1,h)是把含噪信号和滤波器的单位脉冲响应作卷积，其值就是含噪信号经过滤波器的输出。

程序中分别对原始标准心电信号、含噪心电信号的时频域以及滤波器的幅频特性和滤波后心电信号作了图。

（如图2-1、图2-2、图2-3）图2-1图2-2⎩⎨⎧≤≤=-πωωωωαωω c c j j d e e H ,0,)()](sin[)sin()]([)(αωωω-=⎤⎡==n n e H IFT n h c c j图 2-3从图中我们可以看出人为添加上去的50Hz 工频被较好的抑制了。

而且由于滤波器比较简单，故运行速度较快能跟上采样率从而实现实时运行。

这是此类滤波器的优点。

但也因为它简单，所以其幅频响应特性很软，过渡带较大。

这样不仅把要消除的噪声去掉了，一些在过渡带中的有用的信号频率成分也被抑制了。

使心电信号的主要波形QRS 波有较大的削峰，信号衰减较大。

这从滤波后心电信号时域图里可以看出，信号的幅度被衰减了。

2.利用窗函数设计FIR 带阻陷波器理论：FIR 滤波器设计就是在给定技术指标下，求满足这一指标的滤波器单位脉冲响应h(n)。

我们先以设计低通滤波器来说明，而其它滤波器可以通过不同低通滤波器组合得到。

首先选取理想低通滤波器，它在通带上具有单位增益和线性相位，在阻带上具有零响应即式中，ωc 为截止频率，a 为取样时移。

由DFT 时移性质因子e -j αω意味着h d (n)向右移位a 个取样间隔。

得⎩⎨⎧-≤≤=值其他的n N n n h n h d ,0;10),()(2/)1(-=N α10),()()(-≤≤=N n n n h n h d ω这是一个以a 为中心的偶对称无限长非因果序列，故还不是可实现的FIR 滤波器。

解决的方法是截取(即加窗)。

即将h d (n)变为有限长因果序列，截取它的一段来代替它，即将n=0到n=N-1的一段序列截取下来作为 h(n)。

为了得到线性相位，h(n)必须是偶对称的，所以时移a 应取h(n)长度的一半。

即这种截取相当于加矩形窗处理。

即则这样，线性相位因果FIR 滤波器，就是用理想滤波器的单位脉冲响应h d (n)和窗函数ω(n)相乘求得。

加窗处理实际上对h d (n)作了修改，因此h(n)的频响特性H(e jw ) 自然不同于理想滤波器的频响H d (e jw )。

窗函数对信号的影响主要表现在：1）在截止频率ωc 处产生过渡带，过渡带宽为w p -w s2)在过渡带两旁的通带和阻带都产生波动（称为吉布斯现象）。

不同的窗函数对信号的影响也不同，对许多应用来说，哈明窗是较好的选择。

对于带阻滤波器来说就是用一个全通滤波器减去一个带通滤波器构成。

以下是用哈明窗函数构 成的FIR 滤波器,程序如下：%FIR 带阻滤波器(hamming 窗)clear %清除内存中的变量和函数clc %清屏幕load ECG_ xinhao x x1 fs N %调用心电信号参数f=fs/N*(0:N/2-1); %设置频谱分辨率⎩⎨⎧-≤≤=值其它的n N n n ,0;10,1)(ωwc = (ws+wp)/2; %确定截止频率k=0:N-1; %设置窗口长度m=k-(N-1)/2+eps;hz=sin(pi*m)./(pi*m)-(sin(wc(:,2)*m)./(pi*m)-sin(wc(:,1)*m)./(pi*m)); %理想带阻滤波器单位脉冲响应w_ham = (hamming(N))'; %选择Hamming窗h = hz.* w_ham; %求设计的滤波器的h[H,w]=freqz(h,1); %求h的频响mag=abs(H); %求幅频特性db=20*log10((mag+eps)/max(mag)); %化幅频特性以dB为单位y=conv(x1,h); %求含噪信号经过滤波器的输出Y=fft(y,N); %以512点取样点数对输出信号进行快速傅立叶变换t1=(0:length(y)-1)/fs; %设置输出信号取样长度figure(1);subplot(211); plot(k,h);grid;axis([0 N-1 -0.1 0.3]);title('带阻滤波器单位脉冲响应');subplot(212);plot(w*fs/(2*pi),db);grid;axis([0 fs/2 -100 20]);title('带阻滤波器幅频特性');figure(2);subplot(211);plot(t1,y);grid,axis([1.81 4.98 -0.7 0.7]);title('滤波后心电图时域图');subplot(212);plot(f,abs(Y(1:N/2)));grid;title('滤波后心电图频域图');在此程序中，wp =[2*pi*48/fs 2*pi*52/fs]; ws = [2*pi*49/fs 2*pi*51/fs];是设置通阻带边缘，通带边缘为48~52Hz，阻带边缘为49~51Hz。