dsp课程设计报告(基于窗函数 设计法线性相位型FIR)

1、数 字 信 号 处 理课 程 设 计 名 称 基于窗函数设计法线性相位型 FIR数字低通滤波器设计 姓 名 班 级 学 号 指 导 教 师 目录 一 内容概述 3二 题目要求细则 4三 设计原理 5四 设计过程 6五 实验结果及分析 7六 实验心得 14七 参考文献 15八 源代码 16数字信号处理课程设计任务书题目7基于窗函数设计法线性相位型FIR数字低通滤波器设计主要内容1、根据ALPF指标，采用三种窗函数，设计线性相位型FIR滤波器；2、生成一个时间序列，对之滤波，并分析三种滤波器的滤波效果。设计要求1、根据ALPF指标，生成物理可实现的线性相位FIR滤波器之冲击响应函数；2、采用三种窗

2、函数，设计线性相位型FIR滤波器；3、生成一个时间序列，并利用已设计的线性相位型FIR滤波器对之滤波；4、对比分析三种窗函数法设计的数字滤波器性能。主要仪器设备1、计算机1台，安装MATLAB软件主要参考文献美维纳.K.恩格尔，约翰.G.普罗科斯著，刘树棠译.数字信号处理使用MATLABM.西安：西安交通大学出版社，2002. 飞思科技产品研发中心编著.MATLAB7辅助信号处理技术与应用M.北京：电子工业出版社，2005.课程设计进度安排（起止时间、工作内容）课程设计共设16个设计题目，每班3至4人为1组，1人1套设备，每组选作不同的题目，4个班共分4批。完整课程设计共20学时，为期1周，具

3、体进度如下：5学时 学习题目相关知识，掌握实现原理；5学时 用MATLAB语言实现题目要求；5学时 进一步完善功能，现场检查、答辩；5学时 完成并提交课程设计报告。课程设计开始日期课程设计完成日期课程设计实验室名称地 点资料下载地址一【内容概述】 FIR数字滤波器设计最简单的方法是窗函数法，通常也称之为傅里叶级数法。它是在属于进行的，因而必须由理想滤波器的频率响应推导出其单位冲击响应，再设计一个FIR数字滤波器的单位冲击响应去逼近。 窗函数设计法的基本原理是用一定宽度窗函数截取无限脉冲响应序列获得有限长的脉冲响应序列,先通过傅里叶逆变换获得理想滤波器的单位脉冲响应hd(n)。再由性能指标确定窗

4、函数W(n)和窗口长度N。最后求得实际滤波器的单位脉冲响应h(n)。 窗函数在设计FIR数字滤波器中有很重要的作用，正确的选择窗函数可以提高所涉及的数字滤波器的性能，或者在满足设计要求的情况下，减小FIR数字滤波器的阶数。 二【题目要求细则】基于窗函数设计法线性相位型FIR数字低通滤波器设计目的： 熟练掌握MATLAB工具软件在工程设计中的使用； 熟练掌握窗函数法设计线性相位型FIR滤波器的方法。要求： 根据给定ALPF幅频特性指标，生成与之相对应的模拟滤波器h(t)； 根据h(t)取得与之想对应的数字滤波器冲击响应函数h(n)； 选择适当的窗函数，构建线性相位型FIR的冲击响应函数hd(n)

5、； 设计与之相对应的DLPF，对比分析DLPF幅频特性是否符合要求； 产生一有干扰频率的时域序列（借助FFT分析说明其有干扰），使之通过所设计的DLPF，对滤波输出结果作出分析，说明输出结果。步骤： 根据给定设计指标（通带截止频率fb1 000 Hz，通带最大衰减p=1dB，阻带截止频率fc1 200 Hz，阻带最小衰减s=100dB，令T=2）生成与之相对应的模拟滤波器h(t)； 根据h(t)取得与之想对应的数字滤波器冲击响应函数h(n)； 选择适当的窗函数，构建线性相位型FIR的冲击响应函数hd(n)； 设计与之相对应的DLPF，给出窗函数及所设计滤波器的幅度特性，并对比理想冲激响应和实际

