大学毕业设计基于matlab的医学信号频谱分析系统的分析与设计

大学毕业设计基于matlab的医学信号频谱分析系统的分析与设计内容摘要：

矢量 h(n)仿真，结果如图 32。 图 33（ a） 为 h(n)的仿真幅频特性， 参数为：中心频率 f0=50 Hz，带内衰减 50 dB， Fpass2 Fpass1 2Hz。 (b)为 h(n)仿真的相频特性曲线，满足线性关系，保证了滤波中几乎不会产生附加频率。 9 （ a）幅频特 性 (b) 相频特性 图 33 陷波器冲击响应矢量仿真 应用 把滤波器应用于心电信号的预处理中。 FDATool 生成 h(n)后，用 Target菜单 Export to Code Composer Studio to(IDE)选项输出陷波器参数到文件。 在滤波器设计工程中，打开该文件，把陷波器的参数引入。 打开212 格式的心电数据文件 ，进行滤波 [12]。 结果如图 34。 其中 (a)为滤波前的 ECG 信号波形， (b)是滤波后的 ECG 信号波形，证明该滤波器是有效的。 10 图 34 心电信号 50hz 滤波前后对比 结 果分析 由图 34 知，运用 FDATool 工具箱来设计陷波器，不但操作简单，便于修改参数， 重新调整滤波设置， 而且滤波 的效果明显，能够用于所有的医电信号的 50hz 工频干扰的滤除。 167。 小波消噪 一 、 小波函数的选取 在确定了降噪的大致步骤以后，我们首先要确定小波函数。 在小波分解中， 小波函数将成为在各个尺度上构成其分解对象的基本单元。 因而，小波函数的性质直接决定了其逼近基线漂移的能力 [13]。 另外分解尺度 α在小波分解过程中也非常重要。 例如：基线漂移由信号的低频部 分表征，本身不具有很大的突变。 随小波分解尺度的增加，时间分辨率降低，逼近信号所含高频分量减少。 但是基线漂移的低频成分也不是单一频率的，当分解水平超过某一程度时，会丢失大量基线漂移的信号，因此，分层水 11 平 α的选择很重要。 分解水平过小，容易损失有用的心电信息；分解水平过大，则不能完全去除基线漂移 [14]。 通过对比几种典型小波统一分解水平下对 ECG信号分解后得到的逼近信号及重构结果，采用 db5小波优势较大，是能较好地满足要求的小波函数之一。 二 、 用小波变换法降噪 无论用哪种小波函数进行降噪，其步骤大致分为以下三步 [15]： 步骤 1：一维信号的小波分解。 选择一个小波函数，并且确定分解层次，进行分解。 其中，如何选取小波函数是关键。 步骤 2：小波分解高频系数的阈值量化。 对各个分解尺度下的高频系数进行一维小波重构。 如何选取阈值以及进行阈值量化是关键。 对于基线漂移，由于它处于小波分解最高层次的低频层，因此，可以直接将该层细数置零，去除基漂。 步骤 3：一维小波重构。 根据小波分解的最底层低频系数和各层高频系数进行一维小波重构。 利用小波变换方法对信号进行降噪处理，大致可分为以下三种方法 [16]： (1)默认阈值消噪。 (2)给定阈 值消噪。 在实际应用过程中，阈值往往可以通过经验公式获得， 而且这种阈值比默认阈值的可信度高。 (3)强制消噪。 将小波分解结构中的高频系数或想要去掉的频率水平的系数全部置零。 即滤掉所有不需要的部分，然后对信号进行重构。 方法简单，消噪后信号平滑，但容易失去有用成分。 需要注意的是运用阈值法消噪效果强烈依赖于所选择的小波，如果用与信号形状相近的小波去噪，则会得到较好的效果。 另外，阈值的选取也很大程度上影响了降噪的效果。 而且，阈值的选择方法很多，有些是基于原始信号的，而有些则基于样本估计。 有些是经验值，而有些则完全 是计算所得。 (1)从原始信号确定各级阈值：小波变换中，对各层系数降噪所需的阈值一般是根据原信号的信噪比 (噪声强度 )来获取的。 (2)基于样本估计的阈值的选取：在数学上有一 个常用的标准就是在最坏情况下方差最小的约束下的样本估计。 对信号做 无偏似然估计，然后根据最坏情况下降噪信号与原始信号方差最小的原则确定一个统一的阈值，然后截去超出这个阈值的系数 [17]。 12 三 、算法实现 及仿真结果 本设计选取默认阈值消噪，利用 sym8 小波 4 层分解来实现。 其来滤波结果如图 35 所示。 图 35 心电信号小波滤波前后 四 、结果分析 小波阈值消噪方法原理简单，能得到原始信号的最优估计且广泛的信号自适应性，在保持信号的奇异性的同时能有效的去除噪声。 本 设计 以心电信号为例验证了小波阈值消噪技术在信号预处理中的应用，取得了较为满意的结果。 对于其他的生物医学信号诸如心音、胃电等等也可用相同方式预处理。 167。 巴特沃斯滤波器 滤波 相关原理： 设计数字滤波器的任务就是寻求一个因果稳定的线性时不变系统，并使系统函数 H(z)具有指定的频率特性。 对于图 36 所示的 LSI 系统，其时域输入输出关系是： 13      nhnxny  （ 23） h ( n )x ( n ) y ( n ) 图 36 LSI 系统 对于各种类型滤波器，只要设计出滤波器的冲激响应函数，就可以有上式得到系统的输出了。 数字滤波器从实现的网络结构或者从单位冲激响应分类，可以分成无限长单位冲激响应 （ IIR） 数字滤波器和有限长单位冲激响应 （ FIR） 数字滤波器 [18]。 数字滤波器频率响应的三个参数 ： 幅度平方响应 ， 相位响应 及 群时延响应。 （ 1） IIR 数字滤波器： IIR 数字滤波器的系统函数为 的有理分数，即      jnyainxbny MjjMii   00 （ 24） 其中 bi 和 aj 是滤波器的系数，且 aj 中至少有一个非零。 与上述差分方程对应的传递函数为：      NNMMzaazbbbzXzYzH...z1...z11110 （ 25） IIR 数字滤波器的逼近问题就是求解滤波器的系数和 ，使得在规定的物理意义上逼近所要求的特性的问题。 如果是在 s 平面上逼近，就得到模拟滤波器，如果是在 z 平面上逼近，则得到数字滤波器。 （ 2） FIR 数字滤波器： 假设所希望的数字滤波器的频率响应为  jwd eH ，它是频域的周期函数，周期为 2π ， 那么 与  jwd eH 相对应的傅里叶系数为     dweeHnh jnwjwdd  2 1 （ 26） 以 nhd 为冲激响应的数字滤波器将具有频域响应  jwd eH。 14 设 FIR 的单位脉冲响应 nhd 为实 数，长度为 N，则其 z 变换和频率响应分别为按频域采样定理 FIR 数字滤波器的传输函数 H(z)和单位脉冲响应nhd 可由它的 N 个频域采样值 H(k)唯一确定 [19]。 设计内容 用 butter 函数求系统函数分子与分母系数的几种形式 ： （ 1）、 [b,a]=butter(N,wc,39。 high39。 ): 设计 N 阶高通滤波器， wc 为它的 3dB 边缘频率。 （ 2）、 [b,a]=butter(N,wc)： 当 wc 为具有两个元素的矢量 wc=[w1,w2]时，它设计 2N 阶带通滤波器。 （ 3）、 [b,a]=butter(N,wc,39。 stop39。 ): 若 wc=[w1,w2]，则它设计 2N 阶带阻滤波器。 为了设计任意的选项巴特沃斯滤波器，必须知道阶数 N 和 3dB 边缘矢量wc。 这可以直接利用信号处理工具箱中的 butter 函数来计算。 对于不同类型的滤波器，参数 wp 和 ws 有一些限制：对于低通滤波器，wpws；对于高通滤波器， wpws。 对于带通滤波器， wp 和 ws 分别为具有两个元素的矢量， wp=[wp1,wp2]和 ws=[ws1,ws2]，并且 ws1wp1wp2ws2。 对于带阻滤波器 wp1ws1ws2wp2。 本设计以心电信号为例，用 butter 函数来设计带通滤波器，设 N=5，Wn=[ 100],其幅频和相频响应如图 36 所示： 图 36 butter 带通滤波器的幅频和相频响应 仿真结果 15 由 butter 函数可以得到滤波器的传递函数的分子分母系数 [b,a],再利用滤波函数 filter（ b,a,x） 就可以实现滤波。 其仿真结果如图 37 所示。 图 37 心电信号巴特沃斯滤波前后 结果分析 由图 36 知，用 butter 函数设计带通滤波器的效果良好，在 100Hz 处有明显的转折， 对大于 100Hz 的频率分量有明显的衰减作用，从而达到预期的滤波效果。 由图 37 知，滤波后的波形 没有细微的波动，说明将心电频率低于 的强干扰有效地滤除。 另外用该种方法设计滤波器不但操作简便，设置简单，且具有普遍的通用性。 它既可以设计低通，高通滤波器，还能设计带通，带阻滤波器。 16 第三章 信 号频谱分析 任意一个信号都具有时域与频域特性，信号的频谱完全代表了信号，因而研究信号的频谱就等于研究信号本身。 通常从频域角度对信号进行分析与处理，容易对信号的特性获得深入的了解。 因此，信号的频谱分析是数 字信号处理技术中的一种较为重要的工具 [20]。 167。 傅里叶变换 傅里叶变换是信号频谱分析中常用的工具，它是建立在以时间为自变量的“信号”与以频率为自变量的“频谱函数”之间的某种关系。 “时间”或“频率”取连续值还是离散值，就形成各种不同形式的的傅里叶变换对。 其中离散傅里叶变换 (DFT)在实际应用中是非常重要的，利用它可以计算信号的频谱、功率谱和线性卷积等。 但是，当 N很大时，直接计算 DFT，即使使用高速的计算机，所花的时间也太多。 为此，普遍采用在 DFT的基础上发展起来的快速计算方法，即快速傅里叶变换 (FFT)，它是 DFT的快速算法。 本文中采用 FFT对 信号进行频谱分析，用以了解复杂信号的频率成分和幅值 [21]。 167。 离散傅立叶变换（ DFT） 离散傅立叶级数变换是周期序列，有 N个独立的数值，所以它的许多特性可以通过有限长序列延 拓来得到。 对于一个长度为 N的有限长序列 x(n)，也即 x(n)只在 0～（ N1） 各点上有非零值，即  其他,010),()( Nnnxnx （ 31） 把序列 x（ n）以 N为周期进行周期沿拓得到周期序列 )(~nx ，则有  其他,010),(~)( Nnnxnx （ 32） 所以，有限长序列 x（ n） 的离散傅立叶变换 （ DFT） 为 17    10 10,)()]([)( Nn knN NkWnxnxD F TkX （ 33） 逆变换为    10 10,)(1)]([)( Nn knN NkWkXNkXID F Tnx （ 34） 应当注意，如果给出的是连续信号 txa ，则首先要根据其最高频率确定采样速率 sf 以及由频率分辨率选择采样点数 N，然后对其进行软件采样 (即计算    nTxnx a ， 10  Nn )，产生对应序列 nx。 对信号 txa ，频率分辨率的选择要以能分辨开其中的三个频率对应的谱线为准则。 对周期序列，最好截取周期的整数倍进行谱分析，否则有可能产生较大的分析误差 [22]。 167。 主程序框图 开 始读 入 长 度 N调 用 信 号 产 生 子 程 序 产 生 信 号调 用 绘 图 子 程 序 （ 函 数 ） 绘 制 时 间 序 列 波 形 图调 用 F F T 子 程 序 （ 函 数 ） 计 算 信 号 的 D F T调 用 绘 制 子 程 序 （ 函 数 ） 绘 制 幅 频 曲 线结 束 图 31 系统主程序框图 18 167。 算法仿真结果 对原始信号进行频谱分析 由采样信号的文本可知，所用的心电信号的采样频率为 360Hz,设 信号长度 n=1000,则对心电信号进行频谱分析的结果如图 32 所示： 图 32 原始信号及其频谱图 对小波滤波后的信号进行频谱分析 为了从频率角度去观察滤波前后信号的频率组成，所以本设计又对上一章所设计的滤波器滤波后的信号运用频谱分析的方法进行分析，来比较滤波前后信号的频谱分布，从而来验证所设计的滤波器的滤波效果。 图 33 是对小波滤波后的信号进行的频谱分 析： 19 图 33 小波滤波后的信号及其频谱图 对巴特沃斯滤波后的信号进行频谱分析 巴特沃斯滤波适用于所有的医电信号的滤波操作，运用 butter 函数来求出所期待的传递函数的分子分母 [b,a],再用 filter（ b,a,x） 来得到滤波后的波形。 对巴特沃斯各种滤波器。