数字信号处理（MATLAB版）课件 第10章 数字滤波器的优化设计和工具设计法

第10章数字滤波器的优化设计和工具设计法10.1数字滤波器的优化设计方法MATLAB提供了强大的数字滤波器优化设计功能。即预先确定一种最佳设计准则，然后直接求得在该准则下滤波器系统的传递函数[b,a]。这种最优设计方法可方便地用于任意幅频特性要求的多带通复杂滤波器系统的设计。10.1.1数字滤波器的优化设计方法因为数字滤波器在很多场合所要完成的任务与模拟滤波器相同，如作低通、高通、带通及带阻网络等，这时数字滤波也可看作是“模仿”模拟滤波。由于模拟的网络综合理论已经发展得很成熟，产生了许多高效率的设计方法。很多常用的模拟滤波器不仅有简单而严格的设计公式，而且设计参数已表格化，设计起来方便、准确，因此可将这些理论继承下来，作为设计数字滤波器的工具。所以第一种方法，即传统设计法用得较为普遍。但随着计算机技术的发展，最优化设计方法的使用也逐渐增多。最优化设计方法分两步完成：（1）确定一种最优准则，如最小均方误差准则等，使设计出的实际频率响应的幅度特性与所要求的理想频率响应（2）在此最佳准则下，通过迭代运算求滤波器的系数ai、bi。最优设计方法可方便地用于任意幅频特性要求的多带通复杂滤波器系统的设计，可以大大减小滤波器的阶数，从而减小滤波器的体积，并最终降低了滤波器的成本。并可以充分利用MATLAB、FFT等工具进行复杂滤波器的设计，提高设计效率。的均方误差最小，此外还有其他多种误差最小准则。10.1.2最小均方误差准则若以表示逼近误差，则：那么均方误差为：

（10.1.1）均方误差最小准则就是选择一组时域采样值，以使均方误差，这一方法注重的是在整个频率区间内总误差的全局最小，但不能保证局部频率点的性能，有些频率点可能会有较大的误差，对于窗口法FIR滤波器设计，因采用有限项的逼近理想的所以其逼近误差为：

（10.1.2）矩形窗窗口设计法是一个最小均方误差FIR设计，其优点是过渡带较窄，缺点是局部点误差大，或者说误差分布不均匀。10.1.3最大误差最小化准则FIR滤波器的优化设计是按照最大误差最小化准则（最佳一致逼近准则），使所设计的频响与理想频响之间的最大误差，在通带和阻带范围均为最小，而且是等波动逼近的，因此也叫切比雪夫等波纹逼近，以此准则的设计方法又称为切比雪夫最佳一致逼近法。表示为：

其中F是根据要求预先给定的一个频率取值范围，可以是通带，也可以是阻带或两者。在给定频带范围内使频响的最大逼近误差达到最小。优点是可保证局部频率点的性能也是最优的，误差分布均匀，相同指标下，可用最少的阶数达到最佳化。FIR、IIR滤波器的优化设计可使用

雷米兹（Remez）交替算法和Chebyshev逼近理论（Parks-McClellan算法）、Yule-Walker算法、帕德（Pade）逼近法、普罗尼（Prony）算法，等等，每一种算法都有自己的特点和适用范围。

（10.1.3）最佳一致逼近即选择N个频率采样值或时域值，10.2FIR滤波器的优化设计FIR滤波器的优化设计一般使用“最佳一致逼近准则”，有代表性的是“切比雪夫最佳一致逼近”，即“等波纹逼近法”。10.2.1

等波纹逼近法与雷米兹(Remez)交替算法为了简化起见，在优化设计中一般将线性相位

FIR

滤波器的单位脉冲响应h(n)

的对称中心置于

n=0

处，此时，线性相位因子α

=0。当

N

为奇数，且

N=2M+1，则：

（10.2.1）321（10.2.2）就是说，要在通带0

p

内以最大误差

1

逼近

