第3章
CHAPTER 3


数字滤波器设计





主要内容
 数字滤波系统的基本网络结构； 
 数字滤波器的基本概念与分类； 
 IIR型滤波器的设计； 
 FIR型滤波器的设计； 
 有限字长效应分析。
3.1数字滤波系统的基本网络结构
线性时不变系统一个最广泛的应用就是滤波。一般来说，滤波是指改变信号中各个频率分量的相对大小，或者抑制，甚至全部滤除某些频率分量的过程。完成滤波功能的系统称为滤波器，适当地选择或设计系统的频率响应，就可以实现各种不同要求的滤波功能。
3.1.1数字滤波系统的基本概念
一般地，一个线性时不变（LTI）离散系统可用如下差分方程来表示



y(n)=∑N-1i=0bix(n-i)+∑Mj=1ajy(n-j)(3.1.1)


则其系统函数H(z)为


H(z)=Y(z)X(z)=∑N-1i=0biz-i1-∑Mj=1ajz-j (3.1.2)


对于式（3.1.1）或式(3.1.2)表示的系统，如果aj=0，j=1,2,…,M，则


y(n)=∑N-1i=0bix(n-i) (3.1.3)


其系统函数为H(z)=∑N-1i=0biz-i，单位脉冲响应h(n)可表示为


h(n)=bn,0≤n≤N-1


0,其他n(3.1.4)


由于h(n)是有限长的序列，故称该系统为有限脉冲响应（FIR）系统。这种形式的滤波器称为FIR滤波器。FIR系统中不存在输出对输入的反馈支路，没有不为零的极点。若不是FIR系统，则由于系统中存在反馈支路，其单位脉冲响应h(n)为无限长序列，故称该系统为无限脉冲响应（IIR）系统，这种形式的滤波器称为IIR滤波器。这两类系统（滤波器）具有不同的特点，其网络结构也各不相同。




3.1.2IIR滤波系统的基本网络结构
IIR网络的特点是信号流图中含有反馈支路，即含有环路，其单位脉冲响应是无限长的。基本网络结构有三种，即直接型、级联型和并联型。其中直接型包括直接Ⅰ型和直接Ⅱ型，它们各自实现所获得的稳定性和计算误差性能各不相同。
1. 直接Ⅰ型
考虑如下N阶差分方程


y(n)=∑N-1i=0bix(n-i)+
∑Mj=1ajy(n-j) (3.1.5)


y(n)由两部分相加而成，其一部分


∑Mj=1ajy(n-j) 


是对y(n)依次延迟反馈M-1个单元的加权和。另一部分


∑N-1i=0bix（n-i）


是对x(n)依次延迟N个单元的加权和。两者都可用一个链式延迟结构来构成，两部分网络分别实现零点和极点，且一共需要N+M-1个延迟单元和相应的乘法器及一个加法器。直接Ⅰ型网络的优点是物理概念清晰，缺点是使用的延迟单元太多。一般使用如下的直接Ⅱ型。
2. 直接Ⅱ型
直接Ⅱ型又称之为典范型（canonic structure），将式(3.1.2)稍做变化，有


H(z)=∑N-1i=0biz-i1-∑Mj=1ajz-j=H1(z)·H2(z)(3.1.6)


式中


H1(z)=11-∑Mj=1ajz-j


其相应的差分方程为


w(n)=x(n)+∑Mj=1ajw(n-j) 


其中w(n)为中间序列。


H2(z)=∑N-1i=0biz-i


其对应的差分方程为


y(n)=∑N-1i=0biw(n-i)


它可由两个链式延迟结构级联而成，第一个实现系统函数的极点，第二个实现系统函数的零点。两行串行延时支路都对时间序列w(n)进行延迟，因此可予以合并，以节省一半的延迟单元。与直接Ⅰ型相比，除了节省了一半延迟单元外，这种结构参与反馈环路的噪声源减少了一半，可以得到较直接Ⅰ型略小一些的误差，但仍然没有从根本上克服其缺点。
下面通过例子来说明直接Ⅱ型网络结构。
［例3.1.1］设IIR数字滤波系统的系统函数H(z)为


H(z)=8-4z-1+11z-2-2z-31-54z-1+34z-2-18z-3


画出该滤波器的直接Ⅱ型结构。
解由H(z)写出差分方程如下


y(n)=54y(n-1)-34y(n-2)+18y(n-3)

+8x(n)-4x(n-1)+11x(n-2)-2x(n-3)


按照差分方程画出如图3.1.1所示的直接型网络结构。


图3.1.1例3.1.1图


3. 级联型
级联型是以系统函数H(z)经因式分解后的零点cr和极点dk为主要依据的数字滤波系统结构形式，用零点、极点表示的H(z)为


H(z)=A∏N-1i=1(1-ciz-1)∏Mj=11-djz-1 (3.1.7)


式中，A是常数，ci和dj分别表示零点和极点。由于多项式的系数是实数，ci和dj是实数或者是共轭成对的复数，将共轭成对的零点（极点）放在一起，形成一个二阶多项式，其系数仍为实数； 再将分子、分母均为实系数的二阶多项式放在一起，形成一个二阶网络Hj(z)。


Hj(z)=β0j+β1jz-1+β2jz-21-α1jz-1-α2jz-2(3.1.8)



式中β0j、β1j、β2j、α1j和α2j均为实数。这样H(z)就分解为一些一阶或二阶数字网络的级联形式


H(z)=H1(z)H2(z)…Hk(z)(3.1.9)


式中Hi(z)表示一个一阶或二阶的数字网络的系统函数，每个Hi(z)的网络结构均采用直接型网络结构。
［例3.1.2］设系统函数为H(z)=8-4z-1+11z-2-2z-31-54z-1+34z-2-18z-3，试画出其级联型网络结构。
解将H(z)的分子和分母进行因式分解，得到


H(z)=(2-0.379z-1)(4-1.24z-1+5.264z-2)(1-0.25z-1)(1-z-1+0.5z-2)


为减少单位延迟的数目，将一阶的分子、分母多项式组成一个一阶网络，二阶的分子、分母多项式组成一个二阶网络，画出结构图如图3.1.2所示。


图3.1.2例3.1.2图


由式(3.1.9)可知，H(z)得各个因子Hi(z)(i=1,2,…,k)可以互换位置，因而可以得到不同的结构形式，这将对误差有不同的影响，这里存在一个优化的问题，恰当地选择组合形式，会显著地降低计算误差。
级联型结构中每一个一阶网络决定一个零点、一个极点，每一个二阶网络决定一对零点、一对极点。在式(3.1.8)中，调整β0j、β1j和β2j三个系数可以改变一对零点的位置，调整α1j和α2j可以改变一对极点的位置。因此，相对直接型结构，调整方便是优点。对于硬件实现来说，还可以用一个二阶环节进行分时复用。此外，级联结构中后面的网络输出不会再流到前面，运算误差的积累相对直接型也小。
4. 并联型
对式(3.1.2)作另外一种展开，则得到IIR并联型结构。


  H(z)=H1(z)+H2(z)+…+Hk(z) (3.1.10)


式中，Hi(z)通常为一阶网络或二阶网络，网络系统均为实数，其输出Y(z)表示为


Y(z)=H1(z)X(z)+H2(z)X(z)+…+Hk(z)X(z)  (3.1.11)


式(3.1.11)表明将x(n)送入每个二阶（包括一阶）网络后，将所有输出加起来得到输出y(n)。
［例3.1.3］试画出例题3.1.2中的H(z)的并联型结构。
解将例3.1.2中H(z)展成部分分式形式


H(z)=16+81-0.5z-1+-16+20z-11-z-1+0.5z-2 


将每一部分用直接型结构实现，其并联型网络结构如图3.1.3所示。


图3.1.3例3.1.3图



在这种并联型结构中，每一个一阶网络决定一个实数极点，每一个二阶网络决定一对共轭极点，因此调整极点位置方便，但调整零点位置不如级联型方便。由于基本网络并联，可同时对输入信号进行运算，因此并联型结构与直接型和级联型比较，其运算速度最高。对于许多高速数字信号处理系统来说，这种并联思想已经延拓成了“并行”思想，将一个优秀的信号处理算法进行有效的并行分解，使各并行支路的处理速率在性价比高的处理芯片的工作容许范围内也是数字信号处理发展的一个主要方向之一。
3.1.3FIR滤波系统的基本网络结构
FIR系统可由下面的系统函数或差分方程来表示


H(z)=∑N-1n=0h(n)z-n(3.1.12)

y(n)=∑N-1i=0h(i)x(n-i)(3.1.13)


式(3.1.12)中的系统的脉冲响应h(n)与式(3.1.2)中的bi是直接相关的，即bi=h(i)，其主要实现结构包括直接型与级联型这两种。
1. 直接型
按照H(z)或者差分方程直接画出结构图如图3.1.4所示。这种结构称为直接型网络结构或者称为横向卷积型结构。


图3.1.4FIR直接型网络结构


2. 级联型
通常h(n)为实数，H(z)的零点有两种可能： 即为实数或共轭对复数，每一对共轭零点可以合成一个二阶系统，这样级联型网络结构就是由一阶或二阶因子构成的级联结构，其中每一个因式都用直接型实现。

［例3.1.4］设FIR网络系统函数H(z)如下式


H(z)=1-0.4142z-1-0.4142z-2+z-3


画出H(z)的直接型结构和级联型结构。

解将H(z)进行因式分解，得到


H(z)=(1+z-1)(1-1.4142z-1+z-2)


其级联型结构和直接型结构如图3.1.5所示。


图3.1.5例3.1.4图



级联型结构每一个一阶因子控制一个零点，每一个二阶因子控制一对共轭零点，因此调整零点位置比直接型方便，但H(z)中的系数比直接型多，因而需要更多的乘法器。
3. 广义线性相位FIR系统网络结构
对于长度为N的FIR系统单位脉冲响应h(n)，传输函数为


H(ejω)=∑N-1n=0h(n)e-jωn=H(ω)ejθ(ω)(3.1.14)


式中，H(ω)称为幅度特性，θ(ω)称为相位特性。H(ω)为ω的实函数，可能为负数。H(ejω)线性相位是指θ(ω)是ω的线性函数，即


θ(ω)=-τω，τ为常数 (3.1.15)


如果θ(ω)满足下式


θ(ω)=θ0-τω，θ0是起始相位(3.1.16)


严格地说，此时θ(ω)不具有线性相位，但以上两种情况都满足群时延是一个常数，即


dθ(ω)dω=-τ  (3.1.17)


也称这种情况为广义线性相位。一般地，满足式(3.1.15)为第一类线性相位； 满足式(3.1.16)为第二类线性相位。
可以证明，线性相位FIR系统的单位脉冲响应h(n)应满足下面条件： h(n)为实序列，且满足h(n)=±h(N-1-n)，N为长度，即，h(n)关于N-12偶对称或奇对称。

分析h(n)=h(N-1-n)情况：


N为偶数时


H(z)=∑N-1n=0h(n)z-n=∑N/2-1n=0h(n)z-n+∑N-1n=N/2h(n)z-n（令m=N-1-n）

=∑N/2-1n=0h(n)z-n+∑N/2-1m=0h(N-1-m)z-(N-1-m)

=∑N/2-1n=0h(n)［z-n+z-(N-1-n)］(3.1.18)


N为奇数时

H(z)=∑N-32n=0h(n)［z-n+z-(N-1-n)］+hN-12z-N-12  (3.1.19)

分析h(n)=-h(N-1-n)情况：

N为偶数时


 H(z)=∑N/2-1n=0h(n)［z-n-z-(N-1-n)］(3.1.20)


N为奇数时


H(z)=∑N-32n=0h(n)［z-n-z-(N-1-n)］+hN-12z-N-12


但由于  h(n)=-h(N-1-n)，有hN-12=0，因此


H(z)=∑N/2-3n=0h(n)［z-n-z-(N-1-n)］ (3.1.21)


第一类的N为偶数、N为奇数两种情况的网络结构如图3.1.6所示，第二类的网络结构如图3.1.7所示。



图3.1.6第一类线性相位网络结构





图3.1.7第二类线性相位网络结构




图3.1.8例3.1.5图

［例3.1.5］已知H(z)=1.918(1-3.5z-1+7.75z-2-7.75z-3+3.5z-4-z-5)，画出该FIR滤波器的线性相位结构。
解由第二类线性相位结构可作出如图3.1.8网络结构。 
3.1.4线性相位FIR滤波器零点分布特点
第一类和第二类线性相位的系统函数综合起来满足下式


H(z)=±z-(N-1)H(z-1)  (3.1.22)


线性相位零点分布共有如下4种情况： 
(1) 如果z=z1是H（z）的不在单位圆上零点，其倒数z-1i也必然是其零点； 又因为h（n）是实序列，H（z）的零点必定共轭成对，因此z*1和（z-1i）*也是其零点。这样，线性相位FIR滤波器零点分布特点是零点是互为倒数的共轭对，确定其中一个z1，另外三个零点z*1、（z*1）-1和（z1）-1
也就确定了。
（2） 如果零点是实数，则只有两个零点，即图中z2和z-12情况。


图3.1.9线性相位FIR滤波器

零点分布

（3） 如果零点是纯虚数且在单位圆上，则是图中z3和z3情况。
（4） 如果零点在单位圆上且是实数即为1或-1，则没有其他
零点与该零点对应，如图3.1.9中z4情况。
因此，FIR系统的转移函数H(z)用级联形式的结构实现时，可分别用一阶、二阶、四阶子系统级联而成。每一个子系统都是线性相位的。由此级联而成的整个系统也必定保持线性相位特性。
线性相位FIR系统的另一个特点是，即使其滤波器系数h(n)只经过极为粗糙的量化处理，其线性相位特性还能得到保持，这个优势在实际应用中是不言而喻的。
［例3.1.6］已知一个FIR系统的转移函数为


H(z)=1-2.05z-1+3.2025z-2-1.05z-3-1.05z-4+3.2025z-5-2.05z-6+z-7


分析其零点分布，画出用级联形式实现的网络结构。
解由转移函数可知，N=8，且h(n)为实序列，且偶对称，故为线性相位系统，共有7个零点，为7阶系统，因而必存在一个一阶系统，即z=1或z=-1为系统的零点。而最高阶z-7的系数为+1，所以z=-1为其零点。H(z)中包含1+z-1项，则


H1(z)=H(z)1+z-1=1-3.05z-1+6.2525z-2-7.3025z-3+6.2525z-4-3.05z-5+z-6


此时，N=7，为6阶系统，至少有一个二阶系统，设H2(z)=1+az-1+z-2，余下的四阶子系统为H3(z)=1+bz-1+cz-2+bz-3+z-4，且满足H1(z)=H2(z)·H3(z)，代入等式并展开，可得a=-1，b=-2.05，c=4.2025。因此，H(z)=(1+z-1)(1-z-1+z-2)(1-2.05z-1+3.2025z-2-2.05z-3+z-4)。

系统的全部零点为： z1=-1，z2=ejπ3，z2=e-jπ3，z3=0.8ejπ3，z3=0.8e-jπ3，z-13=1.25e-jπ3，(z-13)=1.25ejπ3。
系统网络图如图3.1.10所示。


图3.1.10例3.1.6图

该例题的MATLAB演示程序如下： 

%广义线性相位系统演示程序

b=［1 -2.05 3.2025 -1.05 -1.05 3.2025 -2.05 1］;

a=［1］;

figure(1)

zplane(b,a);

figure(2)

OMEGA=-pi:pi/100:pi;

H=freqz(b,a,OMEGA);

subplot(211),plot(OMEGA,abs(H));

xlabel('＼omega');ylabel('|H(e ^{j＼omega})|');

subplot(212),plot(OMEGA,180/pi*unwrap(angle(H)));

xlabel('＼omega');ylabel('arg［H(e ^{j＼omega})］');

结果如图3.1.11所示。



图3.1.11例3.1.6的MATLAB图

3.1.5数字滤波系统的MATLAB实现
对于一个给定输入和输出关系的数字滤波系统，其系统函数可通过多种算法实现，而不同的算法对应的网络结构也各有不同。从网络结构可以清晰地看到滤波系统的运算步骤、加乘法运算次数和存储单元的数量，这对于数字滤波器软、硬件的实现至关重要。数字滤波系统的网络结构也是数字滤波器设计的一项非常重要的内容，关系到数字滤波器的稳定性、运算速度以及系统的成本和体积等许多重要的性能。下面是两个数字滤波系统的MATLAB实现例子。
［例3.1.7］用直接型实现系统函数为H(z)=1-3z-1+11z-2+27z-3+18z-41+16z-1+12z-2+2z-3-4z-4-z-5的IIR数字滤波器，求单位脉冲响应和单位阶跃信号的输出。
解程序清单如下： 

b=［1,-3,11,27,18］;a=［16,12,2,-4,-1］;

N=25;

h=impz(b,a,N);                 %直接型单位脉冲响应

x=［ones(1,5),zeros(1,N-5)］;       %单位阶跃信号

y=filter(b,a,x);                  %直接型输出信号

subplot(1,2,1);stem(h);title('直接型h(n)');

subplot(1,2,2);stem(y);title('直接型y(n)');

结果如图3.1.12所示。


图3.1.12例3.1.7的MATLAB图


［例3.1.8］FIR滤波器的系统函数为H(z)=0.2n,0≤n≤5


0,其他，试用直接型实现。
解程序清单如下： 

n=0:5;

b=0.2.^n;

N=30;

delta=impseq(0,0,N);

h=filter(b,1,delta);

x=［ones(1,5),zeros(1,N-5)］;

y=filter(b,1,x);

subplot(1,2,1);stem(h);title('直接型h(n)');

subplot(1,2,2);stem(y);title('直接型y(n)');

%--------------------------------------------------------------------

%单位脉冲响应δ(n-n_0 )的生成函数impseq.m

function［x,n］=impseq(n0,ns,nf);

n=［ns:nf］;x=［(n-n0)==0］;

结果如图3.1.13所示。


图3.1.13例3.1.8的MATLAB图



3.2数字滤波器的基本概念
最常用的滤波器就是频率选择性滤波器。所谓频率选择滤波，就是让一个或一组频率范围内的信号尽可能无失真地通过，且衰减或者完全抑制其余滤波范围的信号。先来看一个例子。

假设原始信号为xo(t)=sin (2π×80t)+2sin(2π×140t)，由于某些原因，信号被另一个频率的信号xN(t)=sin (2π×300t)干扰，实际的信号变为x(t)=xo(t)+xN(t)。为了减小干扰信号的影响，恢复原始信号，需要对这个受到干扰的信号进行处理。由于干扰信号的频率成分高于原信号，可让这个受干扰信号通过一个系统，使其低频部分通过这个系统，而高频部分被限制通过。这样的系统就是称为滤波器。
取采样频率为fs=1000Hz，采用MATLAB演示这个过程。

%频率选择性滤波器功能演示

Fs=1000; t=0:1/Fs:2;%取2秒长度的信号

x0=sin(2*pi*80*t)+2*sin(2*pi*140*t);  %原始信号

xN=sin(2*pi*300*t);  %噪声信号

x=x0+xN;  %受污损信号x(t)

figure(1);

subplot(211);plot(t,x0);axis(［0,0.25,-4,4］);

xlabel('t');ylabel('x_{0}(t)');title('原始信号x_{0}(t)');

subplot(212);plot(t,x);axis(［0,0.25,-4,4］);

xlabel('t');ylabel('x(t)');title('受污损信号x(t)');

%设计一个特定的频率选择性滤波器

n=100;  %滤波器长度取100

f=［0 0.13 0.15 0.17 0.19 0.25 0.27 0.29 0.31 1］;

m=［0 0 1 1 0 0 1 1 0 0］;

b=firls(n,f,m); 

［H,W］=freqz(b,1,512,2);figure(2);

plot(W,20*log10(abs(H)));grid

xlabel('归一化频率');ylabel('滤波器的对数幅频特性');

x1=filter(b,1,x);  %完成滤波功能

