实验6 数字滤波器设计

实验6 数字滤波器设计
设计传输函数G （z ），使其频率响应)(jw e G 逼近于给定的频率响应指标，这个过程称为数字滤波器的设计。

在得到G （z ）之后，就以一种适当的滤波器结构的形式实现它。

在前面的实验练习中，已经考虑过了有限冲激响应和无限冲激响应传输函数的实现。

在这个实验练习中，我们将学会设计满足某个给定幅度或增益响应的无限冲激响应或有限冲激响应数字滤波器。

1． 无限冲激响应滤波器设计
最常用于设计无限冲激响应滤波器的方法是基于原型模拟传输函数的双线性变换。

模拟传输函数通常是下列类型之一：巴特沃斯型，切比雪夫I ，切比雪夫II 和椭圆传输函数。

这些滤波器类型之间的差别可通过考虑模拟低通滤波器来说明。

（1）无限冲激响应滤波器的阶数的估计
滤波器设计过程中的第一步是，选择接近所使用的滤波器的类型，然后由滤波器指标来估计 传输函数的阶数。

用来估计巴特沃斯型滤波器的阶数的MATLAB 命令是： [N,Wn]=buttord(Wp,Ws,Rp,Rs)
其中输入参数是归一化通带边界频率wp,归一化阻带频率ws 、单位为Db 的通带波纹RP 和单位为Db 的J 最小阻带衰减纹RS 。

由于抽样频率被假定为2HZ ，WP 和WS 均必须是0和1之间的一个数。

输出数据是满足指标的最低阶数N 和归一化截止频率Wn 。

若RP=3（Db ）,则WN=WP 。

Buttord 也可用于估计高通，带通和带阻的巴特沃斯型滤波器的阶数。

（2）无限冲激响应滤波器设计
在选择了滤波器类型并估计了其阶数后，下一步是确定滤波器的传输函数，为此，MATLAB 对所有的四种类型都提供了函数。

（实验指导书，92） [num,den]=butter(N,Wn)
其中输入参数N 和Wn 可以使用函数buttord 确定，输出是向量num 和den ，它们分别是以z -1的升幂排列的传输函数的分子和分母多项式的系数。

若Wn 是一个标量，butter 返回一个N 阶的低通传输函数，基Wn 是一个双元素向量，它返回一个2N 阶的带通传输函数。

设计N 阶巴特沃斯型数字高通滤波器的MATLAB 命令是：
[num,den]=butter(N,Wn,’high ’)
然而，若Wn 是一个双元素向量，命令
[num,den]=butter(N,Wn,’stop ’)
返回一个2N 阶巴特沃斯型带阻滤波器的传输函数。

程序6.1 说明巴特沃斯型带阻滤波器的设计 %巴特沃斯型带阻滤波器的设计
Ws=[0.4 0.6];Wp=[0.3 0.7];Rp=0.4;Rs=50; %估计滤波器阶数
[N1,Wn1]=buttord(Wp,Ws,Rp,Rs); %设计滤波器
[num,den]=butter(N1,Wn1,'stop'); %显示传输函数
disp('分子系数是');disp(num);
disp('分母系数是');disp(den);
分子系数是
Columns 1 through 8
0.0330 0.0000 0.2972 0.0000 1.1889 0.0000 2.7741 0.0000 Columns 9 through 16
4.1611 0.0001 4.1611 0.0000 2.7741 0.0000 1.1889 0.0000 Columns 17 through 19
0.2972 0.0000 0.0330
分母系数是
Columns 1 through 8
1.0000 0.0000
2.6621 0.0000 4.1451 0.0001 4.1273 0.0001 Columns 9 through 16
2.8977 0.0000 1.4381 0.0000 0.5027 0.0000 0.1178 0.0000 Columns 17 through 19
0.0167 0.0000 0.0011
低通巴特沃斯模拟滤波器设计
要求指标
通带截止频率：fp=3400hz,通带最大衰减：Rp=3Db;
阻带截止频率：fs=4000hz,阻带最大衰减：Rs=40db.
Matlab程序如下
fp=3400;fs=4000;Rp=3;Rs=40;%输入滤波器的指标
[N,fc]=buttord(fp,fs,Rp,Rs,'s');%计算阶数N和3DB截止频率fc
[B,A]=butter(N,fc,'s');%设计低通巴特沃斯模拟滤波器/传输函数的系数
[hf,f]=freqs(B,A,1024); %计算模拟滤波器的频率响应
plot(f,20*log(abs(hf)/abs(hf(1))));
xlabel('f/hz');ylabel('H(db)');
axis([0,4000,-40,5]);
line([0,4000],[-3,-3]);%点坐标，
line([3400,3400],[-90,5]);
这个29阶的滤波器频率特性如图所示：图中幅频曲线在3400HZ处为-3DB，4000处的衰减-40DB，满足设计指标。

line 创建线对象
line([0,4000],[-3,-3]);
line([3400,3400],[-90,5])
>>
>>
2． 有限冲激响应滤波器设计
有限冲激响应滤波器设计的最简单的方法是，对期望的理想频率响应进行离散时间傅里叶逆变换，得到无限长冲激响应，对所得的双无限系数最单地截尾为有限系数。

然而，简单截尾得到的有限冲激响应滤波器的幅度响应中会有一个振荡行为，它通常称为吉布斯现象。

吉布斯现象可用一个合适的有限长窗口函数对无限长冲激响应系数加窗来减小。

MA TLAB 函数fir1和fir2可以用来设计加窗的有限冲激响应数字滤波器。

两个函数均可以产生一 个线性相位设计。

函数fir1可用于设计常规的低通、高通、带通和带阻线性相位有限冲激响应滤波器，对于抽样频率为2HZ 的情况，
命令b=fir1(N,Wn) 在向量b 中返回以z -1的升幂排列的N 阶低通或带通滤波器的冲激响应系数。

详细说明参见实验教程95页。

对于低通设计，归一化截止频率由标量Wn 给定，它是在0和1之间的一个数，对于带通的设计，Wn 是包含指定通带边界的一个双元素向量[Wn1，Wn2]，
其中0<Wn1<Wn2<1。

命令b=fir1(N, Wn,’high ’) , 其中N 为一个偶数，用于设计高通滤波器。

命令b=fir1(N, Wn,’stop ’) , 其中Wn 是一个双元素向量，用于设计带阻有限冲激响应滤波器，若没有指定，用海明窗作为默认值。

命令b=fir1(N, Wn,taper) ,使用在向量taper 中给定的长度为N+1的窗系数。

然而，窗函数必须事先用一个适当的MATLAB 函数来产生，如blackman, hamming, hanning, chebwin 或kaiser.
用到的命令为
taper=blackman(N) taper= hamming (N) taper= hanning (N) taper= chebwin (N) taper= kaiser (N,beta).
最简单元的窗函数是矩形窗，从阻带衰减的观点来看也是性能最差的一种，设计一个线性相位FIR 低通滤波器。

要求通带截止频率rad w c 4/π=，
单位冲激响应h(n)的长度N=21。

MATLAB 程序如下
>> N=21; wc=pi/4;
n=0:N-1;r=(N-1)/2;
hdn=sin(wc*(n-r))/pi./(n-r);%计算理想低通单位冲激响应hdn,具体参见教材P141
subplot(2,1,1);stem(n,hdn);
xlabel('time n');ylabel('hd(n)');title('impulse response');
if rem(N,2)~=0 hdn(r+1)=wc/pi;end %rem()求余数
wn=boxcar(N);%矩形窗
hn=hdn.*wn'; %加窗
hw=fft(hn,512); %幅频特性
w=2*[0:511]/512;
subplot(2,1,2);plot(w,20*log10(abs(hw)));
xlabel('w/pi');ylabel('manitude(db)');
title('frequency response');
程序运行的图像如下：
作业：
设计一个线性相位FIR低通滤波器。

要求通带截止频率wc=pi/7，单位冲激响应h(n)的长度N=29,分别使用海明窗及汉明窗进行设计。