1，在阻带内以最大误差

2

逼近

0。如希望逼近一个低通滤波器，这里M、和固定为某个值。在这种情况下有：等波动逼近的低通滤波器幅频特性，如图

10-2-1

所示。图

10-2-1

等波动逼近的低通滤波器幅频特性p r322优化设计滤波器需要确定五个参数：

M

、

r

、

c

、

1

、

2。其中，

1

为通带波动，

1

为阻带波动。根据式（7.1.15）推导过程，可知：Rp要同时确定上述五个参数较困难，常用的两种逼近方法有：（1）给定M

、

1

、

2

，以

c

和

r

为变量，这是“非线性最优法”。缺点：边界频率不能精确确定。

（2）给定

M

、

r

、

c

，以

1

、

2

为变量，通过迭代运算，使逼近误差

1

和

2

最小，并确定h(n)

。这种方法是所谓“切比雪夫最佳一致逼近（等波纹逼近法）”。特点：能准确地指定通带和阻带边界频率。如果要求通带和阻带具有不同的逼近精度，就要对误差函数进行加权。定义逼近误差函数：E(

)

为在希望的滤波器通带和阻带内算出的误差值，W

(e

j

)

为加权函数，在不同频带可取不同的值：k

应当等于比值

1

/

2

，在这种情况下，设计过程要求|

E(

)

|在区间0

p

和

r

的最大值为最小，它等效于求最小

2。根据数学上多项式逼近连续函数的理论，用三角多项式逼近连续函数，在一定条件下存在最佳逼近的三角多项式，而且可以证明这个多项式是唯一的。这一最佳逼近定理通常称作“交替定理”。在逼近过程中，可以固定

k、M

、

r

、

c

，而允许改变

2

，

1

k

2

，按照交替定理，首先估计出(M+2)个误差函数的极值频率点：Mii dj

n

1式中ρ表示峰值误差。一般仅需求解出ρ，接着便可用三角多项式找到一组新的极值频率点，并求出新的峰值误差ρ。依此反复进行，直到前、后两次ρ值不变化为止，最小的ρ即为所求的

2

。这一算法通常称作“雷米兹(Remez)交替算法”。32510.2.2使用

firpm()和

firpmord()函数优化设计FIR

滤波器前面介绍了

FIR

数字滤波器的两种逼近设计方法，即窗口法（时域逼近法）和频率采样法（频域逼近法），用这两种方法设计出的滤波器的频率特性都是在不同意义上对给定理想频率特性H

d

(e

)

的逼近。窗口法和频率采样法都是先给出逼近方法、所需变量，然后再讨论其逼近特性。如果反过来要求在某种准则下设计滤波器各参数，以获取最优的结果，这就引出了最优化设计的概念，最优化设计需要大量的计算，所以一般需要依靠计算机进行辅助设计。最优化设计的前提是最优准则的确定，在

FIR滤波器最优化设计中，常用的准则有：最小均方误差准则。最大误差最小化准则。firpm()

和

firpmord()

函数，

在以前版本是

remezord()

和

remez()

函数，

使用上述准则和Parks-McClellan

算法进行

FIR

滤波器的优化设计。配合使用

firpm()和

firpmord()函数，可以方便地优化设计FIR

滤波器。1.firpm()函数Parks-McClellan算法实际是使用雷米兹（Remez）交替算法和

Chebyshev

逼近理论，设计一个具有理想与实际频率响应之间最佳匹配的滤波器。firpm()函数可设计

1、2、3

和

4

型

FIR

线性相位滤波器，是所需的频率响应与实际频率响应之间的“最大误差最小化”意义上的优化滤波器。这种滤波器在频率响应中表现为等波纹，有时也称为等波纹滤波器。firpm()函数的使用方法如下：326（1）b=firpm(n,f,a)：返回

n

阶

FIR

滤波器的系数向量

b，行向量

b

具有

n+1

个元素。f

向量是归一化频率，范围是

0～1，1