figure(3);

subplot(211);plot(t,x1);axis(［0,0.25,-4,4］);

xlabel('t');ylabel('x_{1}(t)');title('滤波输出信号x{1}(t)');

subplot(212);plot(t,x1-x0);  %误差信号

xlabel('t');ylabel('x_{1}(t)-x_{0}(t)');title('误差信号x_{e}(t)');

结果如图3.2.1所示。




图3.2.1频率选择性滤波器功能演示



另一类广泛应用的类型是频率成形滤波器，例如信号锐化滤波器、信号平滑滤波器、频率补偿滤波器等，主要的目的是改变信号的频谱形状，进而还原信号的频域特性。在信号处理系统中，常用于系统滤波建模处理。
下面首先设计一个6阶巴特沃斯低通滤波器作为目标响应滤波器，并采用文献［27］中的PRONY模型对目标响应进行拟合建模，拟合模型设置也为6阶，对目标滤波器的脉冲响应进行拟合，在未知滤波器具体设置参数时实现对低通滤波器还原，其脉冲响应与幅度曲线基本一致。这种方法一般可用于还原未知设备的系统函数，通过输入信号和输出信号即可实现对系统设备的整体滤波器建模。

%频率成形滤波器应用演示

d = designfilt('lowpassiir','NumeratorOrder',6,'DenominatorOrder',6, ...

'HalfPowerFrequency',0.2,'DesignMethod','butter'); %设计目标低通滤波器

h = filter(d,［1 zeros(1,49)］);%获取目标滤波器的脉冲响应

bord = 6;aord = 6; %设置拟合参数

［b,a］ = prony(h,bord,aord);%基于PRONY模型对目标响应拟合

figure(1);

subplot(2,1,1);plot(1:50,h,'-b');title('目标脉冲响应');

subplot(2,1,2);plot(1:50,impz(b,a,length(h)),'-r');title('基于PRONY模型拟合的脉冲响应');

［H,W］=freqz(d,512,2);

figure(2);

subplot(2,1,1);plot(W,20*log10(abs(H)),'-b');grid

legend('目标滤波器的幅频特性');

subplot(2,1,2);plot(W,20*log10(abs(freqz(b,a,512,2))),'-r');grid

xlabel('归一化频率');legend('拟合滤波器的幅频特性');axis(［0 1 -400 200］);

如图3.2.2所示的功能演示。


图3.2.2滤波器建模功能演示






图3.2.2（续）


滤波器在通信系统中有一个专用名词——均衡器，利用频域均衡可以有效地补偿传输信道色散引起的频率失真。滤波器均衡和滤波器建模的原理相似，主要应用滤波器的频率成形功能，但滤波器均衡的应用目的与滤波器建模刚好相反，均衡本身是对系统的特性进行改变，进而实现改善修复系统的传输性能。
理想情况下，为了保证音乐信号无损传输，扬声器的幅度频率响应曲线应为一条水平的直线，但实际情况却如图3.2.3所示，参照滤波器建模的方法，大家可以思考如何采用滤波器均衡针对此曲线进行修复？具体的滤波器设计是什么形式？如何进行设计？


图3.2.3扬声器幅度响应示例


除此之外，数字滤波器在语音处理、图像处理、通信、音乐等各个方面都有着极为广泛的应用。随着信号处理技术的发展，一些现代滤波器如卡尔曼滤波器，维纳滤波器等都在实际应用中发挥了重要的作用。本节讨论的重点是频率选择性滤波器设计的相关内容，对频率形成滤波器仅以数字微分为例进行简要介绍，而现代滤波器的有关知识要在后续课程里才讨论得到。

3.2.1频率选择性滤波器
频率选择性滤波器主要分成4种，即低通（LP）、高通（HP）、带通（BP）、带阻（BS）滤波器，每一种又可分成模拟滤波器（AF）和数字滤波器（DF）两种形式，图3.2.4和图3.2.5分别给出了AF和DF的4种滤波器理想幅频响应。



图3.2.4理想模拟滤波器幅度特性


 

图3.2.5理想数字滤波器幅度特性


根据处理信号的性质不同，而选用模拟或数字滤波器，本课程主要学习数字滤波器的设计方法，但与模拟滤波器设计密切相关。应该注意的是，这两类滤波器有不同的特点，数字滤波器的幅度特性以2π为周期，在ω=±(2K+1)π，K=0,1,2,…周围具有高频特性。研究数字滤波器只需研究一个周期的特性即可，一般考虑［0,2π］或者［-π,π］。
对于数字低通滤波器，其频率响应特性可表达为


 H(ejω)=1,|ω|≤ωc


0,ωc<|ω|≤π（3.2.1）


其脉冲响应h(n)为


h(n)=ωcπsinc(ωcn)


很显然，该系统是非因果的，在实践中它是不可实现。根据PaleyWiener定理，h(n)能量有限且对n<0，h(n)=0，即系统为因果系统，则有需要满足  


∫π-πln H(ejω)dω<∞


因此，在一些频率点上，幅度函数H(ejω)可以为0，但是不能在任何有限的频带上均为0（积分变为无限）。只能按某些准则来设计滤波器，使之尽可能地逼近理想滤波器特性。以工程的角度上来说，衡量这种逼近效果好坏的标准是该滤波器一系列的技术指标。需要在严格的技术指标和实现的复杂度之间寻找一个良好的折中。
3.2.2滤波器的技术指标
数字滤波器的传输函数H(ejω)用式（3.2.2）表示


H(ejω)=H(ejω)ejQ(ω)（3.2.2）


式中，H(ejω)称为幅频特性； Q(ω)称为相频特性。
幅频特性表示信号通过该滤波器后各频率成分的衰减情况，而相频特性反映各频率成分通过滤波器后在时间上的延时情况。因此，一般选频的技术要求仅由幅频特性给出，只有当对输出波形有要求，才需要考虑相频特性的技术指标。


图3.2.6低通滤波器的技术要求

可实现的低通滤波器的幅度特性如图3.2.6所示，ωp和ωs分别称为通带截止频率和阻带截止频率。通带频率范围为0≤ω≤ωp，在通带中要求(1-δ1)<H(ejω)≤1。阻带频带范围为ωs≤ω≤π，在阻带中要求H(ejω)≤δ2，不一定衰减到零。从ωp到ωs称为过渡带，一般为单调下降的。通带内允许的最大衰减用αp表示，阻带内允许的最小衰减用αs表示，αp和αs用dB分别定义为



αp=20lgH(ej0)H(ejωp)dB （3.2.3）

αs=20lgH(ej0)H(ejωs)dB（3.2.4）


如将H(ejω)归一化为1，式（3.2.3）和式（3.2.4）表示为


αp=-20lgH(ejωp)（3.2.5）

αs=-20lgH(ejωs) （3.2.6）


当幅度下降到22≈0.707时，对应的频点ω=ωc，此时αp=3dB，称ωc为3dB通带截止频率，它是滤波器设计的重要参数之一。
高通、带通和带阻滤波器的技术指标可类似给出。设计滤波器时应该根据指标参数及对滤波特性的要求，选择适合的滤波器类型和设计方法进行设计。
3.2.3数字滤波器的设计方法
数字滤波器的设计大致包括3个步骤： 
（1） 给出所需要的滤波器的技术指标； 
（2） 设计一个H（z）使其逼近所需要的技术指标； 
（3） 实现所设计的H（z）。
通常，频率选择性滤波器可利用IIR无限脉冲滤波器和FIR有限脉冲滤波器来设计，但两种不同形式的滤波器的设计不同。IIR滤波器的设计方法主要是借助于模拟滤波器转换的设计方法进行的，因为IIR滤波器的单位脉冲响应h(n)为无限长序列，无法由h(n)确定网络结构，而其系统函数H(z)有限，所以设计结果是滤波器系统函数H(z)。FIR的设计方法主要是建立在理想滤波器频率特性作某种近似的基础上的，近似方法有窗函数法、频率采样法，优化设计方法等，其单位脉冲响应h(n)就是直接型网络结构的乘法因子，所以设计的结果是h(n)。
而对于线性相位滤波器，通常采用FIR滤波器，其单位脉冲响应满足一定条件，可以证明其相位特性在整个频带中是严格线性的，这是模拟滤波器无法达到的。当然，也可以采用IIR滤波器，但必须使用全通网络对其线性相位特性进行相位校正，这样增加了设计与实现的复杂性。
FIR与IIR的主要性能对比如表3.2.1所示。


表3.2.1FIR与IIR的主要性能对比表



FIR滤波器
IIR滤波器
无循环反馈
有循环反馈
常稳定，极点只有零点
必须考虑稳定性
截止特性差，阶数N较大
截止特性好
可获得广义线性相位特性
只能与全通相位校正网络级联，以获得近似线性相位特性
不必考虑量化误差的扩大
必须考虑量化误差的扩大
设计方式灵活
主要用于设计具有片段常数特性的滤波器

3.3IIR型滤波器的设计
IIR数字滤波器设计的最通用的方法是借助于模拟滤波器的设计方法。其设计步骤如下： 
（1） 按一定规则将给出的数字滤波器的技术指标转换为模拟低通滤波器的技术指标。
（2） 根据转换后的技术指标设计模拟低通滤波器H（s）。
（3） 再按一定规则将H（s）转换成H（z）。
若所设计的数字滤波器是低通的，那么上述设计工作可以结束； 若所设计的是高通、带通或带阻滤波器，那么还有下面的步骤。
（4） 将高通，带通或带阻数字滤波器的技术指标转化为低通模拟滤波器的技术指标，然后按上述步骤（2）设计出低通H（s），再将H（s）转换为所需的H（z）。
由上述步骤可知，在设计数字滤波器时，先设计模拟低通滤波器，再通过频率变换将其转换成希望滤波器的类型，模拟滤波器的设计方法已经相当成熟，有着大量的现存图表结果可以查阅，而且MATLAB软件包中也包含着许多功能强大的设计调用函数。下面将先在理论上介绍模拟低通滤波器的设计方法。
3.3.1模拟低通滤波器
模拟低通滤波器的设计指标包括αp、Ωp、αs和Ωs，其中Ωp和Ωs分别为通带截止频率和阻带截止频率，要设计一个低通滤波器Ha(s)为


Ha(s)=d0+d1s+…+dN-1sN-1+dNsNc0+c1s+…+cN-1sN-1+cNsN 


使其对数幅频响应20lgH(jΩ)在Ωp，Ωs出分别达到αp，αs的要求。


αp=20lgHa(j0)Ha(jΩp)（3.3.1）

αs=20lgHa(j0)Ha(jΩs)  （3.3.2）


在Ω=0处幅度归一化到1，即Ha(j0)=1，αp和αs表示为


αp=-20lgHa(jΩp)（3.3.3） 

αs=-20lgHa(jΩs)（3.3.4）


如图3.3.1所示，图中Ωc为3dB截止频率。


图3.3.1低通滤波器的幅频特性


由于一般滤波器的单位脉冲响应为实数，所以


Ha(jΩ)2=Ha(s)Ha(-s)|s=jΩ=Ha(jΩ)H*a(jΩ)（3.3.5）


只要求出Ha(jΩ)2，就可以求出所需要的Ha(s)，因此幅度平方函数在模拟滤波器的设计中起了很重要的作用。而Ha(s)必须是稳定的，所以极点必须落在S平面的左平面，相应的Ha(-s)的极点落在右半平面。
3.3.2巴特沃斯低通滤波器的设计
Buttlerworth（巴特沃斯）低通滤波器的设计方法可按以下3个步骤来进行。
1. 确定阶数N
将实际频率Ω归一化，得到归一化幅频平方特性


Ha(jΩ)2=11+ΩΩc2N（3.3.6）


式中，N为滤波器的阶数。Ωc为3dB截止频率，即当Ωc=Ω时，Ha(jΩ)=12； 
当Ω=0时，Ha(jΩ)=1。即在Ω=0处无衰减； 
当0<Ω<Ωc时，Ha(jΩ)随着Ω缓慢减小，N越大，减小得越慢，即在通带内Ha(jΩ)越平坦； 
当Ω=Ωc时，Ha(jΩ)=1/2； 
当Ω>Ωc时，Ha(jΩ)随着Ω迅速衰减，N越大，衰减得速度越快，即过渡带越窄。
幅频特性与Ω和N得关系如图3.3.2所示。所以阶数N的大小影响着幅度特性衰减速度，因此它由技术指标αp、Ωp、αs和Ωs确定。将Ω=Ωp和Ω=Ωs分别代入式（3.3.6），再将所得的幅频平方函数带入式（3.3.3），可得


1+ΩpΩc2N=10αp/10  （3.3.7）

1+ΩsΩc2N=10αs/10（3.3.8）




图3.3.2巴特沃斯幅度特性和N的关系


根据式（3.3.7）和式（3.3.8）可得


ΩpΩsN=100.1αp-1100.1αs-1  


令λsp=Ωs/Ωp，ksp=100.1αp-1100.1αs-1，可求得


N=-lgksplgλsp （3.3.9）


当N有小数时，取大于N的最小整数。
2. 求归一化传输函数
G（p）
将幅度平方函数Ha(jΩ)2写成s的函数


Ha(s)Ha(-s)=11+sjΩc2N（3.3.10） 


可知幅度平方函数有2N个极点


sk=(-1)12N(jΩc)ej2π2Nk

=Ωcejπ12+2k+12N，k=0,1,…,2N-1（3.3.11）




图3.3.3三阶巴特沃斯滤波器

极点分布

2N个极点等间隔分布在半径为Ωc的圆上，称该圆为巴特沃斯圆。为了形成稳定的滤波器，2N个极点中只有取s平面左半平面的N个极点构成Ha(s)，其表达式为


Ha(s)=ΩNc∏N-1k=0(s-sk)  （3.3.12）


如N=3，有6个极点，如图3.3.3所示，它们分别为


s0=Ωcej23π，s1=-Ωc，s2=Ωce-j23π

s3=Ωce-j13π，s4=Ωc，s5=Ωcej13π


取S平面左半平面的极点s0、s1、s2组成Ha(s)


Ha(s)=Ω3c(s+Ωc)(s-Ωcej23π)(s-Ωce-j23π)


对3dB截止频率Ωc归一化，Ha(s)可改写为


Ha(s)=1∏N-1k=0sΩc-skΩc（3.3.13） 


式中，s/Ωc=jΩ/Ωc。
令λ=Ω/Ωc，p=jλ，p称为归一化拉氏复变量，得巴特沃斯传输函数为


G(p)=1∏N-1k=0(p-pk)（3.3.14）


式中，pk为归一化极点，表示为


pk=ejπ12+2k+12N，k=0,1,…,N-1（3.3.15）  


因为N已求得，所以也可以通过查表3.3.1来确定G(p)。
3. 求出Ωc，将G（p）去归一化，得到实际的滤波器传输函数Ha（s）
根据式（3.3.7）和式（3.3.8）可得


Ωc=Ωp(100.1αp-1)-12N（3.3.16）

Ωc=Ωs(100.1αs-1)-12N （3.3.17）


将p=s/Ωc带入G(p)中得到Ha(s)。
表3.3.1列出了常见的巴特沃斯归一化低通滤波器参数，设计时可直接使用这些结果。


表3.3.1巴特沃斯归一化低通滤波器参数




 极点位置



阶数N

P0,N-1
P1,N-2
P2,N-3
P3,N-4
P4,N-5

1
-1.0000

2
-0.7071±j0.7071

3
-0.5000±j0.8660
-1.0000
4
-0.3827±j0.9239
-0.9239±j0.3827
5
-0.3090±j0.9511
-0.8090±j0.5878
-1.0000
6
-0.2588±j0.9659
-0.7071±j0.7071
-0.9659±j0.2588
7
-0.2225±j0.9749
-0.6235±j0.7818
-0.9010±j0.4339
-1.0000
8
0.1951±j0.9808
0.5556±j0.8315
-0.8315±j0.5556
-0.9808±j0.1951
9
-0.1736±j0.9848
-0.5000±j0.8660
-0.7660±j0.6428
-0.9397±j0.3420
-1.0000


分母

多项式


系数

阶数N



B(p)=pN+bN-1 pN-1+bN-2pN-2+…+b1p+b0

b0
b1
b2
b3
b4
b5
b6
b7
b8


1
1.0000
2
1.0000
1.4142
3
1.0000
2.0000
2.0000
4
1.0000
2.6131
3.4142
2.613

5
1.0000
3.2361
5.2361
5.2361
3.2361
6
1.0000
3.8637
7.4641
9.1416
7.4641
3.8637
7
1.0000
4.4940
10.0978
14.5918
14.5918
10.0978
4.4940
8
1.0000
5.1258
13.1371
21.8462
25.6884
21.864
13.1371
5.1258
9
1.0000
5.7588
16.5817
31.1634
41.9864
41.9864
31.1634
16.5817
5.7588


分母因式



阶数N



B(p)=B1(p)B2(p) B3(p)B4(p) B5(p)


1
(p+1)
2
(p2+1.4142p+1)
3
(p2+p+1)(p+1)
4
(p2+0.7654p+1)(p2+1.8478p+1) 
5
(p2+0.6180p+1)(p2+1.6180p+1)(p+1)
6
(p2+0.5176p+1)(p2+1.4142p+1)(p2+1.9319p+1)
7
(p2+0.4450p+1)(p2+1.2470p+1)(p2+1.8019p+1)(p+1)
8
(p2+0.3902p+1)(p2+1.1111p+1)(p2+1.6629p+1)(p2+1.9616p+1)
9
(p2+0.3473p+1)(p2+p+1)(p2+1.5321p+1)(p2+1.8794p+1)(p+1)


［例3.3.1］已知通带截止频率fp=5kHz，通带最大衰减αp=2dB,阻带截止频率fs=12kHz, 阻带最小衰减αs=30dB。请按照以上指标设计巴特沃斯低通滤波器。
解（1）  确定阶数N


ksp=100.1αp-1100.1αs-1=0.0242

λsp=2πfs2πfp=2.4

N=-lg0.0242lg2.4=4.25，取N=5


（2） 由式（3.3.15），其极点为


p0=ej35π，p1=ej45π，p2=ejπ，
p3=ej65π，p4=ej75π


由式（3.3.14），归一化传输函数为


G(p)=1∏4k=0(p-pk)


上式分母可以展开成为五阶多项式，或者将共轭极点放在一起，形成因式分解形式。查表，由N=5，得到
极点： -0.3090±j0.9511，-0.8090±j0.5878，-1.0000


G(p)=1p5+b4p4+b3p3+b2p2+b1p+b0


式中b0=1.0000，b1=3.2361，b2=5.2361，b3=5.2361，b4=3.2361。


G(p)=1(p2+0.6180p+1)(p2+1.6180p+1)(p+1)


（3） 为将G(p)去归一化，先求3dB截止频率Ωc。
根据式（3.3.16），得到


Ωc=Ωp(100.1αp-1)-12N=2π·5.2755krad/s


将Ωc代入（3.2.17）式，得到


Ωs=Ωc(100.1αs-1)12N=2π·10.2525krad/s


此时算出得Ωs比题目中给出的小，因此，过渡带小于要求，或者说在Ωs=12krad/s时衰减大于30dB，所以说阻带指标有富余量。
将p=s/Ωc带入G(p)中得到


Ha(s)=Ω5cs5+b4Ωcs4+b3Ω2cs3+b2Ω3cs2+b1Ω4cs+b0Ω5c


实际上，上例中的设计也可以直接利用MATLAB来完成。

%Butterworth模拟低通滤波器的设计演示

Wp=2*pi*5000;Ws=2*pi*12000;Ap=2;As=30; %滤波器的技术参数

［n,Wn］=butterd(Wp,Ws,Ap,As,'s');	%获得滤波器参数

［b,a］=butter(n,Wn,'s');	%设计butterworth LPF

figure(1);

freqs(b,a,20000);  %画出设计出的LPF频率响应曲线

figure(2);

Wc=2*pi*10.2525*1000;b0=1.0000;b1=3.2361;