6、冲击响应，对比分析DLPF幅频特性是否符合要求； 试说明不同窗函数对所设计数字滤波器性能的影响； 产生一有干扰频率的时域序列（借助FFT分析说明其有干扰），使之通过所设计的DLPF，对滤波输出结果作出分析，说明输出结果。三【设计原理】 窗函数设计法的基本原理是用一定宽度窗函数截取无限脉冲响应序列获得有限长的脉冲响应序列,先通过傅里叶逆变换获得理想滤波器的单位脉冲响应hd(n)。再由性能指标确定窗函数W(n)和窗口长度N。最后求得实际滤波器的单位脉冲响应h(n)。 FIR数字滤波器设计最简单的方法是窗函数法，通常也称之为傅里叶级数法。它是在属于进行的，因而必须由理想滤波器的频率响应推导出其单位冲

7、击响应，再设计一个FIR数字滤波器的单位冲击响应去逼近。 假设理想滤波器的频率响应的表示为：由于滤波器的频率响应与单位冲击响应序列是傅里叶变换对，据此：求得理想滤波器的单位冲击响应后，可得到该离散滤波器的系统传递函数：注意，该为无限长序列，因此是物理不可实现的。为了使系统变为物理可实现的，且使实际的FIR滤波器频率响应尽可能逼急理想滤波器的频率响应，采用窗函数将无线脉冲响应截取一段好来近似表示，可得：，从而有。式中N表示窗口长度，这样就是物理可实现的系统。四【设计过程】先产生一个理想低通，分别用三种不同的窗函数与之时域内点乘后，得到可物理实现的低通滤波器h(n)，观察它的幅度相频特性。产生一个

8、干扰序列，加上噪声后得到一个有毛刺的正弦波，通过三种加不同窗的滤波器，得出三个虑完波的时域序列和频谱。比较三个虑后的时域频域波形，比较出三种不同滤波器的优劣。五【实验结果及分析】1.设计参数：通带截止频率fb1 000 Hz，通带最大衰减p=1dB，阻带截止频率fc1 200 Hz，阻带最小衰减s=100dB，令T=22.不同窗形状对比：Blackman窗、汉宁窗、海明窗、凯塞窗： 分析：若要根据指标设计出理想的数字滤波器，只有布莱克曼窗和凯塞窗可以接近设计要求由于除凯塞窗之外的窗函数全部是参数固定窗函数，所以若要想达到给定设计指标阻带最小衰减s=100dB，只有选择凯塞窗，估算值，才有可能达

9、到要求。由图可知，凯塞窗和布莱克曼窗比较接近实验要求，汉宁窗和海明窗较差，所以matlab默认窗函数海明窗在此设计指标不太合适。3.理想单位脉冲响应对比：布莱克曼： 汉宁：海明： 凯塞：4.加窗后截取实际单位脉冲响应对比：布莱克曼： 汉宁：海明： 凯塞：分析：用不同的窗函数去截取序列会产生不同的误差，也就是通常说的吉布斯效应，会引起过渡带加宽以及通带和阻带内的波动，尤其使阻带衰减减小。由图可知，布莱克曼窗和凯塞窗图中最大正、负峰值对应频率间距较小，而海明窗和汉宁窗较大，所以选择后两者会增大截断效应5.不同窗函数幅度响应对比：布莱克曼： 汉宁：海宁：凯塞窗：分析：要想获得理想的阻带衰减，需要使用

10、窗谱最大旁瓣相对幅度较小的窗函数。由图可知，汉宁窗的阻带波纹最大，所以其阻带衰减就最小。而布莱克曼窗和凯塞窗阻带波纹相对较小，所以其阻带衰减较明显。而凯塞窗可以通过调整窗长度以及值可以有效控制过渡带宽度，获得良好阻带衰减。所以选择不同的窗函数可以减少带内波动以增大阻带衰减。6.滤波器性能分析：（ap:通带最大衰减；as:阻带最小衰减）Blackman窗：Ap1 = 0.0027As1 =73汉宁窗： Ap2 =0.0711As2 =44海明窗： Ap3 =0.0282As3 =51 凯塞窗： Ap4 =0.0056As4 =108例如：布莱克曼窗指数：通过分别对系统加Blackman窗、汉宁窗

