ECG信号分析与处理系统设计.doc

ECG信号分析与处理系统设计

*****************
实践教学
*******************

某某理工大学
计算机与通信学院

2015年春季学期

信号处理课程设计

题 目： ECG信号分析与处理系统设计
专业班级： 通信工程
姓 名：
学 号：
指导教师：
成 绩：

摘要
系统的研究心电信号处理对疾病的早期预测及家庭医疗保健具有十分重要的意义，一直是生物医学工程领域的研究热点。心血管疾病是人类生命的最主要威胁之一，而心电(Electrocardiogram)，ECG信号是诊断心血管疾病的主要依据，心电信号是心脏电生理活动在体表的表现，提供了心脏功能等生理状况的有重要价值的临床医学信息，是临床心脏病诊断的基础。因此，设计心电信号处理系统具有重要意义。本论文综合运用数字信号处理的理论知识对心电信号进行分析与处理，实现ECG信号的频谱分析，基线漂移检测等，设计滤波器实现心电信号的滤波，滤去高频和低频干扰，实现ECG信号的增强。同时使用陷波器对50Hz的工频干扰进一步滤除，得到比较纯净的心电信号。
关键词: 心电信号，工频干扰，基线漂移

目 录
摘要 I
一、前言 1
二、心电信号 2
2.1 原始心电信号分析 2
2.2 心电信号中的噪声 3
2.3 系统总体设计框图 4
三、设计原理及方法 5
3.1 数字滤波器简介 5
3.2 IIR滤波器的设计原理 5
3.3 IIR滤波器的设计 5
3.3.1 IIR数字低通滤波器设计过程 5
3.3.2 IIR数字带通滤波器设计过程 9
3.4 FIR滤波器 10
3.4.1 FIR滤波器的设计 11
3.4.2 FIR数字低通滤波器设计过程 11
3.5 陷波器 13
3.5.1陷波器的基本原理及作用 13
3.5.2双T法设计陷波器 13
四、MATLAB简述 15
五、总结 16
参考文献 17
附录 18

一、前言
心电图(ECG)是用来捕捉心脏在一段时间内情况的反映，它通过外部电极连接到皮肤转换成电信号来采集。心脏外面形成的每个细胞膜都有一个关联电荷，它在每次心跳期间去极化。它以微小电信号的形式出现在皮肤上，可以通过心电图探测到并放大显示。早在1900年Willem Einthoven就发明了第一台实用的心电图。该系统很笨重，需要很多人去操纵它。病人需要把他的胳膊和腿放到含有电解液的大型电极中。今天的心电监护设备结构紧凑，携带方便，这样病人走动时也可以带着。家用十二导联心电图可以装在口袋里。目前，心电信号的采集与处理在医用方面也有了重要地位。同时，在处理信号时也存在着诸多问题。
在信号采集时，身体的任一微小运动都会产生“基线漂移”，这是一种低频干扰，同时，由于肌电的存在又产生了高频的肌电噪声，由于空间电磁场的存在又使心电信号中混有50Hz的工频干扰。这些噪声不去除，就会影响下一步的信号处理。综合运用数字信号处理的理论知识进行生物医学信号分析与处理，实现ECG信号的频谱分析，基线漂移检测等，设计滤波器实现心电信号的滤波，滤去高频和低频干扰，实现ECG信号的增强。信号处理是一项巨大的挑战，因为实际的信号为0.5MV，它处在一个300mv偏移量环境里。其他因素如交流电的干扰，外科设备的射频干扰，手术植入的的设备如起搏器和生理检测系统也会影响精度。心电图里噪声的主要来源是基线漂移（低频噪声）电力线干扰（来自电力线的50Hz或60Hz噪声）肌肉噪声（这种噪声是很难被清除，因为它是在同一地区的实际信号。它通常在软件里纠正。）其他干扰（例如，来自其他设备的射频噪声）。信号采集以后，存在许多软件算法来去除噪声。基线漂移也是目前存在的比较突出的问题，它是一种存在于心电图系统的低频噪声。是由于点击，呼吸和身体运动的偏置电压造成的。这可能会在分析心电图波形是造成问题。这种噪声可以通过使用硬件实现高通滤波。
本文主要介绍了关于几种噪声去除方法以及相应滤波器的介绍。其中着重介绍了IIR滤波器和FIR滤波器。