b2=5.2361;b3=5.2361;b4=3.2361;  %讨论求出的导数值

b=［Wc.^5］;

a=［1 b4*Wc b3*(Wc.^2) b2*(Wc.^3) b1*(Wc^4)(Wc^5)］;

freqs(b,a,20000);  %画出讨论求出的LPF频率响应特性

结果如图3.3.4所示。



图3.3.4例3.3.1图


显然，由MATLAB直接设计出来的巴特沃斯LPF与讨论计算结果一致的。
巴特沃斯滤波器的特点是具有通带内最平坦的幅度特性。所以又称为“最平”幅频响应滤波器。如果在通带边缘满足指标，则在通带内必有富余量。一种更为有效的办法是将指标的精度要求均匀地分布在通带内，或者均匀地分布在阻带内，或同是分布在通带与阻带内，这时就可设计出较低的滤波器，这种具有等波纹特性的三种精度均匀分布情况分别对应着Chebyshev（切比雪夫）Ⅰ型，Chebyshev Ⅱ型的椭圆滤波器显然相比Butterworth而言，椭圆滤波器的阶数较低，Chebyshev Ⅰ和Chebyshev Ⅱ型次之。在这里仅以例3.3.1中的模拟低通指标为例，介绍一下这三种滤波器的MATLAB设计方法。

%Chebyshev I,II型，椭圆滤波器设计演示

Wp=2*pi*5000; 	Ws=2*pi*12000;Ap=2;  As=30;  %LPF技术指标

% Chebyshev I型LPF设计

［n1,Wn1］=cheb1ord(Wp,Ws,Ap,As,'s');  %获取滤波器阶数

［b1,a1］=cheby1(n1,Ap,Wn1,'s');  %设计Chebyshev Ⅰ型LPF

figure(1);

freqs(b1,a1,20000);  %画出Chebyshev Ⅰ型LPF频率响应特性曲线

% Chebyshev Ⅱ型LPF设计

figure(2);

［n2,Wn2］=cheb2ord(Wp,Ws,Ap,As,'s');  %获取滤波器阶数

［b2,a2］=cheby2(n2,As,Wn2,'s');  %设计Chebyshev Ⅱ型LPF

freqs(b2,a2,20000);  %画出Chebyshev Ⅱ型LPF频率响应特性曲线

%椭圆LPF设计

figure(3);

［n3,Wn3］=ellipord(Wp,Ws,Ap,As,'s');

［b3,a3］=ellip(n3,Ap,As,Wn3,'s');

freqs(b3,a3,20000);  %画出椭圆LPF频率响应特性曲线

程序运行结果如图3.3.5所示。


图3.3.5Chebyshev Ⅰ型、Chebyshev Ⅱ型及椭圆LPF频率响应特性曲线



图3.3.5（续）

3.3.3模拟滤波器的频率转换——模拟高通、带通及

带阻滤波器的设计

以上较为详细地了解了巴特沃斯低通滤波器的设计方法。模拟高通、带通及带阻可以利用模拟滤波器转换得到。它们的设计方法是先将要设计的滤波器的指标通过某种频率转换关系转成模拟低通滤波器底指标，并依据这些指标设计出低通转移函数，然后再依据频率转换关系变成所要设计的滤波器的转移函数。
1. 模拟高通滤波器的设计
设低通滤波器G(p)和高通滤波器H(q)的幅频特性如图3.3.6所示，其中p=jλ，q=jη。图中λp、λs分别称为低通滤波器的归一化通带截止频率和归一化阻带截止频率，ηp、ηs分别称为高通滤波器的归一化下限频率和归一化阻带上限频率。


图3.3.6低通与高通滤波器的幅频特性


由于|G(jλ)|和|H(jη)|都是频率的偶函数，可以把|G(jλ)|右边曲线和|H(jη)|曲线对应其来，低通λ从∞经过λs和λp到0时，高通的η则从0经过ηs和ηp到∞，因此λ和η之间的关系为


λ=1/η （3.3.18）


从低通滤波器转换到模拟高通滤波器的设计步骤如下： 
（1） 确定高通滤波器的技术指标，例如，通带下限频率Ω′p，阻带上限频率Ω′s，通带最大衰减αp，阻带最小衰减αs。
（2） 确定相应低通滤波器的设计指标，根据式（3.3.18），将高通滤波器的边界频率转换成低通滤波器的边界频率。各项设计指标如下： 
① 低通滤波器通带截止频率Ωp=1/Ω′p； 
② 低通滤波器阻带截止频率Ωs=1/Ω′s； 
③ 通带最大衰减仍为αp，阻带最小衰减仍为αs。
（3） 设计归一化低通滤波器G（p）。
（4） 求模拟高通的H（s）。将G（p）按照式（3.3.18），转换成归一化高通H（q）=G（p）p=1q，为去归一，将q=s/Ωc代入H（q）中得可H（s），即


H(s)=H(q)q=sΩc
=G(p)p=Ωcs （3.3.19）


上式就是由归一化低通直接转换成模拟高通的转换公式。
［例3.3.2］设计高通滤波器，fp=200Hz，fs=100Hz，幅度特性单调下降，fp处最大衰减为αp=3dB，阻带最小衰减αs=15dB。
解(1) 高通技术要求


fp=200Hz，αp=3dB，fs=100Hz，αs=15dB


归一化频率


ηp=fpfc=1，ηs=fsfc=0.5


（2） 低通技术要求


λp=1，λs=1ηs=2，αp=3dB，αs=15dB


（3） 设计归一化低通G(p)。采用巴特沃斯滤波器，故


ksp=100.1αp-1100.1αs-1=0.18

λsp=λsλp=2

N=-lgksplgλsp=2.47，取N=3

G(p)=1p3+2p2+2p+1


（4） 求模拟高通H(s)


 H(s)=G(p)p=Ωcs=s3s3+2Ωcs2+2Ω2cs+Ω3c


式中Ωc=2πfp。

MATLAB演示程序如下： 

%模拟高通——高通演示

%直接采用MATLAB设计该高通模拟滤波器

Wp=2*pi*200;Ws=2*pi*100;Ap=2;As=15;  %HPF性能指标

［n,Wn］=butterd(Wp,Ws,Ap,As,'s');	  %获得滤波器参数

［b,a］=butter(n,Wn,'high','s');	 %设计butterworth HPF

figure(1);

freqs(b,a,20000);	 %画出设计出的HPF频率响应曲线

Wc=2*pi*200;

b1=［1	0	0	0］;a1=［1	2*Wc	2*Wc.^2		Wc.^3］;

figure(2);

freqs(b1,a1,2000); %画出理论求出的HPF频率响应特性曲

运行结果如图3.3.7所示。


图3.3.7例3.3.2图


2. 模拟带通滤波器的设计
低通滤波器与带通滤波器的幅度特性如图3.3.8所示。


图3.3.8带通与低通滤波器的幅度特性


从低通滤波器转换到模拟带通滤波器的设计步骤如下： 
（1） 确定模拟带通滤波器的技术指标（见图3.3.8）。
带通上限频率Ωu，带通下限频率Ωl，下阻带上限频率Ωs1，上阻带下限频率Ωs2，通带中心频率Ω0，Ω20=ΩlΩu，通带宽度B=Ωu-Ωl。与以上边界频率对应的归一化边界频率如下


ηs1=Ωs1B, ηs2=Ωs2B,

ηl=ΩlB，ηu=ΩuB，η20=ηlηu


将带通和低通的幅度特性对应起来，可以到λ和η的对应关系如表3.3.2所示。


表3.3.2λ和η的对应关系



λ
-∞
-λs
-λp
0
λp
λs
 ∞
η
0
ηs1
ηl
η0
ηu
ηs2

∞
λ和η的对应关系为


λ=η2-η20η  （3.3.20）


还需要确定的技术指标有： 通带最大衰减αp，阻带最小衰减αs。
（2） 确定归一化模拟低通技术要求


λp=1，λs=η2s2-η20ηs2，-λs=η2s1-η20ηs1


λs与-λs的绝对值可能不相等，一般取绝对值较小的λs，这样保证在较大的λs处能满足要求。通带最大衰减仍为αp，阻带最小衰减亦为αs。
（3） 设计归一化模拟低通G（p）。
（4） 直接将G（p）转换成带通滤波器H（s）。

由于p=jλ
代入式（3.3.20）可得


p=jη2-η20η


将q=jη代入上式可得


p=q2+η20q


去归一化，将q=s/B代入上式



p=s2+ΩlΩus(Ωu-Ωl) （3.3.21）


因此

H(s)=G(p)p=s2+ΩlΩus(Ωu-Ωl)（3.3.22）

3. 模拟带阻滤波器的设计
低通与带通滤波器得幅频特性如图3.3.9所示。


图3.3.9低通与带阻滤波器的幅频特性


下面介绍从低通滤波器转换到模拟带阻滤波器的设计步骤。
（1） 确定模拟带阻滤波器的设计要求。
上通带截止频率Ωu，下通带截止频率Ωl，阻带下限频率Ωs1，阻带上限频率Ωs2，阻带中心频率Ω20=ΩlΩu，阻带宽度B=Ωu-Ωl。
以上边界频率对应的归一化边界频率如下： 


ηs1=Ωs1B,ηs2=Ωs2B,

ηl=ΩlB，ηu=ΩuB，η20=ηlηu


将带阻和低通的幅度特性对应起来，可以得到λ和η的对应关系如表3.3.3所示。


表3.3.3λ和η的对应关系



λ
-∞
-λs
-λp
0
0
λp
λs
∞ 
η
η0
ηs2
ηu
∞
0
ηl
ηs1
η0


得到λ和η的关系为


 λ=ηη2-η20  （3.3.23）  


还需要确定的技术指标有： 通带最大衰减αp，阻带最小衰减αs。
（2） 确定归一化模拟低通技术要求


λp=1，-λs=ηs1η2s1-η20，λs=ηs2η2s2-η20


取λs与-λs绝对值较小的λs，这样保证在较大的λs处能满足要求； 通带最大衰减仍为αp，阻带最小衰减亦为αs。

（3） 设计归一化模拟低通G(p)。
（4） 直接将G(p)转换成带阻滤波器H(s)。
可以得到


p=s(Ωu-Ωl)s2+ΩlΩu （3.3.24）


因此  


H(s)=G(p)p=s(Ωu-Ωl)s2+ΩlΩu （3.3.25）


3.3.4模拟与数字滤波器的转换方法
要得到数字滤波器，还需要对上述方法设计的模拟滤波器进行转换，本节将介绍两种方法。
1. 脉冲响应不变法
脉冲响应不变法本质上是一种时域逼近方法。利用数字滤波器的脉冲响应
h(n)与模拟滤波器的脉冲响应在采样点上的值ha（t）t=nT相等，即


h(n)=ha（t）t=nT（3.3.26）


来得到变换关系的。如果给定的模拟滤波器为Ha(s)，则其拉
普拉斯反变换为ha(t)，对其采样得到h(n)，再对h(n)进行Z变换即可得到H(z)。
 以具有单阶极点的系统Ha(s)为例，如果它分母多项式的阶次高于分子多项式的阶次，则可以表达为


Ha(s)=∑Ni=1Ais-si（3.3.27）


则

ha(t)=∑Ni=1Aiesitu(t)（3.3.28）


ha(t)采样有


h(n)=ha(nT)=∑Ni=1AiesinTu(nT) （3.3.29）


则

H(z)=∑Ni=1Ai1-esiTz-1 （3.3.30）

极点情况更复杂的情形，请读者查阅相关的书籍。
比较式（3.3.27）与式（3.3.30），可以看出Z平面与S平面的映射关系： 当s=si时，有z=esiT，映射z=esT将S左半平面映射成Z平面的单位圆。因此，如果模拟滤波器是稳定的，数字滤波器也是稳定的，数字滤波器保持了模拟滤波器的时域瞬态特性。但由于该方法进行了采样，必须要满足采样定理，如果模拟滤波器的频响不是真正带限，用这种脉冲响应不变法设计的数字滤波器的频谱要发生混叠，系统将失真，在使用时要特别注意。由于高通滤波器和带阻滤波器的频带都不是带限的，因此，不能将这种方法应用于高通和带阻滤波器的设计。
［例3.3.3］利用脉冲响应不变法，将系统函数为


Ha(s)=1s2+s+1


的模拟滤波器转换成数字IIR滤波器，并分别取T=0.05和T=0.3利用MATLAB来演示频谱的混叠现象。
解将Ha(s)进行分解，得


Ha(s)=13js--1+3j2-13js--1-3j2


模拟滤波器在s1,2=-1±3j2处有一对共轭极点，根据脉冲响应不变法，相应的数字滤波器的系统函数为


H(z)=13j1-e-1+3j2Tz-1-13j1-e-1-3j2Tz-1

=e-12Tsin 32Tz-11-2e-12Tz-1cos 32T+e-Tz-2

=23e-12Tsin 32Tz-11-e-T2-32jTz-11-e-T2+32jTz-1




%脉冲不变法频域混叠演示

b=［0 0 1］;a=［1 1 1］;

figure(1);

freqs(b,a,10000);

%分析T=0.05和T=0.3情况下的混叠效果

［b1,a1］=impinvar(b,a,1/0.05);

figure(2);

freqz(b1,a1);

［b2,a2］=impinvar(b,a,1/0.3);

figure(3);

freqz(b2,a2);

结果如图3.3.10所示。


图3.3.10例3.3.3脉冲不变法频域混叠




图3.3.10（续）


2. 双线性变换法
由于频谱带宽的严格限制，脉冲响应不变法在应用中受到限制，而双线性变换法则是采用非线性频率压缩，将整个频率轴上的频率范围压缩到±π/T之间，再用z=esT转换到Z平面上，避免了频谱混叠。图3.3.11显示了这一变换过程。首先用反正切变换


Ω1=2Tarctan ΩT2 （3.3.31）




图3.3.11双线性变换法的映射关系

实现频率压缩，当Ω从-∞经过0变化到+∞时，Ω1则从-π/T经过0变化到π/T，实现了S平面上整个虚轴到±π/T之间的变换。令经过非线性频率压缩后的系统函数用Ha(s1)，s1=jΩ1表示，有


s=jΩ=j2Ttan Ω1T2=2T1-e-s1T1+e-s1T（3.3.32）



再通过z=es1T转换到Z平面上，得到


s=2T1-z-11+z-1（3.3.33） 

z=2T+s2T-s（3.3.34）


由于从S平面到S1平面具有非线性频率压缩的功能，因此不可能产生频率混叠现象，这是双线性变换的最大优点。另外，从S1平面转换到Z平面仍然采用标准转换关系z=es1T，S1平面±π/T之间水平带的左半部分映射到Z平面单位圆内部，虚轴映射单位圆。这样，Ha(s)因果稳定，转换成H(z)也时因果稳定的。
由于s=jΩ，z=ejω，代入式（3.3.33），可得


jΩ=2T1-e-jω1+e-jω

Ω=2Ttan ω2 （3.3.35）


图3.3.12显示了这种频率变换关系。可见S平面上Ω与Z平面的ω成非线性正切关系，如图3.3.11所示。


图3.3.12双线性变换法的频率

变换关系

在ω=0附近接近线性关系； 当ω增加时，Ω增加得越来越快； 当ω趋近π时，Ω趋近于∞。正是由于这种非线性关系，消除了频率混叠现象。
显然，要想利用双线性变换法来设计数字滤波器，首先必须利用式（3.3.35）进行数字模拟频率的预畸变处理。
双线性变换法避免了脉冲响应不变法所带来的混叠问题，但是却引入了频率轴非线性压缩，只有当这种压缩在容许范围内进行补偿时，
如在滤波器具有近似理想的分段恒定幅度响应特性的情况下，使用双线性变换法设计离散时间滤波器才会有效。而且又无法避免这种频率轴的失真所引起的滤波器相位响应的畸变。例如，双线性变换法不能将连续时间微分器映射成离散时间微分器。
双线性变换法可由简单的代数式（3.3.33）将Ha(s)直接转换成H(z)，但当阶数稍高时，将H(z)整理成需要的形式，也不是一个简单的工作。为简化设计，已将模拟滤波器各系数和经双线性变换法得到的数字滤波器的各系数之间关系，列成表格供设计时使用。
设


Ha(s)=A0+A1s+A2s2+…+AkskB0+B1s+B2s2+…+Bksk

 H(z)=Ha(s)s=2T1-z-11+z-1

 H(z)=a0+a1z-1+a2z-2+…+akz-k1+b1z-1+b2z-2+…+bkz-k


系数Ak、Bk和ak、bk之间关系见表3.3.4。


表3.3.4系数关系表



k=1
AB0+B1C
a1(A0+A1C)/A
a0(A0-A1C)/A
b1(B0-B1C)/A

k=2
AB0+B1C+B2C2
a0(A0+A1C+A2C2)/A
a1(2A0-2A2C2)/A
a2(A0-A1C+A2C2)/A
b1(2B0-2B2C2)/A
b2(B0-B1C+B2C2)/A

k=3
AB0+B1C+B2C2+B3C3
a0(A0+A1C+A2C2+A3C3)/A
a1(3A0+A1C-A2C2-3A3C3)/A
a2(3A0-A1C-A2C2+3A3C3)/A
a3(A0-A1C+A2C2-A3C3)/A
b1(3B0+B1C-B2C2-3B3C3)/A
b2(3B0-B1C-B2C2+3B3C3)/A
b3(B0+B1C+B2C2+B3C3)/A

k=4
AB0+B1C+B2C2+B3C3+B4C4
a0(A0+A1C+A2C2+A3C3+A4C4)/A
a1(4A0+2A1C-2A3C3-4A4C4)/A
a2(6A0-2A2C2+6A4C4)/A
a3(4A0-2A1C+2A3C3-4A4C4)/A
a4(A0-A1C+A2C2-A3C3+A4C4)/A
b1(4B0+2B1C-2B3C3-4B4C4)/A
b2(6B0-2B2C2+6B4C4)/A
b3(4B0-2B1C+2B3C3-4B4C4)/A
b4(B0-B1C+B2C2-B3C3+B4C4)/A

续表


k=5
AB0+B1C+B2C2+B3C3+B4C4+B5C5
a0(A0+A1C+A2C2+A3C3+A4C4+A5C5)/A
a1(5A0+3A1C+A2C2-A3C3-3A4C4-5A5C5)/A
a2(10A0+2A1C-2A2C2-2A3C3+2A4C4+10A5C5)/A
a3(10A0-2A1C-2A2C2+2A3C3+2A4C4-10A5C5)/A
a4(5A0-3A1C+A2C2+A3C3-3A4C4+5A5C5)/A
a5(A0-A1C+A2C2-A3C3+A4C4-A5C5)/A
b1(5B0+3B1C+B2C2-B3C3-3B4C4-5B5C5)
b2(10B0+2B1C-2B2C2-2B3C3+2B4C4+10B5C5)
b3(10B0-2B1C-2B2C2+2B3C3+2B4C4-10B5C5)
b4(5B0-3B1C+B2C2+B3C3-3B4C4+5B5C5)
b5(B0-B1C+B2C2-B3C3+B4C4-B5C5)

［例3.3.4］用双线性不变法将图3.3.12所示的RC低通滤波器换成数字滤波器。


图3.3.13例3.3.4的RC

低通滤波器



解先按照图3.3.13所示，该滤波器的传输函数Ha(s)为


Ha(s)=aa+s，a=1RC


数字滤波器的系统函数H（z）为


H(z)=Ha(s)s=2T1-z-11+z-1
=a1(1+z-1)1+a2z-1

a1=aTaT+2，a2=aT-2aT+2


［例3.3.5］用双线性变换法设计一个Butterworth低通数字滤波器，频率f在500Hz处衰减3dB，在750Hz处衰减至少15dB，采样频率为2000Hz，确定系统函数H(z)。
解（1） 求数字指标ωp和ωs


ωp=Ω′pT=2π×5002000=π2

ωs=Ω′sT=2π×7502000=0.75π


