冲激响应滤波器课程设计

IIR 数字滤波器设计（Ⅰ）

数字滤波器简介

数字滤波器是对数字信号处理实现滤波的线性时不变系统。数字滤波实质上是一种运算过程，实现对信号的运算处理。输入数字信号(数字序列) 通过特定的运算转变为输出的数字序列，因此，数字滤波器本质上是一个完成特定运算的数字计算过程，也可以理解为是一台计算机。描述离散系统输出与输入关系的卷积和差分方程只是给数字信号滤波器提供运算规则，使其按照这个规则完成对输入数据的处理。时域离散系统的频域特性：Y (e jw ) =X (e jw ) H (e jw ) ， 其 中

Y (e

jw

) , X (e

jw

) 分别是数字滤波器的输出序列和输入序列的频域特性(或称为频

谱特性) ，H (e jw ) 是数字滤波器的单位取样响应的频谱，又称为数字滤波器的频域响应。可以看出， 输入序列的频谱X (e jw ) 经过滤波后变为X (e jw ) H (e jw ) ， 因此，只要按照输入信号频谱的特点和处理信号的目的， 适当选择H (e jw ) ， 使得滤波后的X (e jw ) H (e jw ) 满足设计的要求。

IIR 滤波器简介

数字滤波器在数字信号处理的各种应用中发挥着十分重要的作用。它是通过对采样数据信号进行数学运算处理来达到滤波的目的。其中无限冲击响应数字滤波器也称IIR 是采用对离散采样数据作差分方程运算来进行滤波。IIR DF的优点在于可以利用模拟滤波器设计的结果，然后采用双线性变换法或冲激响应不变法将模拟滤波器转换成数字滤波器，而模拟滤波器的设计方便简单、有大量的图表可查。IIR DF的设计还可利用二阶子系统的串联来有效实现。其主要缺点在于Nyquist 间隔内无法实现线性相位, 若需线性相位，比如在图象处理和数据传输中都要求信道具有线性相位特性，则需要用全通网络进行相位校正。

IIR 滤波器的传递函数模型

IIR 数字滤波器的基本网络结构可用差分方程、单位脉冲响应及系统函数来描述，其差分方程一般表示为

M

N

i

i

y (n ) =

∑b x (n -i ) -∑a y (n -i ) 式(2.1)

i =0

i =1

式中x (n ) 和y (n ) 分别表示输入和输出信号序列，a i 和b i 是滤波器系数。设式(2.1)中输入信号x (n ) 和输出信号y (n ) 在n =0以前处于零起始状态，则有

M M

-i i

H (z ) =

Y (z ) X (z )

∑b z

=

i =0

N i =1

∑b z

i

-i

=

-i

i =1N

式(2.2)

-i

1＋∑a i z ∑a z

i

i =1

对于线性时不变(LTI)系统来说，上式中a i 和b i 均为常数。

[2]

IIR 滤波器的状态方程模型

状态方程是描述系统的一种常用的方式。在系统有不可见的状态变量时可使用滤波器状态方程模型。连续LTI 系统的状态方程可以写成

x (t ) =A x (t ) +B u (t ) 式(2.3)

.

y (t ) =Cx (t ) +Du (t ) 式(2.4)

其中A ，B ，C ，D 分别为常数矩阵，在MA TLAB 中，一般情况下，系统的状态方程可以简记为(A ，B ，C ，D ), 如果D ≡0，则系统的状态方程模型可以简记为(A ，B ，C ) 。

对于离散系统来说，状态方程模型可以写成

x [(i +1) T ]=Ax (iT ) +Bu (iT ) 式(2.5)

y [(i +1) T ]=Cx (i +1) T +Du (i +1) T 式(2.6)

这样的离散时间系统的状态方程模型在MATLAB 中写成(A ，B ，C ，D ) 或(A ，B ，C )

[3]

。

IIR 滤波器的零极点模型

零极点模型实际是传递函数的另一种表达，它的原理是分别对原系统传递函数的分子和分母进行分解因式，以获得系统的零极点表示形式，对单输入单输出系统来说，零极点模型可以简单的表示为：

m

∑(s +z )

i

G (s ) =K

i =1

n

=K

i

(s +z 1)(s +z 2) (s +z m ) (s +p 1)(s +p 2) (s +p n )

式(2.7)

∑(s +p )

