第3章数字滤波器的设计与实现
滤波器是用来对输入信号进行滤波的硬件或软件。如果滤波器的输入、输出都是离散信号，则该滤波器的脉冲响应也必然是离散的，这样的滤波器定义为数字滤波器。数字滤波器的功能，就是把输入序列通过一定的运算变换成输出序列。
数字滤波器在数字信号处理的各种应用中发挥着十分重要的作用，它是通过对抽样数据进行数学运算处理来达到频域滤波的目的。数学运算通常有两种实现方式，一种是频域法，即利用FFT快速算法对输入信号进行离散傅里叶变换，分析其频谱，然后根据所希望的频率特性进行滤波，再利用IFFT快速算法恢复出时域信号，这种方法具有较好的频率选择特性和灵活性，并且由于信号频率与所希望的频谱特性是简单的相乘关系，它比计算等价的时域卷积要快得多。另一种方法是时域法，这种方法是对离散抽样数据做差分数学运算来达到滤波目的。
一般用两种方法实现数字滤波器： 一是采用通用计算机，利用计算机的存储器、运算器和控制器把滤波器所要完成的运算编成程序通过计算机执行，也就是采用计算机软件实现； 二是设计专用的数字处理硬件。
数字滤波器用硬件实现的基本部件包括延迟器、乘法器和加法器； 用软件实现时，它只是一段线性卷积。软件实现的优点是系统函数具有可变性，仅依赖于算法结构，并且易于获得较理想的滤波性能。所以软件滤波在滤波器的使用中起到越来越重要的作用。

3.1IIR滤波器的基本结构

MATLAB的信号处理工具箱的两个基本组成就是滤波器的设计与实现部分以及谱分析部分。工具箱提供了丰富而简单的设计来实现FIR和IIR的方法，使原来烦琐的程序设计简化成函数的调用。本章首先介绍IIR数字滤波器的设计方法和相关的工具箱函数。
一个N阶无限长单位脉冲响应(IIR)数字滤波器的系统函数可以表示为

H(z)=Y(z)X(z)=∑Mk=0bkz－k1－∑Nk=1akz－k


对应的差分方程如下，其中ak和bk是滤波器系数。

y(n)=∑Mk=0bkx(n－k)+∑Nk=1aky(n－k)


实现同一个系统函数H(z)，可以用不同的结构形式，主要的结构形式有直接Ⅰ型、直接Ⅱ型、级联型和并联型等［1］。
1. 直接Ⅰ型
从差分方程可以看出，y(n)由两部分构成，第一部分∑Mk=0bkx(n－k)表示将输入信号进行延时，组成M阶的延时网络，实现了系统的零点； 第二部分∑Nk=1aky(n－k)表示将输出信号进行延时，组成N阶的延时网络，实现了系统的极点，信号流图如图3.1所示。

2. 直接Ⅱ型
还可以将系统函数H(z)改写为如下形式

H(z)=∑Mk=0bkz－k11－∑Nk=1akz－k=H1(z)H2(z)

将H1(z)和H2(z)互换，合并相同的延时支路，得到直接Ⅱ型结构，如图3.2所示。



图3.1直接Ⅰ型结构




图3.2直接Ⅱ型结构




对于N阶差分方程，直接Ⅰ型结构需要M+N个延时单元，而直接Ⅱ型只需要N个延时单元(一般情况下N≥M)，在软件实现时可以节省存储单元，在硬件实现时可以节省寄存器，因此直接Ⅱ型比直接Ⅰ型具有优势。
3. 级联型
一个N阶系统函数可以用它的零、极点表示，因为系统的系数均为实数，因此零极点为实数或者是共轭复数。可以将整个系统函数完全分解为实系数二阶因子的形式［1］，即

H(z)=A∏Nck=11+b1kz－1+b2kz－21－a1kz－1－a2kz－2

式中，Nc=floor(N/2)表示下取整，如floor(7/2)=3。图3.3所示为一个四阶IIR数字滤波器的级联结构。


图3.3四阶IIR数字滤波器的级联结构


级联结构的每一个基本节只是关系到数字滤波器的某一对极点和零点，因此调整对应的系数不会影响其他零极点。因而，级联结构的优点是便于准确地实现数字滤波器的零极点，也便于调整整个数字滤波器的性能，受参数量化影响较小，因此在实际应用中多选用这种级联结构。
4. 并联型
将因式分解后的系统函数展开成部分分式形式，并将实数极点成对组合，则系统函数H(z)可写成

H(z)=∑M－Nk=0Gkz－k+∑Npk=1e0k+e1kz－11－a1kz－1－a2kz－2

同级联型一样，Np=floor(N/2)。图3.4给出了M=N=4时的并联型结构。


图3.4四阶IIR数字滤波器的并联型结构


在并联型结构中，改变系数可以调整极点的位置，但不能像级联型那样可以直接控制零点。在运算误差方面，由于并联型各基本环节的误差互不影响，因此误差比级联型的稍小些。
例3.1自定义函数tf2par和cplxcomp来自陈怀琛.数字信号处理教程——MATLAB释义与实现（第3版）.北京： 电子工业出版社，2013。
已知IIR滤波器的系统函数如下所示，请给出该滤波器直接Ⅱ型、级联型和并联型结构。

H(z)=6+1.2z－1－0.72z－2+1.728z－38－10.4z－1+7.28z－2－2.352z－3

解： 根据系统函数可知系统差分方程为

y(n)=1.3y(n－1)－0.91y(n－2)+0.294y(n－3)+0.75x(n)+

0.15x(n－1)－0.09x(n－2)+0.216x(n－3)

由差分方程可以直接画出滤波器直接Ⅱ型结构如图3.5所示。


图3.5直接Ⅱ型结构



函数tf2sos实现直接型到级联型的转换，函数名tf2sos表示transfer function to secondorder sections，意思是把传递函数表示成二阶节相乘的形式，MATLAB源代码如下（见chp3sec1_1.m）。



% chp3sec1_1.m IIR滤波器直接型转级联型

% 函数tf2sos: transfer function to second-order sections

clear;

clc;

close all;

b=［6,1.2,-0.72,1.728］;%传递函数分子

a=［8,-10.4,7.28,-2.352］; %传递函数分母

［sos,g］=tf2sos(b,a)




程序运行结果如下，g表示级联型结构的增益，sos存储的是级联型结构的系数，每一行表示一个二阶节信息，分子（零点）系数在前，分母（极点）系数在后。



sos =

1.00000.800001.0000-0.60000

1.0000-0.60000.36001.0000-0.70000.4900

g =

0.7500




将结果写成系统函数的形式为

H(z)=34·1+0.8z－11－0.6z－1·1－0.6z－1+0.36z－21－0.7z－1+0.49z－2

级联型滤波器结构如图3.6所示。



图3.6级联型结构



自定义函数tf2par实现直接型到级联型的转换，MATLAB源代码如下（见chp3sec1_2.m）。



% chp3sec1_2.m IIR滤波器直接型转并联型

clear;

clc;

close all;

b=［6,1.2,-0.72,1.728］;%传递函数分子

a=［8,-10.4,7.28,-2.352］; %传递函数分母

［C,B,A］ = tf2par(b,a)




程序运行结果如下，C表示并联型结构的常数项，B表示各部分的分子（零点）系数，A表示对应的分母（极点）系数。



C =

-0.7347

B =

0.01960.2322

1.46510

A =

1.0000-0.70000.4900

1.0000-0.60000




将结果写成系统函数的形式为

H(z)=－0.7347+1.46511－0.6z－1+0.0196z－1+0.23221－0.7z－1+0.49z－2

并联型滤波器结构如图3.7所示。


图3.7并联型结构


在自定义函数tf2par源代码如下。



function ［C,B,A］ = tf2par(b,a);

%  C = 当length(b) >= length(a)时的多项式直通部分

%  B = 包含各bk的K×2维实系数矩阵

%  A = 包含各ak的K×3维实系数矩阵

%  b = 直接型的分子多项式系数

%  a = 直接型的分母多项式系数








%

M = length(b); N = length(a);

［r1,p1,C］ = residuez(b,a);% 先求系统的单根p1,对应的留数r1及直接项C

p = cplxpair(p1,1e-9);	    % 用配对函数cplxpair由p1找共轭复根p，1e-9为容差

I = cplxcomp(p1,p);	        % 找p1变为p时的排序变化，MATLAB无此子程序，列于后面

r = r1(I);　　% 让r1的排序同样变化成为r,以保持与极点对应

% 变为二阶子系统

K = floor(N/2); B = zeros(K,2); A = zeros(K,3);% 二阶子系统变量的初始化

if K*2 == N; 　　%  N 为偶, A(z)的次数为奇,有一个因子是一阶的

for i=1:2:N-2

pi = p(i:i+1,:); 　　% 取出一对极点

ri = r(i:i+1,:);　　% 取出一对对应的留数

［Bi,Ai］ = residuez(ri,pi,［］);　　% 二个极点留数转为二阶子系统分子分母系数

B(fix((i+1)/2),:) = real(Bi);　　% 取Bi的实部，放入系数矩阵B的相应行中

A(fix((i+1)/2),:) = real(Ai); 　　% 取Ai的实部，放入系数矩阵A的相应行中

end

［Bi,Ai］ = residuez(r(N-1),p(N-1),［］);　　% 处理实单根

B(K,:) = ［real(Bi) 0］; A(K,:) = ［real(Ai) 0］;

else　　%  N 为奇, A(z)的次数为偶,所有因子都是二阶的

for i=1:2:N-1

pi = p(i:i+1,:); 　　% 取出一对极点

ri = r(i:i+1,:);　　% 取出一对对应的留数

［Bi,Ai］ = residuez(ri,pi,［］); 　　% 两个极点留数转为二阶子系统分子分母系数

B(fix((i+1)/2),:) = real(Bi); 　　% 取Bi的实部，放入系数矩阵B的相应行中

A(fix((i+1)/2),:) = real(Ai); 　　% 取Ai的实部，放入系数矩阵A的相应行中

end

end




tf2par还用到了自定义函数cplxcomp，该函数用于计算复数极点p1变为p2后留数的新序号，源代码如下。



function I = cplxcomp(p1,p2)

%  计算复数极点p1变为p2后留数的新序号

%  本程序必须用在cplxpair程序之后以便重新确定频率极点向量

%  及其相应的留数向量:

% p2 = cplxpair(p1)

%

I=［］;    % 设一个空的矩阵

for j=1:1:length(p2)                % 逐项检查改变排序后的向量p2

for i=1:1:length(p1)            % 把该项与p1中各项比较

if (abs(p1(i)-p2(j)) < 0.0001)  % 看与哪一项相等

I=［I,i］;                     % 把此项在p1中的序号放入I

end

end

end

I=I';                       % 最后的I表示了p2中各元素在p1中的位置





3.2IIR滤波器的设计

设计典型IIR数字滤波器的步骤如下： 
(1) 按一定规则将给出的数字滤波器的技术指标转换为模拟低通滤波器的技术指标； 
(2) 根据转换后的技术指标使用滤波器阶数选择函数，确定最小阶数N和固定频率Wn； 
(3) 运用最小阶数N产生模拟滤波器原型； 
(4) 运用固定频率Wn把模拟低通滤波器原型转换成模拟低通、高通、带通、带阻滤波器； 
(5) 运用冲激响应不变法或双线性变换法把模拟滤波器转换成数字滤波器。
另外，MATLAB信号处理工具箱提供了几个直接设计IIR数字滤波器的函数，它们把以上几个步骤集成一个整体，这为设计通用滤波器带来了极大的方便。这些函数应与数字滤波器的阶数选择函数配合使用，为特定的滤波器的设计返回所需的阶数N和固有频率Wn。
对于一般条件下使用的滤波器，运用典型设计的方法即可满足滤波器性能的要求。但是，对一些有特殊要求的滤波器，如对于滤波器的过渡带的宽度有特殊要求，典型设计的方法就不宜使用。
3.2.1模拟原型滤波器设计