（2） 求Ωp，Ωs
  利用频率预畸变公式有


Ωp=2Ttanωp2=2Ttanπ4=2T=Ωc

Ωs=2Ttanωs2=2Ttan0.75π2=1T×4.828


（3） 确定滤波器阶数



λsp=ΩsΩp=4.828×1T2T=2.414

 ksp=100.1αp-1100.1αs-1=100.1×3-1100.1×15-1=0.1803

N=-lgksplgλsp=-lg0.1803lg2.414=1.944，N=2 


（4） 确定系统函数


G(p)=1p2+2p+1，Ωp=Ωc=2T

p=sΩcs=2T
1-z-11+z-1=1Ωc×2T×1-z-11+z-1=1-z-11+z-1

H(z)=G(p)p=1-z-11+z-1=11-z-11+z-12+21-z-11+z-1+1

=(1+z-1)2(2+2)+(2-2)z-2



可以利用MATLAB来完成滤波器的设计，分别对比采用双线性法和脉冲响应不变法的区别。

%数字低通双线性变换与脉冲不变法程序

Wp=500*2*pi;Ws=150*2*pi;Ap=3;As=15;Fs=2000;

［n,Wn］=butterd(Wp,Ws,Ap,As,'s');  %选择Butterworth滤波器参数

［b,a］=butter(n,Wn,'s');%采用双线性变换法进行离散化处理

［bn1,an1］=bilinear(b,a,2000);%双线性变换

［H1,W］=freqz(bn1,an1);

%采用脉冲响应不变法进行离散化处理

［bn2,an2］=impinvar(b,a,2000);  %脉冲响应不变法

［H2,W］=freqz(bn2,an2);

plot(W,abs(H1),'-.',W,abs(H2),'-');

grid;xlabel('频率');ylabel('幅度响应');

legend('双线性变法','脉冲响应不变法');

结果如图3.3.14所示。


图3.3.14采用双线性法和脉冲响应不变法的区别

［例3.3.6］设计低通数字滤波器，要求通带内频率低于0.2πrad时，允许幅度误差在1dB之内，频率在0.3π～π之间的阻带衰减大于10dB，试用巴特沃斯型模拟滤波器，采用双线性法设计。
解ωp=0.2πrad,αp=1dB
ωs=0.3πrad,αs=10dB

（1） 频率预畸变


Ωp=2Ttanωp2=2Ttan0.1π=0.649/T(rad/s)

Ωs=2Ttanωs2=2Ttan0.15π=1.019/T(rad/s)

αs=10dB


（2） 确定滤波器阶数


ksp=100.1αp-1100.1αs-1=100.1-1101-1=0.1696

λsp=ΩsΩp=1.019×1T0.649T=1.5682

N=-lgksplgλsp=-lg0.1696lg1.5682=3.9435，取N=4


（3） 查表求归一化低通滤波器函数


G(p)=1p4+2.613p3+3.4142p2+2.6131p+1


（4） 求滤波器系统函数H(z)


Ωc=Ωp(100.1αp-1)-12N=0.649(100.1-1)-18/T=0.7743/T（rad/s）

s=2T1-z-11+z-1，p=sΩc=sT0.7743=2(1-z-1)0.7743(1+z-1)

H（z）=G（p）p=sΩc=H(z)

=0.8329×10-2+0.3331×10-1z-1+0.4997×10-1z-2+0.3331z-3+0.8329z-41-2.0872z-1+1.8948z-2-0.8119z-3+0.1375z-4


可以用MATLAB来验证该设计结果是否正确。 

%低通数字滤波器

b=［0.008329 0.03331 0.04997 0.3331 0.8329］;

a=［1 -2.0872 1.8948 -0.8119 0.1375］;

freqz(b,a,1000);

结果如图3.3.15所示。


图3.3.15例3.3.6图


［例3.3.7］设计一个数字高通滤波器，要求通带截止频率ωp=0.8πrad时，通带衰减不大于3dB，阻带截止频率ωs=0.5π，阻带衰减不小于10dB。试采用巴特沃斯型滤波器进行设计。
解（1） 已知数字高通滤波器指标： ωp=0.8πrad，αp=3dB，ωs=0.5πrad，αs=10dB。
  （2） 由于设计的是高通数字滤波器，所以采用双线性变换法，进行预畸变校正确定相应的模拟高通滤波器指标（为了计算方便，取T=2s）


Ωp=2Ttanωp2=tan0.4π=3.0777（rad/s）

Ωs=2Ttanωs2=tan0.25π=1（rad/s）


（3） 将高通滤波器指标转换成模拟低通指标。由于Ωp=Ωc


ηp=ΩpΩc=1，ηs=ΩsΩc=13.0777=0.3249

λp=1ηp=1，λs=1ηs=3.0777 


（4） 设计归一化低通G（p）


ksp=100.1αp-1100.1αs-1=100.3-110-1=0.333


λsp=λsλp=3.0777，N=-lgksplgλsp=0.97，取N=1


可得G(p)=1p+1
（5） 频率变换，求模拟高通Ha(s)


Ha(s)=G(p)p=Ωcs=ss+Ωc=ss+3.0777



（6） 用双线性变换法将Ha(s)转换成H(z)


H（z）=Ha(s)s=1-z-11+z-1=1-z-14.0777+2.0777z-1


上例利用MATLAB来进行设计： 

%数字高通滤波器设计程序

Wp=0.8;Ws=0.5;Ap=3;As=10;

［N,Wn］=butterd(Wp,Ws,Ap,As); %选择Butterworth滤波器参数

［B,A］=butter(N,Wn,'high');%设计数字高通滤波器

freqz(B,A);


结果如图3.3.16所示。


图3.3.16例3.3.7图



［例3.3.8］设计一个数字带通滤波器，通带范围为0.25πrad到0.45πrad，通带内最大衰减为3dB，0.15πrad以下和0.55πrad以上为阻带，阻带内最小衰减为18dB，采用巴特沃斯方法设计。
解（1） 确定带通滤波器技术指标


ωu=0.45π（rad），ωl=0.25π（rad）

ωs2=0.55π（rad），ωs1=0.15π（rad）


通带内最大衰减αp=3dB，阻带内最小衰减αs=15dB。
（2） 确定相应模拟滤波器技术指标。为计算简单，设T=2s。


Ωu=2Ttanωu2=tan0.225π=0.8541（rad/s）

Ωl=2Ttanωl2=tan0.125π=0.4142（rad/s）

Ωs2=2Ttanωs22=tan0.275π=1.1708（rad/s）

Ωs1=2Ttanωs12=tan0.075π=0.2401（rad/s）



通带中心频率


Ω0=ΩuΩl=0.5948（rad/s）


带宽


B=Ωu-Ωl=0.4399（rad/s）


将以上边界频率对B归一化，得到相应归一化带通边界频率


ηu=ΩuB=1.9416，ηl=ΩlB=0.9416

ηs2=Ωs2B=2.6615，ηs1=Ωs1B=0.5458

η0=ηuηl=1.3521


（3） 由归一化带通指标确定相应模拟归一化低通技术指标。
归一化阻带截止频率为


λs=η2s2-η20ηs2=1.9746,-λs=η2s1-η20ηs1=-2.8037


取λs=1.9746，归一化通带截止频率为λp=1，αp=3dB，αs=18dB。
（4） 设计模拟归一化低通G(p)


ksp=100.1αp-1100.1αs-1=100.3-1101.8-1=0.1266

λsp=λsλp=1.9746，N=-lgksplgλsp=-lg0.1266lg1.9746=3.04


取N=3，因为3.04很接近3，所以取N=3基本满足要求，且系统简单。查表可得归一化低通系统函数G(p)=1p3+2p2+2p+1。
（5） 频率变换，将G(p)转换成模拟带通Ha（s）


Ha(s)=G(p)p=s2+Ω20sB

=B3s3(s2+Ω20)3+(s2+Ω20)2sB+2(s2+Ω20)s2B2+s3B3

=0.085s3s6+0.8798s5+1.4484s4+0.7076s3+0.5124s2+0.1101s+0.0443 


（6） 用双线性变换公式将Ha(s)转换成


H(z)=Ha(s)s=2T1-z-11+z-1

=［0.0181+1.7764×10-15z-1-0.0543z-2-4.4409z-3+0.0543z-4-

2.7756×10-15z-5-0.0181z-6］［1-2.272z-1+3.5151z-2-3.2685z-3+

2.3129z-4-0.9628z-5+0.278z-6］-1




［例3.3.9］设计一个数字带阻滤波器，通带下限频率为ωl=0.19πrad，阻带下限截止频率ωs1=0.198πrad，阻带上截止频率ωs2=0.202πrad，通带上限频率ωu=0.21πrad，阻带最小衰减αs=13dB，ωl和ωu处衰减αp=3dB。采用巴特沃斯型。
解（1） 确定带阻滤波器技术指标


ωu=0.21πradωl=0.19πradαp=3dB

ωs2=0.202πradωs1=0.198πradαs=13dB


（2） 确定相应模拟带阻滤波器技术指标（只能用双线性变换法）。
设T=1，则有


Ωu=2tanωu2=0.685（rad/s）

Ωl=2tanωl2=0.615（rad/s）

Ωs2=2tanωs12=0.657（rad/s）

Ωs1=2tanωs22=0.643（rad/s）


阻带中心频率平方为


Ω20=ΩuΩl=0.421（rad/s）


阻带带宽为


B=Ωu-Ωl=0.07（rad/s）


将以上边界频率对B归一化


ηu=ΩuB=9.786，ηl=ΩlB=8.786

ηs2=Ωs2B=9.386，ηs1=Ωs1B=9.186

η20=ηuηl=85.98


（3） 归一化模拟低通滤波器的技术指标


λp=1,αp=3dB,αs=13dB



由于


λs=ηs2η2s2-η20=4.434,-λs=ηs1η2s1-η20=-5.75


因此，取λs=4.434。
（4） 设计巴特沃思模拟低通滤波器


ksp=100.1αp-1100.1αs-1=0.229

λsp=λsλp=4.434，N=-lgksplgλsp=0.99


取N=1，得归一化巴特沃思模拟低通滤波器传递函数 


  G(p)=1p+1


（5） 将G(p)转换成模拟阻带Ha(s)


 Ha(s)=G(p)p=sBs2+Ω20


（6） 将Ha(s)通过双线性变换，得到数字阻带滤波器H(z)


s=2·1-z-11+z-1

H(z)=H(s)s=21-z-11+z-1


建立p与z的关系以简化运算


p=sBs2+Ω20s=2·1-z-11+z-1=2（1-z2）B4（1-z-1）2+Ω20（1+z-1）2


所以


H（z）=G（p）p=2（1-z2）B4（1-z-1）2+Ω20（1+z-1）2=0.969（1-1.619z-1+z-2）1-1.569z-1+0.939z-2


下面的例题中，直接使用MATLAB来设计。
［例3.3.10］设Fs=200Hz，试设计一低通滤波器，要求： 通带截止频率为30Hz，通带波纹1dB，阻带截止频率为50Hz，阻带衰减50dB。
MATLAB程序实现如下： 

%低通滤波器

Wp=30;Ws=50;Ap=1;As=50;Fs=200;

［n,Wn］=buttord(Wp/(Fs/2),Ws/(Fs/2),Ap,As);  %获得最小阶数n

nb=n;na=n+10;

figure(1);

Maxflat(nb,na,Wp/(Fs/2),'plots');  %获得幅值响应，群延迟及零极点图

［b,a］=maxflat(nb,na,Wn);

figure(2);

freqz(b,a);

结果如图3.3.17所示。


图3.3.17例3.3.10图



图3.3.17（续）

［例3.3.11］用双线性变化法设计一个带通数字滤波器，通带频率为20～30Hz，在通带内的最大衰减为0.5dB，在阻带内的最小衰减为50dB，采样频率为150Hz，其MATLAB程序如下： 

% 带通椭圆数字滤波器设计

Wp1=20*2*pi; Wp2=30*2*pi;Ap=0.5;As=50;Fs=150;

Bw=Wp2-Wp1;

W0=sqrt(Wp2*Wp1);

［z,p,k］=ellipap(7,Ap,As); %采用7阶椭圆数字滤波器来处理

［A,B,C,D］=zp2ss(z,p,k);

［At,Bt,Ct,Dt］=lp2bp(A,B,C,D,W0,Bw);  %模拟频率变换

［At1,Bt1,Ct1,Dt1］=bilinear(At,Bt,Ct,Dt,Fs);  %双线性变换离散化

［b,a］=ss2tf(At1,Bt1,Ct1,Dt1);  %状态变换

freqz(b,a);

结果如图3.3.18所示。


图3.3.18例3.3.11图

［例3.3.12］设数字采样频率为1000Hz，设计一带阻滤波器，要求阻带范围为50～300Hz，通带上限大于250Hz，下限小于100Hz，通带内波纹小于1dB，阻带要求50dB，要求利用最小的阶来实现。
MATLAB程序实现如下： 

%带阻滤波器

Wp1=100;Wp2=250;Ws1=50;Ws2=300;Ap=1;As=50;Fs=1000;  %指标描述

Wp=［Wp1 Wp2］;Ws=［Ws1  Ws2］;

［n,Wn］=butterd(Wp/(Fs/2), Ws/(Fs/2),Ap,As);  %获得最小阶数n

［b,a］=butter(n,Wn,'stop');  %直接设计带阻滤波器

freqz(b,a,512,1000);

运行结果如图3.3.19所示。


图3.3.19例3.3.12图


［例3.3.13］基于MATLAB的语音信号滤波器的设计与实现。
语音很容易受到噪声的污染。语音信号处理的主要目的就是削弱信号中的多余内容，滤出混杂的噪声和干扰。在本例中，通过导入一段语音信号，并添加8kHz的高频余弦噪声来模拟含噪信号。通过分析含噪信号的频谱，设计了合适的低通滤波器，以达到去噪的目的。MATLAB代码如下：

 ［x1,fs］=audioread('voice.wav'); %读取语音信号的数据，赋给变量x1, fs=48000

N=length(x1);

k=0:N-1;

y1=fft(x1,N);  %对信号做L个点FFT变换

f=fs*(0:N-1)/N;



%给原始的语音信号加上一个高频余弦噪声，频率为8kHz。画出加噪后的语音信号时域和频谱图，与

%原始信号对比，可以很明显地看出区别

t=［0:1/fs:(N-1)/fs］;        %将所加噪声信号的点数调整到与原始信号相同，构造采样时间点(模拟时间)

Au=0.02;                    %噪声幅度

d=［Au*cos(2*pi*8000*t); Au*cos(2*pi*8000*t)］'; %噪声为8kHz的余弦信号(模拟时间)

x2=x1+d;



figure(1);

subplot(2,1,1);

plot(x1);        %做原始语音信号的时域图形

title('原始语音信号');xlabel('采样点 n');ylabel('幅值');

subplot(2,1,2);

plot(x2);        %做加噪后语音信号的时域图形

title('加噪后的信号');xlabel('采样点 n');ylabel('幅值');



y2=fft(x2,N);     %对信号做L个点FFT变换

figure(2);

subplot(2,1,1);

plot(f(1:20000),abs(y1(1:20000)));

title('原始语音信号频谱');

xlabel('频率/Hz');ylabel('幅值');

subplot(2,1,2);

plot(f(1:20000),abs(y2(1:20000)));

title('加噪声语音信号频谱');

xlabel('频率/Hz');ylabel('幅值');



wp=2*pi*4000;  %通带边界角频率

ws=2*pi*5000;  %阻带边界角频率

Rp=1;  %通带最大衰减

Rs=15;    %阻带最小衰减

［NN,Wn］=butterd(wp,ws,Rp,Rs,'s');   %选择滤波器的最小阶数

［Z,P,K］=buttap(NN);   %创建butterworth模拟滤波器

［Bap,Aap］=zp2tf(Z,P,K);

［b,a］=lp2lp(Bap,Aap,Wn);

［bz,az］=bilinear(b,a,fs);    %用双线性变换法实现模拟滤波器到数字滤波器的转换



［H,W］=freqz(bz,az);   %绘制频率响应曲线

figure(3);

plot(W*fs/(2*pi),abs(H));grid;

xlabel('频率/Hz');ylabel('频率响应幅度');

title('Butterworth');



x3=filter(bz,az,x2); %滤波

figure(4);

subplot(2,1,1);

plot(t,x2);  %画出滤波前的时域图

title('滤波前的时域波形');

xlabel('采样时刻 t');ylabel('幅值');

subplot(2,1,2);

plot(t,x3);   %画出滤波后的时域图

title('滤波后的时域波形');

xlabel('采样时刻 t');ylabel('幅值');



y2=fft(x2,N);  %对信号做L个点FFT变换

figure(5);  

subplot(2,1,1);

plot(f(1:20000),abs(y2(1:20000)));

title('加躁语音信号频谱');

xlabel('频率/Hz');ylabel('幅值');

y2=fft(x3,N);  %对信号做L个点FFT变换

subplot(2,1,2);

plot(f(1:20000),abs(y2(1:20000)));

title('去噪语音信号频谱');

xlabel('频率/Hz');ylabel('幅值')

如图3.3.20所示为原始语音信号与加噪后的语音信号时域波形。


图3.3.20原始语音信号与加噪后的语音信号时域波形


如图3.3.21所示为原始语音信号频谱与加噪语音信号频谱。
如图3.3.22所示为巴特沃斯低通滤波器的幅频特性。
如图3.3.23所示为滤波前和滤波后的时域信号。
如图3.3.24所示为滤波前和滤波后的频谱。



图3.3.21原始语音信号频谱与加噪语音信号频谱




图3.3.22巴特沃斯低通滤波器的幅频特性





图3.3.23滤波前和滤波后的时域信号




图3.3.24滤波前和滤波后的频谱


3.4FIR型滤波器的设计
3.4.1线性相位FIR滤波器及其特点


IIR数字滤波器可用较少的阶数达到所要求的幅度特性，且实现时所需的运算次数及存储单元都很少，十分适合对于相位特性没有严格要求的场合，如果对其相位特性有要求，就必须加上相位校正网络，因为所设计的滤波器相位特性一般是非线性的，因此使滤波器变得复杂，成本高。
而FIR型滤波器的系统函数H(z)为


H(z)=h(0)+h(1)z-1+…+h(N-1)z-(N-1)=∑N-1n=0h(n)z-n


