实验六用窗函数法设计FIR数字滤波器.doc

实验六 用窗函数法设计FIR数字滤波器
一. 实验目的
(1) 掌握用窗函数法设计FIR数字滤波器的原理与方法。
(2) 熟悉线性相位FIR数字滤波器的特性。
(3) 了解各种窗函数对滤波特性的影响。
二. 实验内容和要求
(1) 复习用窗函数法设计FIR数字滤波器一节内容，阅读本实验原理，掌握设计步骤。
(2) 用升余弦窗设计一线性相位低通FIR数字滤波器，截止频率。窗口长度N =15，33。要求在两种窗口长度情况下，分别求出，打印出相应的幅频特性和相频特性曲线，观察3dB带宽和20dB带宽。总结窗口长度N 对滤波器特性的影响。
设计低通FIR数字滤波器时，一般以理想低通滤波特性为逼近函数，即

其中

(3) ，，用四种窗函数设计线性相位低通滤波器，绘制相应的幅频特性曲线，观察3dB带宽和20dB带宽以及阻带最小衰减，比较四种窗函数对滤波器特性的影响。
三. 实验方法、步骤及结果测试
如果所希望的滤波器的理想的频率响应函数为，则其对应的单位脉冲响应为
（4.1）
窗函数设计法的基本原理是用有限长单位脉冲响应序列逼近。由于往往是无限长序列，而且是非因果的，所以用窗函数将截断，并进行加权处理，得到：
（4.2）
就作为实际设计的FIR数字滤波器的单位脉冲响应序列，其频率响应函数为
（4.3）
式中，N为所选窗函数的长度。
我们知道，用窗函数法设计的滤波器性能取决于窗函数的类型及窗口长度N的取值。设计过程中，要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度N 。各种类型的窗函数可达到的阻带最小衰减和过渡带宽度见表4.1。
表1 各种窗函数的基本参数
窗函数
旁瓣峰值幅度/dB
过渡带宽
阻带最小衰减/dB
矩形窗
-13
4π/N
-12
三角形窗
-25
8π/N
-25
汉宁窗
-31
8π/N
-44
哈明窗
-41
8π/N
-53
不莱克曼窗
-57
12π/N
-74
凯塞窗(α=7.865)
57
10π/N
-80
这样选定窗函数类型和长度N之后，求出单位脉冲响应，并按照式（4.3）求出。是否满足要求，要进行演算。一般在尾部加零使长度满足2的整数次幂，以便用FFT计算。如果要观察细节，补零点数增多即可。如果不满足要求，则要重新选择窗函数类型和长度N ，再次验算，直至满足要求。
如果要求线性相位特性，则还必须满足

根据上式中的正、负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。要根据所设计的滤波特性正确选择其中一类，例如，要设计线性相位低通特性，可以选择这一类，而不能选择这一类。
主程序框图如图4.1所示。其中幅度特性要求用dB表示。
开始
读入窗口长度N
计算hd(n)
调用窗函数子程序求w(n)
调用子程序（函数）计算H(k)=DFT[h(n)]
调用绘图子程序（函数）绘制H(k)幅度相位曲线
结束
图6-1 主程序框图
计算h(n)= hd(n) w(n)

设

画图时，用打印幅度特性。第k点对应的频率。为使曲线包络更接近的幅度特性曲线，DFT变换区间要选大些。例如窗口长度N=33时，可通过在末尾补零的方法，使长度变为64，再进行64点DFT，则可以得到更精确的幅度衰减特性曲线。

