第五次DSP实验

1、第五次DSP实验FIR滤波一、 实验目的1、熟悉 FIR 滤波器设计的基本方法。2、掌握用窗函数设计 FIR 数字滤波器的原理及方法，熟悉相应的计算机高级语言编程。3、熟悉线性相位 FIR 滤波器的幅频特性和相位特性。4、了解各种不同窗函数对滤波器性能的响应。二、实验原理与方法1.FIR 滤波器的设计在前面的实验中， 我们介绍了 IIR 滤波器的设计方法并实践了其中的双线性变换法， IIR具有许多诱人的特性； 但与此同时， 也具有一些缺点。 例如： 若想利用快速傅立叶变换技术进行快速卷积实现滤波器，则要求单位脉冲响应是有限列长的。此外， IIR 滤波器的优异幅度响应， 一般是以相位的非线性为代

2、价的， 非线性相位会引起频率色散。 FIR 滤波器具有严格的相位特性， 这对于语音信号处理和数据传输是和重要的。 目 前 FIR 滤波器的设计方法主要有三种： 窗函数法、 频率取样法和切比雪夫等波纹逼近的最优化设计方法。 常用的是窗函数法和切比雪夫等波纹逼近的最优化设计方法。 本实验项目 中的窗函数法比较简单， 可应用现成的窗函数公式， 在技术指标要求不高的时候是比较灵活方便的。 它是从时域出发， 用一个窗函数截取理想的 hd(n) 得到h(n)，以有限长序列 h(n)近似理想的hd(n) ；如果从频域出发， 用理想的h(ejw) 在单位圆上等角度取样得到 H(k)，根据 h(k)得到 H(z

3、)将逼近理想的 ，这就是频率取样法。2.窗函数设计法同其他的数字滤波器的设计方法一样， 用窗函数设计滤波器也是首先要对滤波提出性能指标。一般是给定一个理想的频率响应 ， 使所设计的 FIR 滤波器的频率响应去逼近所要求的理性的滤波器的响应 。 窗函数法设计的任务在于寻找一个可实现（有限长单位脉冲响应）的传递函数 三、实验内容1. 傅立叶反变换数值计算子程序：代码如下：function xa = myifft(x,m,n) na = 0:n-1; xa=x * ( exp( 1i * pi * 2 /length(m) ) ) . (m' * na)/length(m);end2.窗函数

4、产生子程序function r = window(n) r = ones(1,n-1);endfunction r = hamming(n) r = ones(1,n-1); nr = 0:n-2; r = r.*(0.54-0.46*cos(2*pi*nr/(n-1); %stem(nr,r);endfunction r = hanning(n) r = ones(1,n-1); nr = 0:n-2; r = r.*(0.5-0.5*cos(2*pi*nr/(n-1); %stem(nr,r);endfunction r = blackman(n) r = ones(1,n-1); nr

5、= 0:n-2; r = r.*(0.42-0.5*cos(2*pi*nr/(n-1)+0.08*cos(4*pi*nr/(n-1); %stem(nr,r);endfunction r = kaiser(n,b) r = ones(1,n-1); nr = 0:n-2; tr = b*sqrt(1-(2*nr/(n-1)-1).2); r = r.*besseli(0,tr)/besseli(0,b); %stem(nr,r);end3.主函数function xn = main(hd,n,m) nhd = 0:m-1; hd1 = zeros(1,m); hd1(1:m/2) = hd(1

6、:m/2); hd1(end:-1:m/2+2) = hd1(2:m/2); if mod(n,2) = 0 hd1(1:m/2) = hd1(1:m/2) .* exp(-1i*pi*nhd(1:m/2)/m*(n-1); hd1(m/2+1) = 0; hd1(m/2+2:m) = hd1(m/2+2:m) .* (-exp(-1i*pi*nhd(m/2+2:m)/m*(n-1); xn = real(myifft(hd1,nhd,n); else hd1 = hd1 .* exp(-1i*pi*nhd/m*(n-1); xn = real(myifft(hd1,nhd,n); enden

7、d3上机部分1. 用窗函数法设计一个长度 N 等于 8 的线性相位 FIR 滤波器。其理想的幅频特性为 分别用矩形窗、Hanning 窗、Hamming 窗、Blackman 窗、Kaiser 窗（ =8.5）设计该滤波器，比较设计结果代码如下：function test() n =8; m = 10*n; nhd = 0:m-1; hd = zeros(1,m); hd(1:floor(m*0.2)+1) = 1; subplot(2, 2, 1); stem(nhd(1:m/2)/m*2,abs(hd(1:m/2); xlabel( ' frequency in pi units&

8、#39; ); title( '滤波器频域设计' ); ylabel( 'Hd' ); xn = main(hd,n,m); nx = 0:n-2; xn = xn(1:n-1).*window(n); subplot(2, 2, 2); stem(nx,xn); xlabel( ' n' ); title( '滤波器时域设计' ); ylabel( 'x' ); k = 0:1:1000; %置初值 X, w = dtft ( xn, nx, k ); %调用dtft函数，完成变换 magX = abs( X )

9、; angX = angle( X ); %取模和相位 subplot( 2, 2, 3 ); plot( w(1:501) / pi, magX(1:501) ); grid %分别作出幅度，相位，实部，虚部对的归一化w的图 %这里只做出w取值为0,pi的图像，并对w做归一化 xlabel( ' frequency in pi units' ); title( 'Magnitude part' ); ylabel( 'Magnitude' ); subplot( 2, 2, 4 ); plot( w(1:501) / pi, angX(1:50

10、1) ); grid xlabel( ' frequency in pi units' ); title( 'Angle part' ); ylabel( 'Radians' ); end实验结果：窗函数类型h(0)h(1)h(2)h(3)h(4)h(5)h(6)矩形窗-0.0898 -0.01250.19770.38430.38430.1977-0.0125Hanning窗0 -0.00240.12090.36520.36520.1209-0.0024Hamming窗-0.0072 -0.00320.12700.36680.36680.1270

11、-0.0032Blackman窗0.0000 -0.00110.09080.35370.35370.0908-0.0011Kaiser(=8.5)-0.0001 -0.00120.09180.35410.35410.0918-0.0012可以看到，实验结果和参考答案是基本一样的，所以认为实验程序调试成功2.用 Hanning 窗设计一个线性相位带通滤波器，其长度 N=15上下边带截至频率分别 ， ，求 h(n)，绘制它的幅频和相位特性曲线，观察它的实际3dB和20dB带宽。如果N=45，重复这一设计，观察幅频和相位特性的变化，注意长度 N 变化对结果的影响。代码如下：function ex1(

12、) n =15; m = 20*n; nhd = 0:m-1; hd = zeros(1,m); hd(floor(m*0.15)+1:floor(m*0.25)+1) = 1; subplot(2, 2, 1); stem(nhd(1:m/2)/m*2,abs(hd(1:m/2); xlabel( ' frequency in pi units' ); title( '滤波器频域设计' ); ylabel( 'Hd' ); xn = main(hd,n,m); nx = 0:n-2; xn = xn(1:n-1).*hanning(n); su

13、bplot(2, 2, 2); stem(nx,xn); xlabel( ' n' ); title( '滤波器时域设计' ); ylabel( 'x' ); k = 0:1:1000; %置初值 X, w = dtft ( xn, nx, k ); %调用dtft函数，完成变换 magX = abs( X ); angX = angle( X ); %取模和相位 subplot( 2, 2, 3 ); plot( w(1:501) / pi, magX(1:501) ); grid %分别作出幅度，相位，实部，虚部对的归一化w的图 %这里只做出

14、w取值为0,pi的图像，并对w做归一化 xlabel( ' frequency in pi units' ); title( 'Magnitude part' ); ylabel( 'Magnitude' ); subplot( 2, 2, 4 ); plot( w(1:501) / pi, angX(1:501) ); grid xlabel( ' frequency in pi units' ); title( 'Angle part' ); ylabel( 'Radians' ); end运行

15、结果如下： 这里n=15，如图可以得到一个带通的滤波器，可以测得其通带为0.288pi0.512pi，这要大于0.3pi0.5pi的设计要求，而且可以看到，其过渡带较宽，通带增益衰减快，因此性能比较一般。相位有时会有+-pi的跳变。 这里n=45，如图可以得到一个带通的滤波器，可以测得其通带为0.318pi0.48pi，这要小于0.3pi0.5pi的设计要求，但是其过渡带较窄，通带增益平稳，因此性能比较好。这是由于函数阶次高，长度长，频谱泄露小，但是可以看到为了达到要求的带宽，设计师应放松对截止频率的要求。相位会有+-pi的跳变。而去变化频繁。3.改用矩形窗和 Blackman窗，设计步骤（1

16、）中的带通滤波器，观察并记录窗函数对滤波器幅频和相位特性的影响，比较这三种窗函数的特点实验结果如下：矩形窗，n=15矩形窗，n=45Blackman窗，n=15Blackman窗，n=45综合比较这几个函数，可以发现，随着n的增大，将导致系统函数幅频响应变窄，过渡带变陡，阻带变平，更趋于0，而通带性状变好，变平，更趋于1，通带变窄，已经不能满足设计要求。而相频响应的变化将变得更加剧烈。对比两个不同的窗函数，可以发现blackman窗产生的系统函数具有更好的性能。而矩形窗函数产生的系统函数的幅频响应在通带和阻带内均有波纹，不管n=15还是n=45，通带均过窄，不满足设计要求。相频响应也不是线性的

17、。3.用 Kaiser 窗设计一个专用的线性相位滤波器。N=40，理想的幅频特性。当 值分别 4,6,8 时，设计相应的滤波器，比较它们的幅频和相位特性，观察并分析值不同的时候对结果有什么影响。代码如下：function ex3() n =40; m = 20*n; nhd = 0:m-1; hd = zeros(1,m); hd(floor(m*0.1)+1:floor(m*0.2)+1) = 1; hd(floor(m*0.3)+1:floor(m*0.4)+1) = 1; subplot(2, 2, 1); stem(nhd(1:m/2)/m*2,abs(hd(1:m/2); xlabe

18、l( ' frequency in pi units' ); title( '滤波器频域设计' ); ylabel( 'Hd' ); xn = main(hd,n,m); nx = 0:n-2; xn = xn(1:n-1).*kaiser(n,4); subplot(2, 2, 2); stem(nx,xn); xlabel( ' n' ); title( '滤波器时域设计' ); ylabel( 'x' ); k = 0:1:1000; %置初值 X, w = dtft ( xn, nx, k

19、); %调用dtft函数，完成变换 magX = abs( X ); angX = angle( X ); %取模和相位 subplot( 2, 2, 3 ); plot( w(1:501) / pi, magX(1:501) ); grid %分别作出幅度，相位，实部，虚部对的归一化w的图 %这里只做出w取值为0,pi的图像，并对w做归一化 xlabel( ' frequency in pi units' ); title( 'Magnitude part' ); ylabel( 'Magnitude' ); subplot( 2, 2, 4

20、); plot( w(1:501) / pi, angX(1:501) ); grid xlabel( ' frequency in pi units' ); title( 'Angle part' ); ylabel( 'Radians' ); end结果分析如下：Kaiser窗，b=4Kaiser窗，b=6Kaiser窗，b=8Kaiser窗，b=4，6，8对比可以发现，随着b的增大，过渡带的范围也将增大，导致阻带和通带的性能变差。综合来看，b=4时有最好的性能。，但其实它们都没办法保证应有的带宽。思考题1. 对不同的窗函数和不同的n，3db截止频率都有不同，根据上面的实验，可以发现，随着n的增大，将导致3db截止频率向通带中心频率靠近，也就是说通带变窄。它显然不等于理想频率响应的截止频率。2.由1中的讨论可知，n的大小将直接影响设计通带边缘截至频率和阻带临界频率，所以可以通过选定合适的n已经窗函数以达到设计要求。但这些关系实