i =1

式中z i ，i =1,2,3…，m 和p i ，i =1,2,3…，分别称为系统的零点和极点，它们既可以是实数也可以是复数，而K 称为系统的增益。在MATLAB 中简记为(Z,P,K)。但为了保证滤波器的稳定性，系统的所有极点都位于左半S 平面。如果稳定系统所有的零点都位于左半S 平面，则称该系统为最小相位系统，否则称为非最小相位系统。如果系统中的某个零点值恰好等于其中一个极点的值，则它们

之间可以对消以直接获得一个完全等效的低价系统。

IIR 数字滤波器的实现方法

冲激响应不变法

冲激响应不变法的变换原理是使数字滤波器的单位冲激响应序列h (n ) 模仿模拟滤波器的单位冲激响应序列h a (t ) 。将模拟滤波器的单位冲激响应加以等简隔抽样，使h (n ) 正好等于h a (t ) 的抽样值，即满足

h (n ) =h a (nT )

[5]

式(3.1)

其中T 是抽样周期。

如果令H a (s ) 是h a (t ) 的拉普拉斯变换，H (z ) 为h (n ) 的Z 变换，根据抽样序列的Z 变换和模拟信号的拉普拉斯变换的关系，可得

H (z ) |z =e sT =

1T

∞

∑

k =-∞

H a (s -j

2πT

k ) 式(3.2)

可以看出，利用冲激响应不变法将模拟滤波器变换数字滤波器，实际上是先对模拟滤波器的系统函数H a (s ) 做周期延拓，再经过z =e sT 的映射变换，从而得到数字滤波器的系统函数H (z ) 。假设在S 平面上，S 在j Ω轴上取值，Z 在Z 平面内的单位圆周e j ω上取值，则可得到数字滤波器的频率响应H (e j ω) 和模拟滤波器的频率响应H (j Ω) 间的对应关系：

H (e

j ω

) =

1T

∞

∑

k =-∞

H a (j

ω-2πk

T

)

式(3.3)

如图3.1所示，在S 平面上每一条宽为2π/T 的横条都将重叠地映射到整个Z 平面上，而每一条的左半边映射到Z 平面单位圆内，右半边映射到Z 平面单位圆外，S 平面虚轴映射Z 平面单位圆上，虚轴上每一段长为2π/T 的线段都映射到Z 平面单位圆上一周。由于S 平面每一横条都要重叠地映射到Z 平面上，这正好反映H (z ) 和H a (s ) 的周期延拓函数的变换关系z =e sT ，故有冲激响应不变法并不是简单的从S 平面映射到Z 平面。

图3.1 冲激响应不变法的映射关系

因为不是简单的一一映射关系，且对于任何一个实际的模拟滤波器，它的频率响应是不可能真正带限的。因而将不可避免的出现频率的重叠，即混叠失真。数字滤波器的频率响应不能重现模拟滤波器的频率响应，只有当模拟滤波器的频率响应在超过重叠频率后的衰减很大时，混叠失真才会很小，此时才能够满足设计要求。

综上所述，冲激响应不变法具有以下特点：

(1) 模拟频率和数字频率的转换是线性的，并保持了模拟滤波器的时域瞬态特性。

(2) 当模拟滤波器的频率响应不是严格限带时，则用冲激不变法设计出的数字滤波器在频率出现混叠失真。

(3) 由于(2)而使这种设计方法受到限制，即当G (j ) 不严格限带或g(t)变化不稳定。

巴特沃斯低通滤波器的设计 巴特沃斯滤波器的特性

巴特沃斯滤波器的特点是同频带内的频率响应曲线最大限度平坦，没有起伏，而在阻频带则逐渐下降为零。 在振幅的对数对角频率的波特图上,

从某一边界角频率开始, 振幅随着角频率的增加而逐步减少, 趋向负无穷大。

一阶巴特沃斯滤波器的衰减率为每倍频6分贝，每十倍频20分贝。二阶巴特沃斯滤波器的衰减率为每倍频12分贝、 三阶巴特沃斯滤波器的衰减率为每倍频18分贝、如此类推。巴特沃斯滤波器的振幅对角频率单调下降，并且也是唯一的无论阶数，振幅对角频率曲线都保持同样的形状的滤波器。只不过滤波器阶数越高，在阻频带振幅衰减速度越快。其他滤波器高阶的振幅对角频率图和低级数的振幅对角频率有不同的形状。

