第6章 IIR数字滤波器的设计
数字滤波器是数字信号处理技术的重要内容。和模拟滤波器一样，数字滤波器的主要功能是对数字信号进行处理，最常见的处理是保留数字信号中的有用频率成分，去除信号中的无用频率成分。本章首先介绍数字滤波器的原理和分类，使大家对数字滤波器有个整体的认识，然后介绍由模拟滤波器传递函数转换为数字滤波器传递函数的方法----脉冲响应不变法和双线性变换法。在进行MATLAB设计IIR数字滤波器的讲解之前，第三节介绍滤波器特性分析及使用方法。最后大篇幅举例介绍IIR数字滤波器的经典设计、完全设计函数和直接设计方法。
6.1 概 述
6.1.1 数字滤波器的工作原理
数字滤波器是一个离散时间系统，输入x(n)是一个时间序列，输出y(n)也是一个时间序列。如数字滤波器的系统函数为H(z)，其脉冲响应为h(n)，则在时间域内存在下列的关系。
y(n)?x(n)?h(n) (6-1) 在z域内，输入和输出存在下列关系：
Y(z)?H(z)X(z) (6-2) 式中，X(z)、Y(z)分别为输入x(n)和输出y(n)的z变换。
同样在频率域内，输入和输出存在下列关系：
Y?j???X?j??H?j?? (6-3) 式中，H?j??为数字滤波器的频率特性; X?j??和Y?j??分别为x(n)和y(n)的频谱。?为数字角频率，单位rad。通常设计H?j??在某些频段的响应值为1，在某些频段的响应为0。X?j??和H?j??的乘积在频率响应为1的那些频段的值仍为X?j??，即在这些频段的振动可以无阻碍地通过滤波器，这些频带为通带。X?j??和H?j??的乘积在频率响应为0的那些频段的值不管X?j??大小如何均为零，即在这些频段里的振动不能通过滤波器，这些频带称为阻带。
为讨论方便起见，我们在本章中规定，?为数字角频率，单位为弧度(rad)，?表示模拟角频率，单位弧度/秒（rad/s）。数字角频率?在0～?范围内。
一个合适的数字滤波器系统函数H(z)可以根据需要改变输入x(n)的频率特性。经数字滤波器处理后的信号y(n)保留信号x(n)中的有用频率成分，去除无用频率成分。
6.1.2 滤波器的分类：
按时间域特性，数字滤波器可以分为无限冲激（脉冲）响应数字滤波器（Infinite impulse response digital filter,简称IIR滤波器）和有限冲激（脉冲）响应数字滤波器（Finite impulse response digital filter,简称FIR滤波器）两类。
IIR滤波器的传递函数为：
H(z)?Y(z)??h(n)z?n?X(z)n?0??bzrr?0Nk?1M?r (6-4)
1??akz?kh(n)为滤波器的脉冲响应，n=0~?均有值。M和N为分解的分子和分母多项式的系数个数。FIR滤波器的传递函数为：
Y(z)N?1 H(z)???h(n)z?n (6-5)
X(z)n?0该滤波器的脉冲响应h(n)在n=0,1,?,N-1的有限个点(N个点)上有值。注意（6-4）式中分母ak全为零时，H(z)具有全零点形式，IIR滤波器退化为FIR滤波器。
按频率特性来讲，数字滤波器和模拟滤波器一样可分为低通、高通、带通和带阻等。数字滤波器是一个离散时间系统，在频率特性中具有周期性，因此我们讨论的频率范围仅在??0~?的范围内，相应的归一化频率在0~1之间，?和1对应于Nyquist频率。
和模拟滤波器一样，理想数字滤波器的频率特性Hd?j??在通带内必须满足：
?Hd?j???K ? (6-6)
???Hj?????d?式中，K,?均为常数。
和模拟滤波器一样，数字滤波器的设计目的是使滤波器的频率特性达到所给定的性能指标。其性能指标也包括通带波纹Rp(dB)、阻带衰减Rs(dB)、通带边界频率（Hz）、阻带边界频率（Hz）等。
6.1.3 IIR数字滤波器的设计方法
IIR滤波器的设计就是根据滤波器的性能指标要求，设计滤波器的分子和分母多项式系数。它和FIR滤波器相比优点是在满足相同性能指标要求的前提下，IIR滤波器的阶数要明．．．．．．．．．．．．．．．．．．．．．．．．．．．显低于滤波器；但滤波器的相位是非线性的。 ．．．FIR．．．．．．．．IIR．．．．．．．．．．．．．．在第5章我们学习了模拟滤波器的设计。在数字滤波器设计中，首先让我们想到的是能否利用模拟滤波器的设计成果进行数字滤波器设计，答案是肯定的。现在我们所讲的IIR经典设计就是将已设计好的模拟滤波器按一定变换原理转换为数字滤波器。该方法先根据滤波器的技术指标设计出相应的模拟滤波器，然后再将设计好的模拟滤波器变换成数字滤波器。在MATLAB中，经典法设计IIR数字滤波器主要采用以下步骤： 模拟低通滤波器原型 /buttap,cheb1ap, cheb2ap,ellipap, besselap 频率变换/lp2lp,lp2hp, lp2bp,lp2bs 模拟离散化/bilinear,impinvar IIR数字滤波器 图6-1 IIR滤波器的设计步骤
由图6-1可见，IIR滤波器设计利用了模拟滤波器的设计成果。第一步和第二步在前面已详细讨论过。第二步完成后，一个达到期望性能指标的模拟滤波器（低通、高通、带通、带阻等）已经设计出来。第三步模拟离散化的主要任务就是把模拟滤波器变换成数字滤波器，即把模拟滤波器的系统函数H(s)映射为数字滤波器的系统函数H(z)。这样，数字滤波器的设计工作全部完成。
6.2 模拟滤波器到数字滤波器的转换
由前面一节我们清楚了设计IIR滤波器的非常重要的一步是将设计的模拟滤波器转换为数字滤波器。实现模拟滤波器系统传递函数s域至数字滤波器传递函数z域映射的方法有冲激（脉冲）响应不变法和双线性变换法两种方法。下面不加证明地介绍这两种方法设计要点及MATLAB实现。
6.2.1 脉冲响应不变法
所谓脉冲响应不变法就是使数字滤波器的脉冲响应序列等于模拟滤波器的脉冲响．．．．．．．．．．．．．．．．．．．．．．．．h(n)．．．．．．．．．．．．．．．应(t)的采样值，即 a．．h．．．．．．．．．． h(n)?ha?t?t?nT?ha?nT? (6-7) ．．．．．式中，T为采样周期。
因此数字滤波器的系统函数H(z)可由下式求得：
H(z)?Z?h(n)??Z?ha?nT?? (6-8) Z[-]表示对[-]的内容进行Z变换，Z变换的内容请参考相应的数字信号处理教材。
如果已经获得了满足性能指标的模拟滤波器的传递函数Ha?s?，求与之对应的数字滤波器的传递函数H(z)的方法是：
（1） 求模拟滤波器的单位脉冲响应ha?t?。
ha?t??L?1?Ha?s?? (6-9)
式中，L?1?Ha?s??表示对Ha?s?的Laplace逆变换。Laplace变换内容请参考高等数学的积 分变换或信号处理教材。
（2） 求模拟滤波器单位冲激响应ha?t?的采样值，即数字滤波器冲激响应序列h(n)。 （3） 对数字滤波器的冲激响应h(n)进行z变换，得到传递函数H(z)。
由上述方法推论出更直接地由模拟滤波器系统函数Ha?s?求出数字滤波器系统函数H(z)的步骤是：
（1） 利用部分分式展开将模拟滤波器的传递函数H(s)展开成：
Ha(s)?Rk （6-10） ?s?pk?1kN在MATLAB中这步可通过residue函数实现。
若调用residue函数的形式为[R,P,K]=residue(a,b) 形式，则将下式（传递函数形式）：
b(s)b(1)snb?b(2)snb?1???b(nb)s?b(nb?1) （6-11） Ha(s)??a(s)a(1)snb?a(2)sna?1???a(na)s?a(na?1)变换为：
Ha(s)??Rk（6-12） ?K(1)sM?K(2)sM?1???K(M?1)
s?pk?1kN这种形式为极点留数商向量形式，对于本节所讲的特定情况，K为空矩阵。
若为[b,a]=residue(R,P,K)则为上面调用形式的反过程。 （2） 将模拟极点pk变换为数字极点eNpkT即得到数字系统的传递函数
H(z)?其中T为采样间隔。
Rk （6-13） ?pkT?1zk?11?e(3)将（6-13）转换为传递函数形式，在该步骤中，可采用[R,P,K]=residue(b,a). 对于上面的步骤，MATLAB中已经提供了冲激响应不变法设计数字滤波器的函数，调用格式为：
[bz,az]=impinvar(b,a[[,Fs[，Fp])
式中，b,a为模拟滤波器分子和分母多项式系数向量；Fs为采样频率（所滤波数据），单位Hz，缺省时为1Hz。Fp为预畸变频率（Prewarped frequency）,是一个“匹配”频率,在该频率上，频率响应在变换前后和模拟频率可精确匹配。一般设计中可以不考虑。bz,az分别为数字滤波器分子和分母多项式系数向量。前面已提到过，函数输入变量中的[]表示可添加也可略去的内容。下面我们用例子说明如何使用这个函数。 【例6-1】脉冲响应不变法将模拟滤波器Ha(s)?样周期为T=0.1s。
%Samp6_1
b=[3 2];a=[2 3 1];T=0.1; %模拟滤波器分子和分母多项式系数及采样间隔 [bz1,az1]=impinvar(b,a,1/T)
程序输出为：
bz1 = 0.3000 -0.2807
az1 = 2.0000 -3.7121 1.7214
在应用冲激响应不变法设计数字滤波器时要注意它的特点。脉冲响应不变法由z?esT这一基本关系得到数字角频率?和模拟角频率?满足???T线性变换关系，T为采样间隔。这使得j?轴上每隔2?T便映射到z域中的单位圆一周。如果模拟滤波器频率响应是有限带宽的话，通过变换得到的数字滤波器的频率响应非常接近于模拟滤波器的频率响应。由于数字滤波器的频率响应是模拟滤波器频率响应的周期延拓，因此对于高通和带阻滤波器存在混叠效应，会造成频率响应失真，因此这种方法原则上只适用于有限带宽滤波器。对于高通、带阻等滤波器，由于它们高频成分不衰减，势必产生严重的混迭失真。下面所讲的双线性变换法可以弥补这方面的不足。
3s?2变换为数字滤波器H(z)，采22s?3s?16.2.2 双线性变换法
双线性变换法将s平面的整个频率轴映射到z域的一个频率周期中。因此s平面到z平面的映射是非线性的，其单值双线性映射关系为：
Ts21?z?12 (6-14) s?,z??1TT1?z1?s21?式中，T为采样周期。
因此若已知模拟滤波器的传递函数Ha?s?，将（6-14）式的第一式代入Ha?s?即可得到数字滤波器的传递函数：
H(z)?Ha(s)s?21?z?1 (6-15)
T1?z?1在双线性变换中，模拟角频率和数字角频率存在下面关系： ??可见，模拟角频率?和数字角频率?之间的关系是非线性的。 在MATLAB中，函数bilinear采用双线性变换法实现模拟s域至数字z域的映射，直接用于模拟滤波器变换为数字滤波器。其调用方式为： [zd,pd,kd]=bilinear(z,p,k,Fs) [numd,dend]=bilinear(num,den,Fs)
2??Ttan,??2arctan (6-16) T22式中，z,p分别为模拟滤波器零点、极点列向量；k为模拟滤波器的增益；Fs为采样频率，单位Hz。zd,pd,kd为数字滤波器的零极点和增益。num,den分别为模拟滤波器传递函数分子和分母多项式系数向量，模拟滤波器传递函数具有下面的形式：
num(s)num(1)snn???num(nn)s?num(nn?1) (6-17) H(s)??den(s)den(1)snd???den(nd)s?den(nd?1)numd和dend分别为数字滤波器传递函数分子和分母多项式系数向量。 【例6-2】用双线性变换法将模拟滤波器Ha(s)?样周期（间隔）T=0.1s。
%Samp6_2
b=[3 2];a=[2 3 1];T=0.1; %模拟滤波器分子和分母多项式的系数，采样间隔 [bz1,az1]=bilinear(b,a,1/T) %将模拟滤波器传递函数转换为数字滤波器传递函数 程序输出为：
bz1 = 0.0720 0.0046 -0.0674 az1 = 1.0000 -1.8560 0.8606
双线性变换法克服了脉冲响应不变法的频谱混迭问题，其幅值逼近程度好，可适用于高通、带阻等各种类型滤波器的设计。s域和z域对应关系也简单。缺点是频率变换的非线性导致数字滤波器与模拟滤波器在幅度和频率的对应关系上发生畸变。但一般滤波器的幅频响应具有分段常数的特点，即滤波器允许某一频段信号通过，而不允许另外频段的信号通过的特点，故变换后这一特点仍保留，影响不大。由数字边界频率计算模拟边界频率时，不是按线性关系进行的，这就是所谓的“预畸变”。但如果给定预畸变频率为边界频率，经预畸变频率校正则可以保证所要设计的模拟边界频率精确映射在所要求的数字边界频率上。
3s?2变换为数字滤波器H(z)，采22s?3s?16.3 滤波器特性及使用函数
为了形象地看到滤波器的性能，需要绘出滤波器的幅频和相频特性曲线。本节给出MATLAB信号处理工具箱中提供分析滤波器特性的工具函数，包括时域分析（对任意输入的响应、脉冲响应、阶跃响应等）、频域分析（幅值响应、相位响应、零极点位置）和其他特性（如群延迟等）。滤波器时域和频域分析是设计各类滤波器、评价滤波器性能以及滤波器应用的基础，本节将具体介绍这些函数。
6.3.1 滤波器的特性函数
(1) freqz
在上一章我们讲过，对于模拟滤波器，可以用freqs求解滤波器的频率响应。与之对应的函数freqz用于求数字滤波器的频率响应，其调用格式为：
[[h,w]=]freqz(b,a,n[,’whole’]);或[h,f]=freqz(b,a,n[,’whole’],Fs);
式中，b,a为数字滤波器分子和分母多项式的系数，n为复数频率的响应点数，为整数，最好为2的幂，缺省时为512；Fs为采样频率，单位Hz。如果给定该值，则f位置输出为频率Hz，若没有给定，则按角频率（Angular frequency）给定f的频率矢量；’whole’表示返回的频率f或w值包含z平面整个单位圆频率矢量，即0~2?；缺省时，频率f或w值包含z平面上半单位圆（0~?）之间等间距n个点频率矢量。h为复频率响应；w为n点频率向量（单位rad）；f为n点频率向量(Hz)。函数返回值缺省时，绘制幅频响应和相频响应图。该函数适用于下面形式的数字滤波器：
B(z)b(1)?b(2)z?1???b(nb?1)z?nb (6-18) H(z)???1?naA(z)a(1)?a(2)z???a(na?1)z函数freqz输出的频率向量在0~?。为了获得一个滤波器真正的相频特性图，要对相位角进行解缠绕。为此MATLAB提供了一个函数unwrap来解决这个问题，P=unwrap(angle(H))。
(2) impz
impz用于产生数字滤波器的脉冲响应。调用格式为：
[[h,t]=]impz(b,a[,n,Fs])
式中，b,a分别为滤波器分子和分母多项式系数向量；n为采样点数；Fs为采样频率，缺省值为1；h为滤波器单位脉冲响应向量；t为和h对应的时间向量。当函数输出缺省时，绘制滤波器脉冲响应图；当n缺省时，函数自动选择n值。
(3)零极点图
滤波器的零极点位置决定了滤波器稳定性和性能，因此考察滤波器的零极点的位置是分析滤波器特性的重要方面之一。MTALAB信号处理工具箱提供绘制数字滤波器零极点位置图的工具zplane,调用格式为：
zplane(z,p)或zplane(b,a)
式中，z,p为零极点向量(为复数)，b,a为滤波器分子和分母多项式的系数（为实数）。函数在z平面绘出零点和极点。极点用’×’表示，零点用’o’表示。 【例6-3】设计一个5阶Butterworth低通数字滤波器，截止频率为0.2，绘制其零极点图。 %Samp6_3
[z,p,k]=butter(5,0.2); %设计一个Butterworth数字低通滤波器 zplane(z,p); %绘制其零极点图 xlabel('实部');ylabel('虚部');
程序的输出结果为图6-2。可以明显看出该滤波器零点和极点的分布，从而可以进行进一步的分析。由于该滤波器为一个5阶的Butterworth滤波器，故其极点的个数为5，并且极点为两两共轭，分布在单位圆内。由滤波器的分析理论知，只有极点分布在单位圆内才能保证滤波器的稳定。大家可以换用其他滤波器观看其极点在z平面上的分布。
10.80.60.40.2虚部0-0.2-0.4-0.6-0.8-15-1-0.50实部0.51图6-2 例6-3设计的Butterworth滤波器的零极点图
(4)群延迟
信号传输的不失真条件之一为：滤波器相频特性是一条经过原点的直线，即???????td，
td为常数。但一般滤波器不满足这个条件，衡量实际滤波器相位平均延迟的物理量是群延
迟。群延迟定义为信号通过滤波器的延迟随频率变化的函数，即滤波器相频特性图上切线的负斜率：
?g?????d???? （6-19） d?MATLAB信号处理工具箱提供计算群延迟函数grpdelay，调用格式为： [gd,w]=grpdelay(b,a,n[,’whole’]) [gd,f]=grpdelay(b,a,n[,’whole’],Fs) gd=grpdelay(b,a,w) gd=grpdelay(b,a,f,Fs) grpdelay 其中，gd为群延迟；其他各项意义同函数freqz,函数输出项缺省时，绘制群延迟图。’whole’参数表示绘制包括大于Nyquist频率的一个周期的群延迟。
【例6-4】绘制截止频率为0.2的5阶Butterworth低通滤波器的群延迟。
%Samp6_4
[b,a]=butter(5,0.2); %设计Butterworth低通滤波器 grpdelay(b,a,128) %绘制群延迟曲线
xlabel('归一化频率'); ylabel('群延迟/采样数')
9876群延迟/采样数54321000.10.20.30.40.50.6归一化频率0.70.80.91图6-3 例6-4设计的Butterworth的群延迟图
程序运行结果为图6-3。该图说明，某一频率的波具有最大的延迟。
6.3.2 滤波器的使用函数
（1）
filter函数用来实现数字滤波器对数据的滤波，函数的调用格式为：
y=filter(b,a,x)
其中，b,a分别为滤波器传递函数H(z)的分子和分母多项式系数。x为滤波器的输入。y为滤波器的输出。y为与x具有相同大小的向量。 （2）
filtfilt函数实现零相位前向与后向结合的滤波。调用格式为：
y=filtfilt(b,a,x)
式中，b,a分别为滤波器传递函数H(z)的分子和分母多项式系数。x为滤波器的输入，为值向量。y为滤波器的输出。该函数对序列x进行正常的正向滤波后，将滤波后的输出翻转重新用该滤波器进行滤波，第二次滤波后的输出序列的翻转即得到零相位的滤波输出。这样就可以把延迟后的相位校正至零。但该函数只能用于数字滤波器，FIR滤波器或IIR滤波器均能使用。
6.4 经典设计法
6.4.1 IIR数字滤波器经典设计法
IIR数字滤波器经典设计法的一般步骤为：
（1） 根据给定的性能指标和方法不同，首先对设计性能指标中的频率指标，如数字边界
频率进行变换，转换后的模拟频率指标作为模拟滤波器原型设计的性能指标。
（2） 估计模拟滤波器最小阶数和截止频率，利用MATLAB工具函数buttord、cheb1ord、
cheb2ord、ellipord等。
（3） 设计模拟低通滤波器原型。利用MATLAB工具函数buttap、cheb1ap、cheb2ap、ellipap
等。
（4） 由模拟原型低通滤波器经频率变换获得模拟滤波器（低通、高通、带通、带阻等），
利用MATLAB工具函数lp2lp、lp2hp、lp2bp、lp2bs。
（5） 将模拟滤波器离散化获得IIR数字滤波器，利用MATLAB工具函数bilinear或impinvar。
后面的几个步骤在前面的章节中已经论述，这里主要介绍第一步，关于设计性能指标的转换。 设计IIR滤波器时，给出的性能指标通常分数字指标和模拟指标两种。数字性能指标给出通带截止频率?p，阻带起始频率?s，通带波纹Rp，阻带衰减Rs等。数字频率?p和?s的取值范围为0~?，单位弧度。而MATLAB工具函数常采用归一化频率，?p和?s的取值范围为0~1，对应于0～?，此时需进行转换。
模拟性能指标给出通带截止频率?p，阻带起始频率?s，通带波纹Rp，阻带衰减Rs等。模拟频率?p和?s单位为弧度/秒（rad/s）。
MATLAB信号处理工具箱中，设计性能指标的转换应根据不同设计方法进行不同处理。下面就举例介绍这些方法。
6.4.2 冲激（脉冲）响应不变法采用???T来将?p和?s变换为?p和?s
【例6-5】用脉冲响应不变法设计一个Butterworth低通数字滤波器，使其特征逼近一个低通Butterworth模拟滤波器的下列性能指标：通带截止频率?p?2??2000rad/s，通带波纹Rp小于3dB，阻带边界频率为?s?2??3000rad/s，阻带衰减大于15dB，采样频率Fs=10000Hz。假设一个信号x(t)?sin2?f1t?0.5cos2?f2t，其中f1=1000Hz,f2=4000Hz。试将原信号与通过该滤波器的输出信号进行比较。
%Samp6_5
Wp=2000*2*pi;Ws=3000*2*pi; %滤波器截止频率 Rp=3;Rs=15; %通带波纹和阻带衰减 Fs=10000; %采样频率
Nn=128; %调用freqz所用的频率点数
[N,Wn]=buttord(Wp,Ws,Rp,Rs,'s'); %模拟滤波器的最小阶数 [z,p,k]=buttap(N); %设计模拟低通原型Butterworth滤波器 [Bap,Aap]=zp2tf(z,p,k); %将零点极点增益形式转换为传递函数形式 [b,a]=lp2lp(Bap,Aap,Wn); %进行频率转换
[bz,az]=impinvar(b,a,Fs); %运用脉冲响应不变法得到数字滤波器的传递函数 figure(1)
[H,f]=freqz(bz,az,Nn,Fs); %求解数字滤波器的幅频特性和相频特性 subplot(2,1,1),plot(f,20*log10(abs(H))) xlabel('频率/Hz');ylabel('振幅/dB');grid on; subplot(2,1,2),plot(f,180/pi*unwrap(angle(H))) xlabel('频率/Hz');ylabel('相位/^o');grid on; figure(2)
f1=1000;f2=4000; %输入信号的频率 N=100; %数据长度
dt=1/Fs;n=0:N-1;t=n*dt; %采样间隔和时间序列 x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t); %滤波器输入信号 subplot(2,1,1),plot(t,x),title('输入信号') %绘制输入信号 y=filtfilt(bz,az,x); %用函数filtfilt对输入信号进行滤波 y1=filter(bz,az,x); %用filter函数对输入信号滤波
subplot(2,1,2),plot(t,y,t,y1,':'),title('输出信号'),xlabel('时间/s') legend( ' filtfilt ', 'filter') %加图例
程序的运行结果为图6-4和图6-5。由图6-4上图可知，在小于2000Hz处的衰减小于3dB，而在大于3000Hz处衰减大于15dB，满足滤波器的设计指标。由图6-5可见滤波器对含有1000Hz和4000Hz频率成分的信号进行了滤波，滤除了4000Hz的信号。由程序的输出还可以看出，采用filtfilt函数，输出的1000Hz信号（实线）与输入1000Hz的信号相位一致，即经过滤波后并没有改变信号波形形状。而运用filter函数滤波后（虚线）有一些延迟，改变了信号的形状。
100振幅/dB-10-20-30-40050010001500200025003000频率/Hz35004000450050000-100相位/o-200-300-400050010001500200025003000频率/Hz3500400045005000图 6-4 例6-5所设计Butterworth滤波器的频率响应。
上图：幅频特性；下图：相频特性
输入信号1.510.50-0.5-100.0010.0020.0030.0040.0050.0060.0070.0080.009输出信号1.510.50-0.5-100.0010.0020.0030.0040.0050.0060.0070.0080.009时间/s图6-5例6-5设计滤波器的输入和输出信号
0.01 filtfilt filter0.01
输出信号中实线表示运用filtfilt函数的输出，虚线表示运用filter函数的输出
【例6-6】用脉冲响应不变法设计Butterworth低通数字滤波器，要求通带频率为0???0.2?，通带波纹小于1dB，阻带在0.3?????内，幅度衰减大于15dB，采样周期T=0.01s。假设一个信号x(t)?sin2?f1t?0.5cos2?f2t，其中f1=5Hz,f2=30Hz。试将原信号与经过该滤波器的输出信号进行比较。
注意该题给出的通带边界频率和阻带边界频率均为数字频率，因此设计时首先要将其转换为模拟频率。
%Samp6_6
wp=0.2*pi;ws=0.3*pi;Rp=1;Rs=15; %数字滤波器截止频率、通带波纹和阻带衰减 T=0.01;Nn=128; %采样间隔
Wp=wp/T;Ws=ws/T; %得到模拟滤波器的频率—采用脉冲响应不变法的频率转换形式 [N,Wn]=buttord(Wp,Ws,Rp,Rs,'s'); %计算模拟滤波器的最小阶数 [z,p,k]=buttap(N); %设计低通原型数字滤波器
[Bap,Aap]=zp2tf(z,p,k); %零点极点增益形式转换为传递函数形式 [b,a]=lp2lp(Bap,Aap,Wn); %低通滤波器频率转换
[bz,az]=impinvar(b,a,1/T); %脉冲响应不变法设计数字滤波器传递函数 figure(1)
[H,f]=freqz(bz,az,Nn,1/T); %输出幅频响应和相频响应
subplot(2,1,1),plot(f,20*log10(abs(H))); xlabel('频率/Hz');ylabel('振幅/dB');grid on; subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)))
xlabel('频率/Hz');ylabel('相位/^o');grid on; figure(2)
f1=5;f2=30; %输入信号含有的频率 N=100; %数据点数
n=0:N-1;t=n*T; %时间序列
x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t); %输入信号 subplot(2,1,1),plot(t,x),title('输入信号') y=filtfilt(bz,az,x); %对信号进行滤波
subplot(2,1,2),plot(t,y),title('输出信号'),xlabel('时间/s')
程序的运行结果为图6-6和图6-7。该例所要求的通带频率为0???0.2?，而该滤波器的采样间隔为0.01s，采样频率为100Hz，即2?对应于100Hz，则0.2?对应于10Hz，即该滤波器的通带范围为0~10Hz。观看图6-6上图，其通带范围最大衰减小于1dB，与该分析一致。题意要求阻带在0.3?????内，幅度衰减大于15dB。其中0.3?对应于15Hz，幅度衰减大于15dB。完全符合滤波器设计的要求。测试信号中含有5Hz和30Hz的信号，显然，5Hz可以通过该滤波器，而30Hz的信号不能通过该滤波器。从输出信号可以看出，滤波器滤去了30Hz的高频信号，达到了滤波要求。
0-20振幅/dB-40-60-80051015202530频率/Hz354045500相位/o-200-400-600051015202530频率/Hz35404550图6-6 例6-6设计的滤波器的频率特性 上图：幅频特性；下图：相频特性
输入信号210-1-200.10.20.30.40.5输出信号0.60.70.80.91210-1-200.10.20.30.40.5时间/s0.60.70.80.91图6-7例6-6设计滤波器的输入和输出信号
6.4.3 双线性变换法采用??2?tan将?p和?s变换为?p和?s T2【例6-7】用双线性变换法设计一个椭圆低通滤波器，其性能指标同例6-6.
%Samp6_7
wp=0.2*pi;ws=0.3*pi;Rp=1;Rs=15; %数字滤波器截止频率通带波纹和阻带衰减 Fs=100;Ts=1/Fs;Nn=128; %采样频率
Wp=2/Ts*tan(wp/2.);Ws=2/Ts*tan(ws/2.); %按频率转换公式进行转换 [N,Wn]=ellipord(Wp,Ws,Rp,Rs,'s'); %计算模拟滤波器的最小阶数 [z,p,k]=ellipap(N,Rp,Rs); %设计模拟原型滤波器
[Bap,Aap]=zp2tf(z,p,k); %零点极点增益形式转换为传递函数形式 [b,a]=lp2lp(Bap,Aap,Wn); %低通转换为低通滤波器的频率转换
[bz,az]=bilinear(b,a,Fs); %运用双线性变换法得到数字滤波器传递函数 [H,f]=freqz(bz,az,Nn,Fs); %求出频率特性 subplot(2,1,1),plot(f,20*log10(abs(H))); xlabel('频率/Hz');ylabel('振幅/dB');grid on; subplot(2,1,2),plot(f,180/pi*unwrap(angle(H))) xlabel('频率/Hz');ylabel('相位/^o');grid on;
0振幅/dB-20-40-60051015202530频率/Hz354045500-50相位/o-100-150-200-250051015202530频率/Hz35404550图6-8 例6-7设计椭圆低通滤波器的频率特性
上图：幅频特性；下图：相频特性
程序运行结果见图6-8。在10Hz以前，衰减小于1dB，在15Hz以后衰减均大于15dB,即性能指标完全满足滤波器的设计要求。
6.5 IIR滤波器的完全设计函数
以上介绍了IIR滤波器设计原理和基本方法步骤，并给出了一些例子说明如何用MATLAB编程分步实现这些步骤。由这些步骤可知我们必须多次调用MATLAB信号处理工具箱中的基本工具函数。
实际上，MATLAB信号处理工具箱还提供了IIR滤波器设计的完全工具函数，用户只要调用这些工具函数即可一次性完成设计，而不需要调用那些基本工具函数分步实现。 IIR滤波器设计的完全工具函数有butter,cheby1,cheby2,ellip。这些工具函数既可用于设计模拟滤波器，也可用于设计数字滤波器。关于这些函数设计模拟滤波器已在前面章节中述及，这里介绍这些函数在IIR数字滤波器中的应用。在这两类滤波器设计中，这些工具函数调用格式基本相同，只是在频率处理上有所不同。
在MATLAB滤波器设计工具箱中，数字滤波器采用归一化频率，取值为0~1之间，归一化频率1对应的数字角频率为?，对应的真实频率为采样频率的一半。在应用MATLAB工具函数设计数字滤波器时应注意这一点。 数字IIR滤波器的完全设计函数有： [b,a]=butter(n,wn[,'ftype']) [z,p,k]=butter(n,wn[, 'ftype']) [b,a]=cheby1(n,Rp,wn[,'ftype']) [z,p,k]=cheby1(n, Rp,wn[,'ftype']) [b,a]=cheby2(n,Rs,wn[,'ftype']) [z,p,k]=cheby2(n, Rs,wn[,'ftype']) [b,a]=ellip(n,Rp,Rs,wn[,'ftype']) [z,p,k]=ellip(n, Rp,Rs,wn[,'ftype'])
在上面的调用方式中，n为滤波器的阶数，wn为滤波器的截止频率，取值为0~1。需根据采样频率Fs来定，如滤波器的截止频率为Fc（Hz），则wn的计算公式为：
wn?2*Fc （6-20） Fs这样就转换为0~1的归一化频率。其中wp,ws等边界频率都要根据此公式进行转换。 'ftype'滤波器的类型为：
‘high’为高通滤波器，截止频率为wn.
‘stop’为带阻滤波器，截止频率为wn=[w1,w2] (w a,b分别为滤波器传递函数分子和分母多项式系数向量；z,p,k分别为滤波器的零极点和增益。Rp,Rs分别为所设计滤波器的通带波纹和阻带衰减，单位为dB。 设计好的数字滤波器传递函数具有下面形式： B(z)b(1)?b(2)z?1???b(n?1)z?n (6-21) H(z)??A(z)a(1)?a(2)z?1???a(m?1)z?m上述函数采用双线性变换法和频率的预畸变处理将模拟滤波器离散化为数字滤波器，同时保证模拟滤波器和数字滤波器在wn(或w1,w2)处具有相同的幅频响应。 设计时应注意真实频率和MATLAB归一化数字频率之间的转换，即6-20式的应用。在进行IIR数字滤波器设计之前，请大家注意在模拟滤波器中我们用求取最小阶数和截止频率的函数如buttord,cheb1ord,cheb2ord,ellipord，这些函数完全可以用于IIR数字滤波器设计中，只不过在模拟滤波器设计中需加可选项's'，在数字滤波器中则不加该项。另外，输出的截止频率也是归一化频率（归一化为0~1）。 【例6-8】设计一个Butterworth高通数字滤波器，通带边界频率为300Hz,阻带边界频率为 200Hz，通带波纹小于1dB，阻带衰减大于20dB,采样频率为1000Hz。假设一个信号x(t)?sin2?f1t?0.5cos2?f2t，其中f1=10Hz,f2=400Hz。试将原信号与通过该滤波器的输出信号进行比较。 %Samp6_8 Fs=1000; %采样频率 wp=300*2/Fs;ws=200*2/Fs; %根据采样频率将滤波器边界频率进行转换（6-19式） Rp=1;Rs=20; %通带波纹和阻带衰减 Nn=128; %显示滤波器频率特性的数据长度 [N,Wn]=buttord(wp,ws,Rp,Rs); %求得数字滤波器的最小阶数和截止频率（归一化频率） [b,a]=butter(N,Wn,'high'); %设计Butterworth高通数字滤波器 figure(1) [H,f]=freqz(b,a,Nn,Fs); %用Nn点绘出频率特性 subplot(2,1,1),plot(f,20*log10(abs(H))); xlabel('频率/Hz');ylabel('振幅/dB');grid on; subplot(2,1,2),plot(f,180/pi*unwrap(angle(H))) xlabel('频率/Hz');ylabel('相位/^o');grid on; n=0:127; dt=1/Fs;t=n*dt; %时间序列 f1=10;f2=400; %输入信号频率 x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t); %输入信号 figure(2) subplot(2,1,1),plot(t,x); title('输入信号') %绘制输入信号 y=filter(b,a,x); %对输入信号进行滤波 subplot(2,1,2),plot(t,y),title('输出信号') %绘制输出信号 xlabel('时间/s') 0-50振幅/dB-100-150-200050100150200250300频率/Hz3504004505001000相位/o-100-200-300-400050100150200250300频率/Hz350400450500图6-9例6-8所设计的高通Butterworth滤波器的频率特性 上图：幅频特性；下图：相频特性 输入信号210-1-200.020.040.060.080.10.120.14输出信号10.50-0.500.020.040.060.08时间/s0.10.120.14图6-10 例6-8所设计滤波器的输入和输出信号 程序输出结果为图6-9和图6-10。由图6-9上图可以看到所设计滤波器在大于300Hz为通带，其衰减均小于1dB；小于200Hz为阻带，其衰减大于20dB，完全符合设计要求，但相频特性是非线性的（图6-9下图）。由图6-10可以看出，当滤波器输入10Hz和400Hz两种信号后，滤波器滤除了10Hz的信号，使得400Hz的信号通过了滤波器，起到了滤波的效果。 【例6-9】设计一个带通Chebyshev I型数字滤波器，通带为100Hz~200Hz,过渡带宽均为50Hz,通带波纹小于1dB，阻带衰减大于30dB，采样频率Fs=1000Hz。假设一个信号x(t)?sin2?f1t?0.3cos2?f2t?0.1sin2?f3t，其中f1=30Hz,f2=100Hz,f3=270Hz。试将原信号与经过该滤波器滤波后输出信号进行比较。 该滤波器的通带范围为100~200Hz，过渡带宽为50Hz，因此，通带边界频率为100和200Hz，阻带边界频率为50和250Hz。 %Samp6_9 Fs=1000; %采样频率 wp=[100 200]*2/Fs; %通带边界频率（归一化频率）（6-20式） ws=[50 250]*2/Fs; %阻带边界频率（归一化频率）（6-20式） Rp=1;Rs=30;Nn=128; %通带波纹和阻带衰减以及绘制频率特性的数据点数 [N,Wn]=cheb1ord(wp,ws,Rp,Rs);%求得数字滤波器的最小阶数和归一化截止频率 [b,a]=cheby1(N,Rp,Wn); %按最小阶数、通带波纹和截止频率设计数字滤波器 figure(1) [H,f]=freqz(b,a,Nn,Fs); %求得滤波器的频率特性 subplot(2,1,1),plot(f,20*log10(abs(H))); xlabel('频率/Hz');ylabel('振幅/dB');grid on; subplot(2,1,2),plot(f,180/pi*unwrap(angle(H))) xlabel('频率/Hz');ylabel('相位/^o');grid on; figure(2) f1=30;f2=100;f3=270; %输入信号的三种频率成分 N=100; %输入信号的数据点数 dt=1/Fs;n=0:N-1;t=n*dt; %时间序列 x=sin(2*pi*f1*t)+0.3*cos(2*pi*f2*t)+0.1*sin(2*pi*f3*t); %输入信号 subplot(2,1,1),plot(t,x),title('输入信号') %绘制输入信号 y=filtfilt(b,a,x); %对输入信号进行滤波 subplot(2,1,2),plot(t,y) %绘制输出信号 ylim([-0.2 0.3]) title('输出信号'),xlabel('时间/s') 0-100振幅/dB-200-300-400050100150200250300频率/Hz350400450500400200相位/o0-200-400050100150200250300频率/Hz350400450500图6-11 例6-9所设计滤波器的频率特性 上图：幅频特性；下图：相频特性 输入信号210-1-200.010.020.030.040.050.060.070.080.090.1输出信号0.30.20.10-0.1-0.200.010.020.030.040.050.06时间/s0.070.080.090.1 图6-12 例6-9所设计滤波器的输入和输出信号 程序运行结果为图6-11和图6-12。由图6-11上图可看出，在100~200Hz为通带，其衰减小于1dB；在50Hz以下和250Hz以上，其衰减大于30dB，完全符合滤波器的设计要求。图6-11下图表明滤波器的相位为非线性。将30Hz、100Hz和270Hz的合成振动输入滤波器后，可以看到完全滤除了在阻带范围内的30Hz和270Hz的振动，起到了滤波效果。 【例6-10】设计一个带通Chebyshev II型滤波器，设计参数与例6-9相同。假设一个信号x(t)?sin2?f1t?0.8cos2?f2t，其中f1=10Hz,f2=100Hz。试将原信号与通过该滤波器的输出信号进行比较。采样频率为1000Hz。 %Samp6_10 Fs=1000; %采样频率 wp=[100 200]*2/Fs;ws=[50 250]*2/Fs; %通带和阻带边界频率（归一化频率）（6-20式） Rp=1;Rs=30;Nn=128; %通带波纹和阻带衰减以及绘制频率特性的数据点数 [N,Wn]=cheb2ord(wp,ws,Rp,Rs); %求取数字滤波器的最小阶数和归一化截止频率 [b,a]=cheby2(N,Rs,Wn); %按最小阶数截止频率和阻带衰减设计数字滤波器 figure(1) [H,f]=freqz(b,a,Nn,Fs); %按传递函数系数、数据点数和采样频率求频率特性 subplot(2,1,1),plot(f,20*log10(abs(H))); xlabel('频率/Hz');ylabel('振幅/dB');grid on; subplot(2,1,2),plot(f,180/pi*unwrap(angle(H))) xlabel('频率/Hz');ylabel('相位/^o');grid on; figure(2) f1=10;f2=100; %输入信号的频率成分 N=500; %输入信号的数据点数 dt=1/Fs;n=0:N-1;t=n*dt; %时间序列 x=sin(2*pi*f1*t)+0.8*cos(2*pi*f2*t); %输入信号 subplot(2,1,1),plot(t,x),title('输入信号') %绘制输入信号 y=filtfilt(b,a,x); %对输入信号进行滤波，输出为y. subplot(2,1,2),plot(t,y) %绘制输出信号波形 title('输出信号'),xlabel('时间/s') 0-20振幅/dB-40-60-80050100150200250300频率/Hz350400450500400200相位/o0-200-400050100150200250300频率/Hz350400450500图6-13 例6-10所设计滤波器的频率特性 上图：幅频特性；下图：相频特性 输入信号210-1-200.050.10.150.20.25输出信号0.30.350.40.450.510.50-0.5-100.050.10.150.20.25时间/s0.30.350.40.450.5 图6-14例6-10所设计滤波器的输入和输出信号 程序的运行结果为图6-13和图6-14。可以看出，100~200Hz之间为通带，其衰减不大于1dB；在50Hz以下和250Hz以上，衰减大于30dB，完全符合滤波器的设计要求。但在阻带内有振荡，这是Chebyshev II型滤波器的特点。当滤波器输入10Hz和100Hz两种信号的合成时，滤波器可以滤除处于阻带内的10Hz的低频信号，达到滤波的目的。 【例6-11】设计一个带阻椭圆滤波器，阻带频率从100Hz~200Hz,通带波纹小于1dB,阻带衰减为50dB,两边过渡带宽为50Hz,采样频率1000Hz。假设一个信号x(t)?sin2?f1t?0.5cos2?f2t?0.5sin2?f3t，其中f1=50Hz,f2=200Hz,f3=270Hz。试将原信号与通过该滤波器的输出信号进行比较。 %Samp6_11 Fs=1000; %采样频率 ws=[100 200]*2/Fs;wp=[50 250]*2/Fs; %通带和阻带边界频率（归一化频率） Rp=1;Rs=50;Nn=512; %通带波纹和阻带衰减以及绘制频率特性的数据点数 [N,Wn]=ellipord(wp,ws,Rp,Rs); %求取数字滤波器的最小阶数和归一化截止频率 [b,a]=ellip(N,Rp,Rs,Wn,'stop'); %按最小阶数、截止频率、通带波纹和阻带衰减设计滤波器 figure(1) [H,f]=freqz(b,a,Nn,Fs); %按传递函数系数、数据点数和采样频率求得滤波器的频率特性 subplot(2,1,1),plot(f,20*log10(abs(H))); xlabel('频率/Hz');ylabel('振幅/dB');grid on; subplot(2,1,2),plot(f,180/pi*unwrap(angle(H))) xlabel('频率/Hz');ylabel('相位/^o');grid on; figure(2) f1=50;f2=200;f3=270; %输入信号的频率成分 N=100; %输入信号的数据点数 dt=1/Fs;n=0:N-1;t=n*dt;%时间序列 x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t)+0.5*sin(2*pi*f3*t); %输入信号 subplot(2,1,1),plot(t,x),title('输入信号') %绘制输入信号 y=filtfilt(b,a,x); %对输入信号滤波 subplot(2,1,2),plot(t,y) %绘制输出信号 title('输出信号'),xlabel('时间/s') 0-20振幅/dB-40-60-80-100050100150200250300频率/Hz350400450500500相位/o0-500050100150200250300频率/Hz350400450500图6-15 例6-11所设计滤波器的频率特性 上图：幅频特性；下图：相频特性 输入信号210-1-200.010.020.030.040.050.060.070.080.090.1输出信号1.510.50-0.5-100.010.020.030.040.050.06时间/s0.070.080.090.1 图6-16 例6-11所设计滤波器的输入和输出信号 程序运行结果为图6-15和图6-16。可以看到，在100~200Hz为阻带，衰减均大于50dB，在其他频率处为通带，衰减小于1dB，完全符合滤波器的设计要求。但在通带和阻带内有振荡，这是椭圆滤波器的特点。当滤波器输入50Hz、200Hz和270Hz的合成信号时，滤波器可以滤除处于阻带内的200Hz信号，只剩下50Hz和270Hz的合成信号，达到滤波的目的。 6.6 IIR滤波器直接设计 IIR数字滤波器的经典设计法只限于几种标准的低通、高通、带通、带阻滤波器，而对于任意形状或多频带的滤波器的设计是无能为力的。 如果所设计的IIR滤波器幅频特性比较复杂，可采用最小二乘法拟合给定幅频响应。使设计的滤波器幅频特性逼近期望的频率特性，这种方法称为IIR滤波器的直接设计法。 MATLAB信号处理工具箱函数yulewalk采用直接法设计IIR数字滤波器，函数调用格式为： [b,a]=yulewalk(n,f,m) 这里,n为滤波器的阶数；f为给定的频率点向量，为归一化频率，取值范围为0~1,f的第一个频率点必须是0，最后一个频率点必须为1。其中1对应于Nyquist频率。在使用滤波器时，根据数据采样频率确定数字滤波器的通带和阻带在对此信号滤波的频率范围。f向量的频率点必须是递增的；m为和频率向量f对应的理想幅值响应向量，m和f必须是相同维数向量。b,a分别是所设计滤波器的分子和分母多项式系数向量。IIR滤波器的传递函数具有下面的形式： B(z)b(1)?b(2)z?1???b(n?1)z?n (6-22) H(z)???1?mA(z)a(1)?a(2)z???a(m?1)z在定义频率响应时，应避免通带至阻带的过渡段形状过分尖锐。通常需要调整过渡带的斜率来做到这点。 函数yulewalk首先计算给定幅频响应Fourier逆变换和相关系数，再采用修正的yule-walker方程计算滤波器传递函数分母多项式系数。 函数yulewalk采取下面步骤计算分子多项式：（1）计算与分子多项式相应的幅值平方响应的辅助式；（2）由辅助分子和分母多项式计算完全的频率响应。（3）计算滤波器的脉冲响应。（4）采用最小二乘法拟合脉冲响应最终求得滤波器的分子多项式系数。 注意：yulewalk不能用来设计给定相位指标的滤波器。 【例6-12】 用直接法设计一个10阶多频带数字滤波器，幅频响应值如下：f=[0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0], m=[0 0 1 1 0 0 1 1 1 0 0]。假设一个信号x(t)?sin2?f1t?0.5cos2?f2t，其中f1=6Hz,f2=17Hz。试将原信号与通过该滤波器的输出信号进行比较。 %Samp6_12 Order=10; %滤波器的阶数 f=0:0.1:1; %归一化频率点 m=[0 0 1 1 0 0 1 1 1 0 0]; %幅度点 [b,a]=yulewalk(Order,f,m); %设计滤波器 [h,w]=freqz(b,a,128); %计算128个点的频率特性 figure(1) plot(f,m,'b-',w/pi,abs(h),'k:'); %绘制理想滤波器和设计滤波器的幅频特性 xlabel('归一化频率');ylabel('振幅'); title('运用yulewalk方法设计IIR滤波器'); legend('理想特性','实际设计',1); %给定图例 figure(2) Fs=50; %信号采样频率 f1=6;f2=17; %信号的频率成分 N=100; %数据点数 dt=1/Fs;n=0:N-1;t=n*dt; %时间序列 x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t); %输入信号 subplot(2,1,1),plot(t,x),title('输入信号') %绘制输入信号 y=filtfilt(b,a,x); %对信号进行滤波 subplot(2,1,2),plot(t,y) %绘制输出信号 title('输出信号'),xlabel('时间/s') 运用yulewalk方法设计IIR滤波器1.4理想特性实际设计1.210.8振幅0.60.40.2000.10.20.30.40.50.6归一化频率0.70.80.91图6-17例6-12所设计滤波器（虚线）及理想滤波器（实线）的幅频特性 输入信号210-1-200.20.40.60.81输出信号1.21.41.61.82210-1-200.20.40.60.81时间/s1.21.41.61.82图6-18例6-12所设计滤波器的输入和输出信号 程序运行结果见图6-17和图6-18。由图6-17可见，设计滤波器的幅频响应与理想滤波器的频率响应非常接近。当滤波器输入6Hz和17Hz的以50Hz采样频率采样的信号后，输入信号的归一化频率为6/(50/2)=0.24和17/(50/2)=0.68，由图6-17可见，0.24和0.68均在通带范围内，因此这两个信号可以无阻碍地通过滤波器，因此输入信号和输出信号应该相同。图6-18就验证了这一点。 6.7 IIR数字滤波器在地震数据分析中的应用举例 6.7.1 低频干扰的去除 我们知道，地震台观测的数据中除了地震震源激发的地震波外，还包含了如地脉动、海浪干扰等低频数据，为地震波数据分析带来麻烦。我们以湟源台记录的2003年10月25日23时49分甘肃民乐Ms5.8级近震的一个余震的NS分向记录为例说明低频干扰去除的方法。 在原始波形图上，有一个低频的波浪式背景干扰波与这个高频地震波形叠加，故需要设计高通滤波器，把低频的波浪式背景干扰波滤掉。高通滤波器的阻带边界频率为0.2Hz,通带边界频率为0.5Hz,通带波纹为1dB,阻带衰减为30dB,采用Butterworth滤波器。源程序如下： %Appl6_1 load ml031025ns.txt %加载地震波形记录 dt=0.02; %中国数字地震台网的采样间隔为0.02s,采样频率为50Hz x=ml031025ns'; wp=0.5*2*dt; ws=0.2*2*dt; %根据采样频率将滤波器边界频率进行转换 Rp=1;Rs=30; %通带波纹和阻带衰减 Nn=128; %显示滤波器频率特性的数据长度 [N,Wn]=buttord(wp,ws,Rp,Rs); %求得滤波器的最小阶数和截止频率（归一化频率） [b,a]=butter(N,Wn,'high'); %设计Butterworth高通数字滤波器 figure(1) [H,f]=freqz(b,a,Nn,1/dt); %用Nn点求得频率特性 subplot(2,1,1),plot(f,20*log10(abs(H))); xlabel('频率/Hz');ylabel('振幅/dB');grid on; subplot(2,1,2),plot(f,180/pi*unwrap(angle(H))) xlabel('频率/Hz');ylabel('相位/^o');grid on; n=0:length(x)-1; t=n*dt; %时间序列 figure(2) subplot(2,1,1),plot(t,x); title('输入信号') %绘制输入信号 y=filter(b,a,x); %对输入信号进行滤波 subplot(2,1,2),plot(t,y),title('输出信号') %绘制输出信号 xlabel('时间/s') 程序运行结果见图6-19和图6-20。图6-19为所设计滤波器的频率响应，显示出只有小于0.2Hz的信号才能被完全滤除。图6-20为滤波前和滤波后的波形图，可见，完全去除了低频波浪式的干扰，使得地震波显现出来。 0-10振幅/dB-20-30-400510频率/Hz1520250-100相位/o-200-300-4000510频率/Hz152025 图6-19 设计高通滤波器的频率响应 输入信号3020100-10-200102030输出信号20100-10-20-300102030时间/s405060405060图6-20 运用高通滤波器滤波前后波形的比较 6.7.2 爆破干扰的去除 地震台观测的数据中除了地震震源激发的地震波外，还包含了如地方爆破等高频干扰数据。库尔勒台记录的2003年9月16日一远震UD分向的P波中叠加了一个爆破记录。爆破频率相对地震波的频率明显较高，为了去掉该爆破，需要设计低通滤波器，把高频爆破滤掉。由于远震波形数据的地震波频率成分较低。我们选择低于1.5Hz的频率成分即可满足要求,因此我们设计通带边界频率为1.5Hz,通带波纹为1dB,阻带边界频率为2.5Hz,阻带衰减为30dB,采用Butterworth滤波器。源程序如下： %Appl6_2 load kerl030916ud.txt Fs=50; %采样频率 wp=1.5*2/Fs;ws=2.5*2/Fs; %根据采样频率将滤波器边界频率进行转换 Rp=1;Rs=30; %通带波纹和阻带衰减 Nn=128; [N,Wn]=buttord(wp,ws,Rp,Rs); %求得滤波器的最小阶数和截止频率（归一化频率） [b,a]=butter(N,Wn); %设计Butterworth低通数字滤波器 figure(1) [H,f]=freqz(b,a,Nn,Fs); %用Nn点求得频率特性 subplot(2,1,1),plot(f,20*log10(abs(H))); xlabel('频率/Hz');ylabel('振幅/dB');grid on; subplot(2,1,2),plot(f,180/pi*unwrap(angle(H))) xlabel('频率/Hz');ylabel('相位/^o');grid on; dt=1/Fs; %采样间隔 x=kerl030916ud'; %输入信号 t=[0:(length(x)-1)]*dt; %时间序列 figure(2) subplot(2,1,1),plot(t,x); title('输入信号') %绘制输入信号 y=filter(b,a,x); %对输入信号进行滤波 subplot(2,1,2),plot(t,y),title('输出信号') %绘制输出信号 xlabel('时间/s') 程序运行结果见图6-21和图6-22。图6-21为所设计滤波器的频率响应，显示出只有小于1.5Hz的信号才能被完全通过。图6-22为滤波前和滤波后的波形图，在约130秒处有一个高频爆破，滤波后完全去除了高频爆破的干扰，使得地震波分析不会出现失误。 0-100振幅/dB-200-300-400-5000510频率/Hz0-200152025相位/o-400-600-8000510频率/Hz152025图6-21设计低通滤波器的频率响应 210-1-2x 104输入信号0x 104100200300输出信号400500600210-1-20100200300时间/s400500600图6-22 运用低通滤波器滤波前后波形的比较， 上图：原始地震波形，在约130秒处叠加了一个高频爆破；下图：滤波后的结果 6.7.3汽车干扰的排除 汕头台所用数字地震记录仪器为FBS-3A型宽频带地震计（频带范围为0.05Hz-20Hz）和EDAS-C24型数据采集器，系统的采样频率Fs=50Hz。这里以汕头台受到距台站300m处的汽车干扰的波形记录图的数字地震记录资料为例，说明该干扰的去除方法。 我们采用Fourier分析得到了受干扰地震波的频谱图，其程序如下： %Appl6_3 load iir_signaldata.txt; %调入数据文件 Xt=iir_signaldata; %得到原始信号序列 Fs=50; %采样率为50Hz subplot(2,1,1); %时域坐标方框图 N=length(Xt); n=0:N-1; t=n/Fs; plot(t,Xt); %画出时域中的原始波形图 xlabel('时间/s'),title('时间域'); %时域图的X轴坐标标识 ylabel('振幅'); %时域图的Y轴坐标标识 Xf=fft(Xt); %对原始波形数据进行快速Fourier变换 subplot(2,1,2); %频域坐标方框图 mag=abs(Xf); %得到信号序列Fourier变换后的振幅 semilogy([0:N-1]/(N/Fs),mag*2/N); %画出频域中的FFT波形图 xlabel('频率/Hz'),title('频率域'); %频域图的X轴坐标标识 ylabel('振幅'); %频域图的Y轴坐标标识及标题 xlim([0 Fs/2]); %频率轴只画出Nyquist频率之前的频率成分 该程序的输出结果见图6-23。可见，长周期地震信号中叠加了高频汽车干扰，严重影响数据的分析。这在下图频谱图中可以看出，高频干扰集中在10~15Hz。因此我们可以运用低通滤波器将高频干扰滤除。我们设计的通带边界频率为4Hz,波纹为1dB，阻带起始频率为6Hz,阻带衰减为25dB。 时间域4002000振幅-200-400-600-800050时间/s4100150频率域10102振幅10010-210-40510频率/Hz152025图6-23 上图：汕头台受300米汽车干扰得到的地震波形；下图：其Fourier振幅谱 将高频干扰滤除的程序如下： %Appl6_4 Fcp=4;Fcs=6; %设置边界频率为4Hz，6Hz load iir_signaldata.txt; %调入数据文件 Xt=iir_signaldata; %原始波形数据 Fs=50; %采样频率50Hz wp=Fcp*2/Fs;ws=Fcs*2/Fs; %将通带临界频率Fcp和阻带临界频率Fcs转换为归一化频率 Rp=1;Rs=25; %通带衰减和阻带衰减 Nn=128; %绘频谱图所用点数 [N,Wn]=buttord(wp,ws,Rp,Rs); %求滤波器的最小阶数和临界频率（归一化频率） [b,a]=butter(N,Wn); %求Butterworth数字滤波器传递函数的分子和分母多项式系数b， a figure(1); %图形(一) [H,f]=freqz(b,a,Nn,Fs); %滤波器特性图 subplot(2,1,1),plot(f,20*log10(abs(H))); xlabel('频率/Hz');ylabel('振幅/dB');grid on; subplot(2,1,2),plot(f,180/pi*unwrap(angle(H))) xlabel('频率/Hz');ylabel('相位/^o');grid on; n=0:length(Xt)-1; t=n/Fs; %转换为时间序列 figure(2); %图形(二) subplot(2,1,1),plot(t,Xt);title('输入信号'); %绘制输入信号 xlabel('时间/s');ylabel('振幅'); %坐标轴标识 Yt=filter(b,a,Xt); %对输入信号进行滤波 subplot(2,1,2),plot(t,Yt),title('输出信号');%绘制输出信号 xlabel('时间/s');ylabel('振幅'); %坐标轴标识 程序的输出结果为图6-24和图6-25。图6-24为设计的低通滤波器的频率特性。采用低于4Hz的信号可以通过，高于6Hz的信号不能通过的滤波器，通带波纹为1dB，阻带衰减为25dB。由图可见完全符合设计要求。图6-25为受干扰信号和滤波后信号的比较，可见完全去除了高频干扰。 0-100振幅/dB-200-300-4000510频率/Hz1520250-200-400-600-800-1000相位/o0510频率/Hz152025图6-24 低通滤波器的频率响应 上图：幅频特性；下图：相频特性 输入信号4002000振幅-200-400-600-800050时间/s输出信号2001000100150振幅-100-200-300-400-500050时间/s图6-25 滤波后信号与受干扰信号的比较 上图：受干扰信号；下图：滤波后信号 100150 第六章 习题 1. 用脉冲响应不变法和双线性变换法将下列模拟滤波器H(s)转变为数字滤波器H(z)。(1) 2,采样周期为T=2s。 ?s?1??s?3?s2?s?1s?22（3）Ha(s)?,采样周期为T=2s。（4）,采样周期为H(s)?a222s?3s?1?s?1??4Ha(s)?3,采样周期为T=0.5s。(2) Ha(s)?T=1s。 2. 用脉冲响应不变法和双线性变换法设计一个Butterworth低通数字滤波器，要求满足通 带：0???0.2613?，通带波纹小于2dB；阻带为：0.4018?????，阻带衰减大于30dB，采样周期为T=0.1s。 3. 用双线性变换法设计一个三阶Butterworth滤波器:(1)高通，截止频率为Fc=1.5kHz， 采样频率Fs=6kHz；（2）带通，上截止频率f1=100Hz,下截止频率f2=300Hz，采样频率Fs=1KHz。其中通带波纹为3dB，阻带衰减为20dB。 4. 设计一个低通Chebyshev I型滤波器，通带边界频率为1500Hz,通带波纹小于3dB，阻 带边界频率2000Hz，阻带衰减为40dB，采样频率为8000Hz，试绘出其频率特性图并编制一测试该滤波器满足性能的例子。 5. 设计一个带通Chebyshev II型滤波器，通带边界频率为0.4?~0.5?，阻带边界频率 为0.3?和0.6?，通带波纹0.5dB，阻带衰减50dB。试绘制所设计滤波器的脉冲响应和幅频响应图。设计一个检测滤波器性能的例子。绘制其零极点图和群延迟。 6. 用直接法设计一个多频带数字IIR滤波器，其理想幅频响应为：频率[0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1] ?,幅值[0 0 1 1 1 0 0 1 1 0 0]，绘制其理想幅频图并与实际滤波器幅频图相比较。设计一个例子说明此滤波器满足要求。绘制其群延迟图。