实验七---FIR数字滤波器的设计

1、Guangxi University of Science and Technology实验报告实验课程： 数字信号处理 实验内容： FIR数字滤波器的设计 院 （系）： 计算机科学与通信工程学院 专 业： 通信工程 班 级： 学生姓名： 学 号： 指导教师： 2016年6月14日实验七 FIR数字滤波器的设计一、 实验目的：1. 掌握用窗函数法设计FIR滤波器的原理及方法，熟悉响应的计算机编程； 2. 熟悉线性相位FIR滤波器的幅频特性和相频特性； 3. 了解各种不同窗函数对滤波器性能的影响。 二、 实验原理： 线性相位实系数FIR滤波器按其N值奇偶和h(n)的奇偶对称性分为四种： 1、h(

2、n)为偶对称，N为奇数H(ej)的幅值关于=0，2成偶对称。 2、h(n)为偶对称，N为偶数H(ej)的幅值关于=成奇对称，不适合作高通。 3、h(n)为奇对称，N为奇数H(ej)的幅值关于=0，2成奇对称，不适合作高通和低通。 4、h(n)为奇对称，N为偶数H(ej) =0、20，不适合作低通。(一) 窗口法 窗函数法设计线性相位FIR滤波器步骤 确定数字滤波器的性能要求：临界频率k，滤波器单位脉冲响应长度N； 根据性能要求，合理选择单位脉冲响应h(n)的奇偶对称性，从而确定理想频率响应Hd(ej)的幅频特性和相频特性； 求理想单位脉冲响应hd(n)，在实际计算中，可对Hd(ej)按M(M远

3、大于N)点等距离采样，并对其求IDFT得hM(n)，用hM(n)代替hd(n)； 选择适当的窗函数w(n)，根据h(n)= hd(n)w(n)求所需设计的FIR滤波器单位脉冲响应； 求H(ej)，分析其幅频特性，若不满足要求，可适当改变窗函数形式或长度N，重复上述设计过程，以得到满意的结果。 窗函数的傅式变换W(ej)的主瓣决定了H(ej)过渡带宽。W(ej)的旁瓣大小和多少决定了H(ej)在通带和阻带范围内波动幅度，常用的几种窗函数有： 矩形窗 w(n)=RN(n)； Hanning窗 ； Hamming窗 ； Blackmen窗 ； Kaiser窗 。 式中Io(x)为零阶贝塞尔函数。三、

4、 实验内容及步骤：给定数字带通滤波器参数如下：请分别选用矩形窗、巴特列特窗、汉宁窗、海明窗、布拉克曼窗及凯泽窗（=7.865）设计该滤波器。（1）画出满足指标要求的六个窗函数对应的滤波器冲激响应序列；程序如下：Ws1=0.2*pi;Wp1=0.5*pi;A1s=60;R1s=1;Ws2=0.8*pi;Wp2=0.6*pi;A2s=80;R2s=0.8;Wc1=(Ws1+Wp1)/2;Wc2=(Ws2+Wp2)/2;deltaW=min(Wp1-Ws1),(Ws2-Wp2); %-取两个过渡带中的小者N0=ceil(2*3.1*pi/deltaW); %-查表7-3（P342）汉宁窗N=N0+m

5、od(N0+1,2); %-确保N为奇数hdWindow=ideallp(Wc2,N)-ideallp(Wc1,N);%-ideallp（）函数（非系统自有函数）在系统安装目录的WORK子目录中，ideallp.m%以下注释的内容为非系统函数： % function hd = ideallp(wc,N); % % hd = ideallp(wc,N) % % - % % hd = 点0 到 N-1之间的理想脉冲响应 % % wc = 截止频率(弧度) % % N = 理想滤波器的长度 % % % tao = (N-1)/2; % 理想脉冲响应的对称中心位置 % n = 0: (N-1); %