模拟低通滤波器的振幅响应|H(Ω)|和设计规定如图3.8所示。通带为0~Ωp，通带的幅度响应规定为1，允许的波动为δ1。阻带为Ωs~∞，要求的阻带衰减为δ2。Ωp~Ωs称为过渡带，在过渡带内振幅响应不做明确规定。习惯上称频率Ωp为通带截止频率，称Ωs为阻带截止频率。为了方便，有时把通带截止频率归一化为1，这时滤波器成为标准形式。


图3.8模拟低通滤波器的设计规定


指定低通滤波器的性能可用频率响应的幅值的平方表示。
设定δ1=1－11+ε2，δ2=1A2，则

11+ε2≤|H(Ω)|2≤1,|Ω|<Ωp

0≤|H(Ω)|2≤1A2,Ωs≤|Ω|

此时可得

Rp=－10lg11+ε2，ε=10Rp10－1

As=－10lg1A2，A=10As20

1. Butterworth(巴特沃斯)滤波器

Butterworth滤波器的特点是具有通带内最平坦的幅度特性，而且随着频率升高呈单调递减，因此Butterworth滤波器又称为“最平”的幅频响应滤波器，而且Butterworth滤波器也是最简单的滤波器。
Butterworth低通滤波器是将Butterworth函数作为滤波器的传输函数，它的平方幅频响应函数可以表示为

|G(jΩ)|2=11+C2(Ω2)N

上面的表达式中含有两个参数C和N，其中N表示滤波器的阶数。
在设计时，当ap=3dB，一般有C=1，此时可简化为

|G(jΩ)|2=11+ΩΩp2N



MATLAB工具箱中提供了buttap函数来产生低通模拟Butterworth滤波器，调用格式为［z,p,k］=buttap(n)，n表示滤波器阶数，z，p，k分别表示滤波器的零点、极点和增益。下面通过一个例子具体演示。
例3.2绘制Butterworth低通模拟原型滤波器的平方幅度响应曲线，其中滤波器的阶数分别为2，5，10，20。
解： MATLAB源代码如下(见chp3sec2_1.m)，实验结果如图3.9所示。



%Butterworth滤波器 chp3sec2_1.m

clear;

clc;

n=0:0.01:2;

for i=1:4

switch i

case 1

N=2;

case 2

N=5;

case 3

N=10;

case 4

N=20;

end

［z,p,k］=buttap(N);

［b,a］=zp2tf(z,p,k);

［H,w］=freqs(b,a,n);

magH2=(abs(H)).^2;







hold on;

plot(w,magH2);axis(［0 2 0 1］);

end

xlabel('w');ylabel('|H(jw)|^2');

title('Butterworth模拟原型滤波器平方幅度响应');

grid on;






图3.9Butterworth滤波器幅频响应


Butterworth滤波器除了具有平滑单调递减的频率响应优点外，其过渡带的陡峭程度正比于滤波器的阶数，高阶的Butterworth滤波器的频率响应近似于理想的低通滤波器。
2. Chebyshev(切比雪夫)Ⅰ型滤波器
Butterworth滤波器的典型频率特性为： 无论是在通带还是在阻带，幅值都随频率而单调变化，因而如果在通带边缘满足指标，则在通带内肯定会有富余量，也就是会超过指标的要求，因而并不经济。所以更有效的办法是将指标的精度要求均匀分布在通带内，或均匀分布在阻带内，或同时均匀分布在通带与阻带内，这时就可设计出阶数较低的滤波器。这种精度均匀分布的办法可以通过选择具有等波纹特性的逼近函数来完成。
Chebyshev滤波器的幅度特性就是在一个频带中(通带或阻带)具有这种等波纹特性： 一种是在通带中是等波纹的，在阻带中是单调的，称为Chebyshev Ⅰ型； 一种是在通带内是单调的，在阻带内是等波纹的，称为Chebyshev Ⅱ型。
Chebyshev Ⅰ型滤波器的特性函数为

|Hs(jΩ)|2=11+ε2C2NΩΩc

其中，ε为小于1的正数，它是表示通带波纹大小的一个参数，ε越大波纹也越大。Ω/Ωc为Ω对Ωc的归一化频率，Ωc为截止频率，也是滤波器的某一衰减分贝处的通带宽度。CN(x)是N阶的Chebyshev多项式，可以直接查表计算［2］。Chebyshev Ⅰ型滤波器的主要特点是在阻带内达到最大平滑。
MATLAB工具箱中提供了cheb1ap函数来产生低通模拟Chebyshev Ⅰ型滤波器，调用格式为［z,p,k］=cheb1ap(n,Rp)，n表示滤波器阶数，z，p，k分别表示滤波器的零点、极点和增益，Rp表示通带内的最大衰减。下面通过一个例子具体演示。
例3.3产生一个15阶的Chebyshev Ⅰ型低通模拟滤波器原型，绘制其平方幅度响应曲线，要求通带内的最大衰减为0.3dB。
解： MATLAB源代码如下(见chp3sec2_2.m)，实验结果如图3.10所示。



%Chebyshev Ⅰ型 chp3sec2_2.m

clear;

clc;

n=0:0.01:2;

［z,p,k］=cheb1ap(15,0.3);

［b,a］=zp2tf(z,p,k);

［H,w］=freqs(b,a,n);

magH2=(abs(H)).^2;

plot(w,magH2);grid on;

axis(［0 2 0 1.2］);

xlabel('w');ylabel('|H(jw)|^2');

title('Chebyshev Ⅰ型低通模拟滤波器平方幅度响应');






图3.10Chebyshev Ⅰ型滤波器幅频响应


3. Chebyshev Ⅱ型滤波器
Chebyshev Ⅱ型滤波器的特性函数为

|Ha(jΩ)|2=11+ε2CN(Ωst)CNΩstΩ2

其中，Ωst为阻带衰减达到规定数值的最低频率。Chebyshev Ⅱ型滤波器的主要特点就是在通带内达到最大平滑。
MATLAB工具箱中提供了函数cheb2ap来产生低通模拟Chebyshev Ⅱ型滤波器，调用格式为［z,p,k］=cheb2ap(n,Rs)，n表示滤波器阶数，z，p，k分别表示滤波器的零点、极点和增益，Rs表示阻带内的最大衰减。下面通过一个例子进行具体的演示。
例3.4设计一个10阶的Chebyshev Ⅱ型低通模拟滤波器原型，绘制其平方幅度响应曲线，要求阻带内的最大衰减为50dB。
解： MATLAB源代码如下(见chp3sec2_3.m)，实验结果如图3.11所示。



%Chebyshev Ⅱ型 chp3sec2_3.m

clear;

clc;

n=0:0.01:2;

［z,p,k］=cheb2ap(10,50);

［b,a］=zp2tf(z,p,k);

［H,w］=freqs(b,a,n);

magH2=(abs(H)).^2;

plot(w,magH2);grid on;

axis(［0 2 0 1.2］);

xlabel('w');ylabel('|H(jw)|^2');

title('Chebyshev Ⅱ型低通模拟滤波器平方幅度响应');






图3.11Chebyshev Ⅱ型滤波器幅频响应


4. 椭圆滤波器
椭圆滤波器的特性函数为

|H(jΩ)|2=11+ε2U2NΩΩc

其中，UN(x)为N阶雅可比椭圆函数。
MATLAB工具箱中提供了函数ellipap来产生低通模拟椭圆滤波器，调用格式为［z,p,k］=ellipap(n,Rp,Rs)，n表示滤波器阶数，z，p，k分别表示滤波器的零点、极点和增益，Rp表示通带内的最大衰减，Rs表示阻带内的最小衰减。下面通过一个例子具体演示。
例3.5设计一个三阶低通模拟椭圆滤波器，绘制其平方幅度响应曲线，要求通带内的最大衰减为3dB，阻带内的最小衰减为40dB。
解： MATLAB源代码如下(见chp3sec2_4.m)，实验结果如图3.12所示。



%椭圆滤波器 chp3sec2_4.m

clear;

clc;

n=0:0.01:2;

［z,p,k］=ellipap(3,3,40);

［b,a］=zp2tf(z,p,k);

［H,w］=freqs(b,a,n);

magH2=(abs(H)).^2;

plot(w,magH2);grid on;

axis(［0 2 0 1.2］);

xlabel('w');ylabel('|H(jw)|^2');

title('椭圆低通模拟滤波器平方幅度响应');






图3.12椭圆滤波器幅频响应



3.2.2频率变换
归一化模拟原型滤波器设计完成之后，就可以进行典型滤波器设计的第二个步骤，通过频率转换成所需要类型(低通、高通、带通、带阻)的模拟滤波器。也就是说，设计模拟低通、高通、带通及带阻滤波器，要先将其技术指标通过某种频率转换方法转换成模拟低通滤波器的技术指标，并依据这些技术指标设计出低通滤波器的转移函数，然后再依据频率转换关系变成所要设计的滤波器的转移函数。

符号H1(z)表示原型低通数字滤波器的系统函数，Hd(p)表示所希望的任何类型数字滤波器系统函数，如果能找到变换关系

z－1=f(p－1)

就可以将已知的H1(z)变换到Hd(z)，即

Hd(p)=H1(z)z－1=f(p－1)=H1［f－1(p－1)］


为使一个稳定、因果的原型滤波器H1(z)变换成一个稳定、因果和有理的所需类型的滤波器Hd(p)，必须满足两个要求： ①函数f(p－1)是p-1或p的有理函数； ②z平面单位圆必须映射为p平面单位圆。如令z=ejθ，p=ejω，则必须有

e－jθ=f(e－jω)

或者说

|f(e－jω)|=1

θ=－arg［f(e－jω)］


上述方程规定了p平面和z平面之间的频率关系。可以证明满足上述要求的函数f(z－1)的最一般的形式为

f(p－1)=ejθ0∏Nk=1p－1－ak1－a*kp－1

式中,a*k为ak的复共轭，对于实滤波器，ak为实数或共轭成对出现，且

θ0=±π,ejθ0=±1

代入可得

f(p－1)=±∏Nk=1p－1－ak1－a*kp－1


为使滤波器稳定，|a|必须小于1。通过选择适当的N值和常数ak，可以得到一组映射关系。最简单的情况是由模拟低通滤波器转换为数字低通滤波器的关系 

z－1=f(p－1)=p－1－a1－ap－1

代入z=ejθ，p=ejω可得

e－jθ=e－jω－a1－ae－jω

代入上式可解得

ω=arctan(1－a2)sinθ2a+(1+a2)cosθ

将θc和ωc代入上式，可以确定参数a

a=sinθc－ωc2sinθc+ωc2

如果已知θc和ωc，可以求得a，或者已知θc和a，可以求得ωc。
利用上述方法由一个已知的模拟低通滤波器H1(z)得到所需的Hd(p)，H1(z)和Hd(p)的截止频率分别为θc和ωc，用前面的方法求得a，并代入下式，

Hd(p)=H1(z)z－1=p－1－a1－ap－1


用类似的方法，还可以得到低通到高通，低通到带通，低通到带阻的对应变换关系，具体见表3.1，表中θc始终表示低通原型滤波器截止频率［1,2］。


表3.1根据模拟低通滤波器原型设计各类滤波器



滤波器类型变换函数公式相应的设计公式

低通原型→低通z－1=p－1－a1－ap－1a=sinθc－ωc2sinθc+ωc2

ωc： 变换后低通滤波器截止频率
续表


滤波器类型变换函数公式相应的设计公式

低通原型→高通z－1=－p－1+a1+ap－1a=－cosθc－ωc2cosθc+ωc2

ωc： 变换后高通滤波器截止频率
低通原型→带通z－1=－p－2－2akk+1p－1+k－1k+1k－1k+1p－2－2akk+1p－1+1a=cosωc2+ωc12cosωc2－ωc12

k=tanθc2tanωc2+ωc12

ωc1和ωc2： 变换后带通滤波器的上、下截止频率

低通原型→带阻z－1=－p－2－2akk+1p－1+k－1k+1k－1k+1p－2－2akk+1p－1+1a=cosωc2+ωc12cosωc2－ωc12

k=tanθc2tanωc2+ωc12

ωc1和ωc2： 变换后带阻滤波器的上、下截止频率


