课题二基于MATLAB平台的心电信号分析系统设计与仿真

课题二基于MATLABDE的心电信号分析系统的设计与仿真

一、本课题的目的

本设计课题主要研究数字心电信号的初步分析及滤波器的应用。通过完成本课题的设计，拟主要达到以下几个目的：

(1)了解MATLAB软件的特点和使用方法，熟悉基于Simulink的动态建模和仿真的步骤和过程；

(2)了解人体心电信号的时域特征和频谱特征；

(3)进一步了解数字信号的分析方法；

(4)通过应用具体的滤波器进一步加深对滤波器理解；

(5)通过本课题的设计，培养学生运用所学知识分析和解决实际问题的能力。

二、课题任务

设计一个简单的心电信号分析系统。对输入的原始心电信号，进行一定的数字信号处理，进行频谱分析。采用Matlab语言设计，要求分别采用两种方式进行仿真，即直接采用Matlab 语言编程的静态仿真方式、采用Simulink进行动态建模和仿真的方式。根据具体设计要求完成系统的程序编写、调试及功能测试。

(1)对原始数字心电信号进行读取，由数字信号数据绘制出其时域波形。

(2)对数字信号数据做一次线性插值，使其成为均匀数字信号，以便后面的信号分析。

(3)根据心电信号的频域特征（自己查阅相关资料），设计相应的低通和高通滤波器。

(4)编程绘制实现信号处理前后的频谱，做频谱分析，得出相关结论。

(5)对系统进行综合测试，整理数据，撰写设计报告。

三、主要设备和软件

(1)PC机一台。

(2) MA TLAB6.5以上版本软件，一套。

四、设计内容、步骤和要求

4.1必做部分

4.1.1利用Matlab对MIT-BIH数据库提供的数字心电信号进行读取，并还原实际波形

美国麻省理工学院提供的MIT-BIH数据库是一个权威性的国际心电图检测标准库，近年来应用广泛，为我国的医学工程界所重视。MIT-BIH数据库共有48个病例，每个病例数据长30min，总计约有116000多个心拍，包含有正常心拍和各种异常心拍，内容丰富完整。

为了读取简单方便，采用其txt 格式的数据文件作为我们的原心电信号数据。利用Matlab 提供的文件textread 或textscan 函数，读取txt 数据文件中的信号，并且还原实际波形。

4.1.2对原始心电信号做线性插值

由于原始心电信号数据不是通过等间隔采样得到的，也就是说原始的心电数据并不是均匀的，而用Matlab 中提供的数字滤波器处理数据时，要求数据是等间隔的。因此设计的系统首先应对原始心电信号做线性插值处理，使其变为等间隔的数字信号，否则直接处理后会出现偏差，根据心电信号的特点, 把时间分隔成0.001s 。添加的幅值点采用一次线性插值。对二维数据进行插值,相连幅值间数据的插值根据时间进行,运算公式如下:

1--=?i i t t t ，001.0/t N ?=，1--=?i i A A A ，001.01+=-j j t t ，N A A A j j

/1?+=- 其中i t 是第i 个数据时间点，A i 是与之对应的数据，N 是两数据之间需要的插值数，A ?是

需要插值的两点数据差，，

，，，，，，，1321321-==N j arraysize i 1111----==i j i j A A t t ，， 1=j i ，时数组j j A t ,依次排列,即得到了插值后等间隔的新数据。

将插值后的幅值数据另存为txt 文件，以供后续程序调用。

4.1.3根据心电信号的频域特征，设计相应的低通和带通滤波器

一般正常人的心电信号频率在0.7～100HZ 范围内，幅度为V 10μ(胎儿)～5mV (成人)。人体心电信号微弱，信噪比小，因此，在采集心电信号时，易受到仪器、人体活动等因素的影响，而且所采集的心电信号常伴有干扰。采集心电数据时，由于人的说话呼吸，常常会混有约为0.1Hz 到0.25Hz 频段的干扰，对于这些低频干扰，可以让信号通过一个高频滤波器，低截止频率设置为0.25,来滤除低频信号，对于高频信号干扰，可以让信号再通过一个低频滤波器，其中截止频率设置为99Hz 。也可以直接应用带通滤波器设计。

