六用窗函数法设计FIR数字滤波器

来源：网络收集时间：2025-04-26 下载这篇文档手机版

说明：文章内容仅供预览，部分内容可能不全，需要完整文档或者需要复制内容，请下载word后使用。下载word有问题请添加微信号:xuecool-com或QQ：370150219 处理（尽可能给您提供完整文档），感谢您的支持与谅解。

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

一. 实验目的

(1)掌握用窗函数法设计FIR数字滤波器的原理与方法。 (2)熟悉线性相位FIR数字滤波器的特性。 (3)了解各种窗函数对滤波特性的影响。二. 实验内容和要求

(1) 复习用窗函数法设计FIR数字滤波器一节内容，阅读本实验原理，掌握设计步骤。

(2) 用升余弦窗设计一线性相位低通FIR数字滤波器，截止频率

??c?rad。窗口长度N =15，33。要求在两种窗口长度情况下，分别

4求出h?n?，打印出相应的幅频特性和相频特性曲线，观察3dB带宽和20dB带宽。总结窗口长度N 对滤波器特性的影响。

设计低通FIR数字滤波器时，一般以理想低通滤波特性为逼近函数Hej?，即

Hdej?N?1 2?????j???e，???c??? ??0，?c????其中??1hd?n??2?????Hde??j?1ed??2?j?????ce?j??ej??d??csin??c?n?a?? ??n?a?(3) N?33，?c??4，用四种窗函数设计线性相位低通滤波器，绘制相应的幅频特性曲线，观察3dB带宽和20dB带宽以及阻带最小衰减，比较四种窗函数对滤波器特性的影响。

三. 实验方法、步骤及结果测试

如果所希望的滤波器的理想的频率响应函数为Hdej?，则其对应的单位脉冲响应为

hd?n??12???????Hdej?ej?d? （4.1）

??窗函数设计法的基本原理是用有限长单位脉冲响应序列h?n?逼近由于hd?n?往往是无限长序列，而且是非因果的，所以用窗函数??n?hd?n?。

将hd?n?截断，并进行加权处理，得到：

h?n??hd?n???n?

（4.2）

其频率h?n?就作为实际设计的FIR数字滤波器的单位脉冲响应序列，响应函数Hej?为

??He????h?n?ej?n?0N?1j?n （4.3）

式中，N为所选窗函数??n?的长度。

我们知道，用窗函数法设计的滤波器性能取决于窗函数??n?的类型及窗口长度N的取值。设计过程中，要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度N 。各种类型的窗函数可达到的阻带最小衰减和过渡带宽度见表4.1。

表1 各种窗函数的基本参数

窗函数矩形窗三角形窗

旁瓣峰值幅度/dB 过渡带宽阻带最小衰减/dB -13 -25 4π/N 8π/N -12 -25 2

汉宁窗哈明窗不莱克曼窗凯塞窗(α=7.865) -31 -41 -57 -57 8π/N 8π/N 12π/N 10π/N -44 -53 -74 -80 这样选定窗函数类型和长度N之后，求出单位脉冲响应h?n??hd?n????n?，并按照式（4.3）求出Hej?。Hej?是否满足要求，要进行演算。一般在h?n?尾部加零使长度满足2的整数次幂，以便用FFT计算Hej?。如果要观察细节，补零点数增多即可。如果Hej?不满足要求，则要重新选择窗函数类型和长度N ，再次验算，直至满足要求。

如果要求线性相位特性，则h?n?还必须满足

????????h?n???h?N?1?n?

根据上式中的正、负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。要根据所设计的滤波特性正确选择其中一类，例如，要设计线性相位低通特性，可以选择h?n??h?N?1?n?这一类，而不能选择

h?n???h?N?1?n?这一类。

主程序框图如图4.1所示。其中幅度特性要求用dB表示。

开始读入窗口长度N 计算hd(n) 调用窗函数子程序求w(n) 计算h(n)= hd(n) w(n) 调用子程序（函数）计算H(k)=DFT[h(n)] 调用绘图子程序（函数）绘制H(k)幅度相位曲线结束图6-1 主程序框图

设

H(k)?DFT[h(n)]H(k)?HR(k)?jHI(k)

2H(k)?HR(k)?HI2(k)画图时，用20lgH(k)打印幅度特性。第k点对应的频率?k?2?k。N为使曲线包络更接近Hej?的幅度特性曲线，DFT变换区间要选大些。例如窗口长度N=33时，可通过在h?n?末尾补零的方法，使长度变为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=1.5','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]);

求其频率响应 6

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

程序运行结果：

图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) 如果要求用窗函数法设计带通滤波器，而且给定上、下边带截止频率为?1和?2，试求理想带通的单位脉冲响应hd?n?

%实验中的子函数：产生理想低通滤波器单位脉冲响应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);