MATLAB的信号处理工具箱为实现从低通滤波器向低通、高通、带通和带阻滤波器的转换提供了非常方便、简洁的函数。
1. 从低通到低通的转换lp2lp
lp2lp函数可将截止频率为1rad/s的模拟低通滤波器原型转换为截止频率为Wn的低通滤波器。lp2lp函数有两种表示形式： 传递函数形式和状态空间形式，但是输入系统必须为模拟滤波器原型。
［bt,at］=lp2lp(b,a,Wn)可以将传递函数表示的模拟低通滤波器原型转换成低通滤波器，其截止频率为Wn，bt和at表示变换之后的低通滤波器传递函数的系数，b和a表示滤波器原型传递函数的系数，定义如下

H(s)=b(s)a(s)=b(1)sn+…+b(n)s+b(n+1)a(1)sm+…+a(m)s+a(m+1)


［At,Bt,Ct,Dt］=lp2lp(A,B,C,D,Wn)可以将以连续方程表示的低通滤波器原型转换为截止频率为Wn的低通滤波器，At，Bt，Ct，Dt表示变换之后的状态空间参数，A，B，C，D分别表示滤波器原型的状态空间参数，定义如下

x′=Ax+Bu


y′=Cx+Du


2. 从低通到高通的转换lp2hp
lp2hp函数可将截止频率为1rad/s的模拟低通滤波器原型转换为截止频率为Wn的高通滤波器。lp2hp函数有两种表示形式： 传递函数形式和状态空间形式，但是输入系统必须为模拟滤波器原型。
［bt,at］=lp2hp(b,a,Wn)可以将传递函数表示的模拟低通滤波器原型转换成高通滤波器，其截止频率为Wn，bt和at表示变换之后的高通滤波器传递函数的系数，b和a表示滤波器原型传递函数的系数，具体定义同lp2lp函数。
［At,Bt,Ct,Dt］=lp2hp(A,B,C,D,Wn)可以将以连续方程表示的低通滤波器原型转换为截止频率为Wn的高通滤波器，At，Bt，Ct，Dt表示变换之后的状态空间参数，A，B，C，D分别表示滤波器原型的状态空间参数，具体定义同lp2lp函数。
3. 从低通到带通的转换lp2bp
lp2bp函数可将截止频率为1rad/s的模拟低通滤波器原型转换为具有指定带宽Bw和中心频率为Wn的带通滤波器。lp2bp函数有两种表示形式： 传递函数形式和状态空间形式，但是输入系统必须为模拟滤波器原型。
［bt,at］=lp2bp(b,a,Wn,Bw)可以将传递函数表示的模拟低通滤波器原型转换成带通滤波器，其中心频率为Wn，带宽为Bw，bt和at表示变换之后的带通滤波器传递函数的系数，b和a表示滤波器原型传递函数的系数，具体定义同lp2lp函数。
如果要求滤波器的低端截止频率为w1，高端截止频率为w2，则可计算出Wn和Bw为

Wn=w1w2

Bw=w2－w1


［At,Bt,Ct,Dt］=lp2bp(A,B,C,D,Wn,Bw)可以将以连续方程表示的低通滤波器原型转换为中心频率为Wn的带通滤波器，带宽为Bw。At、Bt、Ct、Dt表示变换之后的状态空间参数，A、B、C、D分别表示滤波器原型的状态空间参数，具体定义同lp2lp函数。
4. 从低通到带阻的转换lp2bs
lp2bs函数可将截止频率为1rad/s的模拟低通滤波器原型转换为具有指定带宽Bw和中心频率为Wn的带阻滤波器。lp2bs函数有两种表示形式： 传递函数形式和状态空间形式，但是输入系统必须为模拟滤波器原型。
［bt,at］=lp2bs(b,a,Wn,Bw)可以将传递函数表示的模拟低通滤波器原型转换成带阻滤波器，其中心频率为Wn，带宽为Bw，bt和at表示变换之后的带通滤波器传递函数的系数，b和a表示滤波器原型传递函数的系数，具体定义同lp2lp函数。
如果要求滤波器的低端截止频率为w1，高端截止频率为w2，则Wn和Bw的计算与函数lp2bp一致。
［At,Bt,Ct,Dt］=lp2bs(A,B,C,D,Wn,Bw)可以将以连续方程表示的低通滤波器原型转换为中心频率为Wn的带阻滤波器，带宽为Bw。At、Bt、Ct、Dt表示变换之后的状态空间参数，A、B、C、D分别表示滤波器原型的状态空间参数，具体定义同lp2lp函数。
可以看出，4种转换函数lp2lp，lp2hp，lp2bp和lp2bs的定义和用法都非常相似，也非常便于记忆： l表示“低”(low)，h表示“高”(high)，b表示“带”(band)，p表示“通”(pass)，s表示“阻”(stop)，2(two)就是英文to的谐音表示。下面通过一个例子进行具体演示。
例3.6设计一个三阶模拟低通滤波器，通带内最大衰减为3dB，在阻带内的最大衰减为40dB，截止频率为8πrad，再把它转换为截止频率为50πrad的高通滤波器，绘出它们的频率响应图。
解： MATLAB源代码如下(见chp3sec2_5.m)，实验结果如图3.13和图3.14所示。



%频率变换 chp3sec2_5.m

clear;

clc;

［z,p,k］=ellipap(3,3,40); %模拟低通滤波器

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

［At,Bt,Ct,Dt］=lp2lp(A,B,C,D,8*pi); %低通->低通

［num,den］=ss2tf(At,Bt,Ct,Dt);

freqs(num,den); title('低通转低通');

figure;

［At1,Bt1,Ct1,Dt1］=lp2hp(A,B,C,D,50*pi); %低通->高通

［num1,den1］=ss2tf(At1,Bt1,Ct1,Dt1);

freqs(num1,den1); title('低通转高通');






图3.13三阶椭圆低通滤波器的频率响应




图3.14三阶椭圆高通滤波器的频率响应


在上面这段程序中用到了一次模拟滤波器生成(ellipap)，三次线性系统表示形式的转换(zp2ss和ss2tf)，以及两次滤波器类型的转换(lp2lp和lp2hp)。

3.2.3滤波器的映射
1. 冲激响应不变法

冲激响应不变法的基本原理是从滤波器的冲激响应出发，对具有传递函数H(s)的模拟滤波器的冲激响应h(t)，以周期T抽样所得到的离散序列h(nT)为数字滤波器的冲激响应。
对下式

h(t)=∑∞n=0h(nT)δ(t－nT)

进行拉普拉斯变换，得到

H(s)=∑∞n=0h(nT)e－nTs


由于H(s)是esT的函数，令z=esT可得

H(z)=∑∞n=0h(nT)z－n

这样就得到了从s域到z域的变换。这实际上是拉普拉斯变换到z变换的标准变换，或者称为冲激响应不变的变换。上式是非递归的，但是如果H(s)是s的有理函数，即

H(s)=∑Mk=0cksk∑Nk=1dksk,M<N

则根据冲激响应不变的变换原理，H(z)也是一个有理函数，这样就可以用递归形式实现IIR滤波器。现在，将H(s)用部分分式展开，并设H(s)无重极点，则

H(s)=∑Nk=1Aks－spk

式中,Ak是s=spk的留数，即Ak=lims=spk(s－spk)H(s)，此时模拟滤波器H(s)的冲激响应是

h(t)=∑Nk=1AkespkTt≥0
0t<0

因此，可得数字滤波器的系统函数为

H(z)=∑∞n=0h(nT)z－n=∑Nk=1Ak1－espkTz－1


通过上述过程可以看出： H(z)是H(s)通过下面的对应关系得到的

1s－spk11－espkTz－1


冲激响应不变法的特点为： ①模拟频率与数字频率之间的转换是线性的，并保持了模拟滤波器的时域瞬态特性，这是冲激响应不变法的优点； ②当模拟滤波器频率响应不是严格带限时，用冲激响应不变法设计出的数字滤波器在频域出现混叠现象，这是冲激响应不变法的缺点。故设计的滤波器性能要求比较高时，不宜使用该方法。

MATLAB信号处理工具箱提供了函数impinvar，可采用冲激响应不变法实现模拟滤波器到数字滤波器的转换。［bz,az］=impinvar(b,a,fs)可以将模拟滤波器(b,a)变换为数字滤波器(bz,az)，两者的冲激响应不变。即模拟滤波器的冲激响应按fs抽样后等同于数字滤波器的冲激响应，输入参数fs的默认值为1Hz。下面通过一个例子进行具体的演示。
例3.7设Butterworth低通模拟原型滤波器的通带波纹Rp=1dB，通带上限角频率为ωp=0.2π，下限频率为ωs=0.3π，阻带的最小衰减为15dB。根据该低通模拟滤波器，利用冲激响应不变法设计相应的数字低通滤波器，并画出设计后的数字低通滤波器的特性曲线(包括幅频、相频、群延迟)。

解： 数字滤波器的单位抽样响应h(n)恰好等于h(t)的采样值，即

h(n)=h(t)t=nTs=h(t)∑∞n=0δ(t－nTs)=h(nTs)

其中Ts为采样周期，对应有

H(ejω)=1Ts∑∞k=－∞H(jΩ－jkΩs)

其中ω=ΩTs。

MATLAB源代码如下(见chp3sec2_6.m)，实验结果如图3.15所示。



%冲激响应不变法设计数字滤波器 chp3sec2_6.m

clear;

clc;

wp=0.2*pi;

ws=0.3*pi;

Rp=1;

As=15;

T=1;

Rip=10^(-Rp/20);

Atn=10^(-As/20);

OmgP=wp*T;

OmgS=ws*T;

［N,OmgC］=buttord(OmgP,OmgS,Rp,As,'s');

［cs,ds］=butter(N,OmgC,'s');

［b,a］=impinvar(cs,ds,T);

［dB,mag,pha,grd,w］=freqz_m(b,a);

%绝对幅频响应

subplot(2,2,1);

plot(w/pi,mag),title('幅频特性'),xlabel('w/＼pi'),ylabel('|H(jw)|');

axis(［0 1 0 1.1］);grid on;

set(gca,'XTickMode','manual','XTick',［0 0.2 0.3 0.5 1］);

set(gca,'YTickMode','manual','YTick',［0 Atn Rip 1］);

%用dB表示的幅频响应

subplot(2,2,2);

plot(w/pi,dB),title('幅频特性/dB'),xlabel('w/＼pi'),ylabel('|H(jw)|');

axis(［0 1 -40 5］);grid on;

set(gca,'XTickMode','manual','XTick',［0 0.2 0.3 0.5 1］);

set(gca,'YTickMode','manual','YTick',［-40 -As -Rp 0］);







%相位响应

subplot(2,2,3);

plot(w/pi,pha/pi),title('相频特性'),xlabel('w/＼pi'),ylabel('pha(＼pi)');

axis(［0 1 -1 1］);grid on;

set(gca,'XTickMode','manual','XTick',［0 0.2 0.3 0.5 1］);

%群延迟

subplot(2,2,4);

plot(w/pi,grd),title('群延迟'),xlabel('w/＼pi'),ylabel('Group Delay');

axis(［0 1 0 12］);grid on;

set(gca,'XTickMode','manual','XTick',［0 0.2 0.3 0.5 1］);






图3.15采用冲激响应不变法得到的数字低通滤波器特性曲线


其中，freqz_m为自定义函数，源代码如下。



function［dB,mag,pha,grd,w］=freqz_m(b,a)

%滤波器的幅频响应、相位响应和群延迟

%dB 用dB表示的幅频响应

%mag 幅频响应的绝对值

%pha 相位响应

%grd 群延迟

%b,a 系统传递函数对应的系数

%w 采样频率

［H,w］=freqz(b,a,1000,'whole');

H=(H(1:501))';

w=(w(1:501))';

mag=abs(H);

dB=20*log10((mag+eps)/max(mag));

pha=angle(H);

grd=grpdelay(b,a,w);





2. 双线性变换法

为了克服冲激响应不变法产生的频率混叠现象，需要使s平面与z平面建立一一映射关系，即求出s=f(z)，然后将它代入H(s)就可以求得H(z)，即