FIR系统没有非零极点，因而它不像IIR系统那样容易取得比较好的通带与阻带衰减特性。FIR数字滤波器的阶数(N-1）一般较大，但由于其优越的稳定性，线性相位特性以及良好的系数量化特性，使得FIR数字滤波器在众多领域拥有非常广泛的应用。

线性相位FIR系统的单位脉冲响应h(n)应满足h(n)为实序列，且满足h(n)=±h(N-1-n)。根据h(n)为偶对称或奇对称以及长度N的奇偶，线性相位FIR滤波系统分为4种情况，其幅度和相位特性各有不同。令H(ejω)=Hg(ω)ejθ(ω)，线性相位FIR滤波器的幅度特性与相位特性如表3.4.1所示。
从表3.4.1可以知道，对脉冲响应偶对称，N为奇数的情况1，由于cosωn关于ω=0，π,2π偶对称，Hg(ω)对这些频率呈偶对称。可以实现各种（低通、高通、带通、带阻）滤波器，对脉冲响应偶对称，N为偶数的情况2，由于cos［(n-1/2)ω］对ω=π奇对称，Hg(ω)对ω=π也呈奇对称，且由于ω=π时，cos［(n-1/2)π］=0，Hg(π)=0，因此，不能用这种情况设计ω=π处，Hg(π)≠0的滤波器，例如高通、带阻滤波器。


表3.4.1线性相位FIR滤波器的幅度特性与相位特性一览表



偶对称单位脉冲响应h(n)=h(N-1-n)

情

况

1





相位响应

θ(ω)=-N-12ω



N为奇数



Hg(ω)=∑(N-1)/2n=0a(n)cosnω



情

况

2

N为偶数



Hg(ω)

=∑N/2n=0b(n)cos ωn-12


续表


奇对称单位脉冲响应h(n)=-h(N-1-n)



情

况

3

情

况

4


相位响应

θ(ω)=-N-12ω-π2





N为奇数



Hg(ω)=∑(N-1)/2n=0c(n)sinnω



N为偶数



Hg(ω)=

∑N2n=0dnsin n-12ω




对脉冲响应奇对称，N为奇数的情况3，当ω=0，π,2π时，Hg(ω)=0，且Hg(ω)对这些频率呈奇对称，故它不能用于低通、高通和带阻滤波器设计，只适合用于实现带通滤波器。
对脉冲响应奇对称，N为偶数的情况4，Hg(ω)关于ω=0和2π两点奇对称，关于ω=π偶对称，因此，它不能用来实现低通和带阻滤波器，可以用来实现高通和带通滤波器。
4种FIR数字滤波器的相位特性只取决于h(n)的对称性而与h(n)的值无关，但其幅度特性则取决于h(n)。在设计线性相位FIR滤波器时，在保证h(n)的对称性的条件下，只要考虑幅度尽可能逼近即可。
FIR数字滤波器的设计方法主要包括3种，窗函数设计法、频率采样设计法以及最优化设计方法。本节主要讨论前面两种方法，对于第3种方法，考虑到其重要性，也将略作介绍并利用MATLAB加以解释。
3.4.2利用窗函数法设计FIR滤波器
  设计的滤波器传输函数为Hd(ejω)，其对应的单位脉冲响应为hd(n)，即


Hd(ejω)=∑∞n=-∞hd(n)e-jωn

hd(n)=12π∫π-πHd(ejω)ejωndω




图3.4.1低通滤波器的幅频特性

然而，一般的hd(n)是无限时宽的，且是非因果的，所以是无法实现的。
例如，假定理想低通滤波器的频率特性为Hd(ejω)，幅频特性Hd(ejω)=1，相频特性θ(ω)=0，如图3.4.1所示，那么与其对应的单位脉冲响应为hd(n)


hd(n)=12π∫π-πHd(ejω)ejωndω=12π∫ωc-ωcejωndω

=sin(ωcn)πn



hd(n)是以hd(0)为对称的sinc函数，这样的系统是非因果的，因此是物理不可实现的。为了保证所设计的滤波器具有线性相位，首先将hd(n)右移α=N-12个采样点，得到 


h′d(n)=hd(n-α)  （3.4.1）


则h′d(n)关于N-12对称。其次，为使系统为因果可实现系统，对h′d(n)截取N个点，得


h(n)=h′d(n)wN(n)=hd(n-α)wN(n) （3.4.2）


式中，wN(n)为有限时宽的窗序列。最简单的窗序列为矩形窗序列RN(n)，得

RN(n)=1,0≤n≤N-1



0,其他

这就是窗函数滤波器设计的基本思想。
用窗函数对h′d(n)进行截断，得到有限长序列h(n)，并以h(n)去代替hd(n)，肯定会引起误差，这就是截断效应。对于矩形窗，对式（3.4.2）进行傅里叶变换，根据复卷积定理可得 


H(ejω)=12π∫π-πH′d(ejω)RN(ej(ω-θ))dθ（3.4.3）


式中，H′d(ejω)和RN(ejω)分别是h′d(n)和RN(n)的傅里叶变换，即


RN(ejω)=∑N-1n=0RN(n)e-jωn=∑N-1n=0e-jωn

=e-j12（N-1）ωsin(ωN/2)sin(ω/2)=RN(ω)e-jωα （3.4.4）


式中

RN(ω)=sin(ωN/2)sin(ω/2)（3.4.5）

α=N-12

RN(ω)称为矩形窗的幅度函数。
由于H′d(ejω)=Hd(ejω)e-jωα，Hd(ejω)=Hd(ω)，将式（3.4.4）代入式（3.4.3）中，有


H(ejω)=12π∫π-πHd(ejθ)e-jθαRN(ω-θ)e-j(ω-θ)αdθ

=e-jωα·12π∫π-πHd(θ)RN(ω-θ)dθ


因此，H(ejω)可写成


H(ejω)=H(ω)e-jωα（3.4.6）


式中


H(ω)=12π∫π-πHd(θ)RN(ω-θ)dθ（3.4.7）


式（3.4.6）和式（3.4.7）表明，H(ejω)的相位函数θ(ω)=-αω是线性的，幅度函数H(ω)则为Hd(ω)与RN(ω)的卷积。
图3.4.2表示了Hd(ω)与RN(ω)的卷积过程，其中图3.4.2（e）为卷积形成的波形，它表明对hd(n)加矩形窗处理后，H(ω)和原理想低通Hd(ω)有着明显的区别。区别之一是，在理想特性不连续点ω=ωc附近形成过渡带，过渡带的宽度近似为RN(ω)的值为宽度，即4π/N。区别之二是，在过渡带两侧形成持续时间很长，逐渐衰减的波纹，即通带内增加了波动，在ω=ωc-2π/N处达到最大正峰。而在阻带内产生余振，并在ω=ωc+2π/N处达到最大负峰。显然这种波动直接取决于窗函数的幅度谱。以上两点就是用窗函数直接截断hd(n)而引起的截断相应在频域中的反应，通常称之为吉布斯（Gibbs）效应。


图3.4.2矩形窗对理想低通滤波器的幅度特性的影响



吉布斯效应直接影响滤波器的性能，因为通带内的波动会影响滤波器通带中的平稳性，阻带内的波动则影响阻带最小衰减，因此，减少吉布斯效应也是FIR数字滤波器设计的关键之一。
由式（3.4.5），有


RN(ω)=sin(ωN/2)sin(ω/2)≈sin(ωN/2)ω/2≈Nsinxx（3.4.8）


可见，随着N的加大，虽然主瓣幅度随之加高，但由于旁瓣也同时加高，主、旁瓣相对比例仍然保持不变。同时，N的加大使起伏振荡变密，但不能改善波动幅度，其最大正、负峰总是8.95％。N加大带来的最大好处是H(ω)的过渡带(4π/N)变窄。因此，加大N并不能有效地减少吉布斯效应。
综上所述，调整窗口长度N只能有效地控制过渡带的宽度，但不能减少带内波动和加大阻带衰减，特别是-20lg(8.95%)=21dB的阻带衰减通常是不能满足工程要求的。换句话说，通带、阻带技术指标仅与窗函数的形状有关。为此必须寻找其他形状的窗函数，是其谱函数的主瓣包含更多的能量，旁瓣幅度就相应地减少，通带、阻带波动也因此减少，从而加大阻带衰减，但这一切总是以加宽过渡带为代价的。
下面介绍几种常用的窗函数。设 


h(n)=hd(n)w(n)


式中w(n)表示窗函数。
1. 矩形窗


wR(n)=RN(n)


其频率相应为


WR(ejω)=e-j12（N-1）ωsin(ωN/2)sin(ω/2)（3.4.9)


WR(ejω)主瓣宽度为4π/N，第一副瓣比主瓣低13dB。
2. 三角形窗


wBr(n)=2nN-1,0≤n≤12(N-1)

2-2nN-1,12(N-1)≤n≤N-1


其频率相应为


WBr(ejω)=2NsinN4ωsinω22e-jω+N-12ω（3.4.10）


其主瓣宽度为8π/N，第一副瓣比主瓣低26dB。
3. 升余弦窗——汉宁窗


wHn(n)=0.51-cos2πnN-1×RN(n)


RN(n)的傅里叶变换为


WR(ejω)=e-j12(N-1)WR(ω)


WHn(n)的傅里叶变换为


WHn(ejω)=0.5WR(ω)+0.25WRω-2πN-1+WRω+2πN-1e-j12(N-1)ω

=WHn(ω)e-j12（N-1）ω（3.4.11）


汉宁的幅度函数WHn(ω)由三部分相加而成，使能量更集中在主瓣中，如图3.4.3所示,但代价使主瓣宽度加宽到8π/N。


图3.4.3汉宁窗的幅度特性


4. 改进的升余弦窗——汉明窗


wHm(n)=0.54-0.46cos2πnN-1×RN(n)


其频域函数WHm(ejω)为


WHm(ejω)=0.54WR(ejω)-0.23WRejω-2πN-1-0.23WRejω+2πN-1 （3.4.12）


其幅度函数WHm(ω)为


 WHm(ω)≈0.54WR(ejω)+0.23WRω-2πN-1+0.23WRω+2πN-1


汉明窗的能量更加集中在主瓣中，主瓣的能量占99.96％，第一旁瓣的峰值比主瓣小40dB，但主瓣宽度和汉宁窗相同，仍为8π/N。
可以利用MATLAB来演示这几种常见的窗函数及其频谱特性。

%窗函数及其频谱特性演示，取N=51点

N=51;n=［0:1:(N-1)］;

%矩形窗

W_box=boxcar(N);

［Hbox,W］=freqz(W_box,1);

subplot(421);

stem(n,W_box);xlabel('n');ylabel('矩形窗');

subplot(422);

plot(W/pi,20*log10(abs(Hbox)/abs(Hbox(1))));ylabel('矩形窗频谱');

%三角窗

W_tri=triang(N);

［Htri,W］=freqz(W_tri,1);

subplot(423);

stem(n,W_tri); ;xlabel('n');ylabel('三角窗');

subplot(424);

plot(W/pi,20*log10(abs(Htri)/abs(Htri(1))));ylabel('三角窗频谱');

%汉宁窗

W_han=hanning(N);

［Hhan,W］=freqz(W_han,1);

subplot(425);

stem(n,W_han); ;xlabel('n');ylabel('汉宁窗');

subplot(426);

plot(W/pi,20*log10(abs(Hhan)/abs(Hhan(1))));ylabel('汉宁窗频谱');

%汉明窗

W_ham=hamming(N);

［Hham,W］=freqz(W_ham,1);

subplot(427);

stem(n,W_ham); ;xlabel('n');ylabel('汉明窗');

subplot(428);

plot(W/pi,20*log10(abs(Hham)/abs(Hham(1))));ylabel('汉明窗频谱');

结果如图3.4.4所示。


图3.4.4几种常见的窗函数及其频谱特性

5. 凯泽（Kaiser）窗
这是一种最有用和最优的窗函数，它是在给定阻带衰减下给出一种主瓣宽度意义上的最优结果，这里面就包含着最为陡峭的过渡带，表达式为


w(n)=I0β1-1-2nN-12I0［β］，0≤n≤N-1（3.4.13）


式中，I0(·)是修正的零阶贝塞尔函数I0(x)=1+∑∞k=11k!x2k2，而β是一个控制阻带衰减的重要参数。一般说来β加大，主瓣加宽，旁瓣幅度减小。β=0相对于矩形窗。β的典型取值为4<β<9。凯泽窗是一种适应性较强的窗函数，其设计经验公式如下。
已知通带截止频率ωp，通带衰减αp，阻带截止频率ωs及阻带衰减指标αs，则
标准过渡带带宽=Δf=
ωs-ωp2π
滤波器的阶数N≈αs-7.9514.3612f+1
控制参数β=0.1102(αs-8.7),αs≥50


0.5842(αs-21)0.4+0.07886(αs-21),21<αs<50
下面用MATLAB来分别演示N与β的作用： 

% 参数B对凯泽窗的影响，N=51，B1=4.5，B2=6.75

W1=kaiser(51,4.5);

［H1,W］=freqz(W1,1);

figure(1);

plot(W/pi,20*log10(abs(H1)/abs(H1(1))));

hold on;

W2=kaiser(51,6.75);

［H2,W］=freqz(W2,1);

plot(W/pi,20*log10(abs(H2)/abs(H2(1))),'r--');

xlabel('频率');ylabel('幅度响应');

title('参数B对凯泽窗的影响');legend('B=4.5','B=6.75');hold off;



%N对凯泽窗的影响，B=5.15，N1=51，N2=81

W3=Kaiser(51,5.15);

figure(2);

［H3,W］=freqz(W3,1);

plot(W/pi,20*log10(abs(H3)/abs(H3(1))));

hold on;

W4=kaiser(81,5.15);

［H4,W］=freqz(W4,1);

plot(W/pi,20*log10(abs(H4)/abs(H4(1))),'r--');

xlabel('频率');ylabel('幅度响应');

title(' N的取值对凯泽窗的影响');legend('N=51','N=81');hold off;

结果如图3.4.5所示。



图3.4.5不同β、N的凯萨窗频谱



表3.4.2列出了几种窗函数的基本参数。


表3.4.2窗函数的基本参数



窗函数
旁瓣峰值幅度/dB
过 渡 带 宽
阻带最小衰减/dB

矩形窗
-13
4π/N
-21
三角窗
-25
8π/N
-25
汉宁窗
-31
8π/N
-44
汉明窗
-41
8π/N
-53
凯泽窗

β=7.865
-57
10π/N
-80

用窗函数设计FIR数字滤波器的步骤如下： 
（1） 根据技术要求确定滤波器的频响特性Hd（ejω），确定其对应的单位脉冲响应hd（n）。


hd(n)=12π∫π-πHd(ejω)ejωndω （3.4.14）


如果H(ejω)比较复杂，或者不能用封闭公式表示时，则无法根据上式求出hd(n)。此时，可以对H(ejω)在ω=0到ω=2π范围内采样M点，得到其采样值为H(ej2πk/M)，k=0,1,…,M-1，并用2π/M代替式（3.4.14）的dω，则式（3.4.14）可以近似表达为


hd(n)=1M∑M-1k=0Hd(ej2πk/M)ej2πkn/M （3.4.15）


根据频域采样定理，可得


hM(n)=∑∞r=-∞hd(n+rM)


当M趋向于∞时，可以使窗口内hM(n)有效逼近于hd(n)，实际计算式（3.4.15）可以用Hd(ejω)的M点采样值，进行M点IDFT得到。
如果给定的设计要求为通带、阻带衰减和边界频率时，可选用理想滤波器作为逼近函数，进而对理想滤波器的特性做傅里叶反变换，求出hd(n)。
（2） 根据对过渡带及阻带衰减指标的要求，选择窗函数形式，并估计窗口长度N。
当待求滤波器的过渡带Δω近似等于窗函数主瓣宽度时，过渡带Δω与窗口长度N成反比，即N≈A/Δω，A决定于窗口形式，例如，矩形窗A=4π，汉明窗A=8π等，A参数选择参考表3.4.2。按照过渡带及阻带衰减情况，选择窗函数形式。原则在保证阻带衰减满足要求的情况下，尽量选择主瓣窄的窗函数。
（3） 计算滤波器的单位取样响应h（n）


h(n)=hd(n)w(n)


式中w(n)是上面选择好的窗函数。
如果要求线性相位，则要求hd(n)和w(n)均对(N-1)/2对称。如果要求h(n)对(N-1)/2奇对称，只要保证hd(n)对(N-1)/2奇对称就可以了。
（4） 计算滤波器频率响应Hd（ejω）


Hd(ejω)=∑N-1n=0h(n)e-jωn


上式计算必要时可用FFT算法验证其是否符合指标要求，如不满足要求，根据具体情况予以修正，重复步骤（2）~步骤（4）直至满足要求。
［例3.4.1］试用窗函数法设计一个FIR低通滤波器。已知


Hd(ejω)=e-jωα,0≤|ω|≤π2

0,π2<|ω|≤π，α=6


（1） 求h(n)的长度N； 
（2） 在矩形窗的条件下，求出h(n)的表达式； 
（3） 写出过渡带宽Δω。
解（1） 由α=N-12，可知N=2α+1=13。
(2) N为奇数，且是低通滤波器，故属于第一类广义线性相位FIR滤波器，截止频率ωc=π/2。



hd(n)=12π∫π-πHd(ejω)ejωndω=12π∫ωc-ωce-jωαejωndω

=sin［ωc(n-α)］π(n-α)=sinπ2(n-6)π(n-6)

h(n)=hd(n)RN(n)=sinπ2(n-6)π(n-6)·R13(n)，0≤n≤12


（3） 过渡带宽： Δω=4π/N=4π/13。
上例的MATLAB的演示如下： 

%矩形窗低通FIR设计演示

N=13;a=(N-1)/2;Wc=pi/2;

n=［0:1:(N-1)］;

m=n-a+eps;  %避免被零除

hd=sin(Wc*m)./(pi*m);

［H,W］=freqz(hd,1);

plot(W/pi,20*log10(abs(H)/abs(H (1))));

xlabel('频率');ylabel('对数幅度响应');

title('FIR加矩形窗时的频谱图');

运行结果如图3.4.6所示。


图3.4.6例3.4.1图


［例3.4.2］利用矩形窗、汉宁窗、汉明窗设计线性相位FIR低通滤波器，要求通带截止频率ωc=π/4rad，N=21。求出分别对应的单位脉冲响应，并进行比较。
解（1） 确定逼近滤波器传输函数Hd(ejω)


Hd(ejω)=

e-jωα,0≤|ω|≤π4,

0,π4<|ω|,
其中α=(N-1)/2=10


（2） 求hd(n)


hd(n)=sin［ωc(n-α)］π(n-α)=sinπ4(n-10)π(n-10)


（3） 加窗得到FIR滤波器单位脉冲响应h(n)。
矩形窗： wR(n)=RN(n)


hR(n)=hd(n)wR(n)=sinπ4(n-10)π(n-10)×R21(n)


汉宁窗


wHn(n)=0.51-cos2πnN-1×RN(n),N=21

hHn(n)=hd(n)wHn(n)=sinπ4(n-10)2π(n-10)1-cos2πn20×R21(n)


汉明窗


  wHm(n)=0.54-0.46cos2πnN-1×RN(n)

  hHm(n)=hd(n)wHm(n)

=sinπ4(n-10)π(n-10)0.54-0.46cos2πn20×R21(n)


矩形窗对应的过渡带最窄，但阻带最小衰减只有21dB，汉明窗对应的阻带衰减最大（大于100dB），但过渡带最宽。
本例题的MATLAB演示如下： 

%利用矩形窗、汉宁窗和汉明窗设计线性相位FIR低通滤波器

N=21;a=(N-1)/2;Wc=pi/4;

n=［0:1:(N-1)］;

m=n-a+eps;  %避免被零除

hd=sin(Wc*m)./(pi*m); 

%加矩形窗

［H1,W］=freqz(hd,1);

figure(1);

subplot(211);

stem(n,hd);title('实际脉冲响应');

xlabel('n');ylabel('h(n)');

subplot(212);

plot(W/pi,20*log10(abs(H1)/max (H1)));

xlabel('频率');ylabel('对数幅度响应');

title('加矩形窗时的频谱图');

%汉宁窗

W_han=(hanning(N))';

h2=hd.*W_han;

［H2,W］=freqz(h2,1);

figure(2);

subplot(211);

stem(n,h2);title('汉宁窗实际脉冲响应');

xlabel('n');ylabel('h(n)');

subplot(212);

plot(W/pi,20*log10(abs(H2)/max (H2)));

xlabel('频率');ylabel('对数幅度响应');

title('加汉宁窗时的频谱图');

%汉明窗

W_ham=(hamming(N))';

h3=hd.*W_ham;

［H3,W］=freqz(h3,1);

figure(3);

subplot(211);

stem(n,h3);title('汉明窗实际脉冲响应');

xlabel('n');ylabel('h(n)');

subplot(212);

plot(W/pi,20*log10(abs(H3)/max (H3)));

xlabel('频率');ylabel('对数幅度响应');

title('加汉明窗时的频谱图');

运行结果如图3.4.7所示。


图3.4.7例3.4.2图





图3.4.7（续）


［例3.4.3］用窗函数法设计第一类线性相位FIR高通数字滤波器，3dB截止频率ωc=3π4rad，阻带最小衰减αs=50dB，过渡带宽度Δω=π/16。
解根据设计要求，N必须为奇数（情况1可以设计任何滤波特性）。
（1） 确定逼近理想高通频响函数Hd(ejω)


Hd(ejω)=e-jωα,ωc<|ω|≤π