模拟低通巴特沃斯滤波器以巴特沃斯函数为滤波器的系统函数，它的幅度平方函数表示为：

|H a (j Ω) |=

1+ε(

2

2

1j Ωj Ωc

)

2N

式(3.10)

式中的N 为正整数，表示滤波器的阶数。Ωc 为通带的截止频率，因为当

Ω=Ωc 时：

|H a (j Ω) |=

2

12

(归一化后的巴特沃斯滤波器，ε

=1)

|H a (j 0) ||H a (j Ωc ) |

|H a (j Ω) |=

1，进而有20log 10(

) =3dB

，可3

分贝。巴特沃斯滤波器在通带有最为平坦的幅度特性，即N 阶低通滤波器在Ω=0处幅度平方函数|H a (j Ω) |2的前(2N-1)阶为零。阻带内随着频率的升高而单调下降。滤波器特征完全由阶数N 决定，当N 增加时，通带更加平坦，也更接近理想的低通滤波器的特性。

低通巴特沃斯滤波器的设计过程包括以下两个方面[7]：

(1) 根据给定的通带和阻带指标确定阶数N 。通常在设计时会给定技术指标，包括通带的截止频率Ωc 、通带内最大衰减A P (dB ) 、阻带截止频率Ωs 、阻带内的最小衰减A s (dB ) 。假设给定的，Ω=Ωc 时通带的最大衰减为A p ，Ω=Ωs 时阻带的最小衰减为A S 。其幅度平方函数图如图3.5所示。这时：

通带容限为：|H a (j Ω) |>

11+ε

2

(|Ω|≤Ωc ) ；

阻带容限为：|H a (j Ω) |

11+λ

2

(|Ω|≥Ωc ) 。

其中ε, λ均为与衰减有关的参数。根据式(3.10)可得到：

图3.5 巴特沃斯滤波器幅度平方函数

A p =-10log(

11+ε

2

) =10log(1+ε2

) A 2

s =-10log(

11+λ

2

) =10log(1+λ) 由式(3.11a)、(3.11b)两式可得：

ε=(10

0.1A p

-1)

0.5

，λ=(100.1A s

-1) 0.5 对于同一滤波器，其幅度平方函数相等，于是可得

ε2

(

j Ωs j Ω)

2N

=λ

2

c

进而可得N 为正整数的取值要求：

log λ10(N ≥

)

log j Ωs 10(

j Ω)

c

当给定的参数指标是归一化的，即Ωc =1, ε=1时，上式表示为：式(3.11a)

式(3.11b) 式(3.12)

式(3.13)

N ≥

log 10(λ) log 10(Ωs )

(2) 从幅度平方函数确定系统函数H a (s ) 。令s =j Ω代入 (2.8)式得到：

H a (s ) H a (-s ) =

1+ε(

2

1s j Ωc

)

2N

巴特沃斯滤波器的全部零点都在s =∞处，没有有限零点，只求解上式的分母多项式的根。当设计归一化的巴特沃斯滤波器时，它的极点均匀分布在单位圆上。得到H a (s ) 的极点为：