H(z)=H(s)s=f(z)


为了得到s=f(z)的函数关系，先从模拟滤波器的传递函数H(s)开始

H(s)=∑Mk=0cksk∑Nk=1dksk,M<N

将H(s)展开为部分分式，设无重极点，则

H(s)=∑Nk=1Aks－spk

上式意味着模拟输入x(t)和输出y(t)之间有如下关系 

y′(t)－spy(t)=Ax(t)

用［y(n)－y(n－1)］/T代替y′(t)，［y(n)－y(n－1)］/2代替y(t)，［x(n)－x(n－1)］/2代替x(t)，则可写出对应的差分方程 

1Ty(n)－y(n－1)－sp2y(n)－y(n－1)=A2x(n)－x(n－1)

对两边取z变换可得

H(z)=Y(z)X(z)=A2T1－z－11+z－1－sp

从上式可得到双线性变换的基本关系如下 

s=2T1－z－11+z－1
z=2T+s2T－s



MATLAB信号处理工具箱提供了函数bilinear实现双线性变换。
［zd,pd,kd］=bilinear(z,p,k,fs)把模拟滤波器的零极点模型转换为数字滤波器的零极点模型，其中fs是抽样频率。
［numd,dend］=bilinear(num,den,fs)把模拟滤波器的传递函数模型转换为数字滤波器的传递函数模型，其中fs是抽样频率。
［Ad,Bd,Cd,Dd］=bilinear(A,B,C,D,fs)把模拟滤波器的状态方程模型转换为数字滤波器的状态方程模型，其中fs是抽样频率。
以上三种形式，都可以增加一个参数fp(预畸变频率)，在进行双线性变换之前，对抽样频率进行预畸变，以保证频率冲激响应在双线性变换前后具有良好的单值映射关系。下面通过两个例子进行具体的演示。
例3.8设计一个数字信号处理系统，抽样频率Fs=100Hz，在该系统设计一个Butterworth型高通滤波器，使其在通带内允许的最小衰减为0.5dB，阻带内的最小衰减为40dB，通带上限临界频率为30Hz，阻带下限临界频率为40Hz。
解： MATLAB源代码如下(见chp3sec2_7.m)，实验结果如图3.16所示。



%双线性变换法 chp3sec2_7.m

clear;

clc;

wp=30*2*pi;

ws=40*2*pi;

rp=0.5;

rs=40;

Fs=100;

［N,Wn］=buttord(wp,ws,rp,rs,'s'); %滤波器阶数选择

［z,p,k］=buttap(N); %Butterworth低通滤波器原型

［A,B,C,D］=zp2ss(z,p,k); %零极点->状态方程

［At,Bt,Ct,Dt］=lp2hp(A,B,C,D,Wn); %低通->高通

［num1,den1］=ss2tf(At,Bt,Ct,Dt);

［num2,den2］=bilinear(num1,den1,100); %双线性变换

［H,W］=freqz(num2,den2);

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

grid on;

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






图3.16Butterworth型高通滤波器的频率响应


例3.9试用双线性变换法设计一个Chebyshev Ⅱ型高通滤波器，使其幅频特性逼近一个具有以下技术指标的模拟Chebyshev Ⅱ型高通滤波器： ωs=2π×1kHz，ωp=2π×1.4kHz，在ωs处的最小衰减为1.5dB，在ωp处的最大衰减不超过0.3dB，抽样频率为20kHz。
解： MATLAB源代码如下(见chp3sec2_8.m)，实验结果如图3.17所示。



%双线性变换法 chp3sec2_8.m

clear;

clc;

ws=2*pi*1000; ws1=ws*2*pi;

wp=2*pi*1400; wp1=wp*2*pi;

rp=0.3;rs=15;Fs=20000;







［N,Wn］=cheb2ord(wp1,ws1,rp,rs,'s');

［z,p,k］=cheb2ap(N,rs);

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

［At,Bt,Ct,Dt］=lp2hp(A,B,C,D,Wn);

［At1,Bt1,Ct1,Dt1］=bilinear(At,Bt,Ct,Dt,Fs);

［num,den］=ss2tf(At1,Bt1,Ct1,Dt1);

［H,W］=freqz(num,den);

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

grid on;

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







图3.17Chebyshev Ⅱ型高通滤波器的频率响应


下面再通过一个例子比较一下冲激响应不变法和双线性变换法设计数字滤波器。
例3.10请分别采用冲激响应不变法和双线性变换法设计一个IIR低通滤波器，要求以小于1dB的衰减通过100Hz信号，以大于40dB的衰减抑制150Hz的信号，采样率为1000Hz。
解： 冲激响应不变法是使数字滤波器的单位冲激响应序列h(n)模仿模拟电路滤波器的冲激响应ha(t)。将模拟滤波器的冲激响应加以等间隔抽样，使h(n)正好等于抽样值，即满足： 

h(n)=ha(nT)

其中，T是抽样周期。
双线性变换法是使数字滤波器的频率响应与模拟滤波器的频率响应相似的一种变换方法。先把整个s平面压缩变换到某一中介平面s1的一条横带里（宽度为2π/T），再通过标准变换关系z=es1T将此横带变换到整个z平面上去。
MATLAB源代码如下（见chp3sec2_9.m），实验结果如图3.18所示。



% 冲激响应不变法和双线性变换法的比较 chp3sec2_9.m

clear;

clc;








close all;

%% 模拟滤波器技术指标

Fs=1*1e3;

OmegaP=2*pi*100;

OmegaS=2*pi*150;

Rp=1; %通带最大衰减

As=40; %阻带最小衰减

%% 数字滤波器技术指标

wp=OmegaP/Fs;

ws=OmegaS/Fs;

w0=［wp,ws］; %数字临界频率

%% 模拟滤波器设计(冲激响应不变法)

% 巴特沃斯型

% ［N,OmegaC］=buttord(OmegaP,OmegaS,Rp,As,'s');

% ［b,a］=butter(N,OmegaC,'s');



% 切比雪夫Ⅰ型

% ［N,OmegaC］=cheb1ord(OmegaP,OmegaS,Rp,As,'s');

% ［b,a］=cheby1(N,Rp,OmegaC,'s');



% 切比雪夫Ⅱ型

［N,OmegaC］=cheb2ord(OmegaP,OmegaS,Rp,As,'s');

［b,a］=cheby2(N,As,OmegaC,'s');



%冲激响应不变法 H(s)->H(z)

［bz,az］=impinvar(b,a,Fs);

［H,w］=freqz(bz,az);



% 检验衰减指标

Hx=freqz(bz,az,w0);

dbHx=-20*log10(abs(Hx)/max(abs(H)))



% 数字滤波器指标

［db,mag,pha,grd,w］=freqz_m(bz,az);

plot(Fs*w/pi/2,db);

hold on;

%% 模拟滤波器设计(双线性变换法)

% 频率预畸

OmegaP1=2*Fs*tan(wp/2);

OmegaS1=2*Fs*tan(ws/2);



% 巴特沃斯型

% ［N,OmegaC］=buttord(OmegaP1,OmegaS1,Rp,As,'s');

% ［b,a］=butter(N,OmegaC,'s');



% 切比雪夫Ⅰ型

% ［N,OmegaC］=cheb1ord(OmegaP1,OmegaS1,Rp,As,'s');

% ［b,a］=cheby1(N,Rp,OmegaC,'s');







% 切比雪夫Ⅱ型

［N,OmegaC］=cheb2ord(OmegaP1,OmegaS1,Rp,As,'s');

［b,a］=cheby2(N,As,OmegaC,'s');



% 双线性变换法 H(s)->H(z)

［bz,az］=bilinear(b,a,Fs);

［H,w］=freqz(bz,az);



% 检验衰减指标

Hx=freqz(bz,az,w0);

dbHx=-20*log10(abs(Hx)/max(abs(H)))

% 数字滤波器指标

［db,mag,pha,grd,w］=freqz_m(bz,az);

plot(Fs*w/pi/2,db);



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

grid on;

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






图3.18两种方法设计数字滤波器的比较



从图3.18可以看出，如果模拟原型滤波器采用巴特沃斯型（N=12阶），那么两种变换方法的衰减特性都可以满足设计要求。如果模拟原型滤波器采用切比雪夫Ⅱ型（N=6阶），则双线性变换法得到的衰减特性满足设计要求，而冲激响应不变法未能达到对150Hz频率分量衰减40dB的要求。

这是由这两种方法各自的特点决定的： 冲激响应不变法得到的频率响应是模拟滤波器频率响应的周期延拓，当其频率响应不是带限时就会造成混叠失真； 双线性变换法的s平面与z平面是一一对应的，消除了多值变换性，也就消除了频谱的混叠失真。但模拟频率与数字频率的非线性特性使得频率标度弯曲，不能保持原模拟滤波器的相频特性。因此，在实际使用中，应该根据使用场合的要求选择合适的IIR滤波器设计方法以满足指标要求。

3.3完全滤波器的设计
MATLAB信号处理工具箱提供了几个直接设计IIR数字滤波器的函数，它们把设计典型滤波器的几个步骤集成了一个函数，直接调用就可以设计出滤波器，这就叫作完全滤波器设计,这为我们设计滤波器提供了非常大的方便。这些函数应与数字滤波器阶数选择配合使用，为特定的滤波器的设计返回所需的阶数和固有频率。对于一般条件下使用的滤波器，运用完全设计的方法就可以满足滤波器性能的要求。下面以Butterworth滤波器设计为例介绍完全滤波器的设计。
MATLAB信号处理工具箱提供了函数butter设计Butterworth滤波器，可以设计低通、高通、带通和带阻的数字和模拟滤波器，其特性为使通带内的幅度响应最大限度的平坦，这会损失截止频率处的下降斜度。在期望通带平滑的情况下，可以使用Butterworth滤波器，但在下降斜度大的场合，应使用椭圆和Chebyshev Ⅰ型滤波器。
［b,a］=butter(n,Rp,Wn)可以设计出n阶低通数字Butterworth滤波器，其中截止频率由Wn确定，通带内波纹由Rp确定，滤波器的系统函数为

H(z)=B(z)A(z)=b(1)+b(2)z－1+…+b(n+1)z－n1+a(2)z－1+…+a(n+1)z－n

截止频率是滤波器幅度下降到Rp分贝处的频率，且是归一化之后的范围［0,1］。通带波纹Rp越小，可得到更宽的变换宽度。
［b,a］=butter(n,Wn,Rp,'ftype')可设计高通、低通、带通和带阻滤波器，ftype分别取值为high、low、band和stop即可。
［b,a］=butter(n,Wn,Rp,'s')可以设计出n阶低通、高通、带通和带阻模拟Butterworth滤波器，其中截止频率由Wn确定，通带内波纹由Rp确定，滤波器的系统函数为

H(s)=B(s)A(s)=b(1)sn+b(2)sn－1+…+b(n+1)sn+a(2)sn－1+…+a(n+1)

用butter函数设计模拟滤波器时，其用法与前面讲过的buttap基本相同，注意此时Wn必须以弧度为单位。
除了传递函数形式，butter函数输出还可以有零极点增益形式和状态方程形式。下面通过例题对butter函数进行演示。
例3.11运用完全设计方法，设计一个低通Butterworth数字滤波器，设计指标如下： Wp=30Hz，Ws=35Hz，Rp=0.5dB，Rs=40dB，Fs=100Hz。
解： MATLAB源代码如下(见chp3sec3_1.m)，实验结果如图3.19所示。



%Buttrerworth滤波器设计 chp3sec3_1.m

clear;

clc;

wp=30;ws=35;Fs=100;rp=0.5;rs=40;Fs=100;

［n,Wn］=buttord(wp/(Fs/2),ws/(Fs/2),rp,rs,'z');

［num,den］=butter(n,Wn);

［H,W］=freqz(num,den);

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

grid on;

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






图3.19Butterworth数字低通滤波器幅频响应


