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

装 订 线评阅人老师实验成绩成绩自动化学院本科生实验报告 数字信号处理 课 程 实 验 报 告实验名称 用窗函数法设计 FIR 数字滤波器 一、实验原理、目的与要求1. 实验原理如果所希望的滤波器的理想频率响应函数为，则其对应的单位脉冲响应为：用窗函数w(n)将截断，并进行加权处理，得到：h(n)就作为实际设计的FIR数字滤波器的单位脉冲响应序列，其频率响应函数为：式中，N为所选窗函数w(n)的长度。用窗函数法设计的滤波器性能取决于窗函数w(n)的类型及窗口长度N。因此，在设计过程中，要根据对阻带最小衰减和过度带宽度的要求选择合适的窗函数类型和窗口长度N。选定窗函数了形和长度N后，求出单位脉冲响应h(n)=hd(n)·w(n)，并可以求出。是否满足要求，要进行验算。一般在h(n)尾部加零使长度满足2的整数次幂，以便用FFT计算。如果要观察细节，补零点数增多即可。如果不满足要求，则要重新选择窗函数类型和长度N，再次验算，直至满足要求。如果要求线性相位特性，则h(n)还必须满足：根据上式中的正、负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。要根据所设计的滤波特性正确选择其中一类。例如，要设计线性相位低通特性，可选择h(n)=h(N-1-n)一类，而不能选h(n)=-h(N-1-n)一类。2. 实验目的（1）掌握用窗函数法设计 FIR 数字滤波器的原理和方法。（2）熟悉线性相位 FIR 数字滤波器特性。（3）了解各种窗函数对滤波特性的影响。3. 实验要求（1）简述实验目的及原理。（2）按照实验步骤及要求，比较各种情况下的滤波性能，说明窗口长度 N 和窗函数类型对滤波特性的影响。（3）总结用窗函数法设计 FIR 滤波器的主要特点。（4）简要回答思考题。二、实验仪器设备（标注实验设备名称及设备号）Windows 计算机 台号 22Matlab 软件三、实验内容步骤及结果分析1.用升余弦窗设计一线性相位低通 FIR 数字滤波器，截至频率wc = / 4 rad。窗口长度 N=15，33。要求在两种窗口长度情况下，分别求出 h(n) ，打印出相应的幅频特性和相频特性曲线，观察 3dB 带宽和 20dB 带宽。总结窗口长度 N 对滤波特性的影响。设计低通 FIR 数字滤波器时，一般以理想低通滤波特性为逼近函数，即：在FIR滤波器的窗函数的设计中，首先是由所要求的阻带最小衰减As确定窗函数的形状，也就是类型，再有过渡带宽确定窗函数的窗长N。以上是标准化的设计流程，但此处我们只研究N对于一个固定窗函数的影响。这里取升余弦窗，也就是汉宁窗，它的主瓣宽度8/N，过渡带宽是6.2/N，阻带最小衰减44dB。N取15或者33。画出N=15和33的h(n),增益，幅频和相频曲线如下：由上可见：N由15变成33,8/N减少，过渡带宽减少这可以在幅频图上看到。增益图和相位图上可以看出起伏震荡变密。而且还可以在幅频图上看出，通带也变宽了，表现在低频为1的部分变多了。过渡带宽减小是通过通带截止点变大，阻带截止点变小实现的，这一点可以在幅频图上看出。另一方面也可以说，通带最大衰减（3dB带宽）减小了，阻带最大衰减（20dB带宽）变大了，这表明了滤波器的效果更好了。同时也可以看出，同一窗函数达到阻带的衰减是一样的。2.N=33，c=/4，用四种窗函数设计线性相位低通滤波器。绘制相应的幅频特性曲线，观察3dB和20dB带宽以及阻带最小衰减，比较四种窗函数对滤波器特性的影响。四种代表不同窗（矩型窗、升余弦窗、改进升余弦窗和二阶升余弦窗）的窗函数的幅频特性，如下：把它们绘制在一张图方便比较：四种代表不同窗（矩型窗、升余弦窗、改进升余弦窗和二阶升余弦窗）的窗函数的幅频特性，如下：把它们绘制在一张图方便比较：矩形窗的旁瓣抑制性能差，旁瓣峰值大，在阻带仍有一定幅度值；汉宁窗旁瓣抑制性能稍好；海明窗和布莱克曼窗旁瓣抑制能力最强，在阻带基本没有幅度值。矩形窗的主瓣窄；汉宁窗的主瓣宽度与海明窗基本一致，但布莱克曼窗还是略微低于汉宁窗和海明窗。主瓣的宽度影响频率的分辨率，旁瓣的强度影响干扰程度。在实际工程中，应在两者之间进行权衡，选择合适的窗函数。相同的窗口长度，不同的窗函数阻带衰减不同，过渡带不同，这是选择窗函数的重要标准。四、思考题目（1）如果给定通带截止频率和阻带截止频率以及阻带最小衰减，如何用窗函数法设计线性相位低通滤波器？写出设计步骤。答：通过通带截止频率和阻带截止频率确定所设计的滤波器的过渡带宽。通过所需的阻带最小衰减查表确定所选的窗函数类型。通过所选的窗函数的过渡带宽公式，确定N值。将hd (n)与窗函数相乘得FIR数字滤波器的单位冲激响应h(n)。h(n)=hd(n)·w(n)（2） 如果要求用窗函数法设计带通滤波器，且给定上、下边带截止频率为 w1和w2 ，试求理想带通的单位脉冲响应 hd (n) 。答：理想带通滤波器的频率响应为:利用傅立叶逆变换推导出对应的单位冲激响应:五、总结针对为什么要用窗函数：每次FFT变换只能对有限长度的时域数据进行变换，因此，需要对时域信号进行信号截断。即使是周期信号，如果截断的时间长度不是周期的整数倍（周期截断），那么，截取后的信号将会存在泄漏。为了将这个泄漏误差减少到最小程度（注意我说是的减少，而不是消除），我们需要使用加权函数，也叫窗函数。加窗主要是为了使时域信号似乎更好地满足FFT处理的周期性要求，减少泄漏。但是，窗函数只能减少泄漏，不能消除泄漏。针对总结用窗函数法设计 FIR 滤波器的主要特点：FIR 滤波器目前常用的设计方法有窗函数法和频率采样法,窗函数法是从时域进行设计,而频率采样法是从频域进行设计。窗函数法由于简单、物理意义清晰,因而得到了较为广泛的应用。窗函数法设计的基本思想是：首先根据技术指标要求,选取合适的阶数 N 和窗函数的类型 w(n),使其幅频特性逼近理想滤波器幅频特性。其次,因为理想滤波器的 hd(n)是无限长的,所以需要对 hd(n) 进行截断,数学上称这种方法为窗函数法。简而言之,用窗函数法设计FIR滤波器是在时域进行的,先用傅里叶变换求出理想滤波器单位抽样相应hd(n),然后加时间窗w(n)对其进行截断,以求得FIR 滤波器的单位抽样响应h(n)。六、代码%-绘制长度为15或者33，汉宁窗构建的滤波器-clearclcwin_type = 1,2,3,4; N1=15,33;omegaC = pi/4;bianer = 6;i=2;for j=1:2N=N1(j);n = 0:N-1;M = 2ceil(log2(N)+bianer);hd = ideallp(omegaC, N);wn = choose_win_type(win_type(i), N)'h = wn.*hd; db, mag, pha, grd, w = freqz_m(h, 1, M);Rp = -min(db(1:omegaC/(pi/M)+1);As = -max(db(omegaC/(pi/M):end); figure subplot(2, 2, 1)stem(n, h)axis(0 N ylim)xlabel('n')ylabel('H(n)') subplot(2, 2, 2)plot(w/pi, db)grid on xlabel('omega/pi');ylabel('20log(|H(w)|)/dB') subplot(2, 2, 3)plot(w/pi, mag)grid onxlabel('omega/pi')ylabel('|H(w)|') subplot(2, 2, 4)plot(w/pi, pha)xlabel('omega/pi')ylabel('arg(|H(w)|)')end%-选窗的函数-function w = choose_win_type(win_type, N) switch win_type case 2 w = hanning(N); case 3 w = hamming(N); case 1 w = rectwin(N); case 4 w = blackman(N); endEnd%-hd-function hd = ideallp(wc, N) tau = (N-1)/2; n = 0:N-1; m = n - tau + eps; hd = sin(wc * m)./(pi * m);end%-取增益，幅频，相频的程序块-function db, mag, pha, grd, w = freqz_m(b, a, M) H, w = freqz(b, a, M*2, 'whole'); H = H(1:M)' w = w(1:M)' mag = abs(H); db = 20*log10(mag+eps)/max(mag); pha = angle(H); grd = grpdelay(b, a, w); end%-分别使用4种窗函数，并进行比较-clearclcwin_type = 1,2,3,4; N1=15,33;omegaC = pi/4;bianer = 6;j=2;N=N1(j);n = 0:N-1;M = 2ceil(log2(N)+bianer);hd, tau = ideallp(omegaC, N);for i=1:4wn = choose_win_type(win_type(i), N)'h = wn.*hd; db, mag, pha, grd, w = freqz_m(h, 1, M);Rp = -min(db(1:omegaC/(pi/M)+1);As = -max(db(omegaC/(pi/M):end); subplot(2, 2, i)plot(w/pi, mag)grid onxlabel('omega/pi')ylabel('|H(w)|')endfigurefor i=1:4wn = choose_win_type(win_type(i), N)'h = wn.*hd; db, mag, pha, grd, w = freqz_m(h, 1, M);Rp = -min(db(1:omegaC/(pi/M)+1);As =