二、心电信号
2.1 原始心电信号分析
用load函数将原心电信号导入b = load('D:\Users\Data.txt')，并画出心电信号的时域波形和频谱图（幅频和相频），如图1所示:

图1 原始心电信号的时域波形图及频谱图
心电信号由于受到人体诸多因素的影响，因而有着一般信号所没有的特点：
(1)信号弱。心电信号是体表的电生理信号，一般比较微弱，幅度在10pV～5mV，频率范围在0.05-100Hz以内，而90％的ECG频谱能量集中0.25-35Hz之间，心电信号频率较低，大量的是直流成分，去掉直流，它的主要频率范围是0.05-100Hz，大部分能量集中在0.05-40Hz。
(2)噪声强。由于人体自身信号弱，加之人体又是一个复杂的系统，因此信号容易受到噪声干扰。
(3)随机性强。心电信号不仅是随机的，而且是非平稳的。同时，在心电图检测过程中极易受到各种噪声源的干扰，从而使图像质量变差，使均匀和连续变化的心电数值产生突变，在心电图上形成一些毛刺。使原本很微弱的信号很难和噪声进行分解。
2.2 心电信号中的噪声
人体心电信号是一种弱电信号，信噪比低。一般正常的心电信号频率范围为0.05-100 Hz，而90％的心电信号(ECG)频谱能量集中在0.25-35 Hz之间。采集一种电信号时，会受到各种噪声的干扰，噪声来源通常有下面几种：
(1)工频干扰
50 Hz工频干扰是由人体的分布电容所引起，工频干扰的模型由50 Hz的正弦信号及其谐波组成。幅值通常与ECG峰峰值相当或更强。
(2)电极接触噪声
电极接触噪声是瞬时干扰，来源于电极与肌肤的不良接触，即病人与检侧系统的连接不好。其连接不好可能是瞬时的，如病人的运动和振动导致松动；也可能是检测系统不断的开关、放大器输入端连接不好等。电极接触噪声可抽象为快速、随机变化的阶跃信号，它按指数形式衰减到基线值，包含工频成分。这种瞬态过渡过程可发生一次或多次、其特征值包括初始瞬态的幅值和工频成分的幅值、衰减的时间常数；其持续时间一般的1s左右，幅值可达记录仪的最大值。
(3)人为运动
人为运动是瞬时的(但非阶跃)基线改变，由电极移动中电极与皮肤阻抗改变所引起。人为运动由病人的运动和振动所引起，造成的基线干扰形 状可认为类似周期正弦信号，其峰值幅度和持续时间是变化的，幅值通常为几十毫伏。
(4)肌电干扰(EMG)
肌电干扰来自于人体的肌肉颤动，肌肉运动产生毫伏级电势。EMG基线通常在很小电压范围内。所以一般不明显。肌电干扰可视为瞬时发生的零均值带限噪声，主要能量集中在30-300 Hz范围内。
(5)基线漂移和呼吸时ECG幅值的变化
基线漂移和呼吸时ECG幅值的变化一般由人体呼吸、电极移动等低频干扰所引起，频率小于5 Hz；其变化可视为一个加在心电信号上的与呼吸频率同频率的正弦分量，在O.015-O.3Hz处基线变化变化幅度的为ECG峰峰值的15％。
(6)信号处理中用电设备产生的仪器噪声
心电信号是由人体心脏发出的极其精密、相当复杂并且有规律的微弱信号，外界干扰以及其它因素的存在都会使其变得更为复杂，要准确地对其进行自动检测、存储、分析却是一项十分艰巨的任务。例如，工频干扰信号对心电图的影响会使心电信号的特征点定位变得十分困难。因此，心电信号的监视、分析必须在建立在有效抑制各种干扰、检测出良好的心电信号的基础之上。
(7)共模信号(commonmode signal)
从体表采集到的信号除了人体心脏产生的电信号外，还包含许多与心电无关的电信号。由于体表各个导联均可看到这些信号，故称为共模信号。共模信号强度可以远远大于心电信号，从而干扰心电图分析。
为了抑制基线漂移，设置了0.5Hz高通滤波；由于心电信号属于低频信号，设置了二阶低通巴特沃斯滤波器，消除100 Hz以上的高频成分（带通滤波）；为了消除50 Hz工频干扰，设置50 Hz陷波器。
心电信号频域信号
心电信号时域分析
Load命令读取心电信号
Butterworth高通滤波器去除基线漂移
Butterworth低通滤波器去除工频干扰产生的毛刺
Butterworth带通滤波器去除肌电噪声
陷波器去除50Hz工频干扰
布拉克曼窗低通滤波器去除工频干扰产生的毛刺
2.3 系统总体设计框图