例3.12设数字抽样频率为100Hz，运用完全设计方法设计一个带通滤波器，性能指标为： 通带范围为100~250Hz，阻带上限为300Hz，下限为50Hz，通带内波纹小于3dB，阻带为-30dB，要求利用最小的阶实现。
解： MATLAB源代码如下(见chp3sec3_2.m)，实验结果如图3.20所示。



%Buttrerworth带通滤波器设计 chp3sec3_2.m

clear;

clc;

Wp1=100;Wp2=250;Ws1=50;Ws2=300;

Rp=3;Rs=30;Fs=1000;

Wp=［Wp1 Wp2］;

Ws=［Ws1 Ws2］;

［n,Wn］=buttord(Wp/(Fs/2),Ws/(Fs/2),Rp,Rs);

［b,a］=butter(n,Wn);

freqz(b,a,512,1000);






图3.20Butterworth数字带通滤波器幅频相频响应


例3.13设数字抽样频率为1000Hz，运用完全设计方法设计一个带阻滤波器，性能指标为： 阻带范围为50~300Hz，阻带上限大于250Hz，下限小于100Hz，阻带内波纹小于3dB，阻带为-30dB，要求利用最小的阶实现。
解： MATLAB源代码如下(见chp3sec3_3.m)，实验结果如图3.21所示。



%Buttrerworth带阻滤波器设计 chp3sec3_3.m

clear;

clc;

Wp1=100;Wp2=250;Ws1=50;Ws2=300;

Rp=3;Rs=30;Fs=1000;

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

［n,Wn］=buttord(Wp/(Fs/2),Ws/(Fs/2),Rp,Rs);

［b,a］=butter(n,Wn,'stop');

freqz(b,a,512,1000);






图3.21Butterworth数字带阻滤波器幅频相频响应


MATLAB信号处理工具箱还提供了函数cheby1来设计Chebyshev Ⅰ型滤波器，函数用法与Butterworth滤波器类似，在此不再介绍。
3.4FIR滤波器的基本结构
FIR滤波器突出的特点是单位脉冲响应h(n)仅有有限个非零值，即h(n)为一个N点长序列，0≤n≤N－1，其系统函数为

H(z)=∑N－1n=0h(n)z－n


FIR滤波器的结构主要是非递归结构，没有输出到输入的反馈，但有些结构，例如频率采样结构也包含有反馈的递归部分。FIR滤波器常见的基本结构包括直接型、级联型、线性相位型等。
1. 直接型
对于因果FIR系统，其系统函数只有有限个零点，这相当于IIR滤波器系统函数H(z)中所有系数ak都为零，因此差分方程为

y(n)=∑Mk=0bkx(n－k)

上式实际上是输入序列x(n)与单位脉冲响应h(n)的线性卷积，其信号流图如图3.22所示。


图3.22直接型FIR数字滤波器


由图3.22可见，这种结构与IIR数字滤波器在所有系数ak都为零时的结构是相同的，故FIR直接型数字滤波器是IIR直接型数字滤波器的一种特殊情况。
2. 级联型
将H(z)分解成实系数二阶因子的乘积形式

H(z)=∑N－1n=0h(n)z－n=∏Nck=1(b0k+b1kz－1+b2kz－2)

式中，Nc表示是N/2的最大整数。其信号流图如图3.23所示。


图3.23级联型FIR数字滤波器


3. 线性相位型

FIR数字滤波器的h(n)为实数，且满足以下任一条件 

h(n)=h(N－1－n)偶对称

h(n)=－h(N－1－n)奇对称

其对称中心在n=(N－1)/2处，则滤波器就具有准确的线性相位。
如果将频率响应表示为如下形式

H(ejω)=H(ω)ejθ(ω)

式中,H(ω)表示幅度函数，θ(ω)表示相位函数。
当h(n)为偶对称时

H(ω)=∑N－1n=0h(n)cosN－12－nω

θ(ω)=－N－12ω

可以看出此时θ(ω)是严格的线性相位，滤波器有(N－1)/2个抽样延时。

当h(n)为奇对称时

H(ω)=∑N－1n=0h(n)sinN－12－nω

θ(ω)=－N－12ω+π2

可以看出此时θ(ω)是严格的线性相位，滤波器有(N－1)/2个抽样延时，同时产生一个π/2的相移，称为正交网络。
存在4种线性相位FIR滤波器，可设计不同类型的数字滤波器［2］，具体情况见表3.2。


表3.24种线性相位FIR滤波器



类型h(n)N高通低通带通带阻

1型偶对称奇数√√√√
2型偶对称偶数×√√×
3型奇对称奇数××√×
4型奇对称偶数√×√×


例3.14已知FIR滤波器的系统函数如下所示，请给出该滤波器直接Ⅱ型、级联型、频率采样型结构。

H(z)=1+16.0625z－4+z－8

解： 根据系统函数可知系统差分方程为

y(n)=x(n)+16.0625x(n－4)+x(n－8)

由差分方程可以直接画出滤波器直接Ⅱ型结构如图3.24所示。


图3.24直接Ⅱ型


与IIR滤波器类似，函数tf2sos实现FIR滤波器直接型到级联型的转换，MATLAB源代码如下（见chp3sec4_1.m）。



% chp3sec4_1.m FIR滤波器直接型转级联型

% 函数tf2sos: transfer function to second-order sections

clear;

clc;

close all;

b=［1,0,0,0,16.0625,0,0,0,1］;%传递函数分子

a=［1］; %传递函数分母，1表示FIR滤波器

［sos,g］=tf2sos(b,a)




程序运行结果如下，g表示级联型结构的增益，sos存储的是级联型结构的系数，具体含义与IIR滤波器类似。



sos =

1.00002.82844.00001.000000

1.0000-2.82844.00001.000000

1.00000.70710.25001.000000

1.0000-0.70710.25001.000000

g =

1




因此，级联形式的传递函数如下所示，其结构见图3.25。

H(z)=(1+2.8284z－1+4z－2)(1－2.8284z－1+4z－2)
(1+0.7071z－1+0.25z－2)(1－0.7071z－1+0.25z－2)



图3.25级联形式



自定义函数tf2fs实现直接型到频率采样型的转换，MATLAB源代码如下（见chp3sec4_2.m）。



% chp3sec4_2.m FIR滤波器直接型转频率采样型

clear;

clc;

close all;

h=［1,0,0,0,16.0625,0,0,0,1］; %FIR滤波器

［C,B,A］=tf2fs(h)




自定义函数tf2ts源代码如下（见tf2ts.m）。



function ［C,B,A］ = tf2fs(h)

% ［C,B,A］ = tf2fs(h)

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

% C = 包含各并行部分增益的行向量

% B = 包含按行排列的分子系数矩阵

% A = 包含按行排列的分母系数矩阵

% h =  FIR滤波器的脉冲响应向量，即FIR滤波器系数向量

%

M = length(h); % 滤波器长度

H = fft(h,M);               % 滤波器的频率特性

magH = abs(H); phaH = angle(H)';    % 滤波器的幅特性和相特性

% 检查M的奇偶性，

if (M == 2*floor(M/2))      % 若M为偶数

L = M/2-1;            % 则A1和C1形式如下

A1 = ［1,-1,0;1,1,0］;

C1 = ［real(H(1)),real(H(L+2))］;










else  % 若M为奇数

L = (M-1)/2;

A1 = ［1,-1,0］;         % 则A1和C1形式如下

C1 = ［real(H(1))］;

end

k = ［1:L］';

% 初始化 B 和 A 数组

B = zeros(L,2); A = ones(L,3);

% 计算分母系数

A(1:L,2) = -2*cos(2*pi*k/M); A = ［A;A1］;

% 计算分子系数

B(1:L,1) = cos(phaH(2:L+1));

B(1:L,2) = -cos(phaH(2:L+1)-(2*pi*k/M));

% 计算增益系数

C = ［2*magH(2:L+1),C1］';




程序运行结果如下，C表示各并行部分增益的行向量，B表示分子系数矩阵，A表示分母系数矩阵。



C =

28.3662

35.1892

30.1250

32.8196

18.0625

B =

-0.93970.9397

0.7660-0.7660

-0.50000.5000

0.1736-0.1736

A =

1.0000-1.53211.0000

1.0000-0.34731.0000

1.00001.00001.0000

1.00001.87941.0000

1.0000-1.00000




因此，频率采样形式的传递函数如下所示，其结构见图3.26。

H(z)=1－z－9928.3662－0.9397+0.9397z－11－1.5321z－1+z－2+35.18920.7660－0.7660z－11－0.3473z－1+z－2
+30.1250－0.5000+0.5000z－11+z－1+z－2+32.81960.1736－0.1736z－11+1.8794z－1+z－2
+18.062511－z－1




图3.26频率采样形式



3.5FIR滤波器的设计
本章主要介绍FIR滤波器的两种设计方法： 窗函数设计法和频率采样设计法。
3.5.1窗函数设计法

设所要求的理想数字滤波器的频率响应为Hd(ejω)，hd(n)是与其对应的单位脉冲响应，因此

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


由于Hd(ejω)是矩形频率特性，故hd(n)一定是无限长非因果序列。而所设计的是FIR滤波器，其单位脉冲响应h(n)必然是有限长的，所以要用有限长的h(n)来逼近无限长的hd(n)，最有效的方法就是截断hd(n)，即用有限长的w(n)来截取，表示为

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

这种方法就称为窗函数设计法。常用的窗函数包括矩形窗、三角窗、汉宁窗(Hanning)、海明窗(Hamming)等［2］。
1. 窗函数法设计数字低通滤波器
理想低通数字滤波器的频率响应Hd(ejω)为

Hd(ejω)=e－jωa，|ω|≤ωc


0，ωc<|ω|≤π

其中，ωc表示截止频率，单位为rad； a表示采样延迟。则理想低通数字滤波器的单位脉冲响应hd(n)为

hd(n)=12π∫π－πHd(ejω)ejωndω=12π∫ωc－ωcejωaejωndω=sin［ωc(n－a)］π(n－a)

其中，hd(n)为无限长非因果序列，关于a偶对称。
为了从hd(n)得到一个FIR数字滤波器，必须同时在两边截取hd(n)，要得到一个因果的线性相位FIR滤波器，h(n)的长度为N，必须有

h(n)=hd(n)，0≤n≤N－1

0，其他


这种截取可以看作是h(n)=hd(n)w(n)，其中w(n)是宽度为N的矩形窗。
h(n)是关于a偶对称的有限因果序列。N为奇数时是1型，N为偶数时是2型。用自定义函数ideal_lp计算理想低通滤波器的单位脉冲响应hd(n)，源代码如下(见ideal_lp.m)： 



function hd=ideal_lp(wc,N)

%计算理想低通滤波器的单位脉冲响应

alpha=(N-1)/2;

n=0:N-1;

m=n-alpha+eps;

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





在频域中，FIR数字滤波器的频率响应H(ejω)为

H(ejω)=12π∫π－πHd(ejθ)W(ej(ω－θ))dθ

因而H(ejω)逼近Hd(ejθ)的好坏，完全取决于窗函数的频率特性W(ejω)。
例3.15根据下面的技术指标，设计一个数字FIR低通滤波器，选择合适的窗函数，确定单位脉冲响应，并绘出所设计的滤波器的幅度响应。

ωp=0.2π,Ap=0.25dB

ωr=0.4π,Ar=50dB


解： 根据窗函数最小阻带衰减的特性，只有海明窗和布莱克曼窗可以提供大于50dB的衰减,故在此选择海明窗，它提供较小的过渡带，因此具有较小的阶数。MATLAB源代码如下(见chp3sec5_1.m)，实验结果如图3.27所示。



%数字低通滤波器的窗函数设计 chp3sec5_1.m

clear;

clc;

wp=0.2*pi;wr=0.4*pi;

tr_width=wr-wp;%过渡带宽

N=ceil(6.6*pi/tr_width)+1;%滤波器的长度，N=奇数为1型，N=偶数为2型

n=0:1:N-1;

wc=(wr+wp)/2;%理想低通滤波器的截止频率

hd=ideal_lp(wc,N);%理想低通滤波器的单位脉冲响应








w_ham=(hamming(N))';%海明窗

h=hd.*w_ham;%截取得到实际单位脉冲响应

