数字滤波器的设计

        设计传递函数, 使其频率响应逼近给定的频率响应指标，这个过程称为数字滤波器设计。得到

传递函数后，就可以一种适当的滤波器结构实现它。

一、滤波器的指标 

        滤波器指标通常按照其幅度响应指定。一般有以下两种，左图为归一化幅度响应，右图为典型的幅度响应。左图中，频率 ωp 和 ωs 分别称为通带边界频率和阻带边界频率。在由 0 ≤ ω ≤ ωp定义的通带中, 幅度的最大值和最小值分别是 1 和  。单位为dB的峰通带波纹是：

由 ωs ≤ |ω| ≤ π 定义的阻带中, 最大波纹由1/A表示; 单位为dB的最小阻带衰减为：

右图是另一种典型的幅度响应指标。在由 0 ≤ ω ≤ ωp 定义的通带中, 幅度值在 1+δp 和 1-δp 间，即幅度以误差 ±δp 逼近1。由 ωs ≤ |ω| ≤ π 定义的阻带中，幅度以误差 δs 接近于0。 δp 和 δs 称为波纹。 

        

  二、IIR滤波器的设计

        最常用的方法基于对原型模拟传递函数的双线性变换。模拟传递函数通常是下列类型之一:

巴特沃兹、切比雪夫1型、切比雪夫2型和椭圆传递函数。

 

（1）第一步：选择所使用的滤波器的类型, 然后由滤波器指标来估计传递函数的阶数。用来估计巴特沃兹滤波器的阶数的MATLAB命令是：

[N, Wn] = buttord( Wp, Ws, Rp, Rs)

其中输入参数是归一化通带边界频率Wp、归一化阻带边界频率Ws、单位为dB的通带波纹Rp和单位为dB的最小阻带衰减Rs。归一化频率Wp和Ws是( 0 , 1 )之间的一个数，1 对应π（弧度/样本）。buttord函数的输出是满足指标的最低阶数N和归一化截止频率Wn（3dB频率），若Rp=3dB,则Wn=Wp。

buttord也可用于估计高通、带通和带阻巴特沃兹滤波器的阶数。对于带通和带阻滤波器设计, Wp和Ws是指定边界频率的双元素向量。

用来估计切比雪夫1型、切比雪夫2型和椭圆滤波器的阶数的MATLAB命令分别是：

[N, Wn] = cheb1ord (Wp,Ws,Rp,Rs)

[N,Wn] = cheb2ord (Wp,Ws,Rp,Rs)

[N,Wn] = ellipord (Wp,Ws,Rp,Rs)

（2）第二步：在选择了滤波器类型并估计了其阶数后，下一步是确定滤波器的传递函数。为此, MATLAB对所有四种类型的滤波器提供了函数。设计巴特沃兹数字低通或带通滤波器的MATLAB命令是

[num, den] = butter(N, Wn)

其中输入参数N和Wn可使用函数buttord确定, 输出是向量num和den,它们分别是以的升幂排列的传递函数的分子、分母多项式的系数。若Wn是一个标量, butter返回一个N阶的低通传递函数。若Wn是一个双元素向量, butter返回一个2N阶的带通传递函数。

 

设计N阶巴特沃兹数字高通滤波器的命令是：

[num, den] = butter(N, Wn,’high’)

设计巴特沃兹带阻滤波器的命令是：(其中Wn是双元素向量)

[num, den] = butter(N, Wn,’stop’)

设计另外三种数字滤波器的MATLAB命令是：

[num, den] = cheby1(N, Rp,Wn)

[num, den] = cheby1(N, Rp,Wn,’filtertype’)

[num, den] = cheby2 (N, Rs,Wn)

[num, den] = cheby2 (N, Rs,Wn,’filtertype’)

[num, den] = ellip (N, Rp,Rs,Wn)

[num, den] = ellip (N, Rp,Rs,Wn,’filtertype’)

 例1：确定一个数字无限冲激响应低通滤波器所有四种类型的最低阶数。指标如下：40kHz的抽样率，4kHz的通带边界频率，8kHz的阻带边界频率，0.5dB的通带波纹，40dB的最小阻带衰减。

clc;clear;close all;
Ft = 40000;
fp = 4000;
fs = 8000;
Rp = 0.5;
Rs = 40;
% buttord ( )函数中 Wp、Ws 是数字频率且对π的归一化
Wp = (2*pi*fp / Ft)/pi;
Ws = (2*pi*fs / Ft)/pi;
% 四个函数的输出是满足指标的最低阶数N和归一化截止频率Wn（3dB频率）
[N1, Wn1] = buttord (Wp,Ws,Rp,Rs);
[N2, Wn2] = cheb1ord (Wp,Ws,Rp,Rs);
[N3, Wn3] = cheb2ord (Wp,Ws,Rp,Rs);
[N4, Wn4] = ellipord (Wp,Ws,Rp,Rs);

例2：设计符合问题6.1中指标的切比雪夫1型低通滤波器，写出所产生的传递函数的准确表达式。

clc;clear;close all;
% 切比雪夫1型低通滤波器
Ws = 0.4; Wp = 0.2; Rp = 0.5; Rs = 40;      
[N1, Wn1] = cheb1ord (Wp,Ws,Rp,Rs);     % 估计滤波器阶数
[num,den] = cheby1(N1,Rp,Wn1,'low');       % 设计滤波器
% 显示传递函数
disp('分子系数是 ');disp(num);
disp('分母系数是 ');disp(den);
[g, w] = gain(num,den);      % 计算增益响应
plot(w/pi,g);       % 绘制增益响应
axis([0 1.1 -60 5]);      grid;
xlabel('\omega/\pi'); ylabel('增益，dB');
title('切比雪夫1型低通滤波器的增益响应');