（1）根据以上指标，设计模拟巴特沃斯（切比雪夫）低通、高通或带通滤波器，画出幅频特性（模拟滤波器幅频特性freqs ）。

（2）根据心电信号频谱范围设计一个3阶以上模拟滤波器对心电信号进行预滤波；

（3）采用直接、级联或并联方式，实现该系统，并画出系统的信号流图；

（4）分析系统的时域特性（阶跃响应、冲击响应等），并用Matlab 绘出相关波形；

（5）用Matlab 分析幅频特性，并绘出相关波形；

（6）分析系统函数零极点与幅频特性的关系。

4.1.4对处理前后的心电信号分别做频谱分析

利用Matlab 软件对处理前后的心电信号编程显示其频谱，分析比较滤波前后的频谱，得出结论。

如果分析频谱，滤波效果不明显，则需变动滤波器参数指标，重新设计滤波器。通过频谱分析，多次试验确定最合适的滤波器。

4.1.5 Simulink仿真

根据前面的设计，进行基于Simulink的动态仿真设计。实现心电信号的分析和处理。给出系统的基于Simulink的动态建模和仿真的系统方框图，同时记录系统的各个输出点的波形和频谱图。

4.2 选作部分

4.2.1减少分析数据的工作量试验

（1）只截取大约2.5s，三个周期左右，大约800个采样数据进行分析；

（2）利用Matlab编程实现对10s原始心电信号数据二次采样并进行结果分析。

4.2.2 simulink仿真设计

（1）将Matlab语言编程进行线性插值后的心电信号结果，直接输入到simulink模块中进行滤波处理及系统仿真；

（2）在仿真环境下，自己设计simulink的线性插值模块，直接输入原始心电信号，进行插值处理再进行系统的功能仿真。

4.2.3 50Hz工频陷波器设计

由于电子设备采集到的信号经常会混有电源线干扰。电源线干扰是以50 Hz为中心的窄带噪声，带宽小于1Hz。设计相应的带阻滤波器滤除电源线干扰，并对处理后的信号做频谱分析。

五、课程设计报告要求

（1）设计报告书包括内容：课程设计题目，课程设计目的和意义，设计方案，详细设计步骤，设计结果（原理图等），测试和仿真结果（图形或数据）及其分析，其它有明确要求的设计内容，结论，参考文献等。

（2）提交课程设计报告时应同时提交相关设计和仿真分析材料（程序、结果等）的电子版。

六、参考文献

[1] 北京迪阳正泰科技发展公司.综合通信实验系统——信号与系统指导书（第二版）. 2006,6

[2] 丁玉美.数字信号处理（第二版）.西安电子科技大学出版社，2001

[3] 吴大正. 信号与线性系统分析（第四版）. 高等教育出版社，2005,8

[4] 谢嘉奎. 电子线路--线性部分(第四版). 高等教育出版社，2003,2

[5] 陈后金. 信号分析与处理实验. 高等教育出版社，2006，8

七、附录——设计原理

1.心电信号的读取

txt格式的数据文件内容及格式如图1-1所示（以100.txt为例）。

图1-1txt格式的心电数据文件

其中文件的第一列为采样时间，第二列是在以MLII这种导联方式所得到的采样数据，第三列是以V5这种导联方式所得到的采样数据，全文件记录了约为10s的心电数据，3600个采样数据，每一行数据之间用Tab符分隔。

由于数据文件中后两列数据是对同一种心电信号进行不同的导联方式所得到的采样数据，所以可以只采用其中的一种采样数据，摒弃另外一种，即可完成对此心电信号的分析。全部的心电文件记录时间约为10s，共计12个左右周期的心电信号。

实际设计心电信号数据文件时应注意：

（1）数据文件的前两行为解释说明文字，不是真正的信号数据，读取信号程序要能够自动忽略前两行文字，只读取真正的数字信号数据（严禁自己手动删除原心电数据文件中的前两行数据，必须通过程序来实现忽略前两行文字的目的）。