对应于

Nyquist

频率。a

是由向量

f

指定的频率点对之间所期望的幅度，如

f(k)、f(k+1)。向量

f、a

的长度必须相同，而且是偶数。（2）b=firpm(n,f,a,w)：在频带内使用向量w

加权，长度是f、a长度的一半。（3）b

=firpm(n,f,a,'ftype')、b

=

firpm(n,f,a,w,'ftype')：ftype

取值为：'hilbert'，用于奇对称线性相位滤波器（3型和

4

型），系数向量

b

为b

k

b

n

2

k

k

1,

,n

1

。这类滤波器包含

Hilbert

转移函数，它有一个整个频带幅度为

1

的期望值。例如，h

=firpm(30,0.1

0.9],1

1],hilbert')，设计一个长度为

31的

FIR

Hilbert

逼近转移函数。'differentiator'，用于

3

类和

4

类滤波器，在非零的幅度带中使用一种特殊的权重技术：误差的权因子是1

/

f

，使低频段的误差远小于高频段。FIR

与众不同之处，有一个振幅与频率成正比的特点，这种滤波器满足最大误差最小化准则。（4）b=firpm(...,lgrid})：使用

lgrid（lgrid

为整数）控制频率网格的密度。有(lgrid*n)/(2*bw)

个频点，bw

是总的频率带间隔[0，1]

的分数，lgrid

默认值是

16，提高该值可以更精确地接近于等波纹滤波器，但需要占用更长的计算时间。{lgrid}

参数必须是一个

1*1

的单元数组。（5）[b,delta]

=

firpm(...)：delta是最大的波纹幅度。【例

10-2-1】优化设计一个

17阶带通滤波器其参数为n=17，f=[0

0.3

0.4

0.6

0.7

1]，

a

=

[0

0

1

1

0

0]。解：程序代码如下：f=[00.30.40.60.71];a=[00110

0];b=firpm(17,f,a);[h,w]=