0,0≤|ω|≤ωc


（2） 求hd(n)


hd(n)=12π∫π-πHd(ejω)ejωndω=12π∫ωc-πe-jωαejωndω+∫πωce-jωαejωndω

=sin ［π(n-α)］-sin ［ωc(n-α)］π(n-α)


式中α=N-12
（3） 选择窗口函数，估算窗函数长度N，根据阻带最小衰减αs=50dB，查表，选择汉明窗，其过渡宽度为8π/N，本题要求过渡带宽度Δω=π/16。由π/16=8π/N，可知N=128，又N必须为奇数，故取N=129。
（4） 加窗计算，h（n）=hd（n）w（n），代入N=129，有

α=N-12=64，ωc=3π/4

 h(n)=sinπ(n-64)-sin3π4(n-64)π(n-64)×0.54-0.46cos2πn128×R129(n)



%汉明窗FIR高通演示

N=129;a=(N-1)/2;Wc=3*pi/4;

n=［0:1:(N-1)］;

m=n-a+eps;  %避免被零除

hd=(sin(pi*m)-sin(Wc*m))./(pi*m);

W_ham=(hamming(N))';

h=hd.*W_ham;

［H,W］=freqz(h,1);

plot(W/pi,20*log10(abs(H)/max(abs(H ))));

xlabel('频率');ylabel('对数幅度响应');

title('FIR高通加汉明窗时的频谱图');

运行结果如图3.4.8所示。


图3.4.8例3.4.3图


［例3.4.4］用长度为41的，β=6.15的凯泽（Kaiser）窗来设计一个数字微分器，理想数字微分器的频率响应为


Hd(ejω)=-jω,0≤ω<π


jω,-π<|ω|<0


解具有线性相位的理想数字微分器的脉冲响应为


hd(n)=12π∫π-πHd(ejω)ejωndω

=cos［π(n-α)］(n-α),n≠α

0,n=α


用MATLAB来完成该设计。

%数字微分器设计

N=41; n=［0:1:(N-1)］;a=(N-1)/2;

hd=(cos(pi*(n-a)))/(n-a);

hd(a+1)=0;

W_kai=(kaiser(41,6.15))';

h=hd.*W_kai;

［H,W］=freqz(h,1);

subplot(211);

stem(n,h);title('数字微分器脉冲响应');

xlabel('n');ylabel('h(n)');

subplot(212);

plot(W/pi,abs(H));

title('数字微分器幅频特性');

xlabel('频率');ylabel('对数幅度响应');

结果如图3.4.9所示。


图3.4.9例3.4.4图


［例3.4.5］带通滤波器的技术指标为： 下阻带边缘ω1s=0.2π，αs=60dB； 下通带边缘ω1p=0.35π，αp=1dB； 上通带边缘ωcp=0.65π，αp=1dB； 上阻带边缘ωcs=0.8π，αs=60dB。试采用合适的窗结构进行设计。
解（1） 显然

Δω1=ω1p-ω1s=0.15π
Δω2=ωcs-ωcp=0.15π

考虑到αs=60dB，可采用凯泽窗进行设计，本题选用Blackman窗。关于Blackman窗，请读者参考相关书籍。
（2） 理想带通数字滤波器的脉冲响应为


hd(n)=sin［ωh(n-α)］-sin［ωl(n-α)］π(n-α)，α=N-12


式中ωh是高通截止频率； ωl是低阻截止频率。可以看出它是由两个低通滤波器相减而得到的。直接用MATLAB来设计，程序如下： 

%带通滤波器

Ws1=0.2*pi;Wp1=0.35*pi;

Wp2=0.65*pi;Ws2=0.8*pi;As=60;

Width=min((Wp1-Ws1),(Ws2-Wp2));  %取过渡带宽

N=ceil(11*pi/Width)+1;  %由带宽公式计算阶数N

n=［0:1:(N-1)］;a=(N-1)/2;

m=n-a+eps;

Wl=( Wp1+Ws1)/2;Wh=(Ws2+Wp2)/2;

hd=(sin(Wh*m)-sin(Wl*m))./(pi*m);

W_bla=(blackman(N))';

h=hd.*W_bla;

［H,W］=freqz(h,1);

subplot(211);

stem(n,h);title('脉冲响应');

xlabel('n');ylabel('h(n)');

subplot(212);

plot(W/pi,20*log10(abs(H)/max (H)));

title('滤波器频响特性');

xlabel('频率');ylabel('对数幅度响应');

运行结果如图3.4.10所示。


图3.4.10例3.4.5图


3.4.3利用频率采样法设计FIR滤波器
FIR数字滤波器的窗函数设计方法的基本思想是使所设计的FIR数字滤波器的单位脉冲响应h(n)逼近所需的hd(n)，因此，它是一种时域设计方法。这种方法的缺点是通带和阻带的截止频率不易控制，且对于复杂频响特性滤波器来说，难以得到hd(n)的闭合表达式。但在工程实际上，更多的是给定频域上的指标，因此采用频域设计法更为方便和直接。频率采样法就是一种频域设计方法，它的基本思想是使所设计的FIR数字滤波器的频率特性
在某些离散的频率点上的值，准确地等于所需的滤波器在这些频率点处的值，在其他频率处的特征按照一定的优化设计则有较好的逼近。因此，频率采样设计的基本设计流程是： 
确定频率特性指标Hd(ejω)频率采样Hd(ej2πk/N)=Hd(k)IDFTh(n)ZT H(z)。
如待设计的滤波器的传输函数用Hd(ejω)表示，对它在ω在0～2π之间等间隔采样N点，得到Hd(k)


Hd（ejω）|ω=
2πNk，k=0,1,…,N-1


再对N点Hd(k)进行IDFT，得到h(n)


h(n)=1N∑N-1k=0Hd(k)ej2πkn/N, n=0,1,…,N-1（3.4.16）


式中，h(n)作为所设计的滤波器的单位脉冲响应，其系统函数H(z)


H(z)=∑N-1n=0h(n)z-n（3.4.17）


以上讨论了频率采样设计的基本方法，但仍存在两个问题，一个是为了保证H(z)具有线性相位，对Hd(k)应加以怎样的约束条件； 另一个是逼近误差的问题以及怎样改进以减少误差？下面将分别解决这两个问题。
1. 用频率采样法设计线性相位滤波器的条件
前面已经指出，具有第一类线性相位的FIR滤波器，其单位脉冲响应h(n)是实序列，且满足条件h(n)=h(N-1-n)，由此推导出其传输函数应满足的条件是


Hd(ejω)=Hg(ω)ejθ(ω)

θ(ω)=-12(N-1)ω （3.4.18）


将Hd(k)写成幅度Hg(k)和相位θ(k)的形式


Hd(k)=Hg(k)ejθ(k)，k=0,1,2,…,N-1（3.4.19）


式中


Hg(k)=Hg(ω)|ω=2πk/N （3.4.20）

θ(k)=θ(ω)|ω=2πk/N=-N-1Nπk（3.4.21）

当N=奇数时，属于情况1，Hg(ω)关于ω=π偶对称，即Hg(ω)=Hg(2π-ω)，将ω=2πNk代入可得：

Hg(k)=Hg(N-k)，N=奇数  （3.4.22）

当N=偶数时，属于情况2，Hg(ω)关于ω=π奇对称。将ω=2πNk代入可得

Hg(k)=-Hg(N-k)，N=偶数  （3.4.23）


式（3.4.21）~式（3.4.23）就是Hd(k)幅度特性必须满足的线性相位的条件。
例如用理想低通作为希望设计的滤波器，截止频率为ωc，采样点数为N，Hg(k)和θ(k)用下面公式计算： 
  N为奇数时


Hg(k)=Hg(N-k)=1，k=0,1,2,…,kc

Hg(k)=0，k=kc+1,kc+2,…,N-kc-1

θ(k)=-N-1Nπk，k=0,1,2,…,N-1


N为偶数时



Hg(k)=1，k=0,1,2,…,kc

Hg(k)=0，k=kc+1,kc+2,…,N-kc-1

Hg(N-k)=-1，k=0,1,…,kc

θ(k)=-N-1Nπk，k=0,1,…,N-1


注意，上面公式中的kc是小于等于Nωc2π的最大整数。另外，对于高通和带阻滤波器，这里N只能取奇数。
对具有第二类线性相位的FIR滤波器可进行类似讨论。θ(k)=π2-N-1Nkπ，当N=奇数时，属于情况3，此时，Hg(k)=-Hg(N-k)； 当N=偶数时，属于情况4，此时，Hg(k)=Hg(N-k)。
2. 频率采样法的设计误差及其改进
如果待设计的滤波器的频率响应为Hd(ejω)，对应的单位脉冲响应为hd(n)


hd(n)=12π∫π-πHd(ejω)ejωndω


根据频域采样定理可知，在频域0～2π之间对Hd(ejω)等间隔采样N点，再利用IDFT得到的h（n）以N为周期，进行周期延拓，再乘以RN(n)，即


h(n)=∑∞r=-∞hd(n+rN)RN(n)


如果Hd(ejω)有间断点，那么相应的单位脉冲响应hd(n)应是无限长的。这样，由于时域混叠，引起所设计的
h(n)和hd(n)有偏差。显然，如果频率的采样点数N越大，所设计出的滤波器越逼近待设计的滤波器Hd(ejω)。
上面是从时域方面分析其设计误差的来源，下面从频域分析。


图3.4.11理想低通滤波器的过渡带

优化示意图

上面已经提出，频率域等间隔采样H(k)，经过IDFT得到h(n)，其Z变换H(z)和H(k)的关系为


 H(z)=1-z-NN∑N-1k=0H(k)1-ej2πNz-1 （3.4.24）


将z=ejω代入式（3.4.24），得


 H(ejω)=∑N-1k=0H(k)Φω-2πNk（3.4.25)


式中

Φ(ω)=1Nsin(ωN/2)sin(ω/2)e-jωN-12 （3.4.26)

式（3.4.26）表明，在采样点ω=2πk/N，k=0,1,…,N-1，Φ(ω-2π/N)=1，因此，采样点处H(ejωk)（ωk=2πk/N）与H(k)相等，逼近误差为0，而在采样点之间的值H(ejωk)，由式（3.4.25）可知，它由有限项的H(k)Φ(ω-2πk/N)之和形成，因而有一定的逼近误差，这种误差的大小取决于理想滤波器响应的形状。理想滤波器响应越陡峭，则逼近误差越大，理想频率特性在非采样点处产生较大的肩峰和波纹。为了减小逼近误差，可以在理想频率响应的边缘加上一些过渡的采样点。如图3.4.11所示，这样大约可以分别获得20～39dB（一个过渡点），40～59dB（二个过渡点），60～70dB（三个过渡点）的阻带增益。但这样处理却会使过渡带加宽。这一点往往限制了频率采样法在FIR滤波器设计中的使用。
因为窄带频率特性的非零值采样点值比较少，频率采样法非常适合于窄带滤波器设计，但是由于存在逼近误差，使得滤波器的截止频率控制比较困难，除非截止频率点正好是采样点。增加采样点数N可提高所需滤波器的性能，但是N太大会使滤波器成本和运算复杂度增加。一般可由通过过渡带Δω-来估算N值，即


Δω-≈(m+1)2πN  （3.4.27)


式中，m为过渡带采样点数目。
［例3.4.6］用频率采样法设计第一类线性相位FIR低通滤波器，要求截止频率ωc=π16,过渡带宽度Δω-=π32，阻带最小衰减αs=30dB。
解（1） 由过渡带Δω-=π/32，αs=30dB，可知过渡带采样点数m=1，故总的频率采样点数为N=2πΔω-(m+1)=4ππ/32=128。
对于第一类线性相位FIR滤波器，N为偶数（也可为奇数），有Hd(ω)=-Hd(π-ω)，即


Hd(ejω)=e-jωN-12=e-j1272ω,0≤ω≤π16

0,π16<ω<31π16

-e-j1272ω,31π16≤ω≤2π



（2） 采样（加入一个过渡采样点）


Hd(k)=Hdej2πNk=e-jN-1Nπk=e-j127128πk,k=0,1,2,3,4

0.3904e-j1271285π,k=5(过渡采样点)

0,k=6~122

-0.3904e-j127128123π,k=123(过渡采样点)

-e-j127128πk,k=124,125,126,127


式中，0.3904（阻带最小衰减达30dB）是程序优化造成的结果。
由式（2.4.9）可知


H(z)=1-z-NN∑N-1k=0Hd(k)1-W-kNz-1


MATLAB演示如下： 

%频率采样法FIR演示

N=128;a=(N-1)/2;

m=0:N-1;W1=(2*pi/N)*m;

Hideal=［ones(1,5),0.3904,zeros(1,117),0.3904,ones(1,4)］;%理想频率域样本

Hdr=［1,1,0,0］;Wd=［0,1/16,1/16,1］; %理想频率响应

k1=0:floor((N-1)/2);k2=(floor((N-1)/2)+1):N-1;

angH=［-a*2*pi*k1/N,a*2*pi*(N-k2)/N］;%采样特性

H=Hideal.*exp(j*angH);%采样点频率性

h=real(ifft(H,N));%实际频率脉冲响应

［H2,w］=freqz(h,1);%获得实际频率特性

subplot(311);

plot(W1(1:65)/pi,Hideal(1:65),'o',Wd,Hdr);%画出理想采样特性及频率点数据

axis(［0,1,-0.2,1.2］);title('频率样本');

xlabel('频率（单位pi）');ylabel('Hideal(k)');

subplot(312);

stem(m,h);title('单位脉冲响应');

xlabel('n');ylabel('h(n)');

subplot(313);

plot(w/pi,20*log10(abs(H2)/abs(H2(1))));

axis(［0,1,-60,10］);grid;

title('幅度响应');xlabel('频率（单位pi）');ylabel('对数幅度（单位dB）');

运行结果如图3.4.12所示。


图3.4.12例3.4.6图


［例3.4.7］利用频率采样法设计一线性相位FIR低通滤波器，给定N=21，通带截止频率ωc=0.15πrad。求出h(n)，为了改善其频率响应采取什么措施？
解（1） 确定逼近滤波器传输函数Hd(ejω)


Hd(ejω)=e-jωα,0≤|ω|≤0.15π,


0,0.15π<|ω|≤π, α=(N-1)/2=10


（2） 采样


Hd(k)=Hdej2πNk=e-jN-1Nπk=e-j2021πk,k=0,1,20

0,2≤k≤19


（3） 求h(n)


h(n)=IDFT［Hd(k)］=1N∑N-1k=0Hd(k)W-knN

=1211+e-j2021πW-n21+e-j2021π·20W-20n21R21(n)

=1211+e-j2π21(n-10)+e-j400π21ej4021πnR21(n)


因为e-j40021π=ej2021π，e-j4021πn=ej4221π-221πn=e-j2π21n，所以


h(n)=1211+e-j2π21(n-10)+e-j2π21(n-10)R21(n)

=1211+2cos2π21(n-10)R21(n)


为了改善阻带衰减和通带波纹，应加过渡带采样点，为了使边界频率更精确，过渡带更窄，应加大采样点数N。
用MATLAB演示其对应的频率特性如下： 

%FIR低通演示

N=21;n=0:1:(N-1);

h=(1+2*cos((n-10)*2*pi/N))/N;

freqz(h,1);

运行结果如图3.4.13所示。


图3.4.13例3.4.7图


［例3.4.8］用频率采样法设计FIR线性相位高通滤波器，截止频率ωp=3π4，采样间隔Δω为π6，设一点过渡H1=0.39。
（1） 求采样点数N。
（2） 该数字滤波器的h(n)、Hg(ω)和θ(ω)各有什么特点？
（3） 求出Hd(k)。
解（1） 采样点数N=2πΔω=2ππ/6=12。
（2） 由于是高通，N为偶数，故只能用第4种情况来设计。h(n)=-h(N-1-n)奇对称，Hg(ω)以ω=0,2π奇对称，以π偶对称，θ(ω)=-N-12ω-π2。
（3） ωc=2πNk，k=ωcN2π=34π×122π=4.5。
 Hg(k)=Hg(N-k)，Hg(k)在5、6、7点上为非零值，因此


θ（ω）=-N-12ω-π2ω=2πkN=-1112πk-π2

Hd(k)=e-j1112πk+π2,k=5,6,7

0.39e-j1112πk+π2,k=4,8

0,k=0～3,9～11


用MATLAB演示的设计过程如下： 

%频率采样法FIR演示

N=12;a=N/2;

m=0:N-1;W1=(2*pi/N)*m;

Hideal=［0,0,0,0,0.39,1,1,1,0.39,0,0,0］;

Hdr=［0,0,1,1］;Wd=［0,0.75,0.75,1］; %理想频率响应

k1=0:N/2;k2=(N/2+1):N-1;

angH=［-a*2*pi*k1/N,a*2*pi*(N-k2)/N］;%采样特性

H=Hideal.*exp(j*angH);  %采样点频率性

h=real(ifft(H,N));  %实际频率脉冲响应

［H2,w］=freqz(h,1);  %获得实际频率特性

subplot(311);

plot(W1(1:7)/pi,Hideal(1:7),'o',Wd,Hdr); 

axis(［0,1,-0.2,1.2］);title('频率样本');

xlabel('频率（单位＼pi）');ylabel('Hideal(k)');

subplot(312);

stem(m,h);title('单位脉冲响应');

xlabel('n');ylabel('h(n)');

subplot(313);

A=max(abs(H2));

plot(w/pi,20*log10(abs(H2)/A));

axis(［0,1,-60,10］);grid;

title('幅度响应');xlabel('频率（单位＼pi）');ylabel('对数幅度响应（单位dB）');

运行结果如图3.4.14所示。


图3.4.14例3.4.8图


［例3.4.9］利用频率采样法设计线性相位FIR带通滤波器，设N=33，理想幅频特性Hg(ω)，如图3.4.15所示。


图3.4.15例3.4.9图


解由图可知理想幅度采样值为

Hg（k）=Hg（ω）ω=2πkN=1，k=7,8,25,26


0,其他


因此Hd(k)=e-j3233πk,k=7,8,25,26
0,其他


h(n)=IDFT［Hd(k)］

=133e-j3233π·7ej2π33π·7n+e-j3233π·8ej2π33·8n+e-j3233π·25ej2π33·25n+e-j3233π·26ej2π33·26nR33(n)

=233cos14π33(n-16)+cos16π33(n-16)R33(n)


用MATLAB画出其对应的频率特性如下： 

%线性相位FIR带通滤波器

N=33;n=0:1:(N-1);

h=(cos((n-16)*14*pi/N)+cos((n-16)*16*pi/N))/N;

freqz(h,1);axis(［0 1 -80 0］);

运行结果如图3.4.16所示。


图3.4.16例3.4.9图


［例3.4.10］用频率采样法设计一个带通滤波器，技术指标如下： 下阻带边缘为0.3π，下通带边缘为0.4π，上通带边缘为0.5π，上阻带边缘为0.6π，αp=0.5dB,αs=50dB,选择适合的滤波器阶数，使通带中有两个优化点。
例题的MATLAB辅助求解源程序如下： 

%频率采样BPF优化设计演示

Ws1=0.3*pi;Ws2=0.6*pi;Wp1=0.4*pi;Wp2=0.5*pi;

Rp=0.5;As=50;  %BPF技术参数

b=min((Wp1-Ws1),(Ws2-Wp2));

N=ceil(2*pi*3/b)+20;  %获得频率采样点数

k=0:N-1,a=(N-1)/2;W1=(2*pi/N);

Wk=(2*pi/N)*k;

T1=0.113;T2=0.605;  %设过渡带采样优化值

Hideal=［zeros(1,ceil(0.3*pi/W1)+1),T1,T2,ones(1,ceil(0.1*pi/W1)+2),T2,T1, zeros(1,ceil(0.8*pi/W1)+2),T1,T2,ones(1,ceil(0.1*pi/W1)+2),T2,T1, zeros(1,ceil(0.3*pi/W1)+1)］; 
%理想频率样本采样值