（2）利用数组函数分别将文件的前两个列分别读入一个一维数组。

（3）最后利用已经转为数值的分别代表心电信号时间和幅值的两个一维数组，图形化还原原始心电信号波形，在此推荐利用plot(x,y)函数对数据做图形化显示。

2.心电信号的线性插值处理

根据上文中提到的插值公式，以此为原理，设计Matlab程序，对心电信号数据做线性

插值处理。插值完以后的数据应该是时间均匀的、以0.001秒为间隔的。

此步骤的实现主要是基于Matlab 中的数组操作函数来实现，建议一定首先熟悉并掌握Matlab 中的所有数组操作函数的作用和操作方法。

其中一种插值方法的思路是：将第一步中读取的心电信号数据的时间数据和幅值数据分别存放在一个一维数组中。然后利用for 循环结构把所有数据依次读取进来。判断时间数据数组中前后两个相邻的数据间隔是否为0.001s ，如果是则判断下一对相邻两个数据；如果间隔大于0.001s 则进行一维插值处理。

注意对时间数据做插值的同时一定不要忘记对幅值数据同样做插值处理，时间数据和幅值数据一定是相互对应的。

3. 滤波器设计

3.1 模拟滤波器设计原理

（1）模拟巴特沃思滤波器原理

巴特沃斯滤波器具有单调下降的幅频特性：在小于截止频率c Ω的范围内，具有最平幅

度的响应，而在c Ω>Ω后，幅频响应迅速下降。

巴特沃思低通滤波器幅度平方函数为： 221()1()a

N c H j Ω=Ω

+Ω (3-1)

式中N 为滤波器阶数，c Ω为3dB 截止角频率。将幅度平方函数写成s 的函数形式：

21()()

1()

a a N c H s H s s

j -=+Ω (3-2) 该幅度平方函数有2N 个等间隔分布在半径为c Ω的圆上的极点121()22 k j N k c s e π++=Ω，

0,1,...21k N =- 为了形成稳定的滤波器，取左半平面的N 个极点构成()a H s ，即：

1

0()()N N a c k k H s s s

-==Ω-∏ (3-3)

为使设计统一，将频率归一化，得到归一化极点121()22k j N k p e

π++=，相应的归一化系统函

数为： 10()1

()N a k k H p p p -==-∏ (3-4)

多项式形式为： 01()1(...)N a H p b b p p =+++ (3-5)

（2）模拟切比雪夫滤波器原理

切比雪夫滤波器的幅频特性具有等波纹特性，有两种形式，在通带内等波纹、阻带单调的是I 型滤波器，在通带内单调、在阻带内等波纹的是Ⅱ滤波器。以I 型滤波器为例。

切比雪夫滤波器的幅度平方函数为：

22221()()1()

a

N p A H j C εΩ=Ω=Ω+Ω (3-6) ε为小于1的正数，表示通带内幅度波动的程度。Ωp 称为通带截止频率。令λ=Ω/Ωp ，称为对Ωp 的归一化频率。C N (x)为N 阶切比雪夫多项式。幅度平方函数的极点是分布在b Ωp 为长半轴，a Ωp 为短半轴的椭圆上的点。同样取s 平面左半平面的极点构成()a H s ：

11

()2

()N

N N a p i i H s s s ε-==Ω??-∏ (3-7) 进行归一化，得到：11()12

()N N a i i H p p p ε-==??-∏ 其中(21)(21)sin[]cos[]22k k p ch jch i N N π

π

ξξ--=-?+? ，1

1()Arsh N ξε

= 3.2 模拟滤波器数字化原理

将模拟滤波器转化为数字滤波器在工程上常用的有脉冲响应不变法和双线性变换法。 脉冲响应不变法是一种时域上的转换方法，它是数字滤波器的单位取样响应在采样点上等于模拟滤波器的单位冲激响应，即：

()()h n h nT a = (3-8)

设模拟滤波器只有单阶极点，其系统函数为：