图2 系统总体设计框图
三、设计原理及方法
3.1 数字滤波器简介
数字滤波器是一种用来过滤时间离散信号的数字系统，通过对抽样数据进行数学处理来达到频域滤波的目的。可以设计系统的频率响应，让它满足一定的要求，从而对通过该系统的信号的某些特定的频率成分进行过滤，这就是滤波器的基本原理。如果系统是一个连续系统，则滤波器称为模拟滤波器。如果系统是一个离散系统，则滤波器称为数字滤波器。信号通过线性系统后，其输出就是输入信号和系统冲激响应的卷积。从频域分析来看，信号通过线性系统后，输出信号的频谱将是输入信号的频谱与系统传递函数的乘积。除非为常数，否则输出信号的频谱将不同于输入信号的频谱，某些频率成分较大的模，因此，中这些频率成分将得到加强，而另外一些频率成分的模很小甚至为零，中这部分频率分量将被削弱或消失。因此，系统的作用相当于对输入信号的频谱进行加权。
3.2 IIR滤波器的设计原理
IIR数字滤波器的设计一般是利用目前已经很成熟的模拟滤波器的设计方法来进行设计，通常采用模拟滤波器原型有butterworth函数、chebyshev函数、bessel函数、椭圆滤波器函数等。
IIR数字滤波器的设计步骤：
(1)按照一定规则把给定的滤波器技术指标转换为模拟低通滤波器的技术指标；
(2)根据模拟滤波器技术指标设计为响应的模拟低通滤波器；
(3)根据脉冲响应不变法和双线性不变法把模拟滤波器转换为数字滤波器；如果要设计的滤波器是高通、带通或带阻滤波器，则首先把它们的技术指标转化为模拟低通滤波器的技术指标，设计为数字低通滤波器，最后通过频率转换的方法来得到所要的滤波器。
3.3 IIR滤波器的设计
3.3.1 IIR数字低通滤波器设计过程
IIR滤波器系统函数的极点可以在单位圆内的任何位置，实现IIR滤波器的阶次较低，所用的存储单元少，效率高，又由于IIR数字滤波器能够保留一些模拟滤波器的优良的特性，因此应用很广。设计数字滤波器的方法主要有基于冲激响应不变法的IIR数字滤波器设计，基于双线性Z变换的IIR数字滤波器设计，数字高通，带通及带阻IIR滤波器设计。我们所使用的方法是基于双线性Z变换的IIR数字滤波器设计。
按照技术要求设计一个模拟滤波器，得到模拟低通滤波器的传输函数H（s），再按一定的转换关系将H(s)转换成数字低通滤波器的系数函数H(z)。这样设计的关键问题就是找到这样的转换关系，将s平面上的H(s)转换成z平面上的H(z)。
(1)巴特沃斯滤波器分母多项式的因式表示，如表1所示：

表1 巴特沃斯滤波器分母多项式的因式表示
(2)巴特沃斯低通滤波器的阶数公式
N=log10((10^(As/10)-1)/(10^(Rp/10)-1))/(2*log10(ws/wp))）
(3)巴特沃斯低通滤波器函数
由巴特沃斯低通滤波器的阶数公式和巴特沃斯滤波器分母多项式的因式表示求出归一化巴特沃斯低通滤波器Has (s )
N=7
则Has(s )=1/((s+1)*(s^2+0.4450s+1)*(s^2+1.247s+1)*(s^2+1.8022s+1))
(4)用于去除工频干扰产生的毛刺的巴特沃斯低通滤波器的频域特性，如图3所示:

