第七章FIR设计

例1：用矩形窗设计设计FIR低通滤波器：N=21，,理想低通截止频率,题目N=21，M=20加矩形窗，相当于限制n取值。中心在：,可以看到矩形窗阻带最小衰减21dB,幅频特性：,Seeword7-1例题,例2：用Hanning窗设计线性相位FIR高通滤波器，N=51，,理想低通截止频率M=N-1,题目N=51，M=50加Hanning窗：中心在：,可以看到汉宁窗阻带最小衰减44dB,幅频特性：,Seeword7-1例题,给定指标窗函数设计法,目的：设计出h(n),以低通为例：已知,未知：,M（或N）窗w(n),设计：,验证其他指标是否满足要求，不满足再重新设计。,我们设计的是第一类线性相位滤波器,高通：,注意：高通，N必为奇数。,带通：,带阻：,注意：带阻，N必为奇数。,我们设计的是第一类线性相位滤波器,窗函数法设计步骤,计算根据阻带衰减指标，确定窗形，并选择过渡带最窄的窗。再根据该窗过渡带与N的关系，确定N。A依不同窗形状确定，可查表7-3。计算出理想滤波器的单位抽样响应，移位并与窗函数相乘。设计后验证滤波器指标，必要时作出调整。,阻带衰减确定窗形过渡带宽确定N,255页表7-3,例3：设计一个数字FIR低通滤波器，其技术指标：,解：衰减50dB的窗查表为海明窗或Blackman窗，选用海明窗，它给出比较小的过渡带，因此有较低的阶。计算N：,see：chap7_expl.m,计算截止频率,Seeword7-1例题例题3仿真图,FIR滤波器设计完毕的标识是：h(n)求出。,再看各项指标是否满足若不满足或改变N或改变窗形因为FIR滤波器不能与要求的幅频特性完全满足,写出h(n)表达式：,注意：1、理想滤波器的单位冲激响应需移位(N-1)/2，不管N是奇数还是偶数，以保证对称，保证线性相位。其线性相位为：用的是第一类线性相位条件，即偶对称：()=-(N-1)/22、窗函数本身均定义在0N-1之间，所以，窗函数直接乘以理想滤波器函数上，不用移位!3、第一类线性相位数字高通滤波器、数字带阻滤波器，N必须为奇数。若计算出N=M（偶数），取N=M+1。4、理想滤波器的频域表示,必须写上相位(不是1):这里,如果hn是偶或奇对称的有限长实数序列,则H(z)具有如下形式:,系统函数的零点特点,Why？,当z=-1时，H(z)=0,当z=1时，H(z)=0,不是四个零点：零点在实轴上，在z=1、z=-1时。,当z=1时，z=-1时，H(z)=0,例：因果线性相位FIR系统，单位脉冲响应hn是实序列，N=7，hn=-h6-n，已知系统函数的一个零点是0.5ej/4，写出系统函数。,四种类型广义线性相位FIR滤波器的应用限制,N奇偶对称，不受限制N偶偶对称，当z=-1时，H(z)=0N奇奇对称，当z=1时，z=-1时，H(z)=0N偶奇对称，当z=1时，H(z)=0,作业：7-47-14（4）并交换p和s设计高通滤波器讲解：7-6,设计步骤按阻带衰减要求选定窗函数计算过渡带，选较小的一个计算N：A与所选窗有关计算两个理想低通滤波器的截止频率写出两个理想截止频率的低通滤波器的N点单位抽样响应理想带通滤波器的单位抽样响应为上述两个低通滤波器之差生成所选窗的窗函数用窗函数与理想带通滤波器相乘，得到实际带通滤波器的单位抽样响应h(n)计算h(n)的幅频响应，用函数freqz(hd,1,w)验证Rp、As是否满足指标，并做调整，记录实际的通带衰减和阻带衰减绘出带通滤波器的单位抽样响应图和幅频特性图Seelab3.m,DTMF,频率的选择是为了避免谐波。,Maxandmaxinsteadofcomprison.,可以用FIR滤波器组检波；其中两个滤波器输出最大；功率检波（平均功率：功率平方和除以点数）；识别出这两个滤波器,检波电路,半波整流,全波整流直流稳压电源,另一种方法用Goertzel算法,要求：按一定采样率仿真0-9双音多频信号，播放拨号音，显示时域波形、频域波形。任意给定某个双音多频信号，通过Goertzel算法（或其他算法）识别出拨号数字，并显示出来。考虑到线路质量，在双音多频信号上加入噪声，再次识别。这里涉及到滤波、DFT、系统响应、比较大小的编程。能做到哪步就做到哪步，仔细在实验报告上写明你的实验过程。,DTMF系统可以用模拟手段或数字手段实现。当用数字手段实现时，可以通过软件来模拟合成DTMF信号。在检测端，可以用8个数字滤波器检测DTMF信号，也可以用其软件实现信号的检测。数字信号的产生可以用查表法实现。电话音频信号在进行数字处理时，其采样周期为125s（采样频率为8kHz），在这段时间内，也可以用计算两个频率的正弦波信号并将它们相加的方法来实现。考虑到DSP芯片的存储容量有限，显然用增加计算复杂度来换取减少存储容量的方法更好些。检测端收到数据后，对数据进行数字处理以检测出接收的信号。这里的检测可用滤波器组实现，也可用FFT算法实现,当频率成分较少时，前者更有效。,2.Goertzel算法N点数据的DFT公式为,k=0,1,N-1,(8-1),式中，WN=e-j2/N。,当用FFT算法实现时，共需进行NlgN次的复数乘法与加法。它的好处是可以得到DFT的所有N点的值,但实际上我们只需知道DFT的8个值即可。Goertzel算法从计算DFT入手，推导出一种线性滤波法，在只需计算DFT的M点且MlgN时更有效。,我们知道，相位因子WNk具有周期性，且WNk。利用这一点，在式（8-1）的两边同时乘以它，得到如下表达式：,(8-2),显然，式（8-2）为一卷积表达式。它是长度为N的有限长输入序列与具有如下单位脉冲响应的滤波器的卷积：,这个滤波器在n=N处的输出就是DFT在频率点k=2k/N的值，就此可以得到该滤波器的系统函数为,令,则可以通过差分方程用迭代方式来计算：,yk(-1)=0,当n=N时，X(k)=yk(n)。,对于每一个k，W-kN为一个固定值，它可以预先计算出来并作为常数项使用,但该值一般为复数，给计算带来不便。通过对滤波器系统函数的变换：,可以推导出如下系统实现的迭代公式：,初始条件是：vk(n-1)=vk(n-2)=0。,这就是Goertzel算法。每输入一个样值，分别用迭代法计算出各个vk(n)，每次计算只需进行两次实数乘和两次实数加。当n=N时，用式（8-3）计算出各yk(n)的值，这些值就是所要求的各点的DFT值。另外，考虑到最后对信号进行判断时需要幅度值|X(k)|或幅度的平方值|X(k)|2，因此对式（8-3）进行变换：,这样，算法中就避免了复数运算，最后对每个值点再进行四次乘法和二次加法运算即可。,3.DTMF信号产生与检测算法的实现在检测信号时，为避免信号前后的0值带来的对DFT计算的影响，同时考虑到DTMF信号的高频和低频分量的最小间隔为73Hz，故选取每200点为一段，分别用Goertzel算法计算每一段的8个频率点的DFT的平方值，然后设定一个阈值对得到的值进行判断（若某一频率点的幅度平方值大于阈值，判为有效，否则为无效），再根据有效频率点的组合给出最后检测到的数字。为了模拟接收到DTMF信号在数据流中的开始位置的检测效果，可以在产生的采样数据前加上不同个数的0。,具体的程序如下（首先输入要产生和检测的信号个数，然后输入加在抽样数据前的0的个数，最后输入单个数字）：functionDTMF%用来模拟产生和检测DTMF信号s=input(请输入检测的数字个数：n);l=input(请输入加在数据前面的0的个数：n);out=zeros(1600*s+1,1);%设置输入样值数组tab=9411336;6971209;6971336;6971447;7701209;7701336;7701477;8521209；8521336;8521477;forI=1:sk=input(请输入数字（09）：n);,forj=1:800 x=tab(k+1,1,)/8000;y=tab(k+1,2)/8000;z=sin(x*j*2*pi)+sin(y*j*2*pi);out(1600*(i-1)+j+1,1)=14000*z;%合成信号采样值out(1600*(i-1)+j+800+1,1)=0endendout=out./28000;wavwrite(out,8000,16,DTMF1);%写入DTMF1.WAV文件w=zeros(8);w(1)=2*cos(2*pi*17/200);w(2)=2*cos(2*pi*19/200);,w(3)=2*cos(2*pi*21/200);w(4)=2*cos(2*pi*23/200);w(5)=2*cos(2*pi*130/200);w(6)=2*cos(2*pi*33/200);w(7)=2*cos(2*pi*36/200);w(8)=2*cos(2*pi*40/200);result=zeros(8*s,8);a=wavread(DTMF);读出DTMF1.WAV文件中的数据fori=1:8*s;b=zeros(8,3);forj=1:200;fork=1:8b(k,3)=w(k)*b(k,2)-b(k,1)+a(i-1)*200+j);,b(k,1)=b(k,2);b(k,2)=b(k,3);endendfork=1:8result(i,k)=b(k,2)2+b(k,1)2-w(k)*b(k,1)*b(k,2)-600;if(result(i,k)0)result(i,k)=0;endendendfori=1:8*s判断数字值a=nnz(result(i,:)if(a=2),q=find(result(i,:);if(q(1)=1)if(q(2)=5)res(i)=1;elseif(q(2)=6)res(i)=2;elseif(q(2)=7)res(i)=3;endelseif(q(1)=2)if(q(2)=5)res(i)=4;elseif(q(2)=6)res(i)=5;,res(i)=6;endelseif(q(1)=3)if(q(2)=5)res(i)=7;elseif(q(2)=6)res(i)=8;elseif(q(2)=7)res(i)=9;endelseif(q(1)=4)if(q(2)=6)res(i)=0;,endendelseres(i)=NaN;endendres=res,只要在MATLAB命令窗口中输入：DTMF可得如下结果：请输入检测的数字个数：4请输入加在数据前面的0的个数：100请输入数字（09）：5请输入数字（09）：7请输入数字（09）：1,请输入数字（09）：4res=Columns1through12NaN555NaNNaNNaNNaNNaN777Columns13through24NaNNaNNaNNaNNaN111NaNNaNNaNNaNColumns25through32444NaNNaNNaNNaN,从结果中可以看到，已准确地检测出了输入的4个数字（NaN表示无信号）。通过上面的实例，可以看到用MATLAB进行模拟的简便性，若用C语言来编写，运算和数据的输入输出将占用较大篇幅。,Matlab产生窗函数的文件1.bartlett(N)（三角窗）2.blackman(N)（布莱克曼窗）3.boxcar(N)（矩形窗）4.hamming(N)（海明窗）5.hanning(N)（汉宁窗）6.chebwin(N)（切比雪夫窗）7.kaiser(N)（凯泽窗）,优点：1.无稳定性问题；2.容易做到线性相位；3.可以设计各种特殊类型的滤波器；4.方法特别简单。,缺点：1.不易控制边缘频率；2.幅频性能不理想；3.较长。,FIRDF设计的窗函数法的特点,其它设计方法？,LPF：,HPF：,BPF：,BSF：,关于窗函数法设计FIR滤波器，错误的说法是（B）（A）增加窗长可以减小过渡带宽；（B）增加窗长可以增加阻带衰减；（C）改变窗形状可以改变阻带衰减；（D）改变窗形状可以改变过渡带宽。关于窗函数法设计FIR滤波器，正确的说法是（B）（A）先根据过渡带宽要求选择窗形状，再根据阻带衰减要求选择窗长；（B）先根据阻带衰减要求选择窗形状，再根据过渡带宽要求选择窗长；（C）先根据过渡带宽要求选择窗长，再根据阻带衰减要求选择窗形状；（D）先根据阻带衰减要求选择窗长，再根据过渡带宽要求选择窗形状。,窗函数法设计线性相位FIR滤波器时，滤波器的阻带衰减由窗形状（窗形状/窗长）决定，过渡带宽由窗形状和窗长（窗形状/窗长）决定。,图中两个序列作为FIR滤波器的单位脉冲响应，正确的说法是（D）,解：第2个序列在时域更平缓地过度到零，所以高频少，作为滤波器的单位脉冲响应，则阻带衰减大，相应地在过度带宽方面付出代价。频响如下图。,FIR滤波器的窗函数法设计步骤：（1）确定FIR滤波器的技术指标,（2）根据技术指标确定理想滤波器的频率响应函数（3）求出理想滤波器的单位冲激响应,（4）根据技术指标要求选择适当的窗函数并确定窗函数的长度N；（5）对理想滤波器加窗得到实际滤波器的单位冲激响应h(n)=w(n)hd(n)；,（6）求实际滤波器的频率响应H(ejw)=FTh(n)并检验设计结果是否符合技术指标要求，若不满足须重新设计,【例】设计线性相位FIR低通滤波器，给定采样频率fs=15kHz，通带截止频率p=2p1.5103rad/s，阻带截止频率s=2p3103rad