6、设定脉冲响应长度 % m = n - tao + eps; % 加一个小数以避免零作除数 % hd = sin(wc*m) ./ (pi*m); % 理想脉冲响应wdWindow1=boxcar(N);%-矩形窗wdWindow2=bartlett(N);wdWindow3=hanning(N);wdWindow4=hamming(N);wdWindow5=blackman(N);wdWindow6=kaiser(N,7.865);hr1=wdWindow1.*hdWindow; %-点乘hr2=wdWindow2.*hdWindow;hr3=wdWindow3.*hdWindow;hr4=w

7、dWindow4.*hdWindow;hr5=wdWindow5.*hdWindow;hr6=wdWindow6.*hdWindow;n=0:N-1; % 矩形窗subplot(4,3,1),stem(n,wdWindow1);xlabel(时间);ylabel(振幅);title(矩形窗);subplot(4,3,4),stem(n,hr1);xlabel(时间);ylabel(振幅);title(矩形窗冲激响应);subplot(4,3,2),stem(n,wdWindow2);% 巴特列特窗xlabel(时间);ylabel(振幅);title(巴特列特窗);subplot(4,3,5)

8、,stem(n,hr2);xlabel(时间);ylabel(振幅);title(巴特列特窗冲激响应);subplot(4,3,3),stem(n,wdWindow3);% 汉宁窗xlabel(时间);ylabel(振幅);title(汉宁窗);subplot(4,3,6),stem(n,hr3);xlabel(时间);ylabel(振幅);title(汉宁窗冲激响应);subplot(4,3,7),stem(n,wdWindow4);% 海明窗xlabel(时间);ylabel(振幅);title(海明窗);subplot(4,3,10),stem(n,hr4);xlabel(时间);yla

9、bel(振幅);title(海明窗冲激响应);subplot(4,3,8),stem(n,wdWindow5);% 布拉克曼窗xlabel(时间);ylabel(振幅);title(布拉克曼窗);subplot(4,3,11),stem(n,hr5);xlabel(时间);ylabel(振幅);title(布拉克曼窗冲激响应);subplot(4,3,9),stem(n,wdWindow6); % 凯泽窗xlabel(时间);ylabel(振幅);title(凯泽窗);subplot(4,3,12),stem(n,hr6);xlabel(时间);ylabel(振幅);title(凯泽窗冲激响应

10、);截图如下：图：7-1（2）在同一坐标图中绘出对应六个窗函数设计出来的滤波器幅频特性曲线（绘制绝对振幅图和相对振幅（分贝）图各一幅），观察幅频特性的变化，逐一验证是否满足设计要求，通过比较说明分贝图相比普通的绝对振幅图在观察滤波器阻、通带衰减或损耗数值上的优势；程序如下：Ws1=0.2*pi;Wp1=0.5*pi;A1s=60;R1s=1;Ws2=0.8*pi;Wp2=0.6*pi;A2s=80;R2s=0.8;Wc1=(Ws1+Wp1)/2;Wc2=(Ws2+Wp2)/2;deltaW=min(Wp1-Ws1),(Ws2-Wp2); %-取两个过渡带中的小者N0=ceil(2*3.1*pi

11、/deltaW); %-查表7-3（P342）汉宁窗N=N0+mod(N0+1,2); %-确保N为奇数hdWindow=ideallp(Wc2,N)-ideallp(Wc1,N);wdWindow1=boxcar(N);%-矩形窗wdWindow2=bartlett(N);wdWindow3=hanning(N);wdWindow4=hamming(N);wdWindow5=blackman(N);wdWindow6=kaiser(N,7.865);hr1=wdWindow1.*hdWindow;%-点乘hr2=wdWindow2.*hdWindow;hr3=wdWindow3.*hdWin