图3 巴特沃斯低通滤波器的相频和幅频特性
(5)经过巴特沃斯低通滤波器器后心电信号的时域波形和频谱图，如图4所示:

图4 经过巴特沃斯低通滤波器后心电信号的时域波形图和频谱图
对比原始信号的时域波形图和频谱图可得通过低通滤波器后的心电信号波形图可以明显看出波形变得平滑，由工频干扰产生的毛刺被低通滤波器成功滤除。
(6)用于去除基线漂移的巴特沃斯高通滤波器的频域特性，如图5所示:

图5 巴特沃斯高通滤波器的频域特性
(7) 经过巴特沃斯高通滤波器器后心电信号的时域波形和频谱图，如图6所示:

图6 经过巴特沃斯高通滤波器器后心电信号的时域波形和频谱图
3.3.2 IIR数字带通滤波器设计过程
根据以上IIR数字滤波器设计方法:
(1)设计模拟低通原型滤波器。用模拟低通滤波器设计方法得到模拟低通滤波器的传输函数Ha(s);借助巴特沃斯(Butterworth)滤波器、切比雪夫(Chebyshev)滤波器、椭圆(Cauer)滤波器、贝塞尔(Bessel)滤波器等。
(2)调用lp2bp函数将模拟低通滤波器转化为模拟带通滤波器。
(3)利用双线性变换法将模拟带通滤波器Ha(s)转换成数字带通滤波器H(z)。
Butterworth模拟低通滤波器原型如图7所示：

图7 Butterworth模拟低通滤波器原型
用于去除肌电噪声的巴特沃斯带通滤波器的频域特性，如图8所示：

图8 巴特沃斯带通滤波器的频域特性

经过巴特沃斯带通滤波器器后心电信号的时域波形和频谱图，如图9所示:

图9 经过巴特沃斯带通滤波器器后心电信号的时域波形和频谱图
由于Butterworth滤波器通带内有最大的平滑特性，信号经过后衰减小，因此我们选用Butterworth带通滤波器滤除基线漂移和呼吸等引起的干扰。
3.4 FIR滤波器
数字滤波器(digital filter)是由数字乘法器、加法器和延时单元组成的一种装置。其功能是对输入离散信号的数字代码进行运算处理，以达到改变信号频谱的目的。由于电子计算机技术和大规模集成电路的发展，数字滤波器已可用计算机软件实现，也可用大规模集成数字硬件实时实现。数字滤波器根据其单位冲激响应函数的时域特性分为两种：无限长冲激响应(IIR) 滤波器和有限长冲激响应(FIR)滤波器。FIR数字滤波器又称有限长单位冲激响应滤波器，它的优点是可以做成具有严格的线性相位，同时又可以具有任意的幅度特性。此外，FIR滤波器的单位抽样响应是有限长的，因而滤波器一定是稳定的。再有，FIR滤波器由于单位冲激响应是有限长的，所以可以用快速傅里叶变换（FFT）算法来实现过滤信号，从而可大大提高运算效率。
在滤波器设计中要对理想滤波器抽样响应进行截断。截断后不可避免的产生了频谱泄漏, 为了尽量减小频谱泄漏, 在设计滤波器时要采用不同的窗函数来满足不同用途的要求各种窗函数的幅频响应都存在明显的主瓣和旁瓣。主瓣宽度和旁瓣的幅值衰减特性决定了窗函数的应用。用于滤波器的窗函数,一般要求窗函数主瓣宽度窄,以获得较好过渡带:旁瓣相对值尽可能小, 以增加通带段的平稳度和增大阻带的衰减。窗函数应满足在0 <Fn<N范围内关于a 对称,在其它区域取零值。由线性系统理论可知，在某种适度条件下，输入到线性系统的一个冲击完全可以表征系统。
当我们处理有限的离散数据时，线形系统的响应（包括对冲击的响应）也是有限的。若线性系统仅是一个空间滤波器，则通过简单地观察它对冲击的响应，我们就可以完全确定该滤波器。通过这种方式确定的滤波器称为有限冲击响应（FIR）滤波器。FIR滤波器是在数字信号处理（DSP）中经常使用的两种基本的滤波器之一，另一个为IIR滤波器。IIR滤波器是无限冲激响应滤波器。
FIR数字滤波器设计的主要方法有：窗函数法设计FIR数字滤波器，频率采样法设计FIR数字滤波器，最优化法设计FIR数字滤波器，此次设计我们采用窗函数法设计FIR数字滤波器。
3.4.1 FIR滤波器的设计
FIR滤波器的设计过程：
(1)给定理想的频率响应函数 Hd(e^jw)及技术指标δ, Δw;
(2)求出理想的单位抽样响应hd(n);
(3)根据阻带衰减选择窗函数w(n);
(4)根据过渡带宽度确定N 值N=A/Δw;
(5)求所设计的FIR滤波器的单位脉冲响应h(n)=hd(n)*w(n);
(6)计算频率响应Hd(e^jw),验算指标是否满足要求。
3.4.2 FIR数字低通滤波器设计过程
(1)布拉克曼窗低通滤波器的频域特性，如图10所示：
图10 布拉克曼窗低通滤波器的相频和幅频特性
(2)经过布拉克曼窗低通滤波器器后心电信号的时域波形和频谱图，如图11所示：