［db,mag,pha,grd,w］=freqz_m(h,［1］);%计算实际滤波器的幅度响应

delta_w=2*pi/1000;

Ap=-(min(db(1:1:wp/delta_w+1)));%实际通带波动

Ar=-round(max(db(wr/delta_w+1:1:501)));%最小阻带衰减

subplot(2,2,1);stem(n,hd);title('理想单位脉冲响应h_d(n)');

subplot(2,2,2);stem(n,w_ham);title('海明窗w(n)');

subplot(2,2,3);stem(n,h);title('单位脉冲响应h(n)');

subplot(2,2,4);plot(w/pi,db);title('幅度响应/dB');axis(［0,1,-100,10］);






图3.271型低通数字滤波器幅度响应


从图3.27的结果可以看出： 滤波器的长度为35，实际通带波动为0.0301dB，最小阻带衰减为55dB，满足设计要求。
2. 窗函数法设计数字高通滤波器

理想高通数字滤波器的频率响应Hd(ejω)为

Hd(ejω)=e－jωa，ωc≤|ω|≤π

0，|ω|<ωc

式中，ωc表示截止频率，单位为rad； a表示采样延迟。则理想高通数字滤波器的单位脉冲响应hd(n)为

hd(n)=12π∫π－πHd(ejω)ejωndω=sin［π(n－a)－ωc(n－a)］π(n－a)

式中，hd(n)为无限长非因果序列，关于a偶对称。
高通数字滤波器的单位脉冲响应为h(n)=hd(n)w(n)，h(n)是关于a偶对称的因果序列，N应为奇数(1型)。用自定义函数ideal_hp计算理想高通滤波器的单位脉冲响应hd(n)，源代码如下(见ideal_hp.m)。



function hd=ideal_hp(wc,N)

%计算理想高通滤波器的单位脉冲响应

alpha=(N-1)/2;

n=0:N-1;

m=n-alpha+eps;

hd=［sin(pi*m)-sin(wc*m)］./(pi*m);





例3.16根据下面的技术指标，设计一个数字FIR高通滤波器，选择合适的窗函数，确定单位脉冲响应，并绘出所设计的滤波器的幅度响应。

ωp=0.6π,Ap=0.25dB

ωr=0.4π,Ar=40dB


解： 根据窗函数最小阻带衰减的特性，只有汉宁窗可以达到44dB的最小阻带衰减，它提供的过渡带宽为6.2π/N，因此具有较小的阶数。MATLAB源代码如下(见chp3sec5_2.m)，实验结果如图3.28所示。



%理想高通数字滤波器  chp3sec5_2.m

clear;

clc;

wp=0.6*pi;wr=0.4*pi;

tr_width=wp-wr;

N=ceil(6.2*pi/tr_width);%应使截取的长度为奇数

n=0:1:N-1;

wc=(wr+wp)/2;%理想高通滤波器的截止频率

hd=ideal_hp(wc,N);%1型理想高通滤波器响应

w_han=(hanning(N))';

h=hd.*w_han;

［db,mag,pha,grd,w］=freqz_m(h,［1］);

delta_w=2*pi/1000;

Ap=-(min(db(wp/delta_w+1:1:501)));%实际通带波动

Ar=-round(max(db(1:1:wr/delta_w+1)));%最小阻带衰减

subplot(221);stem(n,hd);title('理想单位脉冲响应h_d(n)');

subplot(222);stem(n,w_han);title('汉宁窗w(n)');

subplot(223);stem(n,h);title('实际单位脉冲响应h(n)');

subplot(224);plot(w/pi,db);title('幅度响应/dB');






图3.281型高通数字滤波器的幅度响应


结果如图3.28所示，滤波器长度为31，实际通带波动为0.0887dB，最小阻带衰减为44dB，满足设计要求。
3. 窗函数法设计数字带通滤波器
理想带通数字滤波器的频率响应Hd(ejω)为

Hd(ejω)=e－jωa,ωc1≤|ω|≤ωc2

0，其他

式中,ωc1表示下截止频率，ωc2表示上截止频率，单位为rad； a表示采样延迟,则理想带通数字滤波器的单位脉冲响应hd(n)为

hd(n)=12π∫π－πHd(ejω)ejωndω=sin［ωc2(n－a)－ωc1(n－a)］π(n－a)


式中,hd(n)为无限长非因果序列，关于a偶对称。
带通数字滤波器的单位脉冲响应为h(n)=hd(n)w(n)，h(n)是关于a偶对称的因果序列，当N为奇数时是1型，当N为偶数时是2型。用自定义函数ideal_bp计算理想带通滤波器的单位脉冲响应hd(n)，源代码如下(见ideal_bp.m)。



function hd=ideal_bp(wc1,wc2,N)

%计算理想带通滤波器的单位脉冲响应

alpha=(N-1)/2;

n=0:N-1;

m=n-alpha+eps;

hd=［sin(wc2*m)-sin(wc1*m)］./(pi*m);




例3.17根据下面的技术指标，设计一个数字FIR带通滤波器，选择合适的窗函数，确定单位脉冲响应，并绘出所设计的滤波器的幅度响应。

低端阻带边缘ωr1=0.2π,Ar1=60dB

低端通带边缘ωp1=0.4π,Ap1=1dB

高端通带边缘ωp2=0.6π,Ap2=1dB

高端阻带边缘ωr2=0.8π,Ar2=60dB


解： 根据窗函数最小阻带衰减的特性，选择布莱克曼窗可以达到75dB的最小阻带衰减，它提供的过渡带宽为11π/N。MATLAB源代码如下(见chp3sec5_3.m)，实验结果如图3.29所示。



%理想带通数字滤波器  chp3sec5_3.m

clear;

clc;

wr1=0.2*pi;wp1=0.4*pi;

wp2=0.6*pi;wr2=0.8*pi;

tr_width=min((wp1-wr1),(wr2-wp2));

N=ceil(11*pi/tr_width)+1;%长度为奇数1型，长度为偶数2型

n=0:1:N-1;

wc1=(wr1+wp1)/2;wc2=(wr2+wp2)/2;%理想带通滤波器的上、下截止频率

hd=ideal_bp(wc1,wc2,N);%1、2型理想带通滤波器响应

w_blk=(blackman(N))';








h=hd.*w_blk;

［db,mag,pha,grd,w］=freqz_m(h,［1］);

delta_w=2*pi/1000;

Ap=-(min(db(wp1/delta_w+1:1:wp2/delta_w+1)));%实际通带波动

Ar=-round(max(db(wr2/delta_w+1:1:501)));%最小阻带衰减

subplot(2,2,1);stem(n,hd);title('理想单位脉冲响应h_d(n)');

subplot(2,2,2);stem(n,w_blk);title('布莱克曼窗w(n)');

subplot(2,2,3);stem(n,h);title('实际单位脉冲响应h(n)');

subplot(2,2,4);plot(w/pi,db);title('幅度响应/dB');






图3.29带通数字滤波器的幅度响应


结果如图3.29所示，滤波器长度为56，实际通带波动为0.0027dB，最小阻带衰减为73dB，满足设计要求。
4. 窗函数法设计数字带阻滤波器

理想带阻数字滤波器的频率响应Hd(ejω)为

Hd(ejω)=e－jωa，|ω|≤ωc1,π≥|ω|≥ωc2

0，ωc1<|ω|<ωc2

式中，ωc1表示下截止频率，ωc2表示上截止频率，单位为rad； a表示采样延迟，则理想带阻数字滤波器的单位脉冲响应hd(n)为

hd(n)=12π∫π－πHd(ejω)ejωndω=sin［ωc1(n－a)+π(n－a)－ωc2(n－a)］π(n－a)

式中，hd(n)为无限长非因果序列，关于a偶对称。
带阻数字滤波器的单位脉冲响应为h(n)=hd(n)w(n)，h(n)为关于a偶对称的因果序列，N应为奇数(1型)。用自定义函数ideal_bs计算理想带阻滤波器的单位脉冲响应hd(n)，源代码如下(见ideal_bs.m)。



function hd=ideal_bs(wc1,wc2,N)

%计算理想带阻滤波器的单位脉冲响应

alpha=(N-1)/2;

n=0:N-1;

m=n-alpha+eps;

hd=［sin(wc1*m)+sin(pi*m)-sin(wc2*m)］./(pi*m);




例3.18根据下面的技术指标，设计一个数字FIR带阻滤波器，选择合适的窗函数，确定单位脉冲响应，并绘出所设计的滤波器的幅度响应。

低端阻带边缘ωr1=0.4π,Ar1=40dB

低端通带边缘ωp1=0.3π,Ap1=1dB

高端通带边缘ωp2=0.8π,Ap2=1dB

高端阻带边缘ωr2=0.6π,Ar2=40dB


解： 根据窗函数最小阻带衰减的特性，选择汉宁窗可以达到44dB的最小阻带衰减，它提供的过渡带宽为6.2π/N。MATLAB源代码如下(见chp3sec5_4.m)，实验结果如图3.30所示。



%带阻数字滤波器  chp3sec5_4.m

clear;

clc;

wr1=0.4*pi;wp1=0.2*pi;

wp2=0.8*pi;wr2=0.6*pi;

tr_width=min((wr1-wp1),(wp2-wr2));

N=ceil(6.2*pi/tr_width);

n=0:1:N-1;

wc1=(wr1+wp1)/2;wc2=(wr2+wp2)/2; %理想带阻滤波器的上、下截止频率

hd=ideal_bs(wc1,wc2,N);%理想带阻滤波器响应

w_han=(hanning(N))';

h=hd.*w_han;

［db,mag,pha,grd,w］=freqz_m(h,［1］);

delta_w=2*pi/1000;

Ap=-(min(db(1:1:wp1/delta_w+1)));%实际通带波动

Ar=-round(max(db(wr1/delta_w+1:1:wr2/delta_w+1)));%最小阻带衰减

subplot(221);stem(n,hd);title('理想单位脉冲响应h_d(n)');

subplot(222);stem(n,w_han);title('汉宁窗w(n)');

subplot(223);stem(n,h);title('实际单位脉冲响应h(n)');

subplot(224);plot(w/pi,db);title('幅度响应/dB');








图3.30带阻数字滤波器的幅度响应


结果如图3.30所示，滤波器长度为31，实际通带波动为0.0885dB，最小阻带衰减为44dB，满足设计要求。
在前面介绍的四种自定义函数源代码中(如ideal_lp，ideal_hp等)，都出现了eps，目的就是使0/0有输出。在MATLAB中，eps表示数与数之间的最小分辨率，或者可以说是MATLAB中绝对值最小的数，取值一般为2.2204e16。
以例3.16为例，在子函数ideal_hp中有一段代码如下： 



m=n-alpha+eps;

hd=［sin(pi*m)-sin(wc*m)］./(pi*m)




如果没有eps，那么数组m的第16个取值就是0，在计算数组hd的第16个元素时就会出现0/0的情况，此时MATLAB的输出为NaN。如果加上eps，数组m的第16个取值就是eps，再计算hd的第16个元素，取值就为1/2。有兴趣的读者，可以尝试在MATLAB命令行中输入0/0和eps/eps，前者输出为NaN，后者输出为1。
3.5.2频率采样设计法

频率采样设计法是从频域出发，对给定的理想频率响应Hd(ejω)等间隔采样N个点得到

Hd(ejω)ω=2πNk=Hd(k)=H(k)

再对Hd(k)作离散傅里叶逆变换(IDFT)，得到h(n)为

h(n)=1N∑N－1k=0H(k)ej2πNkn,k=0,1,…,N－1

将h(n)作为所设计的滤波器的单位脉冲响应。
频率采样一般有两种方法： 
1. 频率采样法1
第一个采样点在ω=0处




H(k)=Hd(k)=Hd(ejω)ω=2πNk

h(n)=1N∑N－1k=0H(k)ej2πNkn

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

H(ejω)=1Ne－jN－12∑N－1k=0H(k)e－jkπNsinωN2sinω2－πkN


由于

H(k)=∑N－1n=0h(n)e－j2πNkn


当h(n)为实数时，满足H(k)=H*((N－k))N，故

