《基于频域抽样法的fir数字带阻滤波器.docx》由会员分享，可在线阅读，更多相关《基于频域抽样法的fir数字带阻滤波器.docx（16页珍藏版）》请在三一办公上搜索。
1、武汉理工大学数字信号处理课程设计说明书目 录1前言12 设计原理12.1线性相位FIR数字滤波器的特点12.2 频率抽样法设计原理12.3线性相位FIR数字滤波器的约束条件33 带阻滤波器设计过程33.1 带阻滤波器概念33.2 初始条件解析43.3 设计步骤43.4改善滤波器性能的措施54 MATLAB程序仿真55 Simulink仿真106 设计心得13参考文献：131前言 数字滤波器是一种用来过滤时间离散信号的数字系统，通过对抽样数据进行数学处理来达到频域滤波的目的。根据其单位冲激响应函数的时域特性可分为两类：无限冲激响应（IIR）滤波器和有限冲激响应（FIR）滤波器。与IIR滤波器相比
2、，FIR的实现是非递归的，总是稳定的；更重要的是，FIR滤波器在满足幅频响应要求的同时，可以获得严格的线性相位特性。因此，它在高保真的信号处理，如数字音频、图像处理、数据传输、生物医学等领域得到广泛应用。 有限长单位冲激响应(FIR) 数字滤波器具有严格的线性相位，又具有任意的幅频特性。同时FIR 系统只有零点，系统是稳定的,因而容易实现线性相位和允许实现多通道滤波器。只要经过一定的时延，任何非因果有限长序列都能变成因果的有限长序列， 因而总能用因果系统来实现。FIR 滤波器由于单位冲激响应是有限长的，可以用快速傅立叶变换(FFT) 算法来实现过滤信号，从而大大提高运算效率。由于FIR 滤波器
3、具有以上优点，在信号处理和数据传输中得到了广泛的应用。 Matlab 语言是一种用于科学计算的高效率语言。随着Matlab信号处理工具箱(Signal Pro2cessing Toolbox) 的不断完善,使数字滤波器的计算机辅助设计得以实现。2 设计原理2.1线性相位FIR数字滤波器的特点FIR数字滤波器的系统函数无分母，为，系统频率响应可写成：，令，H(w)称为幅度函数，称为相位函数。这与模和幅角的表示法有所不同，H(w)为可正可负的实数，这是为了表达上的方便。如某系统频率响应，如果采用模和幅角的表示法，的变号相当于在相位上加上，从而造成相位曲线的不连贯和表达不方便，而用这种方式则连贯而方
4、便。 线性相位的FIR滤波器是指其相位函数满足线性方程： （是常数）根据群时延的定义，式中表示系统群时延，表示附加相移。线性相位的FIR系统都具有恒群时延特性，因为为常数，但只有0的FIR系统采具有恒相时延特性。 FIR滤波器的极点都在原点上，而h(n)是因果稳定的有限长序列，因此H(z)在有限z平面上是稳定的。线性相位FIR DF的零点有自己的特点：它们必定是互为倒数的共轭对。2.2 频率抽样法设计原理 周期系列的离散傅里叶级数的系数的值和的一个周期的z变换在单位圆的N个均分点上的抽样值相等，这就实现了频域的抽样。通过有限的频率特性取样值去逼近理想滤波特性，然后由有限的频率特性取样值（如系统
5、冲击响应的DFT）取得系统函数 。窗函数设计法是从时域出发，把理想的hd(n)用一定形状的窗函数截取成有限长的h(n) ，以h(n)来近似hd(n)，从而使频率相应函数近似理想频率响应。频率取样法是从频域出发，对理想的频率响应进行等间隔取样，以有限个频率相应去近似理想频率响应即： （式2-1）等间隔取样，并且 相比较而言，窗函数法思路为：；频率抽样法思路为：。对于离散值有插值公式： （式2-2）应用于系统函数得： （式2-3）当采样点数N已知后，便是常数，只要采样值确定，则系统函数就可确定，要求的FIR滤波器就设计出来了。上式形式的FIR滤波器很容易以频率采样型结构实现。对滤波系统的频率特性有
6、： （式2-4）代入得： （式2-5）即其中为频率取样内插函数。重构的频率响应，在采样点上严格等于H(k)，而在采样点之间，频率相应有加权的内插函数延伸叠加而成。2.3线性相位FIR数字滤波器的约束条件若，其中、分别是对幅度函数和相位函数的第个抽样。因为是实数，所以一定满足共轭偶对称式（3-59）： （式2-6）又因为线性性，满足对称性，所以对一般滤波用的第1、2类FIR滤波器，必须满足条件： （式2-7）对于正交网络、微分器（第3、4类FIR滤波器，必须满足条件： （式2-8）综合以上条件，只有当满足式（式2-6），满足式（式2-7）、（式2-8）之一时，才有线性相位。如果理想频响给得不好或
7、采样点位置安排得不恰当，都将得不到线性相位。只有当满足上面的约束条件时，对0，区间上抽取一半频率样点，而其余的一半根据约束条件强行推出。3 带阻滤波器设计过程3.1 带阻滤波器概念 数字带阻滤波器也具有频率响应的周期性，频率变量以数字频率来表示（，为模拟角频率，为抽样时间间隔，为抽样频率），所以数字滤波器设计中必须给出抽样频率。图3.1数字带阻滤波器理想幅度频率响应（只表示了正频率部分），这样的理想频率响应是不可能实现的，原因是频带之间幅度响应是突变的，因而其单位抽样响应是非因果的。图3.1 理想带阻滤波器的幅频特性一般来说，滤波器的性能要求往往以频率响应的幅度特性的允许误差来表征。以低通滤波
8、器为例，如图3.2所示，频率响应有通带、过渡带和带阻三个范围（非理想的）。在通带内，幅度响应以误差为逼近于1，即， （式3-1）在带阻中，幅度响应以误差为逼近于0，即， （式3-2）其中，分别为通带截止频率和阻带截止频率，他们都是数字域频率。为了逼近理想低通滤波器特性，还必须有一个非零宽度的过渡带，在这个过渡带内的频率响应平滑地从通带下降到阻带。图3.2 理想低通滤波器逼近的误差容限3.2 初始条件解析 本次设计带阻数字滤波器，老师给出了以下初始条件： 由以上初始条件可以算出：上边带通带截止频率：；同理可得上边带阻带截止频率：；下边带阻带截止频率：；下边带通带截止频率：。 由以上结果可以算得过
9、渡带。3.3 设计步骤1）选取频率幅度特性的样本点，。样本点数的选取主要考虑过渡带宽度，其中为过渡带点数。为了得到线性相位特性样本点取为实偶对称序列。频率取样法所得FIR滤波器通带和阻带波动主要是由于过渡带的突变引起的。通过在过渡带优化过渡样点可以得到较好的通带和阻带特性，缺点是过渡带加宽，但可以通过增加样点数N来克服。一般一个过渡样点可以使最小阻带衰耗达-40dB，而2个过渡样点可以使最小阻带衰耗达-60dB。因此，为了满足最小阻带衰减，必须取两个过渡样点才能满足题设要求。将和K=2代入中就可以确定抽样点数N=31。2）根据阻带衰耗指标要求，选取频率幅度特性的过渡带样本点，。由1）分析可以知
10、道K=2，并且，为了更好看出效果，我取：，。3）根据傅立叶反变换计算冲击响应系列（FIR滤波器系数）。，。由于是对称的，所以只要计算点。3.4改善滤波器性能的措施 如果给出的理想低通滤波器在通带的频谱等于1而阻带为0，则不论样点N取得如何密，在临界频率处总有两个幅度突变的样点，它们之间的落差为1。于是阻带边缘产生反冲和阻尼振荡，其最大幅度取决于sinc函数，是个固定的值。这样设计出来的滤波器的阻带最小衰耗固定为-20dB，与矩形窗一样。增加采样点数N不能改善阻带最小衰耗。改善阻带衰耗的唯一办法是加宽过渡带。具体方法是：在通、阻带交界处人为地安排一到几个过渡点，其值介于零和1之间，这样可减小样点
11、间的落差，使过渡平缓，反冲减小，阻带最小衰耗增大。经验表明：每多加一个过渡点，过渡带宽增加，最优情况下阻带衰耗可增大2030dB。兼顾过渡带宽和阻带最小衰减。增加采样点，同时在通、阻带交界处安排过渡点。频率采样法特别适用于设计窄带选频滤波器。因为这时只有少数几个非零值的，计算量大为降低。但由于频率抽样点的分布必须符合一定规律，在规定通、阻带截止频率方面不够灵活。比如当截止频率不是整数倍时会产生较大逼近误差，因而该方法的应用不及窗口法普遍。4 MATLAB程序仿真 根据MATLAB相关函数编写如下程序：wsl=0.1067*pi;%低阻带边缘wsh=0.8067*pi;%高阻带边缘wpl=0.3
12、067*pi;%低通带边缘wph=0.6067*pi;%高通带边缘delta=(wpl-wsl);%过度带M=ceil(2*pi*3/delta);%抽样点数al=(M-1)/2;wl=(2*pi/M); %抽样间隔k=0:M-1;T1=0.11; T2=0.59;%过渡带样本点Hrs=ones(1,ceil(0.1067*pi/wl)+1),T2,T1,zeros(1,ceil(0.3*pi/wl),T1,T2,ones(1,ceil(0.3734*pi/wl),T2,T1,zeros(1,ceil(0.3*pi/wl),T1,T2,ones(1,ceil(0.1133*pi/wl)+1);
13、wdl=0 0.1133 0.3133 0.6133 0.8133 1;k1=0:floor(M-1)/2);k2=floor(M-1)/2)+1:M-1;angH=-al*(2*pi)/M*k1,al*(2*pi)/M*(M-k2);H=Hrs.*exp(j*angH);h=real(ifft(H);%傅立叶反变换figure(1);%冲击响应图stem(k,h);title(impulse response);xlabel(n);ylabel(h(n);grid;figure(2);%幅频曲线图Hf=abs(H);w=k*wl/pi;plot(w,Hf,*b-)axis(0 1 -0.1
14、1.1);title(amplitude response);xlabel(frequency in pi units);ylabel(Hr(w);set(gca,xtickmode,manual,xtick,wdl);set(gca,ytickmode,manual,ytick,0 0.11 0.59 1);grid;figure(3);fs=15000;c,f3=freqz(h,1);f3=f3/pi*fs/2;plot(f3,20*log10(abs(c);title(频谱特性);xlabel(频率/HZ);ylabel(衰减/dB);grid; 冲击响应系列如图4.1所示，幅频曲线如图
15、4.2所示，衰减系数和相频曲线如图4.3所示。图4.1冲击响应系列图4.2 幅频曲线图4.3 衰减系数和相频曲线 为了验证我所做滤波器的正确性，继续编程实现验证。在以上程序后面添加以下画图程序：t=(0:100)/fs;x=sin(2*pi*t*750)+sin(2*pi*t*3000)+sin(2*pi*t*6100);q=filter(h,1,x);a,f1=freqz(x);f1=f1/pi*fs/2;b,f2=freqz(q);f2=f2/pi*fs/2;figure(4);subplot(2,1,1);plot(f1,abs(a);title(输入波形频谱图);xlabel(频率);
16、ylabel(幅度)subplot(2,1,2);plot(f2,abs(b);title(输出波形频谱图);xlabel(频率);ylabel(幅度)用figure(4)画第四个图形，即分别输入750Hz、3000Hz、6100Hz验证所做滤波器是否具有滤波功能，用滤波函数filter后，滤波效果如图4.4所示。图4.4 滤波效果图6 设计心得课程设计是培养我们综合运用所学知识,发现、提出、分析和解决实际问题,锻炼实践能力的重要环节。本次课程设计的题目是基于频域抽样法的FIR数字带阻滤波器设计，通过仔细阅读课本相关章节和借阅MATLAB教程书籍，我为实际设计打好了理论基础，在此基础上，通过自
17、己动手设计完成了课程设计要求。通过这次课设，我更进一步理解数字滤波器设计原理，学会了数字滤波器设计的方法和一般步骤，能够独立设计一个数字滤波器，实现了把理论知识转化为实际动手能力的过程。我还从本次课程设计中体会到了MATLAB软件的强大功能，了解到它在各种工程计算中的重要作用，为我以后进一步学习打下了良好的基础。当然这次课程设计也暴露了我的一些问题，比如学习程序设计教程不够快，虽然MATLAB使用的语言和语法都继承于C语言，但还是花了不少时间学习其中的函数，最后才能把课程设计顺利完成。参考文献：1董长虹. MATLAB信号处理与应用.北京：国防工业出版社，20052程佩青.数字信号处理(第2
18、版) M .北京：清华大学出版社，20033王济.MATLAB在振动信号处理中的应用.北京：中国水利水电出版社、知识产权出版社，20064张志涌.精通MATLAB6. 5 版M .北京：北京航空航天大学出版社，20045候正信译. 数字信号处理基础.北京：电子工业出版社，2003程序1wsl=0.1067*pi;%低阻带边缘wsh=0.8067*pi;%高阻带边缘wpl=0.3067*pi;%低通带边缘wph=0.6067*pi;%高通带边缘delta=(wpl-wsl);%过度带N=ceil(2*pi*3/delta)-1;%抽样点数al=(N-1)/2;wl=(2*pi/N); %抽样间隔
19、k=0:N-1;T1=0.11; T2=0.59;%过渡带样本点%建立符幅特性样本序列Hrs=ones(1,ceil(0.1067*pi/wl)+1),T2,T1,zeros(1,ceil(0.3*pi/wl),T1,T2,ones(1,ceil(0.3734*pi/wl),T2,T1,zeros(1,ceil(0.3*pi/wl),T1,T2,ones(1,ceil(0.1067*pi/wl)+1);wdl=0 0.1067 0.3067 0.6067 0.8067 1;k1=0:floor(N-1)/2);k2=floor(N-1)/2)+1:N-1; %k分为两部分k1=0(N-1)/2
20、,k2=(N-1)/2N-1angH=-al*(2*pi)/N*k1,al*(2*pi)/N*(N-k2);%建立相位特性样本序列 求H(k)的相角(k)H=Hrs.*exp(j*angH);%建立频率特性样本序列 求H(k)h=real(ifft(H);%傅立叶反变换求脉冲序列由H(k)傅立叶反变换得到h%冲击响应图subplot(321);stem(k,h);title(冲击响应);xlabel(n);ylabel(h(n);grid;%幅频曲线图Hf=abs(H);w=k*wl/pi;subplot(322);plot(w,Hf,*b-)axis(0 1 -0.1 1.1);title(
21、幅频响应);xlabel(频率/);ylabel(Hr(w);set(gca,xtickmode,manual,xtick,wdl); %设置所要显示数据的位置set(gca,ytickmode,manual,ytick,0 0.11 0.59 1); %给这些数据加标签，用来在指定范围内标注网格线grid;fs=15000;c,f3=freqz(h,1);f3=f3/pi*fs/2;subplot(323);plot(f3,20*log10(abs(c);title(频谱特性);xlabel(频率/HZ);ylabel(衰减 );grid;t=(0:100)/fs;x=sin(2*pi*t*
22、750)+sin(2*pi*t*3000)+sin(2*pi*t*6100);q=filter(h,1,x);a,f1=freqz(x);f1=f1/pi*fs/2;b,f2=freqz(q);f2=f2/pi*fs/2;subplot(324);plot(f1,abs(a);title(输入波形频谱图);xlabel(频率/HZ);ylabel(幅度)subplot(325);plot(f2,abs(b);title(输出波形频谱图);xlabel(频率/HZ);ylabel(幅度)程序2wsl=0.1067*pi;%低阻带边缘wsh=0.8067*pi;%高阻带边缘wpl=0.3067*p
23、i;%低通带边缘wph=0.6067*pi;%高通带边缘delta=(wpl-wsl);%过度带N=ceil(2*pi*3/delta)-1;%抽样点数al=(N-1)/2;wl=(2*pi/N); %抽样间隔k=0:N-1;T1=0.11; T2=0.59;%过渡带样本点%建立符幅特性样本序列Hrs=ones(1,ceil(0.1067*pi/wl)+1),T2,T1,zeros(1,ceil(0.3*pi/wl),T1,T2,ones(1,ceil(0.3734*pi/wl),T2,T1,zeros(1,ceil(0.3*pi/wl),T1,T2,ones(1,ceil(0.1067*pi
24、/wl)+1);wdl=0 0.1067 0.3067 0.6067 0.8067 1;k1=0:floor(N-1)/2);k2=floor(N-1)/2)+1:N-1; %k分为两部分k1=0(N-1)/2,k2=(N-1)/2N-1angH=-al*(2*pi)/N*k1,al*(2*pi)/N*(N-k2);%建立相位特性样本序列 求H(k)的相角(k)H=Hrs.*exp(j*angH);%建立频率特性样本序列 求H(k)h=real(ifft(H);%傅立叶反变换求脉冲序列由H(k)傅立叶反变换得到h%冲激响应图figure(1)subplot(211);stem(k,h);tit
25、le(冲激响应);xlabel(n);ylabel(h(n);grid;%幅频曲线图Hf=abs(H);w=k*wl/pi;subplot(212);plot(w,Hf,*b-)axis(0 1 -0.1 1.1);title(幅频响应);xlabel(频率/);ylabel(Hr(w);set(gca,xtickmode,manual,xtick,wdl); %设置所要显示数据的位置set(gca,ytickmode,manual,ytick,0 0.11 0.59 1); %给这些数据加标签，用来在指定范围内标注网格线grid;%相频响应曲线pha=unwrap(angle(H);figu
26、re(2)subplot(211);plot(w,pha); %相频曲线title(相频响应);xlabel();ylabel();grid;fs=15000;c,f3=freqz(h,1);f3=f3/pi*fs/2;subplot(212);plot(f3,20*log10(abs(c);title(频谱特性);xlabel(频率/HZ);ylabel(衰减 );grid;t=(0:100)/fs;x=sin(2*pi*t*1000)+sin(2*pi*t*3000)+sin(2*pi*t*6000);q=filter(h,1,x);a,f1=freqz(x);f1=f1/pi*fs/2;b,f2=freqz(q);f2=f2/pi*fs/2;figure(3)subplot(211);plot(f1,abs(a);title(输入信号频谱图);xlabel(频率/HZ);ylabel(dB)subplot(212);plot(f2,abs(b);title(输出信号频谱图);xlabel(频率/HZ);ylabel(dB)15