图11 经过布拉克曼窗低通滤波器器后心电信号的时域波形和频谱图
对比原始信号的时域波形图和频谱图可得通过低通滤波器后的心电信号波形图可以明显看出波形变得平滑，由工频干扰产生的毛刺被低通滤波器成功滤除。

3.5 陷波器
3.5.1陷波器的基本原理及作用
陷波器也称带阻滤波器（窄带阻滤波器），它能在保证其他频率的信号不损失的情况下，有效的抑制输入信号中某一频率信息。所以当电路中需要滤除存在的某一特定频率的干扰信号时，就经常用到陷波器。
在我国采用的是50hz频率的交流电，所以在平时需要对信号进行采集处理和分析时，常会存在50hz的工频干扰，对我们的信号处理造成很大干扰，因此50Hz陷波器在日常成产生活中被广泛应用，其技术已基本成熟。
工频陷波器不仅在通信领域里被大量应用，还在自动控制、雷达、声纳、人造卫星、仪器仪表测量及计算机技术等领域有着广泛的应用。
陷波器的设计方法有：文氏法、双T法、反相带通法、模拟电感法。
3.5.2双T法设计陷波器
陷波器的传输函数为
B(1/z) (z-exp(j*2*pi*f0))*(z-exp(-j*2*pi*f0))
H(z) = -------- = --------------------------------------------
A(1/z) (z-a*exp(j*2*pi*f0))*(z-a*exp(-j*2*pi*f0))
其中f0为陷波器要滤除信号的频率，a为与陷波器深度相关的参数，a越大，深度越深。
带阻滤波器的频率特性如图12所示:

图12 带阻滤波器的频率特性
去除50HZ工频干扰的陷波器的设计，如图13所示：

图13 陷波器
经过陷波器的心电信号的时域与频域波形，如图14所示：

