7-2 第七章FIR设计-2

例1：用矩形窗设计设计FIR低通滤波器：N=21，理想低通截止频率题目N=21，M=20加矩形窗，相当于限制n取值。中心在：12N可以看到矩形窗阻带最小衰减21dB幅频特性：Seeword7-1例题例2：用Hanning窗设计线性相位FIR高通滤波器，N=51，理想低通截止频率M=N-1题目N=51，M=50加Hanning窗：中心在：12N可以看到汉宁窗阻带最小衰减44dB幅频特性：Seeword7-1例题给定指标窗函数设计法目的：设计出h(n)以低通为例：已知未知：M（或N）窗w(n)设计：验证其他指标是否满足要求，不满足再重新设计。我们设计的是第一类线性相位滤波器高通：注意：高通，N必为奇数。带通：带阻：注意：带阻，N必为奇数。我们设计的是第一类线性相位滤波器窗函数法设计步骤1.计算2.根据阻带衰减指标，确定窗形，并选择过渡带最窄的窗。3.再根据该窗过渡带与N的关系，确定N。A依不同窗形状确定，可查表7-3。4.计算出理想滤波器的单位抽样响应，移位并与窗函数相乘。5.设计后验证滤波器指标，必要时作出调整。阻带衰减确定窗形过渡带宽确定N255页表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、窗函数本身均定义在0~N-1之间，所以，窗函数直接乘以理想滤波器函数上，不用移位!3、第一类线性相位数字高通滤波器、数字带阻滤波器，N必须为奇数。若计算出N=M（偶数），取N=M+1。4、理想滤波器的频域表示,必须写上相位(不是1):这里,()0,0jacjdceHe范围表示通带情况，如为高通此处范围表示阻带情况，如N-1=2如果h[n]是偶或奇对称的有限长实数序列,则H(z)具有如下形式:)zz11)(zz11)(zz1)(zz1(1*0101*010z5z4z3*z31/z21/z1*1/z1z1z1*z2系统函数的零点特点Why？(1)1()()NHzzHz(1)2(1)0()()[]NnNnnNHzhnzz偶、偶对称：1(1)12(1)201()()()[]2NNnNnnNNHzhzhnzz奇、偶对称：当z=-1时，H(z)=0(1)2(1)0()()[-]NnNnnNHzhnzz偶、奇对称：当z=1时，H(z)=0不是四个零点：零点在实轴上，在z=1、z=-1时。1(1)12(1)201()()()[-]2NNnNnnNNHzhzhnzz奇、奇对称：当z=1时，z=-1时，H(z)=0例：因果线性相位FIR系统，单位脉冲响应h[n]是实序列，N=7，h[n]=-h[6-n]，已知系统函数的一个零点是0.5ejπ/4，写出系统函数。四种类型广义线性相位FIR滤波器的应用限制N为奇数N为偶数低通高通带通带阻低通高通带通带阻h[n]偶对称YYYYYNYNh[n]奇对称NNYNNYYNN奇偶对称，不受限制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设计步骤1.按阻带衰减要求选定窗函数2.计算过渡带，选较小的一个3.计算N：A与所选窗有关4.计算两个理想低通滤波器的截止频率5.写出两个理想截止频率的低通滤波器的N点单位抽样响应6.理想带通滤波器的单位抽样响应为上述两个低通滤波器之差7.生成所选窗的窗函数8.用窗函数与理想带通滤波器相乘，得到实际带通滤波器的单位抽样响应h(n)9.计算h(n)的幅频响应，用函数freqz(hd,[1],w)10.验证Rp、As是否满足指标，并做调整，记录实际的通带衰减和阻带衰减11.绘出带通滤波器的单位抽样响应图和幅频特性图Seelab3.mDTMF频率的选择是为了避免谐波。Maxandmaxinsteadofcomprison.可以用FIR滤波器组检波；其中两个滤波器输出最大；功率检波（平均功率：功率平方和除以点数）；识别出这两个滤波器检波电路半波整流全波整流直流稳压电源……另一种方法用Goertzel算法要求：1.按一定采样率仿真0-9双音多频信号，播放拨号音，显示时域波形、频域波形。2.任意给定某个双音多频信号，通过Goertzel算法（或其他算法）识别出拨号数字，并显示出来。3.考虑到线路质量，在双音多频信号上加入噪声，再次识别。4.这里涉及到滤波、DFT、系统响应、比较大小的编程。5.能做到哪步就做到哪步，仔细在实验报告上写明你的实验过程。DTMF系统可以用模拟手段或数字手段实现。当用数字手段实现时，可以通过软件来模拟合成DTMF信号。在检测端，可以用8个数字滤波器检测DTMF信号，也可以用其软件实现信号的检测。数字信号的产生可以用查表法实现。电话音频信号在进行数字处理时，其采样周期为125μs（采样频率为8kHz），在这段时间内，也可以用计算两个频率的正弦波信号并将它们相加的方法来实现。考虑到DSP芯片的存储容量有限，显然用增加计算复杂度来换取减少存储容量的方法更好些。检测端收到数据后，对数据进行数字处理以检测出接收的信号。这里的检测可用滤波器组实现，也可用FFT算法实现,当频率成分较少时，前者更有效。2.Goertzel算法N点数据的DFT公式为10)()(NnnkNWnxkXk=0,1,…,N-1(8-1)式中，WN=e-j2π/N。当用FFT算法实现时，共需进行NlgN次的复数乘法与加法。它的好处是可以得到DFT的所有N点的值,但实际上我们只需知道DFT的8个值即可。Goertzel算法从计算DFT入手，推导出一种线性滤波法，在只需计算DFT的M点且MlgN时更有效。我们知道，相位因子WNk具有周期性，且WN－kＮ。利用这一点，在式（8-1）的两边同时乘以它，得到如下表达式：10)()()()(NnnNkNkNNWnxkXWkX(8-2)显然，式（8-2）为一卷积表达式。它是长度为N的有限长输入序列与具有如下单位脉冲响应的滤波器的卷积：)()(nuWnhknNk这个滤波器在n=N处的输出就是DFT在频率点ωk=2πk/N的值，就此可以得到该滤波器的系统函数为111)(zWzHkNk令10)()()(NmmnkNkWmxny则可以通过差分方程用迭代方式来计算：)()1()(nxnyWnykkNkyk(-1)=0当n=N时，X(k)=yk(n)。)()(nuWnhknNk10)()()()(NnnNkNkNNWnxkXWkX对于每一个k，W-kN为一个固定值，它可以预先计算出来并作为常数项使用,但该值一般为复数，给计算带来不便。通过对滤波器系统函数的变换：211)/π2cos(211)(zzNkzWzHkNk可以推导出如下系统实现的迭代公式：)1()()()()2()1(π2cos2)(nvWnvnynxnvnvNknvkkNkkkkk初始条件是：vk(n-1)=vk(n-2)=0。111)(zWzHkNk这就是Goertzel算法。每输入一个样值，分别用迭代法计算出各个vk(n)，每次计算只需进行两次实数乘和两次实数加。当n=N时，用式（8-3）计算出各yk(n)的值，这些值就是所要求的各点的DFT值。另外，考虑到最后对信号进行判断时需要幅度值|X(k)|或幅度的平方值|X(k)|2，因此对式（8-3）进行变换：NkNvNvNvWNvNykXkkkkNkk2cos2)1()(|)1()(||)(||)(|22222这样，算法中就避免了复数运算，最后对每个值点再进行四次乘法和二次加法运算即可。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(′请输入数字（0~9）：＼n′);forj=1:800x=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)