1()N

i a i i A H s s s ==-∑

(3-9)

对()a H s 进行拉氏反变换得到()a h t ，对()a h t 进行等间隔采样，得到()()h n h nT a =，对()h n 进行Z 变换，得到数字滤波器系统函数：

11()1i N

i s T i A H z e z -==-∑ (3-10)

这种方法s 和z 的关系是：sT z e =。该方法的优点是频率坐标变换时线性的切数字滤波器的单位脉冲响应完全模仿模拟滤波器的单位冲激响应，时域特性逼近好；缺点是会产生频谱混叠现象，适合低通、带通滤波器的设计，不适合高通、带阻滤波器的设计。

双线性变换法为了克服频谱混叠现象，采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到/T π±之间，再用sT z e =转换到Z 平面上。

这种方法s 和z 的关系是：11(2/)(1/1)s T z z --=-+。该方法克服了频谱混叠现象，但带来了频率坐标变换的非线性：(2/)tan(/2)T ωΩ=，由模拟滤波器系统函数转换为数字滤波器系统函数公式为： 11

211()()|a z z T z H z H s ---=+= (3-11)

3.3数字高通、带通、带阻滤波器的设计

上述滤波器可以借助于模拟滤波器的频率变换设计一个所需类型的模拟滤波器， 再通过双线性变换法将其转换成所需类型的数字滤波器。

首先确定所需类型数字滤波器的技术指标；然后将数字滤波器技术指标按照公式(2/)t a n (/2)

T ωΩ=转换成所需类型滤波器的模拟域技术指标；将所需类型滤波器的模拟域技术指标转换成低通滤波器技术指标；设计归一化模拟低通滤波器；去归一化得到模拟低通滤波器的系统函数；将模拟低通滤波器转换为所需类型的模拟滤波器；最后通过双线性变换法转换成所需类型的数字滤波器。

4.使用Simulink 建模和仿真的过程

启动Matlab 后，在命令窗口中输入命令“simulink ”，打开simulink 模块库窗口（使用命令‘simulink3’可以打开老版本的simulink 模块库界面）。典型的Simulink 模块包括三个部分：输入模块、状态模块、输出模块。

4.1 Simulink 模块库简介

(1)Continuous （连续模块）库

(2)Discrete （离散模块）库

(3)Math （数学模块）库

(4)Sinks （信号输出模块）库：常用模块为Scope （示波器模块）、XYGraph （二维信号显示模块）、Display （显示模块）

(5)Sources （信号源模块）库（如图4-1所示），常见模块有：Constant （输入常数模块）、Signal Generator （信号源发生器模块）。

(6)Signal Generator 用于产生不同的信号波形，其中包括：正弦波、方波、锯齿波信号。

Sources （信号源模块）还包括其它常用模块：Ramp （斜坡输入信号）、Sine Wave （正弦波输入信号）、Step （阶跃输入信号）、Clock （时间信号）、Pulse （脉冲信号）等。

图4-1 Sources （信号源模块）库

4.2 利用Simulink 建立滤波器仿真实例

设某一六阶Butterworth 数字低通滤波器的系统函数为：

)2155.09044.01)(3583.00106.11)(7051.02686.11()1(0007378.0)(21212161-------+-+-+-+=z z z z z z z

z H

（1）采用级联型结构实现该数字低通滤波器；

（2）观察对实际心电图信号的滤波效果。

人体心电图信号在测量过程中往往受到工业高频干扰，所以必须经过低通滤波处理。已知某一实际心电图信号的采样序列如下：

x(n)=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0]，共56点。

解:滤波器的系统函数H(z)可以写成以下形式: 3

12(12)()()121)1A z z H z H z k B z C z k k k --++=

=----=∏

其中，0.09036, 1.26860.7051

11A B C ===-，1.0106,0.35830.9.440.215522333B C B C ==-===，，

H(z)由三个二阶子网络H 1(z)、H 2(z)和H 3(z)级联组成，如图4-2所示。

图4-2 滤波器H(z)的级联组成 仿真过程如下：