图14 经过陷波器的心电信号的时域与频域波形
四、MATLAB简述
MATLAB 是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。
MATLAB是矩阵实验室（Matrix Laboratory）的简称，和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。
MATLAB是矩阵实验室（Matrix Laboratory）之意。除具备卓越的数值计算能力外，它还提供了专业水平的符号计算，文字处理，可视化建模仿真和实时控制等功能。MATLAB的基本数据单位是矩阵，它的指令表达式与数学,工程中常用的形式十分相似,故用MATLAB来解算问题要比用C,FORTRAN等语言完相同的事情简捷得多。
MATLAB包括拥有数百个内部函数的主包和三十几种工具包(Toolbox).工具包又可以分为功能性工具包和学科工具包.功能工具包用来扩充MATLAB的符号计算,可视化建模仿真,文字处理及实时控制等功能.学科工具包是专业性比较强的工具包,控制工具包,信号处理工具包,通信工具包等都属于此类。开放性使MATLAB广受用户欢迎.除内部函数外,所有MATLAB主包文件和各种工具包都是可读可修改的文件,用户通过对源程序的修改或加入自己编写程序构造新的专用工具包。
MATLAB是一个数据分析和处理功能十分强大的工程实用软件，运用它来进行信号的处理相当便捷，文章介绍了在MATLAB环境中对信号进行处理的方法，并对信号进行时域和频域的分析。另外，利用MATLAB环境设计数字滤波器滤除高频及低频成分。给出了设计IIR数字滤波器的方法,并通过用MATLAB语言来实现。根据要求，设计了50HZ陷波器，并用MATLAB语言来实现。

五、总结
这次课程设计虽然遇到了很多问题，很多困难，但是也学到了很多东西。不仅学到了书本上的东西，而且学到了很多课本上没有的东西，很多程序里的东西，特别是程序语法，总是有错误，但是总是不知道错在哪里，在细心的检查下，终于找出了错误和警告，排除困难后，程序编译就通过了，心里终于舒了一口气。还有各种各样问题，通过查网络和请教同学来弄明白，这个过程是痛苦的，有时候有些问题不能马上解决，感到很头痛，真想放弃这个问题，但是坚持下来，并且解决这些问题的时候，真的有种苦尽甘来的感觉。
应用MATLAB进行心电信号的处理是与我们所学课程及专业紧密相连的，有着很强的实践性。做这个课程设计的时候，并不是非常的顺利，我也有遇到很多困难。
刚开始由于对滤波器的滤波原理并不是很了解，于是我又翻出学过的数字信号处理课本，认真研究起各种滤波器了，这才使我明白了大多数滤波器是如何工作地，不再单单只是懂理论，理论与实际相结合是很重要的，只有理论知识是远远不够的，只有把所学的理论知识与实践相结合起来，从理论中得出结论。实验过程中，我感觉到初始心电信号和滤波输出后的心电信号有一定的差别，这说明了信号在处理过程中有损耗。不管对于什么样的课题，其实也是有很多东西可以发掘的，这需要我们在平时多积累，多思考，只有这样，才能取得更大的进步，才能学有所用，学有所长。
通过这次设计，进一步加深了对数字信号处理的了解，让我对它有了更加浓厚的兴趣。通过这次课程设计使我懂得了，平时的理论知识只有通过自己动手做一个课题，从做这个课题的过程中发现问题，解决问题，这个学习的过程，会比我们平时只通过课堂上听讲得到的知识更加生动立体，更让人记忆深刻。在设计的过程中，我发现同学间的互帮互助真的很重要。当我们有问题的时候，大家一起讨论，将自己的观点表达出来，当发现别人的观点与自己的不同的时候，我们通过查阅资料找到最终正确的答案，这个过程是互利互惠的。这也培养了我们以后走上工作岗位后的团队精神，对我们以后的为人处世都有很大帮助。
总的来说，通过这次的课程设计我对心电信号有了全面的认识，对数字信号处理的知识又有了深刻的理解，让我感受到只有在充分理解课本知识的前提下，才能更好的应用这个工具；并且熟练的应用MATLAB也可以很好的加深我对课程的理解，方便我的思维。这次设计使我了解了MATLAB的使用方法，学会分析滤波器的优劣和性能，提高了分析和动手实践能力。同时我相信，进一步加强对MATLAB的学习与研究对我今后的学习将会起到很大的帮助！