下面给出MATLAB主程序：
%实验四，用窗函数法设计FIR数字滤波器
b=1;
close all;
i=0;
while(b);
temp=menu('选择窗函数长度N','N=10','N=15','N=20','N=25','N=30','N=33','N=35','N=40','N=45','N=50','N=55','N=60','N=64');
menu1=[10,15,20,25,30,33,35,40,45,50,55,60,64];
N=menu1(temp);
temp=menu('选择逼近理想低通滤波器截止频率Wc','Wc=pi/4','Wc=pi/2','Wc=3*pi/4','Wc=pi','Wc=0.5','Wc=1.0','Wc=15','Wc=2.0','Wc=2.5','Wc=3.0');
menu2=[pi/4,pi/2,3*pi/4,pi,0.5,1,1.5,2,2.5,3];
w=menu2(temp);
n=[0:(N-1)];
hd=ideal(w,N); %得到理想低通滤波器
k=menu('请选择窗口类型:','boxcar','hamming','hanning','blackman');
if k==1
B=boxcar(N);
string=['Boxcar','N=',num2str(N)];
else if k==2
B=hamming(N);
string=['Hamming','N=',num2str(N)];
else if k==3
B=hanning(N);
string=['Hanning','N=',num2str(N)];
else if k==4
B=blackman(N);
string=['Blackman''N=',num2str(N)];
end
end
end
end
h=hd.*(B)'; %得到FIR数字滤波器
[H,m]=freqz(h,[1],1024,'whole'); %求其频率响应
mag=abs(H); %得到幅值
db=20*log10((mag+eps)/max(mag));
pha=angle(H); %得到相位
i=i+1;
figure(i)
subplot(2,2,1);
n=0:N-1;
stem(n,h,'.');
axis([0,N-1,-0.1,0.3]);
hold on;
n=0:N-1;
x=zeros(N);
plot(n,x,'-');
xlabel('n');
ylabel('h(n)');
title('实际低通滤波器的h(n)');
text((0.3*N),0.27,string);
hold off;
subplot(2,2,2);
plot(m/pi,db);
axis([0,1,-100,0]);
xlabel('w/pi');
ylabel('dB');
title('衰减特性(dB)');
grid;
subplot(2,2,3);
plot(m,pha);
hold on;
n=0:7;
x=zeros(8);
plot(n,x,'-');
title('相频特性');
xlabel('频率(rad)');
ylabel('相位(rad)');
axis([0,3.15,-4,4]);
subplot(2,2,4);
plot(m,mag);
title('频率特性');
xlabel('频率W(rad)');
ylabel('幅值');
axis([0,3.15,0,1.5]);
text(0.9,1.2,string);
b=menu('Do You want To Continue ?','Yes','No');
if b==2
b=0;
end
end
temp=menu('Close All Figure ?','Yes','No');
if temp==1
close all
end
程序运行结果：
运行程序，根据实验内容要求和程序提示选择你要进行的实验参数。三个实验参数选定后，程序运行输出用所选窗函数设计的实际FIR低通数字滤波器的单位脉冲响应h(n)、幅频衰减特性(20lg|H(ejw)|)、相频特性及幅频特性|H(ejw)|的波形，h(n)和|H(ejw)|图中标出了所选窗函数类型及其长度N值。对四种窗函数（N=15和N=33）的程序运行结果如图4-2到图4-9所示，由图可以看出用各种窗函数设计的FIR滤波器的阻带最小衰减及过渡带均与教材中一致。在通带内均为严格相位特性。
图6-2 矩形窗（N=15）

图6-3 矩形窗（N=33）

图6-4 哈明窗（N=15）

图6-5 哈明窗（N=33）

图6-6 汉宁窗（N=15）

图6-7 汉宁窗（N=33）

图6-8 布莱克曼窗（N=15）

图6-9 布莱克曼窗（N=33）

四. 实验报告要求
(1) 简述实验原理及目的。
(2) 按照实验步骤以及要求，比较各种情况下的滤波性能，说明窗口长度N和窗函数类型对滤波特性的影响。
(3) 总结用窗函数法设计FIR滤波器的主要特点。
(4) 简要回答思考题。
五. 思考题
(1) 如果给定通带截止频率和阻带截止频率以及阻带最小衰减，如何用窗函数法设计线性相位低通滤波器，写出设计步骤。

(2) 如果要求用窗函数法设计带通滤波器，而且给定上、下边带截止频率为和，试求理想带通的单位脉冲响应
16

%实验中的子函数：产生理想低通滤波器单位脉冲响应hd(n)
function hd=ideal(w,N);
alpha=(N-1)/2;
n=[0:(N-1)];
m=n-alpha+eps;
hd=sin(w*m)./(pi*m);

相关练习：
例1 用凯塞窗设计一FIR低通滤波器，低通边界频率 ，阻带边界频率 ，阻带衰减 不小于50dB。
解 首先由过渡带宽和阻带衰减 来决定凯塞窗的N和

图4.1给出了以上设计的频率特性，(a) 为N=30直接截取的频率特性(b)为凯塞窗设计的频率特性。凯塞窗设计对应的MATLAB程序为：
wn=kaiser(30,4.55);
nn=[0:1:29];
alfa=(30-1)/2;
hd=sin(0.4*pi*(nn-alfa))./(pi*(nn-alfa));
h=hd.*wn';
[h1,w1]=freqz(h,1);
plot(w1/pi,20*log10(abs(h1)));
axis([0,1,-80,10]);
grid;
xlabel('归一化频率/p')
ylabel('幅度/dB')

例2 利用雷米兹交替算法，设计一个线性相位低通FIR数字滤波器，其指标为：通带边界频率fc=800Hz，阻带边界fr=1000Hz，通带波动 阻带最小衰减At=40dB，采样频率fs=4000Hz。
解
在MATLAB中可以用remezord 和remez两个函数设计，其结果如图4.2，MATLAB程序如下：
fedge=[800 1000];
mval=[1 0];
dev=[0.0559 0.01];
                 fs=4000;
[N,fpts,mag,wt]=remezord(fedge,mval,dev,fs);
b=remez(N,fpts,mag,wt);
[h,w]=freqz(b,1,256);
plot(w*2000/pi,20*log10(abs(h)));
grid; 
                                      xlabel('频率/Hz')
ylabel('幅度/dB')

        函数remezord中的数组fedge为通带和阻带边界频率，数组mval是两个边界处的幅值，而数组dev是通带和阻带的波动，fs是采样频率单位为Hz。