（1） 建立模型

首先打开simulink 模块库窗口，在simulink 模块库窗口中单击菜单项“File/New/Model ”，即可以建立一个新的simulink 模型文件。如图4-3所示。

利用鼠标单击Simulink 模块库窗口中的Continuous 子库，选取传递函数模块Transfer Fcn ，将它拖动到新建模型文件窗口的合适的位置。然后对模型模块进行参数设置和修改，单击右键从快捷菜单中选取“TransferFcn parameters …”修改传递函数参数，在弹出对话框中的传递函数分子系数“Numerator ：”栏填入[1 2 1]，在传递函数分母系数“Denominator ：”栏填入[1 1.2686 -0.7051]，其余参数使用默认值。若需要进一步了解该模块的参数设置说明，可以单击该对话框下方的“help ”按钮。然后确认，就得到了需要仿真的传递函数。如图4-4所示。

采用同样的方法，在Simulink 模块库中的Sources 子库中选取激励信号源，本例选取simin(From Workspace)，并将它拖动到建模文件窗口的合适的位置。然后在Sinks 子库中选取示波器（Scope ）模块作为系统输出波形显示。接下来利用鼠标将三个模块连接起来。

连接方法如下：

模块外部的大于符号“>”分别表示信号的输入输出节点，为了连接两个模块的输入输

出节点，可以将鼠标置于节点处，这时鼠标显示为“十”字形状，拖动鼠标到另一个模块的端口，然后释放鼠标按钮，则可以形成带箭头的连线，箭头方向表示信号的流向。

图4-3 Simulink模块库窗口和新建模型文件窗口

图4-4 修改仿真模型的参数

建立好的系统模型可以存盘为模型文件，扩展名为“mdl”。如图4-5所示。

图4-5完成的建模方框图

（2）参数设置

建立好系统模型后，接着要对输入信号源进行参数设置。输入信号源的参数设置界面如图4-6所示。最后双击示波器模型图标，打开示波器显示窗口。在快捷键设置菜单中设置为自动刻度。

（3）启动仿真

以上工作完成后，可通过建模窗口菜单项“Simulation/Start”启动仿真，也可以单击工具栏上的小三角按钮启动仿真。

图4-6输入信号源的参数设置界面

4.3Signal Processing Blockset模块仿真过程

Signal Processing Blockset提供了包括变换、矩阵运算、FIR、IIR、自适应和多速率滤波、谱分析和实时数据I/O等模块。利用这个工具箱对数字滤波器进行仿真。图4-7为Signal

Processing Blockset模块库，图4-8为DSP Sources模块库，图4-9为Filter Designs模块库。

图4-7 Signal Processing Blockset模块库

图4-8 DSP Sources模块库

图4-9 Filter Designs模块库

通过调用Simulink中的功能模块可以构成数字滤波器的仿真，在仿真过程中，可以双击各功能模块，随时改变参数，获得不同状态下的仿真结果。

（1）Matlab自带的滤波器设计与分析工具(FDAT001)介绍

在Matlab command窗口下，执行FDATool命令，将弹出一个FDATool窗口，Matlab自带的滤波器设计与分析工具(FDAT001)设计滤波器，FDATool是交互式信号分析处理环境，通过这一窗口用户能够方便地进行滤波器设计和分析。FDATool(Filter Design&Analysis T001)是Matlab信号处理工具箱里专用的滤波器设计分析工具，Matlab6.0以上的版本还专门增加了滤波器设计工具箱(Filter Design Toolbox)。FDATool可以设计几乎所有的基本的常规滤波器，包括FIR和IIR的各种设计方法。它操作简单，方便灵活。FDATool界面总共分两大部分，一部分是Design Filter，在界面的下半部，用来设置滤波器的设计参数，另一部分则是特性区，在界面的上半部分，用来显示滤波器的各种特性。Design Filter部分主要分为：Filter Type(滤波器类型)选项，包括Lowpass(低通)、Highpass(高通)、Bandpass(带通)、Bandstop(带阻)和特殊的FIR滤波器。Design Method(设计方法)选项，包括IIR滤波器的Butterworth(巴特沃思)法、Chebyshev Type I(切比雪夫I型)法、Chebyshev Type II(切比雪夫II型)法、Elliptic(椭圆滤波器)法和FIR滤波器的Equiripple法、Least—Squares(最小乘方)法、Window(窗函数)法。Filter Order(滤波器阶数)选项，定义滤波器的阶数，包括Spec ify Order(指定阶数)和