参考文献
[1] 陈天华. 数字图像处理[M]. 北京. 清华大学出版社，2009
[2] 刘卫国. MALTAB程序设计与应用[M]. 北京. 高等教育出版社，2008
[3] 程正兴. 小波分析算法与应用. 西安. 西安交通大学出版社，1998
[4] 程佩青. 数字信号处理. 北京. 清华大学出版社，2007
[5] 苏金明，王永利. MALTAB应用指南[M]. 上册. 北京电子工业出版社，2004
[6] 夏良正. 数字图像处理（修订版）[M]. 南京. 东南大学出版社，1999
[7] 霍红涛，林小竹，何薇. 数字图像处理[M]. 北京. 北京理工大学出版社，2003
[8] 张开滋，郭继鸿，刘海洋. 临床心电信息学[M]. 长沙. 湖南科技出版社，2002
[9] 黄宝晨，朱怡然. 心电图基本知识[J]. 中国乡村医药杂志，2004(第7页)
[10] 许原. 心电图解读心电图如何解读和诊断[J]. 中国临床医生，2004(第5页)
附录
窗函数设计法
窗函数设计法是一种通过截短和计权的方法使无限长非因果序列成为有限长脉冲响应序列的设计方法。通常在设计滤波器之前，应该先根据具体的工程应用确定滤波器的技术指标。在大多数实际应用中，数字滤波器常常被用来实现选频操作，所以指标的形式一般为在频域中以分贝值给出的相对幅度响应和相位响应。常用的窗函数有以下几种:矩形窗，三角窗，汉宁窗，海明窗，布拉克曼窗，切比雪夫窗，巴特里特窗以及凯塞窗。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
b=load('D:\Users\Data2.txt');
t=b(:,1);%时间
c=b(:,2);%幅值
figure(1);
subplot(311);
plot(t,c);
title('原始心电信号的时域波形图');xlabel('时间t/s');ylabel('幅值/A');
n=3600;
m=(fft(c,n));
fs=360;%采样频率
f=fs/n*(0:n-1);%定位坐标
subplot(312);
plot(f,abs(m));
title('原始心电信号的频谱图');
xlabel('频率/HZ');ylabel('幅值/db');axis([0,100,0,150]);
subplot(313);
plot(f,angle(m));
title('原始心电信号的相频图');
xlabel('频率/Hz');ylabel('相角/rad');axis([0,100,-5,5]);
% 低通IIR滤波器;
wp=2*pi*10/fs; %通带数字频率转换成模拟频率
ws=2*pi*15/fs; %通带截至数字频率转换成模拟频率
rp=2; %通带最大衰减
rs=20; %阻带最小衰减
[N,wc]=buttord(wp,ws,rp,rs,'s'); %确定最小阶数N和频率参数Wc
[Bz,Az]=butter(N,wc);%得巴特沃斯归一化低通原型
[H,w]=freqz(Bz,Az); %生成频率响应参数
f1=w/pi*fs/2; %采样频率转换成模拟采样频率
y1=filter(Bz,Az,c);%使用filter函数对信号进行滤波
figure(2);
subplot(211);
plot(f1,angle(H));
xlabel('频率/Hz');ylabel('幅度');
title('低通滤波器相频特性');
subplot(212);
plot(f1,abs(H));
xlabel('频率/Hz');ylabel('幅度');
title('低通滤波器幅频特性');
% 低通IIR滤波后图形
figure(3);
subplot(311);
plot(t,y1);
title('滤波后时域波形');xlabel('时间t/s');ylabel('幅值/A');
subplot(312);
f2=fs/n*(0:997-1);
plot(f2,abs(fft(y1)));
title('滤波后心电信号的频谱图');
xlabel('频率/Hz');ylabel('幅值/db');axis([0,100,0,150]);
subplot(313);
plot(f2,angle(fft(y1)));
title('滤波后心电信号的相频特性');
xlabel('频率/Hz');ylabel('相角/rad');axis([0,100,-5,5]);
%低通FIR滤波器
N=50; %定义窗函数的长度
wc=0.3;
window=blackman(N);%根据N的值产生一个布拉克曼窗window
hn=fir1(N-1,wc,window);%可以指定窗函数向量 window。如果缺省 window参数,则 fir1默认为 hamming窗。
y2=filter(hn,1,b(:,2));%使用filter函数对信号进行滤波
figure(4);
freqz(hn,1);
% 低通滤波器滤波后图形
figure(5);
subplot(311);
plot(t,y2);
title('滤波后时域波形');xlabel('时间t/s');ylabel('幅值/A');
subplot(312);
plot(f2,abs(fft(y2)));
title('滤波后心电信号的频谱图');
xlabel('频率/Hz');ylabel('幅值/db');axis([0,100,0,150]);
subplot(313);
plot(f2,angle(fft(y2)));
title('滤波后心电信号的相频特性');
xlabel('频率/Hz');ylabel('相角/rad');axis([0,100,-5,5]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
b=load('D:\Users\Data2.txt');
t=b(:,1);%时间
c=b(:,2);%幅值
%高通滤波器，去除基线漂移
fs=40;
wp=1*2/fs;
ws=0.25*2/fs;
rp=1;
rs=60;
Nn=256;
[N,wn]=buttord(wp,ws,rp,rs)
[Bz,Az]=butter(N,wn,'high')
[H,w]=freqz(Bz,Az); %生成频率响应参数
f1=w/pi*fs/2; %采样频率转换成模拟采样频率
y1=filter(Bz,Az,c);%使用filter函数对信号进行滤波
figure(2);
subplot(211);
plot(f1,angle(H));
xlabel('频率/Hz');ylabel('幅度');
title('高通滤波器相频特性');
subplot(212);
plot(f1,abs(H));
xlabel('频率/Hz');ylabel('幅度');
title('高通滤波器幅频特性')% 低通IIR滤波后图形
figure(3);
subplot(311);
plot(t,y1);
title('滤波后时域波形');xlabel('时间t/s');ylabel('幅值/A');
subplot(312);
n=40;
f2=fs/n*(0:997-1);
plot(f2,abs(fft(y1)));
title('滤波后心电信号的频谱图');
xlabel('频率/Hz');ylabel('幅值/db');axis([0,350,0,50]);
subplot(313);
plot(f2,angle(fft(y1)));
title('滤波后心电信号的相频特性');
xlabel('频率/Hz');ylabel('相角/rad');axis([0,350,-5,5]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
b=load('D:\Users\Data2.txt');
t=b(:,1);%时间
c=b(:,2);%幅值
%带通滤波器，滤除肌电噪声
fs=1000;
n=2;
wn=100*2/fs;
m=256;
[b,a]=butter(n,wn);
freqz(b,a,m,fs);
[H,w]=freqz(b,a); %生成频率响应参数
f1=w/pi*fs/2; %采样频率转换成模拟采样频率
y1=filter(b,a,c);%使用filter函数对信号进行滤波
figure(3);
subplot(311);
plot(t,y1);
title('滤波后时域波形');xlabel('时间t/s');ylabel('幅值/A');
subplot(312);
f2=0:996;
plot(f2,abs(fft(y1)));
title('滤波后心电信号的频谱图');
xlabel('频率/Hz');ylabel('幅值/db');axis([0,300,0,200]);
subplot(313);
plot(f2,angle(fft(y1)));
title('滤波后心电信号的相频特性');
xlabel('频率/Hz');ylabel('相角/rad');axis([0,300,-5,5]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%模拟高通滤波器的设计
m=0:0.5:10;
N=30;
[z,p,k]=buttap(N);
[b,a]=zp2tf(z,p,k);
[bt,at]=lp2hp(b,a,0.25*2*pi); %由低通原型滤波器转换为截止频率为100Hz的高通滤波器
[Ht,w]=freqs(bt,at,m);
figure;
subplot(311);
hold on;
plot(w,abs(Ht));
xlabel('w/pi');
ylabel('|H(jw)|^2');
title('模拟高通滤波器');
grid on;
%模拟高通滤波器性能
pha=angle(Ht); %输出系统的相频特性
subplot(312);
plot(w,20*log10(pha));
grid;
xlabel('w/wc');ylabel('相位/dB');
title('模拟高通滤波器相频特性');
mag=abs(Ht); %输出系统的幅频特性
subplot(313);
plot(w,20*log10(mag));
grid;
xlabel('w/wc');ylabel('幅度/dB');title(' 模拟高通滤波器幅频特性')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%