Hdr=［0,0,1,1,0,0］;Wd=［0,0.3,0.4,0.5,0.6,1］; 

k1=0:floor(( N-1)/2);k2=floor((N-1)/2)+1:(N-1);

angH=［-a*2*pi*k1/N, a*2*pi*(N-k2)/N］;  %相位样本

H=Hideal.*exp(j*angH);  %频率采样数据

h=real(ifft(H,N));  %实际频率脉冲响应

［H1,w］=freqz(h,1); 

figure(1);

plot(Wk(1:ceil(N/2))/pi,Hideal(1:ceil(N/2)),'o',Wd,Hdr); 

xlabel('频率（单位＼pi）');ylabel('幅度响应');

title('频域样本数据');grid;axis(［0,1,-0.2,1.2］);

figure(2);

stem(k,h);title('实际单位脉冲响应');

xlabel('n');ylabel('h(n)');

figure(3);

plot(w/pi,20*log(abs(H1))); title('幅度响应'); grid; 

xlabel('频率（单位＼pi）');ylabel('对数幅度响应 （单位dB）');

axis(［0,1,-120,10］);

运行结果如图3.4.17所示。


图3.4.17例3.4.10图




图3.4.17（续）

3.4.4FIR滤波器的最优等波纹设计法
对于FIR设计中的窗函数法和频率采样法，它们在设计中不能将边缘频率ωs和ωp精确给定，而且它们的逼近误差在频带区间上不是均匀分布的。在靠近频带边缘处逼近误差大，在远离频带边缘处误差小，按照IIR滤波器设计时的切比雪夫逼近原理可知，如果我们能将这种逼近误差均匀铺开，显然能得到阶数较窗函数法和频率采样法要小的FIR线性相位滤波器。
最优等波纹逼近方法实际上是一种加权最大误差最小化的设计方法。设所设计的线性相位滤波器为 


H(ejω)=H(ω)ej(N-1)ω2（3.4.28)


具有所希望特性的滤波器Hd(ejω)为



Hd(ejω)=Hd(ω)ej(N-1)ω2（3.4.29)


则衡量Hd(ω)与H(ω)的逼近误差的评价函数为


J=max ω|W(ω)［Hd(ω)-H(ω)］|，0≤ω≤π（3.4.30)


式中，W(ω)为误差加权函数。W(ω)的引入能同时优化边带波纹δ1和阻带δ2，这一点与窗函数法和频率采样法是显然不同的。
过渡带的范围由设计者自行选择，若过渡带过宽，能得到的滤波器的截止特性就比较缓慢，反之，若过渡带过窄，滤波器的波纹就会比较大。
目前应用普遍的方法是由Parks 和McClellan 提出来的 Remez（瑞米兹）交换算法。在此对其不作过多的解释，仅借助MATLAB来演示这种滤波器的设计过程及主要特性。
［例 3.4.11］利用Remez算法来设计一个等波纹FIR低通滤波器，技术特性如下


ωp=0.2π,αp=0.25dB,ωs=0.3π,αs=50dB


MATLAB程序如下： 

%Remez算法FIR LPF设计演示

Wp=0.2*pi;Ws=0.3*pi;Rp=0.25;As=50;

a_w=2*pi/1000;

Wsi=Ws/a_w+1;

a1=(10^(Rp/20)-1)/(10^(Rp/20)+1);%参数转换

a2=(1+a1)*(10^(-As/20));  %参数转换  

aH=max(a1,a2);aL=min(a1,a2);

weights=［a2/a1,1］;

af=(Ws-Wp)/(2*pi);

N=ceil((-20*log10(sqrt(a1*a2))-13)/(14.6*af)+1);  %根据经验公式估计N

%可以计算出N=43

f=［0 Wp/pi Ws/pi 1］;  %频率向量

m=［1 1 0 0］;  %对应频点期望幅度响应向量

h1=remez(N-1,f,m,weights);%

［H1,w］=freqz(h1,1);

db=20*log10(abs(H1));

Asd=-max(db(Wsi:1:501));

%可以计算出Asd=34.017，显然不合要求，按N=N+1进行循环调用，可知当N=47时满足要求

N=47;n=0:1:N-1;

h2=remez(N-1,f,m,weights);

［H2,w］=freqz(h2,1);

figure(1);

subplot(211);

stem(n,h2);title('单位脉冲响应');xlabel('n');ylabel('h(n)');

subplot(212);

B=max(abs(H2));

plot(w/pi,20*log10(abs(H2)/B));title('幅度响应');grid;

xlabel('频率（单位＼pi）');ylabel('对数幅度响应 (单位 dB)');

figure(2);

L=(N-1)/2;

b=［h2(L+1) 2*h2(L:-1:1)］;

n1=［0:1:L］;

w1=［0:1:500］'*pi/500;

Hr=cos(w1*n1)*b';

subplot(211);

plot(w1/pi,Hr);grid;  %画出振幅响应曲线

xlabel('频率（单位＼pi）');ylabel('幅度');

title('振幅响应');

subplot(212);  %画出逼近误差等波纹曲线

H3=(Hr(1:101))'-ones(1,101);

H31=［H3 zeros(1,400)］;

H4=(Hr(151:501))';H41=［zeros(1,150) H4］;

plot(w1/pi,H31+H41);grid;

xlabel('频率（单位＼pi）');ylabel('幅度');

title('逼近误差等波纹曲线');axis(［0 1 -0.015 0.015］);

运行结果如图3.4.18所示。



图3.4.18例3.4.11图

请读者观察一下图3.4.18中的逼近误差等波纹交错点的个数，看看它和N之间的关系。
3.5有限字长效应
前面所讨论的离散时间信号和系统都未涉及系数的精度问题，即无论是序列还是滤波器的系数，都是以无限精度的数据来表示的。但是一个实际的数字信号处理系统不论是用计算机的软件来实现还是专用的数字硬件来实现，它们的输入、输出、中间结果以及滤波系数等都必须存储在有限字长的存储器中，另外对于输入信号也要进行量化，无限精度的模拟信号经过模数变换后只能取有限多个可能值。这样在信号处理时所得到的结果和理论值之间存在误差，这就是有限字长效应。
和有限字长效应有关的误差如下： 
（1） 模/数（A/D）转换器将输入信号变成一组离散值时产生的量化误差。
（2） 用有限位二进制数来表示数字系统各参数时产生的量化误差。
（3） 按所需的算法进行运算时，为限制位数扩展而进行尾数处理以及防止溢出而压缩信号大小所产生的误差。
（4） 溢出振荡器产生的误差。
（5） 固定信号输入时产生的极限环振荡带来的误差。
上述各种量化误差与数字系统的结构形式、数的表示方法、数位长短以及采用的运算方法有关。当然可以采用位数多的通用计算机或字长较长的数字硬件系统，以减小有限字长效应，但付出的代价是成本昂贵、设备复杂、运行时间长等。通常在一般的通用计算机上做信号处理，由于字长较长，所以有限字长效应影响不大，可以不予考虑。而用专用的硬件做信号处理时，考虑到成本、复杂性以及运算速度等因素，字长受到较大的限制，因此在精度和造价之间作合理的折中。
专用机大多数采用定点制表数和运算，具有快速、经济等优点，本节主要讨论定点制情况下有限字长效应，由于运算中多用补码，所以讨论中侧重补码运算。
数字系统的有限长效应是系统实现时必须考虑的实际问题，但由于其复杂性及随机性因而有些误差在理论上尚无固定的解释。本节对一些常用的数字系统中的有限字长所引起的误差或现象作原理上的分析，了解它们的本质，供读者实际应用时参考。
3.5.1数的表示方法对量化的影响
1. 数制和数的表示方法

1） 定点制
通常定点制表示的数限制为±1。对于一个数x用定点数示为
x=β0β1β2…βb，其中β0是符号位，当x≥0时β0=0，当x<0
时β0=1。以b表示数据位的位数，则表示寄存器的位长为b+1位。定点制的运算结果的绝对值不能大于1，否则出现“溢出”错误。为使整个运算的绝对值不超过1，需要引入比例因子，以减少或限制输入信号序列的动态范围。定点制加法运算可能溢出，乘法运算不会溢出，但字长会增加一倍，每次要对结果截尾，引入截尾误差。解决的办法是对信号序列及运算结果做归一化处理。但是滤波器系数不能随便归一化，应预留小数位。
定点制的缺点是动态范围小，可能产生溢出。
2） 浮点制
浮点制是将一个数表示为尾数和指数两部分，即x=±M·2c其中c称为阶码，决定数的范围大小； M称为尾数，决定数的精度。浮点数无论进行加法还是乘法都会使尾数部分加长，需要进行尾数处理。浮点知道的特点是在运算过程中小数点的位置时浮动的。
3） 数的三种表示方法
二进制数有三种表示方法： 原码、反码和补码。正数的三种表示方法是相同的，且符号位0； 而负数的三种表示方法，它们的符号位都是1，但是尾数的二进制码各不相同。
对于原码有

x=

|x|,x≥0


1+|x|,x<0


对于反码有

x=

|x|,x≥0


2-|x|-ab,x<0


对于补码有

x=

|x|,x≥0


2-|x|,x<0


原码的优点是乘法运算方便，缺点是加减运算复杂。采用补码，加法运算方便，但是乘法运算复杂。反码和补码相似，加法运算方便，但是运算规则稍有不同。现在已有将补码乘法做成专用集成芯片，使补码获得广泛的应用。
2. 尾部量化方式及量化误差
量化方式有截尾和舍入两种方式。截尾处理是保留b位码，抛掉余下的尾数； 而舍入处理是按最接近的值取b位码。这两种处理所产生的误差是不一样的。通常用量化阶距q表示量化精度，它是最小码位所代表的数值，即q=2-b。
1） 定点制的量化误差
对于定点制的截尾误差有Er=Q［x］-x，分布范围如下： 
 原码、反码： -q<Er<q。
 补码： -q<Er<0。
定点制舍入误差不论是正数还是负数，也不论是原码、反码还是补码，其误差总是在±q/2之间，即正数负数截尾相对误差范围是

-q/2<Er≤q/2


2） 浮点制的量化误差
在浮点制中，截尾或舍入的处理只影响尾数的字长。对于尾数进行截尾或者舍入处理所产生的误差和定点量化完全相同。如以EM表示尾数的量化误差，即EM=［M］Q-M。由于x=M·2c，因此x的量化误差Ex=［x］Q-x=EM·2c，与阶码值有关。所以浮点制的量化误差是和数字本身的大小有关。通常用相对误差ε=Exx来描述。
浮点制舍入相对误差范围为-q<εr≤q。
浮点制截尾相对误差范围如下： 
 原码、反码： -2q<εr≤0。
 补码： -2q<εr≤2q。
需要指出的是，由于截尾误差量化噪声具有直流分量，将影响信号的频谱结构，所以一般采用舍入处理。
［例3.5.1］将下列十进制数分别用b＝4的原码、反码和补码表示

x1=0.4375，x2=0.625，x3=-0.4375，x4=-0.625

解因为 ［0.4375］10=［0.0111］2，而且对于正数，原码＝反码＝补码，所以
［x1］原=［x1］反=［x1］补=0.1111
因为［0.625］10=［0.1010］2，而且对于正数，原码＝反码＝补码，所以
［x2］原=［x2］反=［x2］补=0.1010
因为［0.4375］10=［0.0111］2而且对于负数，原码＝1+|x|，所以
［x3］原=1+0.0111=1.0111
对于负数的反码=2-2-b-|x|，所以
［x3］反=10-0.0001-0.0111=1.0000
对于负数的补码=2-|x|，所以
［x3］补=10-0.0111=1.1001，［x4］原=1.1010，［x4］反=1.0101，［x4］补=1.0110 
［例3.5.2］若以下二进制数码分别是原码、反码、补码时，请算出其对应的十进制。
x1=0.1001，x2=0.1101，x3=1.1011，x4=1.0000

解x1=0.1001，x2=0.1101，它们的符号位为0，所以不管它是原码、反码还是补码，它均表示相同的数。x1=0.1001表示的十进制数是-(2-1+2-4)=-0.5625。x2=0.1101所表示的十进制数是-(2-1+2-2+2-4)=-0.8125。x3=1.1011，x4=1.0000，它们的符号位为1，它们若是原码，则x3=1.1011所表示的十进制数是-(2-1+2-3+2-4)=-0.6875； x4=1.0000
所表示的十进制数是-0。它们若是反码，则x3=1.1011所对应的原码是1.0100，它表示的十进制数是-(2-2)=-0.25； x4=1.0000所对应的原码是1.1111，它表示的十进制数是-(2-1+2-2+2-3+2-4)=-0.9375。它们若是补码，则x3=1.1011所对应的原码是1.0101，它表示的十进制数是-(2-2+2-4)=-0.3125； x4=1.0000
表示的十进制数是-1。
3.5.2A/D转换的量化效应
A/D转换实际上可分为如图3.5.1所示的两部分： 采样和量化。采样实现时域上的离散化，而量化则把一个无限精度的模拟信号转换成数字信号。量化


图3.5.1A/D转换器功能原理图