|H(k)|=|H(N－k)|

θ(k)=－θ(N－k)


当N为奇数时，有线性相位约束条件 


θ(k)=－2πNkN－12k=0,1,…,N－12

2πN(N－k)N－12k=N+12,…,N－11型

θ(k)=π2－2πNkN－12k=0,1,…,N－12
－π2+2πN(N－k)N－12k=N+12,…,N－13型


当N为偶数时，有线性相位约束条件 

θ(k)=－2πNkN－12k=0,1,…,N2－1

2πN(N－k)N－12k=N2+1,…,N－1
0k=N22型

θ(k)=π2－2πNkN－12k=0,1,…,N2－1
－π2+2πN(N－k)N－12k=N2+1,…,N－1
0k=N24型


例3.19利用频率采样法1，根据下面的技术指标，设计2型FIR低通滤波器，确定单位脉冲响应，并绘出所设计的滤波器的幅度响应。

ωp=0.3πAp=1dB

ωr=0.4πAr=40dB


解： 选择N=60，则ωp在k=9处，ωr在k=12处。T1=0.7，T2=0.2。MATLAB源代码如下(见chp3sec5_5.m)，实验结果如图3.31所示。



% 频率采样法设计低通滤波器(2型)(N=60) chp3sec5_5.m

clear;

clc;

close all;



N=60; %滤波器阶数(该数值计算由过渡带宽等决定)



% wd表示频率采样点，MATLAB规定一定是从0到1

% 0表示起始频率，1表示归一化的折叠频率，也就是pi

% 0.35表示低通的cutoff频率，对应两个值，表示在0.35处从0变为了1







wd=［0,0.35,0.35,1］; %频率采样点

mhi=［1,1,0,0］; % 频率采样点处对应的幅值(理想值)

wd2=［0,0.35,0.35,1.65,1.65,2］; %为了显示0-2*pi的频率采样点，

mhi2=［1,1,0,0,1,1］; % 频率采样点处对应的幅值(理想值)



m=0:N-1;

w1=(2*pi/N)*m; % N阶滤波器的频域采样点index



% ///////////////////////////加过渡带抽样///////////////////////////

% 理想振幅响应采样加了两个点的过渡带抽样

T1=0.7;T2=0.2; %查表得到

Hrs=［ones(1,10),T1,T2,zeros(1,37),T2,T1,ones(1,9)］;%理想振幅响应采样

h=fir2(N-1,w1(1:1+N/2)/pi,Hrs(1:31)); %取wd1的一半，而且保证一定是0开头，1结尾，
%对π归一化

% //////////////////////////////////////////////////////////////////



［db,mag,pha,grd,w］=freqz_m(h,［1］);

［Hr,ww,a,L］=hr_type2(h);%实际振幅响应(2型振幅响应是奇对称的,这就体现了相位信息)



subplot(221);plot(w1/pi,Hrs,'r.',wd2,mhi2);legend('抽样值','理想值');

title('频率样本');axis(［0 2 -2 2］);grid on;

subplot(222);stem(m,h,'b.');

title('脉冲响应h(n)'); grid on;

subplot(223);plot(w1/pi,Hrs,'r.',ww/pi,Hr);

title('振幅响应');legend('抽样值','实际值');grid on;

subplot(224);plot(w/pi,db);

title('幅频响应/dB');axis(［0 1 -150 10］); grid on;







图3.31加过渡点的频率采样设计法设计低通滤波器(2型，N=60)


从幅度响应曲线可以看出，N=60时最小的阻带衰减为45dB，达到设计指标要求。此处使用自定义函数hr_type2计算2型滤波器的振幅响应。为便于查找，在此给出计算1型、2型、3型和4型滤波器振幅响应的源代码。



%计算1型滤波器设计的振幅响应

function［Hr,w,a,L］=hr_type1(h)

%Hr=振幅响应

%a=1型滤波器的系数

%L=Hr的阶次

%h=1型滤波器的单位脉冲响应

M=length(h);L=(M-1)/2;

a=［h(L+1) 2*h(L:-1:1)］;

n=［0:1:L］;

w=［0:1:500］'*2*pi/500;

Hr=cos(w*n)*a';



%计算2型滤波器设计的振幅响应

function［Hr,w,b,L］=hr_type2(h)

%Hr=振幅响应

%b=2型滤波器的系数

%L=Hr的阶次

%h=2型滤波器的单位脉冲响应

M=length(h);L=M/2;

b=2*［h(L:-1:1)］;

n=［1:1:L］;n=n-0.5;

w=［0:1:500］'*2*pi/500;

Hr=cos(w*n)*b';



%计算3型滤波器设计的振幅响应

function［Hr,w,c,L］=hr_type3(h)

%Hr=振幅响应

%c=3型滤波器的系数

%L=Hr的阶次

%h=3型滤波器的单位脉冲响应

M=length(h);L=(M-1)/2;

c=［2*h(L+1:-1:1)］;

n=［0:1:L］;

w=［0:1:500］'*2*pi/500;

Hr=sin(w*n)*c';



%计算4型滤波器设计的振幅响应

function［Hr,w,d,L］=hr_type4(h)

%Hr=振幅响应

%c=4型滤波器的系数

%L=Hr的阶次

%h=4型滤波器的单位脉冲响应

M=length(h);L=M/2;

d=2*［h(L:-1:1)］;

n=［1:1:L］;n=n-0.5;

w=［0:1:500］'*2*pi/500;

Hr=sin(w*n)*d';





2. 频率采样法2
第一个采样点在ω=π/N处

H(k)=Hd(k)=Hd(ejω)ω=2πNk+πN

h(n)=1N∑N－1k=0H(k)ej2πNkn+πN=ejπNIDFT［H(k)］

H(z)=1+z－NN∑N－1k=0H(k)1－W－k+12Nz－1

H(ejω)=cosωN2Ne－jN－12ω∑N－1k=0H(k)e－jπNk+121jsinω2－πkNk+12


由于

H(k)=∑N－1n=0h(n)e－j2πNk+12n



当h(n)为实数时，满足H(k)=H*((N－1－k))N，故

|H(k)|=|H(N－1－k)|

θ(k)=－θ(N－1－k)


当N为奇数时，有线性相位约束条件

θ(k)=－2πNk+12N－12，k=0,1,…,N－32
0，k=N－12
2πNN－k－12N－12，k=N+12,…,N－11型

θ(k)=π2－2πNk+12N－12，k=0,1,…,N－32
0，k=N－12
－π2+2πNN－k－12N－12，k=N+12,…,N－13型


当N为偶数时，有线性相位约束条件

θ(k)=－2πNkN－12，k=0,1,…,N2－1

2πN(N－k)N－12，k=N2+1,…,N－1
0，k=N22型

θ(k)=π2－2πNkN－12，k=0,1,…,N2－1
－π2+2πN(N－k)N－12，k=N2+1,…,N－1
0，k=N24型


例3.20利用频率采样法2，根据下面的技术指标，设计1型FIR低通滤波器，确定单位脉冲响应，并绘出所设计的滤波器的幅度响应。

ωp=0.25π，Ap=2dB

ωr=0.35π，Ar=20dB


解： 选择N=41，则ωp在k=5附近，ωr在k=7附近。MATLAB源代码如下(见chp3sec5_6.m)，实验结果如图3.32所示。



% 频率采样法设计低通滤波器(1型)(N=41) chp3sec5_6.m

% 不加过渡带抽样

clear;

clc;

close all;



N=41; %滤波器阶数(该数值计算由过渡带宽等决定)



% wd表示频率采样点，MATLAB规定一定是从0到1

% 0表示起始频率，1表示归一化的折叠频率，也就是pi

% 0.35表示低通的cutoff频率，对应两个值，表示在0.35处从0变为了1

wd=［0,0.3,0.3,1］; %频率采样点

mhi=［1,1,0,0］; % 频率采样点处对应的幅值(理想值)

wd2=［0,0.3,0.3,1.7,1.7,2］; %为了显示0-2*pi的频率采样点，

mhi2=［1,1,0,0,1,1］; % 频率采样点处对应的幅值(理想值)



m=0:N-1;

w1=(2*pi/N)*m; % N阶滤波器的频域采样点index



% ///////////////////////////不加过渡带抽样/////////////////////////////////

% 直接代入wd和mhi

h=fir2(N-1,wd,mhi);

T1=1;T2=0; %1和0相当于在过渡带没有加抽样点

Hrs=［ones(1,6),T1,T2,zeros(1,26),T2,T1,ones(1,5)］;%理想振幅响应采样

% ////////////////////////////////////////////////////////////////////////



［db,mag,pha,grd,w］=freqz_m(h,［1］);

［Hr,ww,a,L］=hr_type1(h);%实际振幅响应(1型振幅响应是偶对称的,这就体现了相位信息)



subplot(221);plot(w1/pi,Hrs,'r.',wd2,mhi2);legend('抽样值','理想值');

title('频率样本');axis(［0 2 -2 2］);grid on;

subplot(222);stem(m,h,'b.');

title('脉冲响应h(n)'); grid on;

subplot(223);plot(w1/pi,Hrs,'r.',ww/pi,Hr);

title('振幅响应');legend('抽样值','实际值');axis(［0 2 -2 2］);grid on;

subplot(224);plot(w/pi,db);

title('幅频响应/dB');axis(［0 1 -100 10］); grid on;







图3.32加过渡点的频率采样设计法设计低通滤波器(1型，N=41)


从幅度响应曲线可以看出，N=41时最小的阻带衰减为24dB，达到设计指标要求。
例3.21分别用频率采样法和窗函数法设计一个阶数为21的FIR低通滤波器，要求的技术指标如下： 


ωp=0.3π，Ap=5dB

ωr=0.4π，Ar=40dB

ωs=2.1π


解： 对于频率采样法

kp=N·ωpωs=3

kp=N·ωrωs=4

则

|H(k)|={1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1}

θ(k)=－20π21k，k=0,1,…,10
20π21(21－k)，k=11,12,…,20

H(k)=|H(k)|ejθ(k)

从而得到FIR滤波器的系数为

h(n)=IDFT［H(k)］


对于窗函数法，根据窗函数最小阻带衰减特性，汉宁窗可以提供44dB的衰减，满足设计要求。MATLAB源代码如下(见chp3sec5_7.m)，实验结果如图3.33所示。





%分别利用频率采样法和窗函数法设计FIR低通滤波器 chp3sec5_7.m

clear;

clc;

N=21;

%频率采样法

Hrs=［ones(1,4),zeros(1,13),ones(1,4)］;

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

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

alpha=(N-1)/2;

angH=［-alpha*2*pi*k1/N,alpha*2*pi*(N-k2)/N］;

H=Hrs.*exp(j*angH);

h1=real(ifft(H));

freqz(h1);

hold on;

%窗函数法

wp=0.3*pi; wr=0.4*pi;

tr_width=wr-wp;

n=0:N-1;

wc=(wp+wr)/2;

m=n-alpha+eps;

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

w_hann=hanning(N);

h2=hd.*w_hann';

freqz(h2);







图3.33频率采样法和窗函数法实现比较


