实验4-基于MATLAB的FIR数字滤波器设计.doc

实验4 基于MATLAB的FIR数字滤波器设计 作者： 日期：19 个人收集整理 勿做商业用途实验4 基于MATLAB的FIR数字滤波器设计实验目的:加深对数字滤波器的常用指标和设计过程的理解。实验原理：低通滤波器的常用指标：通带边缘频率，阻带边缘频率 ,通带起伏，通带峰值起伏，阻带起伏，最小阻带衰减。数字滤波器有IIR和FIR两种类型，它们的特点和设计方法不同。 在MATLAB中，可以用b=fir1（N,Wn，ftype,taper) 等函数辅助设计FIR数字滤波器。N代表滤波器阶数；Wn代表滤波器的截止频率(归一化频率），当设计带通和带阻滤波器时,Wn为双元素相量;ftype代表滤波器类型，如high高通，stop带阻等;taper为窗函数类型，默认为海明窗，窗系数需要实现用窗函数blackman， hamming,hanning chebwin, kaiser产生。例1 用凯塞窗设计一FIR低通滤波器，通带边界频率 ,阻带边界频率 ,阻带衰减 不小于50dB。 解 首先由过渡带宽和阻带衰减 来决定凯塞窗的N和 ， , 上图给出了以上设计的频率特性，(a） 为N=30直接截取的频率特性(b）为凯塞窗设计的频率特性.凯塞窗设计对应的MATLAB程序为： wn=kaiser（30,4.55); nn=0:1:29; alfa=（30-1）/2; hd=sin(0.4*pi（nn-alfa)）./(pi*（nnalfa); h=hd.wn; h1，w1=freqz(h，1); 或者：b = fir1（29,0。4，kaiser（30，4.55)）； h1，w1=freqz(b,1)； plot(w1/pi,20log10(abs(h1）)）; axis（0，1,80，10）; grid;xlabel（归一化频率/p) ;ylabel（幅度/dB） ；还可以使用n,Wn，beta,ftype = kaiserord（f，a,dev)函数来估计滤波器阶数等，得到凯塞窗滤波器：fcuts = 0.3 0。5; %归一化频率omega/pimags = 1 0;devs = 0。05 10（-2.5)；n,Wn，beta，ftype = kaiserord（fcuts,mags,devs)； 计算出凯塞窗N，beta的值hh = fir1（n,Wn,ftype，kaiser（n+1，beta)，noscale); freqz（hh）;实际中，一般调用MATLAB信号处理工具箱函数remezord来计算等波纹滤波器阶数N和加权函数W（），调用函数remez可进行等波纹滤波器的设计,直接求出滤波器系数。函数remezord中的数组fedge为通带和阻带边界频率，数组mval是两个边界处的幅值，而数组dev是通带和阻带的波动，fs是采样频率单位为Hz。例2 利用雷米兹交替算法设计等波纹滤波器，设计一个线性相位低通FIR数字滤波器，其指标为：通带边界频率fc=800Hz，阻带边界fr=1000Hz，通带波动 阻带最小衰减At=40dB，采样频率fs=4000Hz。 解 在MATLAB中可以用remezord 和remez两个函数设计，其结果如图2，MATLAB程序如下: fedge=800 1000； mval=1 0; dev=0。0559 0。01； fs=4000; N，fpts，mag,wt=remezord（fedge，mval,dev，fs）； b=remez（N,fpts，mag，wt）; h,w=freqz(b,1,256）； plot(w2000/pi，20log10（abs(h）)）; grid; xlabel（频率/Hz） ； ylabel(幅度/dB); 一、实验内容： 利用MATLAB编程设计一个数字带通滤波器，指标要求如下：通带边缘频率：，,通带峰值起伏:。阻带边缘频率：，最小阻带衰减： 。 分别用窗函数法和等波纹滤波器法设计两种FIR数字滤波器。实验要求：给出FIR数字滤波器的冲激响应，绘出它们的幅度和相位频响曲线，讨论它们各自的实现形式和特点。1-1）用窗函数法实现：调用函数n，wn，bta,ftype=kaiserord（f，a，dev,fs）参数：f=0.3 0.45 0。65 0.8为对应数字频率， a=0 1 0为由f指定的各个频带上的幅值向量，一般只有0和1表示；和f长度关系为（2a的长度)-2=(f的长度）devs=0。01 0.1087 0.01用于指定各个频带输出滤波器的频率响应与其期望幅值之间的最大输出误差或偏差,长度与a相等，计算公式：阻带衰减误差=，通带波动衰减误差=fs缺省值为2HZ n,wn,bta，ftype=kaiserord(0。3 0.45 0。65 0。8,0 1 0，0。01 0.1087 0。01)；用kaiserord函数估计出滤波器阶数n和beta参数 h1=fir1（n,wn,ftype，kaiser（n+1，bta），noscale); hh1,w1=freqz(h1，1，256）; figure（1） subplot（2，1,1) plot(w1/pi,20*log10(abs(hh1）) grid xlabel（归一化频率w）;ylabel（幅度/db）； subplot(2,1，2) plot(w1/pi，angle（hh1） grid xlabel（归一化频率w);ylabel(相位/rad）;h1 = Columns 1 through 8 0。0041 0.0055 -0.0091 0.0018 0。0056 -0.0000 0。0391 0。0152 Columns 9 through 16 0.0381 0.0077 -0.0293 0.0940 0.0907 0.2630 0。0517 0.3500 Columns 17 through 24 0。0517 0。2630 0。0907 0。0940 0.0293 0。0077 0。0381 -0.0152 Columns 25 through 31 0。0391 -0。0000 -0。0056 0。0018 0.0091 0。0055 0.0041 图41如果直接用freqz（h1,1,256）;画图得： 12）用等波纹法设计：调用函数n，fpts，mag,wt=remezord（f，a，dev）f=0。3 0。45 0.65 0。8a=0 1 0dev=0.01 0.1087 0.01其含义同函数n，wn,bta,ftype=kaiserord（f，a，dev,fs）中的参数相同。 n,fpts，mag，wt=remezord（0.3 0.45 0。65 0.8,0 1 0，0.01 0.1087 0.01);%用remezord函数估算出remez函数要用到的阶n、归一化频带边缘矢量fpts、频带内幅值响应矢量mag及加权矢量w，使remez函数设计出的滤波器满足f、a及dev指定的性能要求。 h2=remez(n，fpts,mag，wt）；设计出等波纹滤波器 hh2，w2=freqz（h2,1,256); figure（2） subplot(2,1,1) plot(w2/pi，20log10(abs（hh2）) grid xlabel(归一化频率w）;ylabel(幅度/db）； subplot(2,1,2） plot（w2/pi,angle（hh2) grid xlabel（归一化频率w)；ylabel(相位/rad）; h2h2 = Columns 1 through 9 -0.0013 0。0092 -0.0255 0。0642 0.1177 0。0922 0。2466 0。0466 0。3116 Columns 10 through 17 -0。0466 0.2466 0。0922 0。1177 -0。0642 -0。0255 0。0092 -0.0013 图4-2用freqz（h2，1，256）；直接得图：二、对课本作业9.23画图。21)用汉宁窗实现:1）手动计算 n=0.001：58。001; hd=sin（0.18125pi*(n-29)）。/(pi(n-29）； win=0.5+0。5*cos(2pi(n29)/58）; h1=2cos（pi(n29）/2).hd。*win； hh1，w1=freqz（h1,1，256）； figure（1） subplot（2，1，1） plot(w1,20*log10（abs(hh1)）) grid xlabel（数字频率w/rad);ylabel（幅度/db）; subplot(2,1，2) plot（w1，angle(hh1) grid xlabel(数字频率w/rad)；ylabel(相位/rad）; 图4322）用自带的fir1函数： n=59; wn=3/8 5/8； h2=fir1(n，wn，bandpass，hann(n+1)； hh2,w2=freqz（h2，1,256）; figure(2) subplot(2，1,1) plot(w2,20log10(abs（hh2）) xlabel(数字频率w/rad);ylabel(幅度/db); grid subplot（2，1，2） plot（w2，angle(hh2)） xlabel(数字频率w/rad);ylabel（相角/rad）; grid 图4-42-3）用等波纹法设计：调用函数n,fpts，mag，wt=remezord(f，a，dev）f=0.2625 0。375 0。625 0.7375a=0 1 0dev=0.01 0。1087 0.01 n,fpts,mag，wt=remezord（0。2625 0.375 0.625 0。7375，0 1 0,0。01 0。1087 0。01）； h2=remez(n,fpts，mag,wt）; hh2,w2=freqz（h2，1,256); figure（2） subplot(2,1，1) plot(w2/pi,20*log10(abs（hh2）) grid xlabel（归一化频率w）;ylabel(幅度/db); subplot(2,1,2) plot（w2/pi,angle(hh2）) grid xlabel（归一化频率w）；ylabel（相位/rad); h2画图为：另外带通滤波器还可以用低通和高通级联的方法实现,因为步骤过多，这里不作讨论。三、实验分析：1）总结:FIR滤波器实现一般采用窗函数法和等纹波设计法.窗函数法还包含两个分支,一种是用公式先手动算出N值和其他对应得窗函数参数值，再代入窗函数和fir1实现，一种是用函数*rord估算出N和相应参数再用fir1实现。不过要注意*rord会低估或高估阶次n,可能会使滤波器达不到指定的性能,这时应稍微增加或降低阶次。如果截止频率在0或Nyquist频率附近,或者设定的dev值较大,则得不到正确结果。2）滤波器实现形式及特点： 由于一般的滤波器在利用窗函数是其通带波纹和阻带波纹不同（一般为第一个阻带波纹最大）因此，在满足第一个阻带衰减旁瓣时，比其频率高的旁瓣，它们的衰减都大大超出要求.而根据阻带衰减与项数的近似关系可得当阻带衰减越大，所需项数越多.2)等波纹设计和窗函数设计不同之处在于实现形式及特点： 窗函数设计是通过最小平方积分办法来设计的即该滤波器的误差为： 即要求最小方法来设计滤波器，这样的滤波器更忠实于理想滤波器（即滤波系数更接近于理想滤波器。证明如下： = = = = 因此,幅度频谱差值越小，实际滤波器就越接近理想滤波器。 而等波纹滤波器是通过最大加权误差最小化来实现，其误差为： 要求该误差最小来实现滤波器，得出来的滤波系数较窗函数设计相差较远。以下通过对题目中的h1及h2作比较。% sigsum是用来对数组各元素进行求和function y=sigsum（n1,n2，n，x）;y=0； for i=n1+1min（n）:n2+1-min（n) y=y+x（i）;end n=0。001:30.001； h=2*cos（0。55pi（n-15).*sin（0。175pi*（n15)）./（pi*（n15)）; delta1=h-h1; n=0.001:16.001; h=2cos（0。55pi*（n-15).sin（0。175pi*（n-15)./(pi*（n15); delta2=hh2; y1=sigsum（0,30,0：30,（abs(delta1）.2)）/31; y2=sigsum（0,16,0:16，（abs（delta2）.2)/17； y1y1 = 1.9099e-004 y2y2 = 0.0278由此得到用窗函数实现的滤波系数比用等波纹滤波器系数的每一项更接近于理想滤波器（y1为用窗函数实现的与理想滤波器的差值，y2为用等波纹滤波器实现的与理想滤波器的差值)；对比图4-1及图4-2可见在幅度频谱上等波纹滤波器阻带边缘比用窗函数实现的更平滑（理想滤波器为垂直下降的）。从设计的角度考虑，由于窗函数设计法都是通过已有的窗函数对理想滤波器的改造，因此，可以用手算的办法方便的设计滤波器.而等波纹滤波器，其实现是通过大量的迭代运算来实现，这样的方法一般只能通过软件来设计。项数的问题由于等波纹滤波器能较平均的分布误差，因此对于相同的阻带衰减，其所需的滤波系数比窗函数的要少.图41图42从设计的角度考虑,由于窗函数设计法都是通过已有的窗函数对理想滤波器的改造，因此，可以用手算的办法方便的设计滤波器。而等波纹滤波器，其实现是通过大量的迭代运算来实现，这样的方法一般只能通过软件来设计.项数的问题由于等波纹滤波器能较平均的分布误差,因此对于相同的阻带衰减，其所需的滤波系数比窗函数的要少。3)画相角时采用不同的单位其形状不一样：（如下图）上面第一个图是用角度为单位画出来的，下面的图是用rad单位画出来的。从图形可以观察到在0。3到0。8数字频率间两个图都是严格的线性相位,至于下面的图为什么在这个区间会有跳变是因为rad的区间只有-pipi,当相位由pi继续增加时只能跳到pi而不能大于pi，而角度表示则可以连续增大。疑问:两个图中都可以看出,在小于0。3和大于0.8的区间并不严格按照线性相位，这可能是软件仿真的问题和设计的滤波器与理想滤波器之间的误差（通带波纹和阻带截频）造成的。4)调用firl或者reme函数时，用scale（缺省方式）对滤波器进行归一化，即滤波器通带中心频率处的响应幅值为0db。用noscale不对滤波器归一化。