12、dow;hr4=wdWindow4.*hdWindow;hr5=wdWindow5.*hdWindow;hr6=wdWindow6.*hdWindow;H1,W=freqz(hr1,1);H2,W=freqz(hr2,1);H3,W=freqz(hr3,1);H4,W=freqz(hr4,1);H5,W=freqz(hr5,1);H6,W=freqz(hr6,1);subplot(1,2,1),plot(W/pi,abs(H1),W/pi,abs(H2),:,W/pi,abs(H3),-.,W/pi,abs(H4),W/pi,abs(H5),W/pi,abs(H6),-.);xlabel(om

13、ega/pi);ylabel(绝对振幅);title(带通滤波器幅频特性（1）);legend(矩形窗,巴特列特窗,汉宁窗,海明窗,布拉克曼窗,凯泽窗);subplot(1,2,2),plot(W/pi,20*log10(abs(H1),W/pi,20*log10(abs(H2),:,W/pi,20*log10(abs(H3),-.,W/pi,20*log10(abs(H4),W/pi,20*log10(abs(H5),W/pi,20*log10(abs(H6),-.);xlabel(omega/pi);ylabel(相对振幅);title(带通滤波器幅频特性（2）);legend(矩形窗,巴

14、特列特窗,汉宁窗,海明窗,布拉克曼窗,凯泽窗);for step=1:1:512end截图如下：图：7-2（3）以凯泽窗为例分析：在已经满足设计指标要求的情况下，当N取值增大时，滤波器幅频特性的相应变化；又，固定N值，当在4到9间变化时，滤波器幅频特性的相应变化；程序如下：Ws1=0.2*pi;Wp1=0.5*pi;A1s=60;R1s=1;Ws2=0.8*pi;Wp2=0.6*pi;A2s=80;R2s=0.8;Wc1=(Ws1+Wp1)/2;Wc2=(Ws2+Wp2)/2;deltaW=min(Wp1-Ws1),(Ws2-Wp2); %-取两个过渡带中的小者N0=ceil(2*3.1*pi

15、/deltaW); %-查表7-3（P342）汉宁窗N=N0+mod(N0+1,2); %-确保N为奇数hdWindow5=ideallp(Wc2,N)-ideallp(Wc1,N);hdWindow6=ideallp(Wc2,N+18)-ideallp(Wc1,N+18);%-ideallp（）函数（非系统自有函数）在系统安装目录的WORK子目录中，ideallp.mwdWindow5=kaiser(N,5.658);wdWindow6=kaiser(N,8.960);hr5=wdWindow5.*hdWindow;hr6=wdWindow6.*hdWindow;H5,W=freqz(hr5

16、,1);H6,W=freqz(hr6,1);subplot(1,2,1),plot(W/pi,abs(H5),W/pi,abs(H6),:);xlabel(omega/pi);ylabel(绝对振幅);title(带通滤波器幅频特性（1）);legend(=5.658,=8.960);subplot(1,2,2),plot(W/pi,20*log10(abs(H5),W/pi,20*log10(abs(H6),:);xlabel(omega/pi);ylabel(相对振幅);title(带通滤波器幅频特性（2）);legend(=5.658,=8.960);截图如下;图：7-3（4）若阻带最小

17、衰减要求均提高到60dB,分析要用什么窗函数才能实现。答：如果阻带衰减60dB,则需要用布莱克曼窗，因为布莱克曼窗最小衰减为74dB，只有此窗符合技术要求。海明窗最小衰减小于60dB，故不选此窗。实验总流程图：图：7-4四、实验总结与分析：通过实验，认识了不同窗函数的特性曲线，基本掌握了用窗函数法设计FIR滤波器的原理及方法，也熟悉了线性相位FIR滤波器的幅频特性和相频特性。进一步知道了各种窗函数的区别。以及相应函数的用法：（1）ceil:向正无穷方向取最近整数，即向数轴右边取最近整数；（2）freqz:H,F = freqz(B,A,N,Fs) 其中B/A 提供滤波器系数，N 表示选取单位圆的上半圆等间距的N个点作为频响输出；Fs 为采样频率，该参数可以省略；H 为N个点处的频率响应复值输出向量，其模即为频响幅值曲线幅值20log10(abs(H)DB，其幅角angle(H)即为频响相位曲线相位值。F 为与第N点处对应的频率值f（Hz），如果Fs 参数省略时，则频率值w为rad/sample，w = 2*pi*f/Fs。（3）legend:添加图例的标注。（4）mod:是一个求余函数，其格式为：mod(nExp1,nExp2)，即是两个数值表达式作除法运算后的余数。（5）矩形窗:调用格式：w=boxcar(