附：求增益响应的函数 gain 

function [g,w] = gain(num,den)
% 在单位圆上半部分的256个等分点，
% 计算传递函数的增益响应（dB） 
% 分子系数在向量 num 中
% 分母系数在向量 den 中
% 频率值以向量w返回
% 增益数值以向量g返回
w = 0:pi/255:pi;
H = freqz(num,den,w);
g = 20*log10(abs(H));

 

根据 MATLAB 的输出数据可以得到传递函数分子、分母对应的系数，即可得到传递函数：

 三、FIR滤波器的设计

（1）第一步：估计传递函数阶数。对于FIR数字滤波器的设计，可以由下列数字滤波器指标来直接估计最小滤波器阶数N：归一化通带边界角频率ωp, 归一化阻带边界角频率ωs, 峰通带波纹δp, 峰阻带波纹δs。一个较简单的逼近式为凯泽公式：

其中Δω=|ωs-ωp|是过渡带的宽度。

在许多滤波器设计问题中，估算的阶数可能得到不满足给定指标的滤波器。在这种情况下, 需要逐渐增加N值, 直到满足给定指标。

函数kaiord可根据指标估计线性相位低通FIR滤波器的阶数

function N = kaiord(Fp, Fs, dp, ds, FT)
% 使用凯泽公式计算线性相位有限冲激响应多频带滤波器的阶数
% dp是通带波纹
% ds是阻带波纹
% Fp是单位为Hz的通带边界
% Fs是单位为Hz的阻带边界
% FT是单位为Hz的抽样频率
% N是估计出的有限冲激响应滤波器的阶数
if nargin == 4,
	FT = 2;
end
if length(Fp) > 1,
	TBW = min(abs(Fp(1) - Fs(1)), abs(Fp(2) - Fs(2)));
else
	TBW = abs(Fp - Fs);
end
num = -20*log10(sqrt(dp*ds))-13;
den = 14.6*TBW/FT;
N = ceil(num/den);

（2）FIR滤波器设计的最简单方法是, 对期望的理想频率响应进行离散时间傅里叶逆变换, 得到双边无限长冲激响应, 对所得的双边无限系数简单地截尾为有限系数。然而, 简单截尾得到的有限冲激响应滤波器的幅度响应中会有一个振荡行为, 它通常称为吉布斯现象。吉布斯现象可用一个合适的有限长窗函数对双边无限长冲激响应系数加窗来减小。

MATLAB函数fir1可以用来设计加窗的有限冲激响应数字滤波器，并产生一个线性相位设计。fir1可用于设计常规的低通、高通、带通和带阻线性相位FIR滤波器。命令

b=fir1(N,Wn)

在向量b中返回以的升幂排列的N阶低通或带通滤波器的冲激响应系数。对于低通设计, 归一化截止频率Wn是在0和1之间的一个数。

例3：观察吉布斯效应：使用函数sinc，以产生截止频率在ωc＝0.4π处、长度分别为81，61，41和21的四个低通滤波器的冲激响应系数，然后计算并画出它们的幅度响应。可使用冒号“:”运算符从长度为81的滤波器的冲激响应系数中抽出较短长度滤波器的冲激响应系数。

clc;clear;close all;
wc=0.4*pi; 
n = -40:40;     % length of hn is 81 
% sinc 产生的函数是 sin(πx)/(πx)
hn_81 = 0.4 * sinc(wc*n/pi);
subplot(2,2,1);     stem(n,hn_81);
title('序列长度为81的冲激响应系数的幅度响应');
hn_61 = hn_81(11:71);
subplot(2,2,2);     stem(hn_61);
title('序列长度为61的冲激响应系数的幅度响应');
hn_41 = hn_81(21:61);
subplot(2,2,3);     stem(hn_41);
title('序列长度为41的冲激响应系数的幅度响应');
hn_21 = hn_81(31:51);
subplot(2,2,4);     stem(hn_21);
title('序列长度为21的冲激响应系数的幅度响应');

例4：使用函数kaiord，估计具有以下指标的线性相位低通FIR滤波器的阶数：通带边界为2kHz，阻带边界为2.5kHz，通带波纹δp=0.005，阻带波纹δs=0.005，抽样率为10kHz。

clc;clear;close all;
dp = 0.005;
ds = 0.005;
Fp = 2000;
Fs = 2500;
FT = 10000;
N = kaiord(Fp, Fs, dp, ds, FT);

可以通过改变 阻带边界、通带波纹δp、阻带波纹δs、抽样率FT等来考虑阶数受哪些因素的影响

例4：使用函数fir1，设计一个线性相位FIR低通滤波器，使用 例3 的系数和指标显示滤波器系数，并画出其增益和相位响应。

clc;clear;close all;
dp = 0.005;
ds = 0.005;
Fp = 2000;
Fs = 2500;
FT = 10000;
N = kaiord(Fp, Fs, dp, ds, FT);

% 显示滤波器系数，并画出其增益和相位响应
Wp=2*Fp/FT;     % Wp对采样频率的一半(FT/2)归一化
Ws=2*Fs/FT;
Wn=(Wp+Ws)/2;
h = fir1(N,Wn) % % 在向量h中返回以的z^(-1)升幂排列的N阶低通或带通滤波器的冲激响应系数
w = 0:pi/255:pi;
H = freqz(h,1,w);   % 求频响
g = 20*log10(abs(H)); % 增益数值以向量g返回
subplot(2,1,1)
plot(w/pi,abs(H));  grid;
title('FIR低通滤波器的增益');
subplot(2,1,2)
plot(w/pi,angle(H));  grid;
title('FIR低通滤波器的相位相应');