基于Matlab的IIRButterworth低通数字滤波器设计

问题：信号滤波是一个非常重要的信号处理手段，大量的信号处理系统中，信号滤波的质量将直接影响到系统整体性能。现在有家医院向你提出一个问题，她们的心电图信号记录仪器由于受到了比较严重的电源干扰(50Hz附近)，记录到的心电图ECG波形有较大的畸变，请问如何设计一个干扰滤除系统，让记录下的ECG数据继续可以进行疾病诊断。(ECG的波形信息非常重要，设计滤波器需要注意这一点)基于Matlab的IIRButterworth低通数字滤波器设计0引言心电信号是心脏电活动在人体体表的表现，它一般比较微弱，其频率为0.05Hz~100Hz（能量集中在0.05Hz~44Hz），幅度为几百微伏至几毫伏。因此，在心电信号的采集过程中，极易受到内、外界环境的干扰，其中50Hz的工频干扰尤为突出。心电图信号的干扰对心电图数据分析和压缩有一定影响。如何消除50Hz工频干扰，成为处理心电信号的首要任务。为了消除ECG数据中的主要干扰，保留有用信号或者从干扰中“抽取”有用信号，人们提出了许多方法，诸如，平滑滤波，中心频率固定的带阻滤波器，自适应滤波，低通数字滤波等等。这些信号处理方法各有优缺点，平滑滤波算法简单，处理速度快，滤波效果较好，但存在一定的削峰作用；中心频率固定的带阻滤波器原理简单，能较大程度地抑制50Hz工频干扰，但存在“群延时”现象；自适应滤波器的中心频率能够跟随工频信号的频率幅度变化而自动调节并抵消工频干扰，但需要附加参考信号通道，算法相对复杂，难以用于实时处理。考虑到工频干扰是50Hz的低频信号，假设要处理的ECG信号为0.05HZ~44HZ段信号，且Matlab软件中有一个数字滤波器软件包，可以直接进行调用处理信号，所以本文设计IIR巴特沃斯(Butterworth)低通数字滤波器来实现人体ECG信号的滤波。1数字滤波器介绍数字滤波器是一种用来过滤时间离散信号的数字系统，它可以用软件(计算机程序)或用硬件来实现，而且在两种情况下都可以用来过滤实时信号或非实时信号(记录信号)。尽管数字滤波器这个名称一直到六十年代中期才出现。但是随着科学技术的发展及计算机的更新普及，数字滤波器有着很好的发展前景。有些时候，它可以完全取代模拟滤波器。数字滤波器具有如下的优点：它没有漂移，它能处理低频信号，频率响应特征可做成极接近于理想的特性，它可以做成没有插人损耗和有线性相位特性，可相当简单地获得自适应滤波，滤波器的设计者可以控制数字字长，因而可以精确地控制滤波的精度。其中，它最重要的优点是：随着滤波器参数的改变，很容易改变滤波器的性能。这一特点允许我们只用一种程序滤波器完成多重滤波任务。本文采用双线性变换法设计IIR巴特沃斯低通数字滤波器实现了对工频干扰信号的有效滤波。2巴特沃斯（Butterworth）低通数字滤波器的特性所谓低通滤波器，就是只能让低频分量通过。在低通滤波器的设计中，有许多不同的逼近方法：巴特沃斯逼近，切比雪夫逼近，椭圆逼近方法等。其中，最简单的逼近为巴特沃斯逼近。巴特沃斯低通滤波器的平方幅度响应为2211ncHj其中，n为滤波器的阶数，cw为低通滤波器的截止频率。该滤波器具有一些特殊的性质：①对所有的n,都有当0w时,201Hj；②对所有的n,都有当cww时,20.5cHjw；③2Hjw是的单调递减函数,即不会出现幅度响应的起伏；④当n时,巴特沃斯滤波器趋向于理想的低通滤波器；⑤在0w处平方幅度响应的各级导数均存在且等于0，因此2Hjw在该点上取得最大值，且具有最大平坦特性。图1展示了2阶、4阶、8阶巴特沃斯低通滤波器的幅频特性。可见阶数n越高，其幅频特性越好，低频检测信号保真度越高，过渡带变窄，即衰减加剧，但半功率点不变。图1巴特沃斯低通滤波器的幅频特性3IIR数字滤波器的性质无限长冲激响应(IIR)数字滤波器是数字滤器的一种，数字滤波器根据其冲激响应函数的时域特性,还包括有限长冲激响应(FIR)数字滤波器。IIR数字滤波器的特征是：①具有无限持续时间冲激响应；②需要用递归模型来实现，这可以从其差分方程00MNikikynbxniaynk得出，也可以从其系数函数为：10101MiiNkkbzHzaz得出。数字巴特沃思滤波器属于IIR滤波器，该类滤波器具有特定的性质和设计方法。目前比较成熟的IIR数字滤波器设计方法有两种：1）直接法目前所用的方法主要是：零极点累试法、频域幅度平方误差最小法和时域单位脉冲响应逼近法。直接法的最大优点在于可以设计任意幅频特征的滤波器。2)间接法，目前所用的方法主要是：冲激响应不变法、阶跃响应不变法和双线性法。它们都是借助于已经成熟的现有低通滤波器原型进行设计，即对数字低通数字滤波器，先将数字低通滤波器的技术指标按希望的设计方法转换为模拟低通滤波器的技术指标，再按指定的模拟低通滤波器的类型设计模拟滤波器H(s)，然后，将模拟滤波器的系统函数H(s)从s平面转换到z平面，得到数字低通滤波器的系统函数H(z)；如所设计的数字滤波器为高通、带通或带阻滤波器，则可借助模拟滤波器的频带变换转换为低通模拟滤波器。由于直接法设计巴特沃思滤波器相对复杂，在不需要任意幅频特征的情况下，一般采用问接法，同时由于冲激响应不变法和阶跃响应不变法，从s平面转换到z平面的映射为多值映射，容易造成频谱混叠，故而本文采用不会产生频谱混叠的双线性变换法。4数字巴特沃思滤波器的双线性变换法设计设计IIR滤波器的任务就是寻求一个因果、物理上可实现的系统函数H(z)使其频率响应jwHe满足所希望得到的频域指标，即符合给定的通带截止频率、阻带截止、通带衰减和阻带衰减等指标。双线性变换法的基本步骤如下：(1)根据任务，确定性能指标：在设计之前，首根据工程实际的需要确定巴特沃思滤波器的技术指标：边界频率：pw，sw，cw；通带最大衰减p和阻带最小衰减s；(2)将数字滤波器的技术指标转换成模拟滤波器指标。双线性变换法的转换关系为：22tanTw，w为数字角频率，为模拟角频率，T为取样周期；(3)用模拟滤波器设计方法得到模拟巴特沃斯滤波器的传输函数H(s)，这可借助巴特沃思滤波器相关的严格的设计公式、现成的曲线和图表得到；(4)映射实现：利用双线性变换法关系式11211zsTz将模拟巴特沃思滤波器H(s)转换成数字巴特沃思滤波器H(z)；(5)用有限精度算法实现这个系统函数H(z)，同时要确定该滤波器的结构，因为如果在硬件上实现该滤波器必须考虑计算的复杂度、存储容量限制等因素，而用Matlab实现无此影响；(6)用合适的软、硬件技术实现该滤波器。5数字巴特沃思滤波器的Matlab实现按照传统的方法设计巴特沃思数字滤波器，不论采用采用何种方法都面临着繁杂的数学计算或反复的累加试验，以验证设计结果是否满足指标要求，工作量巨大、且效率极低。而现代分析信号时，使用软件Matlab是一种高效、简便的设计方法。Matlab是一款高性能数值计算和可视化软件，应用及其广泛。Matlab软件的诞生，尤其是数字信号处理工具箱(SignalProcessingToolbox)和滤波器设计工具箱(FilterDesignToolbox)推出后，使数字信号处理系统的分析与设计问题变得简单了,它为数字信号处理系统的设计与仿真提供了一个强有力的工具。在设计数字滤波器时，通常采用Matlab来进行辅助设计和仿真。可方便地分析显示滤波器系统的幅频、相频特性,以及对给定信号经所设计的滤波器进行滤波前和滤波后的幅频特性。由于冲击响应不变法和阶跃响应不变法的不足，下面以双线性变换法设计一个巴特沃思低通滤波器，下面以双线性变换法设计一个巴特沃思低通滤波器。这可以严格按照上述步骤进行，并且有很多文献。正是按照这个步骤进行，但通过仔细研究Matlab软件自带的信号处理工具箱下Matlab\toolbox\signal\signal\butter.m和Matlab\toolbox\signal\signal\buttord.m等文件，可以清楚看到：Matlab默认采用双线性变换法设计，因而无需完全严格按照上述步骤进行，完全可以采用更简洁的方法。所以设计滤波器可以运用Matalb中的程序实现也可以通过它中的Smulink仿真模块来搭建系统实现，具体的做法是：运行Matalb软件，在路径Smulink\SignalProcessingBlockset\Filtering\FilterDesigns\AnalogFilterDesign，可以设置highpass,lowpass,bandpass等。如下图所示。接着就可以看到滤波器设计方法，滤波器类型，阶数及边界频率的参数，如下图所示。本文采用的是数字滤波器的设计，所以选择路径Smulink\SignalProcessingBlockset\Filtering\FilterDesigns\DigtalFilterDesign或者直接在Matlab命令窗口中直接敲入FDATOOL运行即可，即通过FDATool界面设计，则可以相应的设置滤波器的类型及参数。本文采用FDATool界面设计滤波器。下面详细介绍一下FDATool界面。如下图所示。FDATool（FilterDesign&AnalysisTool）是Matlab信号处理工具箱里专用的滤波器设计分析工具，Matlab6.0以上的版本还专门增加了滤波器设计工具箱（FilterDesignToolbox）。FDATool可以设计几乎所有的基本的常规滤波器，包括FIR和IIR的各种设计方法。它操作简单，方便灵活。FilterType（滤波器类型）选项，包括Lowpass（低通）、Highpass（高通）、Bandpass（带通）、Bandstop（带阻）和特殊的FIR滤波器。DesignMethod（设计方法）选项，包括IIR滤波器的Butterworth（巴特沃思）法、ChebyshevTypeI（切比雪夫I型）法、ChebyshevTypeII（切比雪夫II型）法、Elliptic（椭圆滤波器）法和FIR滤波器的Equiripple法、Least-Squares（最小乘方）法、Window（窗函数）法。FilterOrder（滤波器阶数）选项，定义滤波器的阶数，包括SpecifyOrder（指定阶数）和MinimumOrder（最小阶数）。在SpecifyOrder中填入所要设计的滤波器的阶数（N阶滤波器，SpecifyOrder＝N-1），如果选择MinimumOrder则Matlab根据所选择的滤波器类型自动使用最小阶数。FrenquencySpecifications选项，可以详细定义频带的各参数，包括采样频率Fs和频带的截止频率。它的具体选项由FilterType选项和DesignMethod选项决定，例如Bandpass（带通）滤波器需要定义Fstop1（下阻带截止频率）、Fpass1（通带下限截止频率）、Fpass2（通带上限截止频率）、Fstop2（上阻带截止频率），而Lowpass（低通）滤波器只需要定义Fstop1、Fpass1。采用窗函数设计滤波器时，由于过渡带是由窗函数的类型和阶数所决定的，所以只需要定义通带截止频率，而不必定义阻带参数。MagnitudeSpecifications选项，可以定义幅值衰减的情况。例如设计带通滤波器时，可以定义Wstop1（频率Fstop1处的幅值衰减）、Wpass（通带范围内的幅值衰减）、Wstop2（频率Fstop2处的幅值衰减）。当采用窗函数设计时，通带截止频率处的幅值衰减固定为6dB，所以不必定义。WindowSpecifications选项，当选取采用窗函数设计时，该选项可定义，它包含了各种窗函数。在设计过程中，可以对比滤波器幅频相频特性和设计要求，随时调整参数和滤波器类型，以便得到最佳效果。最后设置参数，设计本文所需滤波器，参数选择如下图所示。如果问题中提供了被干扰的ECG信号数据