s pk =-sin(

2k -12N

2k -12N

1, 2…（, 3

）N +为1奇数/2

π) +j cos(

π) k ={1，，2…3

为 偶 数 N

式(3.14)

令归一化的系统函数用H aN ，则在左半平面上的极点组成H aN ，H aN （s ）（s ）（s ）

（s /Ωc ）和H a (s ) 的关系为H a (s ) =H aN 。

H a (s ) 的极点表达式形式可表示为：

H a (s ) =

1

(s -s 1)(s -s 2)(s -s 3) ⋯（s -s pk ）

式(3.15)

(3) 利用MA TLAB 设计巴特沃斯滤波器的实例。

一、设计内容

设计一个数字巴特沃斯低通滤波器，设计指标如下：

ωp =0. 2π R p =1dB ωs =0. 3π A s =15dB

采样时间间隔T =1S 。

二、设计要求

（1）用单位冲激响应不变变换法进行设计。

（2）给出详细的滤波器设计说明书。

（3）给出经过运行是正确的程序清单并加上详细的注释。 （4）画出所设计滤波器的幅度特性和相位特性。

三、设计说明

（1） 把数字频率转换为模拟频率：Ωp =ωp T ，Ωs =ωs T 。 （2） 计算巴特沃斯模拟滤波器的截止频率Ωc 和阶数N 。 （3） 设计巴特沃斯模拟低通滤波器，给出参数b 和a （此处使用了

MATLAB 中的buttap(N)函数。

（4） 把模拟滤波器用单位冲激响应不变变换法转换成数字滤波器

（此处使用了MATLAB 中的residuez 函数）。

（5） 变直接形式为并联形式，并给出结构图。 （6） 画出幅度特性和相位特性。

四、设计步骤框图

数字滤波器的设计步骤如图1所示

五、Matlab 程序代码

%主程序 wp=0.2*pi; ws=0.3*pi; Rp=1;%通带波动; As=15;%阻带衰减 T=1;%采样周期 Omigrp=wp*T; Omigrs=ws*T;

t1=[1,2*0.2588,0.2588^2+0.9659^2];

b=1;

a=conv(conv(t1,t1),t1);

[N,wc]=buttord(wp,ws,Rp,As,'s' ) %计算巴特沃斯数字滤波器的阶数N 和截止频率wc

[B,A]=butter(N,wc,'s' )

[z,p,k]=buttap(N)%模拟低通原型零、极点系数和增益因子

[bz,az]=impinvar(b,a)%用冲激响应不变法将模拟滤波器转化为数字滤波器, 采样频率默认1Hz wz=[0:pi/512:pi];

hz1=freqz(bz,az,wz);%巴特沃斯模拟低通滤波器频率响应 [C,B,A]=dir2par(b,a)%直接型转换成并联型 %绘图

subplot(1,2,1); plot(wz/pi,abs(hz1)/hz1(1)); title(' 幅度响应' )

xlabel('' ); ylabel('|H|');axis([0,1,0,1]);grid; subplot(1,2,2); plot(wz/pi,hz1/pi); title(' 相位响应' ) xlabel('' ); ylabel(' 单位:pi'); axis([0,1,-1,2]);grid;

%直接型转换成并联型子程序

function [C,B,A]=dir2par(b,a); M=length(b); N=length(a);

[r1,p1,C]=residuez(b,a); p=cplxpair(p1,10000000*eps);

I=cplxcomp(p1,p);

r=r1(I);

K=floor(N/2);B=zeros(K,2);A=zeros(K,3);

if K*2==N;

for i=1:2:(N-2)

Brow=r(i:1:(i+1),:);

Arow=p(i:1:(i+1),:);

[Brow,Arow]=residuez(Brow,Arow,[]);

B(fix((i+1)/2),:)=real(Brow);

A(fix((i+1)/2),:)=real(Arow);

end

[Brow,Arow]=residuez(r(N-1),p(N-1),[]);

B(K,:)=[real(Brow) 0];A(K,:)=[real(Arow) 0];

else

for i=1:2:(N-1)

Brow=r(i:1:(i+1),:);

Arow=p(i:1:(i+1),:);

[Brow,Arow]=residuez(Brow,Arow,[]);

B(fix((i+1)/2),:)=real(Brow);

A(fix((i+1)/2),:)=real(Arow);

end

end

%比较两个含同样标量元素但（可能）有不同下标的复数对及其相应留数向量子程序

function I=cplxcomp(p1,p2);

I=[];

for j=1:length(p2)

for i=1:length(p1)

if (abs(p1(i)-p2(j))

I=[I,i];

end

end

End

六、运行结果

Rp = 1

N = 6

wc = 0.7087

B = 0 0 0 0 0 0 0.1266

A = 1.0000 2.7380 3.7484 3.2533 1.8824

z = []

p = -0.2588 + 0.9659i

-0.2588 - 0.9659i

-0.7071 + 0.7071i

-0.7071 - 0.7071i

-0.9659 + 0.2588i

-0.9659 - 0.2588i

k = 1

bz = -0.0000 0.0060 0.6905 0.1266 0.1049 0.1946 0.0625

0.0021

az = 1.0000 -2.6340 4.1006 -3.8164 2.4437 -0.9355 0.2117

C = []

B = 0.4019 0.1523

0.3135 -0.0693

0.2846 0.0829

A = 1.0000 0.5176 0.9999

1.0000 0.5176 0.9999

1.0000 0.5176 0.9999

七、仿真波形