Minimum Order(最小阶数)。在Specify Order中填入所要设计的滤波器的阶数(N阶滤波器，

Specify Order=N.1)，如果选择Minimum Order则Matlab根据所选择的滤波器类型自动使用最小阶数。

心电滤波需要设计低通滤波器、高通滤波器和50Hz陷波器。

（2）低通滤波器设计

首先设计低通滤波器，在Filter Type中选择Lowpass(低通滤波器)；在Design Method选项中选择FIR(FIR滤波器)，接着在Window Specifications选项中选取Equiripple；指定Filter Order项中的Minimum；由于采用Equiripple法设计，只要给出Fstop(阻带截止频率)、Fpass(通带截止频率)，选取Fpass=100Hz、Fstop=102Hz。

设置完以后点击Design Filter即可得到所设计的FIR滤波器。通过菜单选项Analysis可以在特性区看到所设计滤波器的幅频响应、相频响应、零极点配置和滤波器系数等各种特性。设计完成后将结果保存为1．fda文件。在设计过程中，可以对比滤波器幅频相频特性和设计要求，随时调整参数和滤波器类型，低通滤波器的幅频如图4-10所示。

图4-10低通滤波器的参数设定及幅频特性

（3）高通滤波器设计

设计高通滤波器，在Filter Type中选择highpass(高通滤波器)；在Design Method选项中选择IIR(IIR滤波器)，接着在Window Specifications选项中选取Butterworth；指定Filter Order 项中的Minimu；由于采用Butterworth法设计，只要给出Fstop(阻带截止频率)、Fpass(通带截止频率)。选取Fstop=0.25Hz、Fpass=1Hz。设置完以后点击Design Filter即可得到所设计的IIR滤波器。通过菜单选项Analysis可以在特性区看到所设计滤波器的幅频响应、相频响应、零极点配置和滤波器系数等各种特胜。设计完成后将结果保存为2．fda文件。在设计过程中，可以对比滤波器幅频相频特性和设计要求，随时调整参数和滤波器类型，高通滤波器的幅频和相频特性如图4-11示。

图4-11 高通滤波器的参数设定及幅频特性

（4）50Hz陷波器设计

设计50Hz陷波器，在Filter Type中选择Bandstop(带阻滤波器)；在Design Method选项中选择IIR Window(HR滤波器巴特沃思法)，接着在Window Specifications选项中选取Butterworth；指定Filter Order项中的Minimu；由于采用Butterworth法设计，只要给出Fstopl(下阻带截止频率)、Fpassl(通带下限截止频率)、Fpass2(通带上限截止频率)、Fstop2(上阻带截止频率)。选取Fstopl=49Hz、Fpassl=48Hz、Fstop2=51Hz、Fpass2=53Hz。设置完以后点击Design Filter即可得到所设计的IIR滤波器。通过菜单选项Analysis可以在特性区看到所设计滤波器的幅频响应、相频响应、零极点配置和滤波器系数等各种特胜。设计完成后将结果保存为3．fda文件。在设计过程中，可以对比滤波器幅频相频特性和设计要求，随时调整参数和滤波器类型，带阻滤波器的幅频和相频特性如图4-12示。

图4-12带阻滤波器的参数设定及幅频特性

通过Simulink环境下的Digital Filter Design（数字滤波器设计）模块导入FDA Tool中所设计的滤波器文件1.fda、2.fda、3.fda。仿真图如图4-13所示。然后加载插值后的心电信号*.txt,对信号进行滤波，仿真结果如图4-14所示。

图4-13心电信号系统仿真图

图4-14仿真结果