11、、海明窗、凯塞窗得到以上结果。从加窗后系统的幅频响应曲线可以看出，Blackman窗凯塞窗引起的过渡带加宽现象和阻带波动较小，即截断效应较小，而且通过比较可以看出，凯塞窗还具有最小的通带最大衰减和阻带最小衰减，综合比较发现凯塞窗blackman窗的各项性能优于海明窗和汉宁窗，其中凯塞窗性能最好。所以，选择凯塞窗选择不同的值，可以设计出性能最优的FIR滤波器。7产生一个有干扰频率的时域序列，使之通过所设计的DLPF，对滤波输出结果作出分析，说明输出结果。布莱克曼：汉宁窗：海明窗：凯塞窗：分析：由图明显看出时域和频域部分分别经过低通滤波器后高频部分被滤掉，时域部分毛刺被滤掉。8扩展部分：第一种：第

12、二种：分析：由图可知，通过带通滤波器后低频和高频被滤掉七【参考文献】参考文献 1、 高西全，数字信号处理（第三版）西安电子科技大学出版社2、 程佩青，数字信号处理教程（第三版），清华大学出版社3、 MATLAB7.0在数字信号处理中的应用罗军辉机械工业出版社4、 飞思科技产品研发中心编著.MATLAB7辅助信号处理技术与应用M.北京。 电子工业出版社，2005. 5、数字信号处理原理及其MATLAB实现从玉良电子工业出版社 八【实验部分代码】% % 基于窗函数设计法线性相位型FIR数字低通滤波器设计% 布莱克曼 %clc,clearfp=1000; fc=1200;fs=4000; %对信号的