freqz(b,1,512);plot(f,a,w/pi,abs(h));legend('理想滤波器','firpm

设计的滤波器')firpm()函数设计的带通滤波器，如图

10-2-2

所示。2.firpmord()函数firpmord()函数用于估算Parks-McClellan

优化

FIR

滤波器阶次。用法如下：（1）[n,fo,ao,w]

=

firpmord(f,a,dev)：返回最接近的阶次

n、归一化频率带边沿

fo、频带幅度

ao

和权重

w，输入参数为

f、a

和dev。f

是一个频率向量，位于[0，fs/2]，fs是抽样频率。a

是由向量

f

指定的频率点对之间所期望的幅度向量。f

的长度是

a

长度的

2

倍减

2

：length(f)=2*length(a)-2，length(a)=length(f)/2+1。所得到的函数是分段常数。dev

是一个长度与

a

长度相同的向量，指定频率响应与理想幅度之间允许的最大偏差（即：[

p327

]），阻带误差：通带误差：

（2）[n,fo,ao,w]

=firpmord(f,a,dev,fs)：指定抽样频率

fs，默认是

2Hz。【例

10-2-2】利用雷米兹交替算法，设计一个低通

FIR

数字滤波器利用雷米兹交替算法，设计一个线性相位低通

FIR

数字滤波器，其指标为：通带边界频率

fc

800Hz

，阻带边界

f

1000Hz

，通带波动R

0.5dB

，阻带最小衰减R

40dB

，采样频率

f

4000Hz

。r p r s图

10-2-2

firpm()函数设计的带通滤波器p r解：已知，通带最大波动R

0.5dB

，阻带最小衰减R

40dB

，得即dlt1=0.0559、dlt2=0.01，数组dev=[dlt1dlt2]是通带和阻带的波动，fs=4000Hz

是采样频率。fc

800Hz

、fr

1000Hz

，数组

f=[fc

fr]为通带和阻带边界频率，数组

a是两个边界处的幅值：mp=1、mr=0,a=[mp

mr]。在MATLAB

中可以用firpmord()计算出阶次等参数，然后使用

firpm()函数设计，程序如下：grid;title('雷米兹(Remez)交替算法')328fc=800;fr=1000;mp=1;mr=0;dlt1=0.0559;dlt2=0.01;dev=[dlt1

dlt2]; f=[fcfr];a=[mpmr];fs=4000;[N,f0,a0,w]=firpmord(f,a,dev,fs);b=firpm(N,f0,a0,w);[h,w]=freqz(b,1,256);plot(w*2000/pi,20*log10(abs(h)));xlabel('频率/Hz');

ylabel('幅度/dB');雷米兹(Remez)交替算法设计的低通滤波器，如图

10-2-3

所示。同样是设计一个

FIR

低通数字滤波器，综合分析可以看出：

（1）窗函数法在阶数较低时，阻带特性不满足设计要求，只有当滤波器阶数较高时，使用海明窗和凯塞窗基本可以达到阻带衰耗要求；

（2）频率采样法偏离设计指标最明显，阻带衰减最小，而且设计比采用窗函数法复杂。只有适当选取过渡带样点值，才会取得较好的衰耗特性；（3）利用等波纹切比雪夫逼近法则的优化设计可以获得最佳的频率特性和衰耗特性，具有通带和阻带平坦，过渡带窄等优点。330图10-2-3雷米兹(Remez)交替算法设计的低通滤波器使用频率抽样法的fir2()函数和等波纹法的firpm()函数，可以很方便的设计多带滤波器。10.2.3多带滤波器的设计【例10-2-3】设计多带滤波器分别用频率抽样法和等波纹法设计出62阶FIR四通带数字滤波器，并与理想滤波器比较，如图10-2-4所示。

采用最优化设计方法时大大减小了滤波器的阶数，从而减小了滤波器的体积，并最终降低了滤波器的成本。这样使得设计出来的滤波器更为简单经济。在实际应用中，如果需要对某一信号源进行特定的滤波，并要检验滤波效果，应用传统方法实施起来比较繁琐。在Matlab环境下，可先用软件模拟产生信号源，再设计滤波器对其进行滤波。图10-2-4

62阶FIR四通带数字滤波器332根据功率频率响应的加法分解式计算相应的分子多项式。根据分子和分母多项式计算完整的频率响应。使用谱分解技术获得滤波器的脉冲响应。使用最小二乘法拟合该脉冲响应获得分子多项式。10.3.1使用

yulewalk()函数优化设计

IIR

滤波器MATLAB

提供

yulewalk()函数直接设计递归

IIR

滤波器，这个函数在离散时域使用最小二乘法拟合逼近给定的频率特性，近似一种分段线性幅度响应。它使用改进的

Yule-Walker

公式计算分母系数，由给定的频率响应的傅里叶逆变换计算相关的系数。yule-walk按以下步骤计算分子系数：10.3 IIR

滤波器的优化设计Yulewalk()函数用法如下：[b,a]=yulewalk(n,f,m)：该函数返回一个

Yule-Walk

滤波器的系数矩阵[b，a]，包含

n+1

个系数的行向量b和

a；其中

n是

IIR

滤波器的阶数，f

和

m是已知的频率响应，n

阶

IIR

滤波器的幅频特征与给出的向量

f

和

m

匹配：f

是频率点向量，f

元素的值必须在

0

和

1

之间，允许出现相同的频率值。而且必须是升序，以0开始、以

1

结束。1对应于半采样频率（奈奎斯特频率）。m是一个向量，代表幅度值，每点代表在

f

点指定的幅度响应。f

和

m

必须具有相同的长度。plot(f,m)显示该滤波器的形状。所得到的滤波器按

z的降序排列，系统函数可写成：当指定频率响应时，为避免通带到阻带过渡的转折过于锐利，可能需要多次检验过渡区，以得到最佳设计的滤波器边沿。【例

10-3-1】优化法设计

IIR

低通数字滤波器设计一个

10阶低通滤波器，并绘制实际频率响应和所要求的理想频率响应。f=[00.60.61];m=[110

0];解：理想的响应是频率在

0～0.6

之间，幅值为

1；在频率在

0.6～1

之间的幅值为

0。设计程序为：n=10;f=[00.60.61];m=[1 1 00

];[b,a]=yulewalk(n,f,m);[H，w]=freqz(b,a,

128)plot(f,m,

w/pi,

abs(H),'--')legend('理想'

,

'yulewalk

设计')title('理想与

yulewalk

设计比较')如图

10-3-1

所示，是用

yulewalk()函数设计的

Yule-Walk

滤波器幅频响应与理想的幅频响应的比较图。图10-3-1yulewalk滤波器和理想低通滤波器频响对比10.3.2

普罗尼（Prony）算法333根据式（10.1.2）可知，均方误差最小准则就是选择一组时域采样值

a(i)、b(i)，由于

是

a(i)、b(i)的非线性函数，求解这个最小化问题比较困难。普罗尼（Prony）算法通过两步过程求一个近似的最小二乘解。MALAB提供一个

Prony

函数完成普罗尼（Prony）算法：[A,B]

=

prony(hn,M,N)

返回因果系统函数：最小，其中U是预先设定的上限。使均方误差系数多项式的分子向量

A

和分母向量

B，已知单位脉冲响应hd

(n)

用hn

表示。系统函数的分子

A

的长度为

M+1，分母

B

的长度是N+1。

M、N

是分子和分母多项式的阶数，如果hn

的长度小于

M、N

中最大者，用

0补齐。M

为

0

时是全极点系统函数，N

为

0时是全零点系统。【例

10-3-2】普罗尼（Prony）算法设计滤波器已知单位脉冲响应hd

(n)

用hn

表示

hd=[0.9

0.80.6

0.3-0.100.1

0.01

0.003

-0.001]，Prony

函数设计一个

IIR

滤波器，逼近其响应。解：

h(n)

的长度

Nh=10，选

N=6、M=4，M+N=10。hd=[0.90.80.6

0.3-0.1 00.10.010.003

-0.001];N=6;M=4;[A,B]=prony(hd,M,N);hn=impz(A,B,20);tf(A,B,0.01)[H,w]=freqz(hn,1);plot(w/pi,20*log10(abs(H/H(1))));

grid;xlabel('归一化频率

w/

\pi')

;

ylabel('幅度(dB)')

;画出用

Prony

函数设计的滤波器频响曲线，如图

10-3-2

所示。系统函数为：图

10-3-2

用

Prony

函数设计的滤波器频响曲线在

MATLAB

工具箱“Toolbox”中，在滤波器设计“Filter

Design”和信号处理工具箱“Signal

Processing”目录下，提供了一种有用的数字滤波器的设计工具

FDATool（FilterDesign

&

Analysis

Tool），可以快速有效的设计数字滤波器。10.4.1FDATool

简介计算机辅助设计方法是集电路理论、网络图论、数值分析、矩阵运算、元件建模、优化技术、高级计算机语言等多交叉学科于一身的新领域，它把计算机的快速、高精度、大存储容量、严格的逻辑判断和优良的数据处理能力与人的思维创造能力充分结合起来，极大地简化了数字滤波器的设计过程。滤波器设计方法有程序设计法和工具设计法。FDATool

即滤波器设计和分析工具，是

MATLAB

信号处理工具箱里专用的滤波器设计分析工具。DATool

有一个强大的用户界面用于设计和分析滤波器，FDATool

几乎可以设计低通、高通、带通和带阻等所有的常规滤波器，包括

FIR

和

IIR

的各种设计方法。利用

FDATool

设计滤波器，可以随时对比设计要求和滤波器特性调整参数、操作简单、直观方便灵活，极大的减轻了工作量，有利于滤波器设计的最优化。FDATool

可以从

MATLAB

工作空间输入，也可以直接指定滤波器性能系数。FDATool

还提供了滤波器分析工具，如幅度、相位响应和零极点图等。33510.4

使用

FDATool

设计数字滤波器在

MATLAB

命令行窗口输入命令：>>

fdatool如图

10-4-1

所示，FDATool

界面共分两大部分，一部分是“Design

Filter”，在界面的下半部，用来设置滤波器的设计参数；另一部分则是特性区，在界面的上半部分，用来显示滤波器的各种特性。“Design

Filter”部分主要分为

6

部分：Response

Type、Design

Method、Filter

Order、Options、Frenquency

Specifications

和

MagnitudeSpecifications。滤波器响应类型Response

Type

“（滤波器响应类型）”选项，包括：Lowpass（低通）、Highpass（高通）、Bandpass（带通）、Bandstop（带阻）和特殊的

FIR

滤波器。图

10-4-1 FDATool界面滤波器的设计方法DesignMethod（设计方法）选项，包括：IIR

滤波器的

Butterworth（巴特沃思）法；Chebyshev

Type

I（切比雪夫

I

型）法；Chebyshev

Type

II（切比雪夫

II

型）法；Elliptic（椭圆滤波器）法；FIR

滤波器的

Equiripple法；Least-Squares（最小乘方）法；Window（窗函数）法。3.滤波器规格特性（1）FilterOrder（滤波器阶数）选项。（2）Options.可用的选项“Options”取决于选定的滤波器设计方法。只有“FIREquiripple”和“FIRWindow（窗函数设计方法）”有可设置的选项。（3）“FrenquencySpecifications（频率特性）”选项，可以详细定义频带的各参数，包括采样频率Fs和频带的截止频率。（4）“MagnitudeSpecifications（幅度特性）”选项，可以定义幅值衰减的情况。4.结构转换在FDATool中，可以实现滤波器的结构形式转换.10.4.2带通滤波器的设计以带通滤波器的设计为例，介绍FDATool的使用方法，其它滤波器的设计方法与此类似。【例10-4-1】设计一带通数字滤波器参数要求：96阶带通FIR数字滤波器，通带下限截止频率Fc1＝60Hz和通带上限截止频率Fc2＝124Hz，采样频率1000Hz，采用Hamming窗函数设计。设计方法和步骤如下：1.设置参数（1）首先在“滤波器类型（FilterType）”选择框中选择“带通滤波器（Bandpass）”。（2）在“设计方法（DesignMethod）”选项中选择“FIR滤波器窗函数法（FIRWindow）”，然后在“WindowSpecifications”选项中选取“Hamming”。（3）指定“滤波器阶次（FilterOrder）”项中的SpecifyOrder＝95；采样频率1000Hz。（4）由于采用窗函数法设计，只要给出通带下限截止频率Fc1和通带上限截止频率Fc2，选取Fc1＝60Hz，Fc2＝124Hz。2.滤波器系数计算与结果的保存（1）设置完以后点击“DesignFilter”按钮，即可得到所设计的FIR滤波器。（2）设计完成后将结果保存为Mat文件，扩展名为“.fda”。l在“File”菜单选择“SavesessionAs”，打开“SaveFilterDesignFile”对话框。l在文件名文本框中输入要保存的名称，例如“filter1.fda”。3.设计结果观察与调整通过菜单项“Analysis”可以在特性区看到所设计滤波器的幅频响应、相频响应、零极点配置和滤波器系数等各种特性，如图10-4-4所示。在设计过程中，可以对比滤波器幅频相频特性和设计要求，随时调整参数和滤波器类型，以便得到最佳效果。其它类型的FIR滤波器和IIR滤波器也都可以使用FDATool来设计。设计完成后可将输出结果保存为下列文件：1.向工作空间输出为滤波器系数（Coefficients）或对象（Objects）可以将自己设计的滤波器保存为滤波