由图3.33可以看出，由频率采样法得到的FIR滤波器的最小阻带衰减约为-30dB，没有达到设计要求，可以采取增大采样点数的方法提高衰减性能。窗函数法得到的同样阶数的FIR滤波器的最小阻带衰减约为-45dB，达到了设计要求。
由此可以得出结论： 在同样滤波阶数的条件下，窗函数法可以比频率采样法得到更低的阻带衰减。这是因为频率采样法只保证设计的频率响应和要求的频率响应在有限个点处一致，而在其他频率处则没有约束，这就会造成比较大的波动； 通过选择合适的窗，可以有效地抑制吉布斯效应。因此，在同等条件下，窗函数法一般会比频率采样法得到更佳的滤波性能。
FIR滤波器最大的优点就是可以在通带内得到线性相位，窗函数的选择对于结果的优劣起着重要的作用。加窗的主要作用就是消除由于对无限长序列的截断导致的吉布斯效应。窗函数要针对不同的信号和不同的处理目的来选择，才能收到良好的效果。利用频率采样法比较方便，可以更好地满足频域的要求。
3.6IIR和FIR滤波器的比较
至今为止，我们讨论了IIR和FIR两种滤波器的设计方法，但在实际应用时应该如何选择它们呢？下面对这两种滤波器进行比较。
从性能上说，IIR滤波器系统函数的极点可以位于单位圆内的任何位置，因此可以用较低的阶数获得高的选择性，所用的存储单元少，所以经济效率高。但是这个高效率是以相位的非线性为代价的。FIR却可以得到严格的线性相位，然而FIR滤波器系统函数的极点固定在原点，所以只能用较高的阶数达到高的选择性。对于同样的滤波器设计指标，FIR滤波器所需的阶数比IIR滤波器的阶数高5~10倍，结果成本较高，信号延时也较大。
从结构上看，IIR滤波器必须采用递归结构，极点位置必须在单位圆内，否则系统将不稳定。在这种结构中，由于运算过程中对序列进行舍入处理，这种有限字长效应有时会产生寄生振荡。相反，FIR滤波器采用非递归结构，不存在稳定性的问题，运算误差也较小。此外FIR滤波器可以采用快速傅里叶变换算法，在相同阶数的条件下，运算速度快得多。
从设计工具看，IIR滤波器可以借助模拟滤波器的成果，因此，一般都有有效的封闭形式的设计公式可供准确计算，计算工作量比较小，对计算工具的要求不高。FIR滤波器一般没有封闭形式的设计公式。窗口法仅对窗口函数可以给出计算公式，但计算通带阻带衰减无显式表达式。一般FIR滤波器的设计只有计算程序可循，因此对计算工具要求较高。
另外也应看到，IIR滤波器虽然设计简单，但主要是用于设计具有片段常数特性的滤波器，如低通、高通、带通及带阻等，往往脱离不了模拟滤波器的格局。FIR滤波器则要灵活得多，尤其它能易于适应某些特殊的应用，如构成微分器或积分器，或用于巴特沃斯、切比雪夫等逼近不可能达到预定指标的情况。例如，由于某些原因要求三角形振幅响应或一些更复杂的幅频响应，因而有更大的适应性和更广阔的天地。
从上面的简单比较我们可以看到，IIR和FIR滤波器各有所长，所以在实际应用时应该从多方面考虑来加以选择。例如，从使用要求来看，在对相位要求不敏感的场合，如语音通信等，选用IIR较为合适，这样可以充分发挥其经济高效的特点，而对于图像信号处理、数据传输等以波形携带信息的系统，则对线性相位要求较高，采用FIR滤波器较好。当然，在实际应用中还应考虑经济上的要求。
3.7课外知识点： 滤波器设计与分析工具FDATool
设计一个“好”的数字滤波器，需要根据实际应用需求事先给定设计要求。例如选定的是FIR还是IIR类型的滤波器，滤波器的阶数，窗函数类型，以及阻/通带衰减，截止频率等技术指标。当给定这些技术指标以后，设计出来的数字滤波器就是一组系数(Filter Coefficients)。
本章前面已详细介绍利用MATLAB自带的函数设计数字滤波器，在此介绍MATLAB信号处理工具箱中的一种集成工具FDATool(Filter Design & Analysis Tool)，该工具可以用来设计和分析数字滤波器。在MATLAB命令窗口中输入fdatool命令后就会出现如图3.34所示窗口界面，该工具主要有以下6个功能区域。
 区域1 Current Filter Information，该区域简明扼要地给出滤波器的基本信息，如结构、阶数、是否稳定等。
 区域2 Filter Specifications，该区域主要描述滤波器的基本性能，如滤波器的理想特性、幅频特性、相频特性、零极点分布等。
 区域3 Response Type & Design Method，该区域主要给定滤波器的响应特性和设计方法，前者包括低通、高通还是带通滤波器等，后者主要包括是IIR类型还是FIR类型的滤波器，



图3.34FDATool工具界面


以及相应的滤波器类型(如Butterworth、Chebyshev Ⅰ和Chebyshev Ⅱ等)。
 区域4 Filter Order，该区域主要给定滤波器的阶数和相应的选项。
 区域5 Frequency Specifications，该区域主要给定滤波器的频率参数，如通带频率、截止频率、采样率、频率单位等。
 区域6 Magnitude Specifications，该区域主要给定滤波器的幅度参数，如通带波纹(Apass)和阻带衰减(Astop)等。
需要注意的是： 当其中某一个输入参数发生变化时，其他区域的输入参数也会相应地进行更新，每个区域中的输入并不是一成不变的。
FDATool工具很好地集成了数字滤波器设计中现有的算法和技术指标，而且软件的人机交互也非常不错，使用者只需要根据技术要求进行输入即可完成滤波器的设计。下面通过一个具体的例子进行讲解。
例3.22设计一个12阶的低通滤波器，采用汉宁窗(Hanning Window)函数，截止频率为7kHz，采样频率为48kHz。
解： 根据题目给出的要求在FDATool参数设置区域中进行选择和输入即可，如图3.35所示。可以发现，选定滤波器类型为Lowpass后，区域6的幅度特性无法再进行设置(截止频率的衰减为固定值6dB)。图3.36给出了该低通滤波器的理想幅频响应。


图3.35输入滤波器设计的技术指标




图3.36低通数字滤波器的理想特性


参数设置完毕后，单击Design Filter按钮即可自动完成滤波器设计工作。单击Analysis下拉菜单，可以对滤波器的各种特性进行分析和比较，如幅频特性(图3.37)、相频特性(图3.38)、单位阶跃响应(图3.39)和零极点分布(图3.40)等。


图3.37滤波器的幅频特性




图3.38滤波器的相频特性




图3.39滤波器的单位阶跃响应




图3.40滤波器的零极点分布



当根据给定的技术指标设计好数字滤波器后，FDATool的输出就是滤波器的系数，如图3.41所示。
单击File→Export下拉菜单，可以将设计好的滤波器系数导出。如果选择Export To Workspace，则将滤波器系数以数组的形式存储到MATLAB的工作空间Workspace里； 选择Export To Coefficient File(ASCII)，则MATLAB将会生成一个后缀为fcf的文件，该文件相当于滤波器设计报告，详细描述了滤波器的算法、系数和类型等信息； 选择Export To MATFile，将滤波器系数存为后缀为mat的MATLAB数据文件； 选择Export To SPTool，则存储为供信号处理工具软件SPTool(Signal Processing Tool)使用的数据格式,如图3.42所示。




图3.41滤波器的系数




图3.42导出滤波器系数




下面的fcf文件就是12阶低通滤波器的相关信息(Hanning Window，Fc=7kHz，Fs=48kHz)。




%Generated by MATLAB(R) 7.6 and the Signal Processing Toolbox 6.9

%

%Generated on: 18-Jul-2017 14:40:56

%

%Coefficient Format: Decimal

%Discrete-Time FIR Filter (real)    








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

%Filter Structure  : Direct-Form FIR

%Filter Length     : 13             

%Stable            : Yes            

%Linear Phase      : Yes (Type 1)   



Numerator:                           

-0                                   

-0.0042038899084595278               

-0.009890319761490041                

0.020185928072081524                

0.1146397834545644                  

0.23426883725633207                 

0.28999932177394327                 

0.23426883725633207                 

0.1146397834545644                  

0.020185928072081524                

-0.009890319761490041                

-0.0042038899084595278               

-0      




上面介绍了如何使用FDATool工具来设计数字滤波器。可以看出FDATool是一种非常强大和完备的工具，不仅可以按照要求设计数字滤波器，还可以对数字滤波器的各种特性进行分析。在FDATool中需要输入的各种参数和限制条件非常丰富，甚至多于事先给定的技术指标，此时需要设计者进行合理比较和取舍。FDATool工具极大减轻了设计者的劳动强度，让设计者可以把精力集中在技术含量更“高”的工作上面，如根据应用需求、性价比或硬件资源等因素选择滤波器的类型，确定滤波器的技术指标等。

3.8参考文献


［1］李素芝，万建伟.时域离散信号处理［M］.长沙： 国防科技大学出版社,1994.

［2］万建伟，王玲.信号处理仿真技术［M］.长沙： 国防科技大学出版社,2008.


本章用到的MATLAB函数总结



函数名称及调用格式函 数 用 途对应章节数

［sos,g］=tf2sos(b,a)直接型转级联型3.1节（例3.1）
［C,B,A］ = tf2par(b,a)直接型转并联型（自定义函数）3.1节（例3.1）
I = cplxcomp(p1,p2)计算复数极点p1变为p2后留数的新序号（自定义函数）3.1节（例3.1）
［C,B,A］ = tf2fs(h)直接型转频率采样型（自定义函数）3.4节（例3.14）
［z,p,k］ = buttap(n)产生低通模拟Butterworth滤波器3.2节（例3.2，例3.8）
［n,Wn］ =buttord(Wp,Ws,Rp,Rs)计算Butterworth滤波器的阶数和截止频率3.2节（例3.7，例3.8）
3.3节（例3.11，例3.12，例3.13）
［b,a］ = butter(n,Rp,Wn)设计n阶Butterworth完全滤波器3.2节（例3.7）

3.3节（例3.11，例3.12，例3.13）
续表



函数名称及调用格式函 数 用 途对应章节数

［z,p,k］ =cheb1ap(n,Rp)产生低通模拟Chebyshev Ⅰ型滤波器3.2节（例3.3，例3.10）
［z,p,k］ = (n,Rs)产生低通模拟Chebyshev Ⅱ型滤波器3.2节（例3.4，例3.9）
［z,p,k］ =ellipap(n,Rp,Rs)产生低通模拟椭圆型滤波器3.2节（例3.5，例3.6）
［bt,at］ = lp2lp(b,a,Wo)从低通到低通的转换3.2节（例3.6，例3.10）
［bt,at］ = lp2hp(b,a,Wn)从低通到高通的转换3.2节（例3.6，例3.8，例3.9）
［bt,at］=lp2bp(b,a,Wn,Bw)从低通到带通的转换3.2节
［bt,at］=lp2bs(b,a,Wn,Bw)从低通到带阻的转换3.2节
freqs(num,den)模拟滤波器频率响应3.2节（例3.2，例3.4，例3.5，例3.6）
freqz(num,den)数字滤波器频率响应3.2节（例3.7，例3.8，例3.9，例3.10）

3.3节（例3.11，例3.12，例3.13）

3.5节（例3.15，例3.21）
［dB,mag,pha,grd,w］=freqz_m(b,a)计算数字滤波器的幅频响应、相位响应和群延迟（自定义函数）
3.2节（例3.7）
3.5节（例3.15，例3.16，例3.17，例3.18，例3.19，例3.20）
［bz,az］ = impinvar(b,a,fs)冲激响应不变法将模拟滤波器变换为数字滤波器
3.2节（例3.7，例3.10）
［zd,pd,kd］ = bilinear(z,p,k,fs)双线性变换法将模拟滤波器变换为数字滤波器3.2节（例3.8）
epsMATLAB中绝对值最小的数，使0/0有意义3.5节
hd=ideal_lp(wc,N)计算理想低通滤波器的单位脉冲响应（自定义函数）3.5节（例3.15）
hd=ideal_hp(wc,N)计算理想高通滤波器的单位脉冲响应（自定义函数）3.5节（例3.16）
hd=ideal_bp(wc1,wc2,N)计算理想带通滤波器的单位脉冲响应（自定义函数）
3.5节（例3.17）
hd=ideal_bs(wc1,wc2,N)
计算理想带阻滤波器的单位脉冲响应（自定义函数）
3.5节（例3.18）
［Hr,w,a,L］=hr_type1(h)
计算1型滤波器设计的振幅响应（自定义函数）
3.5节（例3.20）
［Hr,w,a,L］=hr_type2(h)
计算2型滤波器设计的振幅响应（自定义函数）
3.5节（例3.19）
［Hr,w,a,L］=hr_type3(h)
计算3型滤波器设计的振幅响应（自定义函数）
3.5节
［Hr,w,a,L］=hr_type4(h)
计算4型滤波器设计的振幅响应（自定义函数）3.5节