13、采样频率为4000ws=2*pi*fc/fs %t通带截至频率wp=2*pi*fp/fs %阻带截至频率tr_width=ws-wp; %求过渡带M=ceil(11*pi/tr_width)+1 %利用blackman窗的过渡带=(11pi/M),求窗宽Mn=0:1:M-1;wc=(ws+wp)/2; %求通带截止频率wcalpha=(M-1)/2;n=0:(M-1);m=n-alpha;h1=sin(wc*m)./(pi*m); %求线性低通滤波器的单位脉冲响应hd=ideal_lp(wc,M); %求系统单位脉冲响应hd(n)w_black=(blackman(M); %根据阶数M求bla

14、ckman窗h=hd.*w_black; %对FIR滤波器低通冲激响应加窗，h为低通滤波器的冲击相应db,mag,pha,w=freqz_m2(h,1) % 计算实际滤波器幅度响应delta_w=2*pi/1000;Ap1=-(min(db(1:1:wp/delta_w+1) %实际通带纹波As1=-round(max(db(ws/delta_w+1:1:501),% 汉宁窗 %clc,clearfb=1000;fc=1200;fs=4000; %对信号的采样频率为4000ws=2*pi*fc/fs; %t通带截至频率wp=2*pi*fb/fs %阻带截至频率tr_width=ws-wp; %

15、求过渡带M=ceil(6.2*pi/tr_width)+1 %利用Hanning窗的过渡带=(6.2pi/M),求窗宽Mn=0:1:M-1;wc=(ws+wp)/2; %求通带截止频率wchd=ideal_lp(wc,M); %求单位脉冲响应hd(n)w_han=(hanning(M);h=hd.*w_han; %求FIR滤波器低通冲激响应h(n)db,mag,pha,w=freqz_m2(h,1); % 计算实际滤波器幅度响应delta_w=2*pi/1000;Ap2=-(min(db(1:1:wp/delta_w+1) %实际通带纹波As2=-round(max(db(ws/delta_w

16、+1:1:501) % 海明窗 %clc,clearfp=1000;fc=1200;fs=4000; %对信号的采样频率为4000ws=2*pi*fc/fs; %t通带截至频率wp=2*pi*fp/fs %阻带截至频率tr_width=ws-wp; %求过渡带M=ceil(6.6*pi/tr_width)+1 %利用Hamming窗的过渡带=(6.6pi/M),求窗宽Mn=0:1:M-1;wc=(ws+wp)/2; %求通带截止频率wchd=ideal_lp(wc,M); %求单位脉冲响应hd(n)w_ham=(hamming(M); %海明窗h=hd.*w_ham; %求FIR滤波器低通冲激

17、响应h(n) 实际单位脉冲响应db,mag,pha,w=freqz_m2(h,1); % 计算实际滤波器幅度响应delta_w=2*pi/1000;Ap3=-(min(db(1:1:wp/delta_w+1) %实际通带纹波As3=-round(max(db(ws/delta_w+1:1:501)% 凯塞窗 %clc,clearfp=1000; fc=1200;fs=4000; %对信号的采样频率为4000ws=2*pi*fc/fs %t通带截至频率wp=2*pi*fp/fs %阻带截至频率tr_width=ws-wp %求过渡带M=ceil(11*pi/tr_width)+1 %利用blac

18、kman窗的过渡带=(11pi/M),求窗宽Mn=0:1:M-1;wc=(ws+wp)/2; %求通带截止频率wchd=ideal_lp(wc,M); %求系统单位脉冲响应hd(n)w_kai=(kaiser(M,10.2256); %根据阶数M求blackman窗 理想的(100-8.7)*0.112 h=hd.*w_kai; %对FIR滤波器低通冲激响应加窗，h为低通滤波器的冲击相应db,mag,pha,w=freqz_m2(h,1); % 计算实际滤波器幅度响应delta_w=2*pi/1000;Ap4=-(min(db(1:1:wp/delta_w+1) %实际通带纹波As4=-rou

19、nd(max(db(ws/delta_w+1:1:501)a=7.9856;L=(100-8)/(2.285*tr_width) %产生干扰序列%产生随机白噪声f1=600;f2=1800;W1=f1*(2*pi);W2=f2*(2*pi);t=0:(1/4000):0.05;m=sin(W1.*t)+cos(W2.*t);N=10000;M=abs(fft(m,N);k=filter(h,1,m);%k=conv(m,h2);K=abs(fft(k,N);figure();subplot(2,2,1),plot(m);ylabel(m);title(输入信号);subplot(2,2,2),

20、stem(M(1:ceil(N/2);ylabel(fft(m);title(输入信号幅频响应);subplot(2,2,3),plot(k);ylabel(k);title(凯塞窗 通过低通滤波器);subplot(2,2,4),stem(K(1:ceil(N/2);ylabel(fft(k);title(通过低通滤波器后幅频响应); 扩展：1. clf,clc,clearwph=0.65*pi;wpl=0.35*pi;wsth=0.83*pi;wstl=0.17*pi;tr_width=wsth-wph;Mham=ceil(6.6*pi/tr_width)+1 %利用Hamming窗的过渡

21、带=(6.6pi/M),求窗宽Mn=0:1:Mham-1;hd=ideal_bp(wpl,wstl,wph,wsth,Mham); %求单位脉冲响应hd(n)w_ham=(hamming(Mham);h=hd.*w_ham; %求FIR滤波器低通冲激响应h(n)H,w=freqz(h,1,1000,whole);db=20*log10(abs(H)/max(abs(H); %求FIR滤波器频响的dB值delta_w=2*pi/1000;rp=-min(db(1:1:(wpl/delta_w+1) %求通带最大衰减as=-max(db(wstl/delta_w+1:1:501) %求阻带最小衰减Horg=fft(hd,1000);db_org=20*log10(abs(Horg)/max(abs(Horg); %求FIR滤波器频响的dB值2. clear all;wpl=0.4*pi;wph=0.6*pi;wsl=0.2*pi;wsh=0.8*pi;tr_width=min(wpl-wsl),(wsh-wph);N=cei