时要进行截尾或舍入处理，会产生误差。A/D转换通常采用定点制的线性量化。
1. 量化效应的统计分析
A/D转换通常采用定点制的线性量化。若A/D转换采用定点补码截尾处理方式，量化就是把信号的理想值x(n)用其量化值表示x^(n)=Q{x(n)}，必然会引入量化误差e(n)，即
x^(n)=
x(n)+e(n）。一般量化噪声的统计分析是基于如下一些假设： 


（1） e(n)是一个平稳的随机序列。
（2） e(n)与信号x(n)是不相关的。
（3） e(n)序列本身的任意两个值之间也是不相关的，即e(n)是白噪声序列。
（4） e(n)在其误差范围内具有均匀等概率的分布特性。
通常情况下，信号越不规则，这种假设就越接近于实际，因此作为一种平均的大概的特性分析来说，这种假设是合适的。
在以上假设下，截尾量化噪声e(n)的均值和方差分别为-q/2和q2/12。舍入量化噪声（误差）e(n)的均值和方差分别为0和q2/12。
由于量化噪声方差和A/D转换器的字长有直接的关系，字长越长，q越小，量化误差σe便越小。量化后的信噪比定义为信号的功率与量化噪声功率之比

SNR=10lgσ2xσ2e=6.02b+10.79+10lgσ2x

从上式可以看出，字长每增加一位，信噪比提高6dB。另外，输出信噪比还和输入信号的功率有关，为提高输出信噪比，可提高输入信号功率，但信号xa(t)的幅度不能超过A/D转换器的动态范围。
2. 量化噪声通过线性系统
量化噪声通过线性系统

y^（n）=x^（n）*h（n）=x（n）*h（n）+e（n）*h（n）

记f(n)=e(n)h(n)，利用随机信号处理的基本理论，有μf=μeH(ej0)

σ2f=σ2e∑∞n=-∞h2(n)=12πσ2e∫π-π|H(ejω)|2dω=σ2e2πj ∮cH(z)H(z-1)z-1dz

［例3.5.3］有一个b=7的A/D转换器，它的输出x^（n）通过线性系统

H(z)=11-0.99z-1

求系统输出端的量化噪声功率。
解由题意已知b=7，故q=2-7，量化噪声的功率为

σ2e=112q2=112·2-14

设量化噪声通过线性系统后的输出噪声功率为σ2f，则有


σ2f=σ2e∑∞n=-∞h2(n)=12πσ2e∫π-π|H(ejω)|2dω=σ2e2πj∮cH(z)H(z-1)z-1dz
 =σ2e2πj∮c11-0.99z-1·11-0.99zdz
=2-141211-0.992=2.556×10-4


3.5.3数字滤波器的有限字长效应
数字滤波器中量化效应表现在两个方面： 运算量化效应和数字滤波器的系数量化效应。
数字滤波器的实现可采用定点运算实现，也可以采用浮点运算来实现，两种方法实现的运算误差是不同的。分别讨论如下： 
1. 数字滤波器定点制运算误差分析
定点制加法运算可能会产生溢出，乘法运算会引起尾数增长，需要进行截尾或者舍入处理，这都将影响滤波器的正常工作。下面分别讨论。
1） 加法运算的溢出问题及其解决办法
定点制中，两个b位的小数相加，仍然是b位小数。因此，加法运算不会产生量化误差，但总和的绝对值可能超过1，将出现溢出问题，由于溢出，符号位变号，总和为正数时，由于溢出将变为负数，而负数则变成正数。因此带来很大的误差。对于IIR滤波器它还使滤波器在最大幅度界限之间振荡。因此要避免溢出。
首先要理解IIR滤波器中的零输入极限环振荡的现象： IIR滤波器是一个反馈系统，当它的所有极点均位于单位圆内时，系统肯定是稳定的。即当输入为零输入时，系统的输出逐渐衰减趋向于零。但在有限字长情况下，数字滤波器零输入时，由于补码的加法溢出非线性特性，这时输出不为零（可能会衰减到某一非零的幅度范围，而后呈现振荡特性），这种现象称为零输入极限环振荡。一阶系统的极限环振荡幅度与量化阶距成正比，因此增加字长将使极限环振荡减弱。
通常消除溢出有以下几种常用的方法。
（1） 用补码饱和输出特性。具体实现时，先判断补码加法结果是否产生溢出，若产生了，再作相应的处理。
（2） 限制输入信号的动态范围。
（3） 为防止溢出产生，采用更好的滤波器结构，用标准的滤波器结构实现。标准滤波器的设计涉及矩阵论的知识，这里只要了解就行。
2） 乘法运算的量化效应
一般用其统计模型分析方法进行分析，了解这种方法的基本假设和基本分析方法： 认为每一个乘法运算都有一个独立的白噪声误差源，再根据线性系统理论，可以分别计算出各噪声源单独作用是在输出端形成的输出噪声，总的输出噪声则是所有输出噪声的线性叠加。
不同的滤波器结构和不同的编排顺序，量化效应是不同的，一般直接型误差大于级联型误差，级联型误差又大于并联型误差。这是因为在直接型的结构中所有舍入误差都要经过全部网络的反馈环节，因此使这些误差在反馈过程中积累起来，致使误差增大。在级联型结构中，每个输入误差只通过其后面的反馈环节，而不通过它前面的反馈环节，因而误差比直接型小。在并联型结构中，误差仅仅通过本通路的反馈环节，与其他并联网络无关，因此累积作用小。同时在级联时，按系数大小顺序或极点距原点距离的大小编排，可望得到最小的误差。
2. 数字滤波器浮点制运算误差分析
浮点制的优点是动态范围大，一般不考虑溢出问题，同时也不考虑零输入极限环振荡。但浮点制运算中，每次加法之后和乘法一样需作尾数处理，因此都会引入量化误差。浮点量化误差是一个乘性误差，与量化信号幅度无关。常以相对误差来表示。它的量化处理可表示为

Q［x(n)］=x(n)［1+ε(n)］=x(n)+x(n)ε(n)=x(n)+e(n)

式中，Q ［x(n)］是量化后的x(n); ε(n)是相对量化误差; e(n)是量化误差。
只要求了解浮点量化相对误差的统计分析方法，对分析结果应该理解： 浮点制中不论信号大小，也不论信号的结构、分布如何，滤波器的输出都具有相同的信噪比，或者说都具有同样的相对精度。这也与定点制不同，在定点制中输出噪声的方差是和信号无关的。因此，信号越大，输出的信噪比也就越大； 信号小，信噪比也就小。但是，定点制中，由于受到溢出的限制，信噪比的提高会受到引入比例因子的限制。
在相同的尾数字长下，浮点制的误差比定点制小。
3. 数字滤波器的系数量化效应
由于滤波器的所有系数都必须以有限位的二进制码的形式存放在寄存器中，所以就必须对理论的系数值进行量化。这样就使实际系数存在误差，使滤波器的零极点位置发生偏离，从而影响到滤波器的性能，严重时会使单位圆内的极点偏离到单位圆外，使系统失去稳定性。
系数量化影响大小不仅和字长有关，同时也和滤波器的结构有密切的关系。对于高阶滤波器来讲，应该尽量避免使用直接型结构，而是应该尽量分解为低阶的级联或者并联系统，这样，在给定的字长情况下，可以使系数量化的影响最小。
下面直接利用MATLAB对T阶椭圆低通滤波器分析其量化前后的频率特性及零极点位置偏移情况。该滤波器的技术指标为ωc=0.5π,通带波纹为1dB，阻带衰减为60dB，量化时对滤波器进行截尾处理。

% IIR量化效应演示 f327.m

［b,a］=ellip(7,0.5,60,0.5);

［H,w］=freqz(b,a,512);%求量化时的特性分析

H1=20*log10(abs(H));

%下面的程序将a,b中的参数进行6位截尾量化处理

m=1,b1=abs(b);

while fix(b1)0

b1=abs(b)/(2^m);

m=m+1

end

c1=fix(b1*2^6);

c=sign(b).*c1.*2^(m-6-1);%b量化后变成c

n=1,b2=abs(a);

while fix(b2)0

b2=abs(a)/(2^n);

n=n+1

end

d1=fix(b2*2^6);

d=sign(a).*d1.*2^(m-6-1);%a量化后变成d

［Hq,w］=freqz(c,d,512);

H2=20*log10(abs(Hq));

subplot(121);

plot(w/pi,H1, 'b',w/pi,H2, 'r:');

grid;

axis(［0 1 -100 10］);

xlabel('归一化频率');ylabel('幅频响应特性');

legend('量化前','量化后');title('相位量化前后幅频特性分析');

subplot(122);

［z1,p1,k1］=tf2zp(b,a);

［z2,p2,k2］=tf2zp(c,d);

zplaneplot(［z1,z2］,［p1,p2］,{'s','+','d','*'}); 

legend('量化前的零点', '量化前的极点', '量化后的零点', '量化后的极点');

title('相位量化前后零极点分析');


结果如图3.5.2所示。


图3.5.2量化前后的频率特性及零极点位置偏移情况

［例3.5.4］为在数字计算机上处理序列，必须将序列的幅度量化成一组离散电平。这种量化过程可用输入序列x(n)通过一


图3.5.3例3.5.4图

个量化器Q［x］表示，Q［x］的输入输出关系如图3.5.3所示。
如果量化间距和输入序列电平的变化相比很小，则可以假设量化器输出y(n)
的形式为y(n)=x(n)+e(n)，e(n)是一个平稳随机过程，它是在［-Δ/2，Δ/2］之间均匀分布的，它在各取样间互不相关，它也与x(n)也独立无关。因此对于所有的m和n有： E［e(n)x(n)］=0。令x(n)是均值为0，方差为σ2x的平稳白噪声随机过程。
(1) 求e（n）的平均值、方差和自相关序列。
（2） 求信号—量化噪声比σ2x/σ2e。
(3) 把量化的信号用一个单位取样响应为h(n)=12［an+(-a)n］u(n)的数字滤波器滤波，这里假定a为实数，且|a|<1，试确定输出端的信噪比。
解(1) 由题知e(n)在［-Δ/2，Δ/2］之间均匀分布，所以有

均值

me=E［e(n)］=∫Δ2-Δ2epen(e,n)de=∫Δ2-Δ2e1Δde=0

方差

σ2e=E［(e(n)-me)2］=∫Δ2-Δ2e2(n)pen(e,n)de=Δ212

因为n≠0时，自相关序列为

φee(m)=E{x(n)x(n+m)}=0

所以

φee(m)=σ2eδ(m)

(2) 量化后的信噪比为

σ2xσ2e=σ2x·12Δ2


(3) 设输出噪声为f(n)，则

f(n)=e(n)h(n)


由于me=0，因此mf=0。因为e(n)是白色的，所以


σ2g=σ2e∑∞n=-∞|h(n)|2=σ2e14∑∞n=-∞［an+(-a)n］2
＝σ2e14∑∞n=-∞［2a2n］2=σ2e11-a4


设输出信号为s(n)，则

s(n)=x(n)h(n)

由于mf=0，因此

ms=mx∑∞n=-∞h(n)=0

因为x(n)是白色的，所以

σ2s=σ2x∑∞n=-∞|h(n)|2=σ2x·11-a4


故通过线性系统后的输出信噪比为

σ2sσ2g=σ2xσ2e=12·σ2xΔ2

输出信噪比和输入信噪比是一样的。
［例3.5.5］一个一阶IIR网络，理想运算的差分方程是y(n)=ay(n-1)+x(n)
，用定点制原码运算，尾数作截尾处理。求证： 只要|a|<1，就不会发生极限振荡环。
证明： 用Q［］代表原码截尾，当输入为0时，差分方程为

w(n)=Q［a·w(n-1)］

在形成极限环时，|w(n)|=|w(n-1)|，即|w(n)|=±|w(n-1)|，现在分别对下面几种情况进行讨论： 
情况1： Q［a·w(n-1)］=w(n-1)。此时a>0

Er=Q［a·w(n-1)］-a·w(n-1)=(1-a)w(n-1)

① 若w（n-1）保持大于0，由于正数的截尾误差为-2-b<Er≤0，即（1-a）w（n-1）≤0，所以得到（1-a）≤0，a≥1。
② 若w（n-1）保持小于0，由于负数的截尾误差为0≤Er≤2-b，即（1-a）≥0，故（1-a）≤0,a≥1。
情况2： Q［a·w(n-1)］=-w(n-1)。此时a<0

Er=Q［a·w(n-1)］-a·w(n-1)=-(1+a)w(n-1)

若某一个w(n-1)>0，由于a·w(n-1)<0。又由于负数的截尾误差为0≤Er≤-2-b，即0≤-(1+a)w(n-1)，所以得到(1+a)≥0，a≥-1。
若某一个w（n-1）＜0，由于a·w（n-1）<0。又由于正数的截尾误差为-2-b<Er≤0，即-(1+a)w(n-1)≤0，所以得到(1+a)≤0，a≤-1。
总结如上情况，如果存在极限环，一定是|a|≥1。对于|a|<1，系统不存在极振荡环。
［例3.5.6］差分方程y(n)-ay(n-1)=x(n)-1ax(n-1)描述了一个全通系统。系统的实现如图3.5.4所示。
这个系统的幅频响应是一个与频率无关的常数。下面比较一下定点运算和浮点算法实现全通网络时，两种算法的舍入效应。将所有定点数看作小数，因此定点数在±1之间。令寄存器字长为b位，不包含符号位。对于浮点数，令t表示尾数的位数，不包含位数的符号位。假设滤波器输入x(n)是一个白色随机过程，它的幅度在±x0间均匀分布，设a为实数，且1/2<a<1。
 (1) 试确定定点实现的输出噪声信号比（信噪比），即求舍入噪声产生的方差与x(n)产生的输出方差之比。
(2) 试确定浮点算法实现的输出信噪比。
解(1) 定点实现的系统如图3.5.5所示。




图3.5.4例3.5.6图（1）




图3.5.5例3.5.6图（2）


一般假设噪声源e1(n)和e2(n)无关，因此

σ21=σ2e1∑∞n=-∞h2(n)+σ2e2=σ2e1(1-a-2)=(1-a-2)2-2b12

为了不发生溢出，需要y(n)≤1，要研究h(n)，整个滤波器的系统函数为


H(z)=1-1az-11-az-1=1-1az-1(1+az-1+a2z-2+…)

因此
h(n)=1,n=0
an-an-2,n≥1
0,n<0


已知1/2<a<1，有



∑∞n=-∞|h(n)|=1+∑∞n=-∞|an+an-2|+1+∑∞n=-∞1a2-1an

＝1+1a2-1a1-a=2+1a



所以为了y(n)≤1，必须有

xmax≤1∑∞n=-∞|h(n)|=a2a+1

现在研究y1(n)≤1的条件，先看系统函数


H(z)=11-az-1∑∞n=0anz-n

∑∞n=-∞|h(n)|=∑∞n=0|an|=11-|a|



因此要求xmax ≤1-|a|。
由于12<a<1，所以

14<a2a+1<13

而a=12时，1-a=12

a=12时，1-a=a2a+1

因此，取

x0=

a2a+1,12<a<12
1-a,12<a<1


由于题目中假设x(n)在±x0之间均匀分布，因此


px=12x0，mx=0，σ2x=∫x0-x0x212x0dx=x203

σ2y=1a2·σ2x=13a2·x20=31(1+2a)2，12<a<12
13a2（1-a）2，12<a<1



可以得到输出中的信号噪声


图3.5.6例3.5.6图（3）


σ2fσ2y=

2-2b-2(1+a-2)（1+2a-2）,
12<a<12
2-2b-21+a2(1-a)2,

12<a<1


(2) 在浮点情况下加法和乘法都引进误差，设误差源为e1(n)，e2(n)，…，e4(n)，如图3.5.6所示，可以写出


e1(n)≈ε1(n)y1(n)

e2(n)≈aε2(n)y1(n-1)

e3(n)≈1aε3(n)y1(n-1)

e4(n)≈ε4(n)y1(n)


由于ε1,2,3,4(n)的均值都为0，而且和信号无关。因此


σ2e1(n)=E［y21(n)］σ2ε1=σ2y1σ2ε1

σ2e2(n)=a2E［y21(n-1)］σ2ε2=a2σ2y2σ2ε2

σ2e3(n)=1a2E［y21(n-1)］σ2ε3=1a2σ2y3σ2ε3

σ2e4(n)=E［y21(n)］σ2ε4=σ2yσ2ε4



由(1)中已知

∑∞n=-∞h21(n)=1(1-a)2，∑∞n=-∞h2(n)=1a2

可得到


σ21=σ2x∑∞n=-∞h21(n)=σ2x(1-a)2

σ2y=σ2x∑∞n=-∞h21(n)=σ2xa2



按题意，尾数长为t位，因此

σ2ε1=σ2ε2=σ2ε3=σ2ε4=
13×2-2t

输出噪声的方差为

σ2f=(σ2ε1+σ2ε2)∑∞n=-∞h2(n)+(σ2ε3+σ2ε4)=2-2tσ2x(1-a2)a2

由此可以得到浮点实现的输出信噪比为

σ2fσ2y=2-2t(1-a)2


3.5.4FFT运算中的有限字长效应
就运算角度而言，FFT同样可以看作是一个系统。因此，有限字长效应的分析和数字滤波器的分析方法基本是相同的。FFT运算主要有三个方面的误差源，它们是输入量化误差，系数量化误差和运算量化误差，除了FFT中的运算是一个复数运算外，和数字滤波器的分析方法是相同的，引入噪声源，采用统计的分析方法。三种噪声源以运算量化误差最重要，另外，考虑到硬件专用FFT器件主要采用定点制，因此只讨论定点制FFT运算量化效应和防溢出措施。FFT运算结构不同时，有限字长效应的分析也会略有不同，但是结果是相似的。这里以时间抽取FFT为例说明运算中的注意事项。
FFT的核心运算是蝶形运算，通过对蝶形运算进行分析，可知


max［|x(k)|］≤2Mmax［|x(n)|］=Nmax［|x(n)|］


所以为了避免溢出，应限制x(n)使max［|x(n)|］≤1N。
但是输入幅度的限制会使输出信噪比降低，输出信噪比与N2成反比，N增加一倍，信噪比（SNR）下降为14倍，或要保持运算精度不变，每增加一级运算，字长也须增加一倍。
输出信噪比的降低是由于输入幅度被限制地过小，这种状况可通过逐级引入比例因子来改善。这种方法使信噪比得以提高，为保证精度不变，N
增加4倍，字长才需要增加一倍。
3.6本章小结
本章详细介绍了数字滤波系统以及IIR和FIR数字滤波器的各种设计方法，讨论了频域特性和量化效应分析的基本现象、基本理论和优缺点，IIR借助于已经非常成熟的各种模拟滤波器的设计图表，再通过频率变换和离散化来得到，在设计方法上来说相对成熟一些。而FIR数字滤波器除了本章介绍的窗函数法、频率采样法、优等波纹逼近法之外，还有许多非常灵活的设计方法。需要指出的是，从实际应用的角度来说，没有哪一种滤波器能够给出最好的性能指标，需要在算法开销和硬件复杂性之间选择一个良好的折中。
实际上，目前市面上流行版本的MATLAB软件包提供了两个功能非常强的互交式集成信号处理开发图形界面工具。FDATool工具（见图3.6.1）以及SPTool（见图3.6.2），它们分别由命令fdatool和sptool进行激活。


图3.6.1FDATool的界面





图3.6.2SPTool的界面


这两个工具给用户的设计使用带来了非常大的方便。在掌握了本章所讲述的基本理论之后再来融会贯通这两个设计工具将是一件非常有意义的事情。

习题
31画出H(z)=(2-0.379z-1)(4-1.24z-1+5.264z-2)(1-0.25z-1)(1-z-1+0.5z-2)级联型网络结构。
32画出 H(z)=(2-3z-1)(4-6z-1+5z-2)(1-7z-1)(1-z-1+8z-2)级联型网络结构。
33已知某三阶数字滤波器的系统函数为

H(z)=3+53z-1+23z-21-13z-11+12z-1+12z-2

试画出其并联型网络结构。
34已知一FIR滤波器的系统函数为 H(z)=(1-0.7z-1+0.5z-2)(1+2z-1)，画出该FIR滤波器的网络结构。
35已知一个FIR系统的转移函数为

 H(z)=1+1.25z-1-2.75z-2-2.75z-3+1.25z-4+z-5

求用级联形式实现的线性相位结构流图，并用MATLAB画出其零点分布及其频率响应曲线。
36给定 |H(jΩ)|2=1/(1+64Ω6)，确定模拟滤波器的系统函数 H(s)。
37模拟低通滤波器的参数为： αp=3dB，αs=25dB，fp=25Hz，fs=50Hz，试用巴特沃斯近似求 H(s)。
38已知 Ha(s)=11+s/Ωc，使用脉冲响应不变法和双线性方法分别设计数字低通滤波器，使得3dB截止频率为 ωc= 0.25πrad。
39用脉冲响应不变法将 H(s)转换为 H(z)，采样周期为 T

H(s)=A(s-s0)m，其中 m为任意整数

310要求设计一个数字低通滤波器，在频率低于ω=0.2613πrad的范围内，低通幅度特性为常数，并且不低于0.75dB，在频率ω=0.4018π和 πrad之间，阻带衰减至少为20dB。试求出满足这些指标的最低阶巴特沃斯滤波器的传递函数H(z)，采用双线性变换。
311试设计一巴特沃斯数字低通滤波器，设计指标为： 在0.3π通带频率范围内，通带幅度波动小于1dB，在0.5π ～ πrad阻带频率范围内，阻带衰减大于12dB。
312用双线性变换法设计数字低通滤波器，等效模拟滤波器指标参数如下： 输入模拟信号 xa(t)的最高频率fd=100Hz； 选用巴特沃斯滤波器，3dB截止频率fc=100Hz，阻带截止频率fs=150Hz，阻带最小衰减αs=20dB。采样频率f=400Hz。
313试设计一个数字高通滤波器，要求通带下限频率 ωp=0.8πrad。阻带上限频率为 ωs=0.44πrad，通带衰减不大于3dB，阻带衰减不小于20dB。
314一个数字系统的采样频率Fs=2000Hz，试设计一个为此系统使用的带通数字滤波器H(z)，希望采用巴特沃斯滤波器，通带范围为300~400Hz，在带边频率处的衰减不大于3dB； 在200Hz以下和500Hz以上衰减不小于18dB。
315一个数字系统的采样频率为1000Hz,已知该系统受到频率为100Hz的噪声干扰，现设计一带阻滤波器 H(z)去掉该噪声。要求3dB
的带边频率为95Hz和105Hz，阻带衰减不小于14dB，阻带的下边和上边频率分别为99Hz和101Hz。
316试用矩形窗口设计法设计一个FIR线性相位低通数字滤波器，已知 ωc=0.5πrad，N=21。画出h(n)和20lgH(ω)H(0)曲线，再计算正、负肩峰值的位置和过渡带宽度。
317试用窗函数法设计一个第一类线性相位FIR数字高通滤波器，已知Hd(ejω)=e-jωα,3π4≤|ω|≤π,Hd(ejω)=0,0≤|ω|<3π4。对于矩形窗，过渡带宽度为Δω=π16。求： 
（1） h（n）的长度。
（2） h（n）的表达式。
（3） α。
318用矩形窗设计线性相位数字低通滤波器，理想滤波器传输函数Hd(ejω)为

 Hd(ejω)=

e-jωα,0≤|ω|≤ωc


0,ωc≤|ω|≤π



（1） 求出相应的理想低通滤波器的单位脉冲响应hd（n）。
（2） 求出用矩形窗函数法设计的FIR滤波器的h(n)表达式。
319用矩形窗设计线性相位高通滤波器，逼近滤波器传输函数Hd(ejω)为


Hd(ejω)＝

e-jωα,ωc≤|ω|<π


0,其他
 

(1) 求出相应于理想高通的单位脉冲响应hd（n）。
(2) 求出矩形窗设计法的h(n)表达式，确定α与N之间的关系。
(3) N的取值有什么限制？为什么？
320使用频率取样设计法（第一种形式取样）设计一个FIR线性相位低通数字滤波器。已知ωc=0.5πrad,N=51。
321用频率采样法设计第一类线性相位FIR低通滤波器，要求通带截止频率ωp=π3，阻带最大衰减25dB，过渡带宽度 Δω=π16，问滤波器长度至少为多少才可能满足要求？
322利用频率采样法设计线性相位FIR低通滤波器，设N＝16，给定希望逼近滤波器的幅度采样为

 Hg(k)＝

1,k=0,1,2,3

0.389，k=4

0,k=5,6,7


323一个IIR网络的差分方程为y(n)＝14y(n-1)+x(n)，当输入序列 x(n)＝12u(n)时：
(1) 试求在无限精度运算下网络输出 y(n),以及 n→∞时的输出稳态值。
(2) 当网络采用 b=4位字长的定点运算时，尾数采取截尾处理，试计算 0≤n≤20以内 21点输出值 y^(n)，并求其稳态响应 y^(∞)。
324在用模型表示数字滤波器中舍入和截尾效应时，把量化变量表示为 y(n)=Q［x(n)］=x(n)+e(n)，式中 Q［  ］表示舍入或截尾操作，e(n)表示量化误差。在适当的假定条件下，可以假设 e(n)是白噪声序列，即E［e(n)e(n+m)］ = σ2xδ(m)。舍入误差的一阶概率分别是如题324（a）图所示的均匀分布，截尾误差是如题324（b）图所示的均匀分布。


题324图


（1） 求输入噪声的均值和方差。
（2） 求截尾噪声的均值和方差。
325某因果LTI系统的系统函数为

H(z)=11-1.04z-1+0.98z-2

该系统稳定吗？若系统按“四舍五入”舍入，所得到的系统是稳定的吗？
326理想离散时间Hilbert变换器是一个 0<ω<π引入-90°相移，而对 -π<ω<0引入+90°相移的系统，对于 0<ω<π和 -π<ω<0，频响幅度为常量（单位1），这类系统也称为90°移相器。
（1） 给出一个理想离散时间Hilbert变换器的理想频率响应 Hd(ejω)的方程，该变换器还包括稳定（非零）群延迟，画出该系统对于 -π<ω<π的相位响应曲线。
（2） 可用哪类（Ⅰ、Ⅱ、Ⅲ、Ⅳ）FIR线性相位系统来逼近（1）中的理想Hilbert变换器？
（3） 假设我们要用窗函数法设计一个逼近理想Hilbert变换器的线性相位系统。若FIR是当n<0和n>M时，hd(n)=0，请用（1）中给出的 Hd(ejω)求理想脉冲响应 hd(n)。
（4） 当M=21时该系统的延迟是多少？若采用矩形窗，请画出在这种情况下的FIR逼近的频率响应之幅度曲线。