FIR滤波器.

FIR滤波器的设计一、窗函数设计法MATLAB提供了几个子程序来实现各种窗函数。w=boxcar(M)数组w中返回M点矩形窗函数w=triang(M)数组w中返回M点Bartlett(三角)窗函数w=hanning(M)数组w中返回M点汉宁窗函数w=hamming(M)数组w中返回M点哈明窗函数w=blackman(M)数组w中返回M点布莱克曼窗函数w=kaiser(M,beta)数组w中返回beta值M点凯泽窗函数设计FIR滤波器还需要一个理想低通脉冲响应hd(n)。下面的子程序用于产生hd(n)。functionhd=ideal_lp(wc,M);%IdealLowPassfiltercomputation%[hd]=ideal_lp(wc,M)%hd=0到M-1之间的理想脉冲响应%wc=截止频率（弧度）%M=理想滤波器的长度alpha=(M-1)/2;n=[0:1:(M-1)];m=n-alpha+eps;hd=sin(wc*m)./(pi*m);根据给定的滤波器技术指标，选择滤波器长度M和窗函数w(n)窗函数名称旁瓣峰值衰减近似过渡带宽精确过渡带宽最小阻带衰减矩形13db4π/M1.8π/M21dB三角25db8π/M6.1π/M25dB汉宁31db8π/M6.2π/M44dB哈明41db8π/M6.6π/M53dB布莱克曼57db12π/M11π/M74dB此外，凯泽窗的设计公式如下：给定ωp、ωs、Rp、As归一化过渡带宽=Δf=(ωs-ωp)/2π滤波器阶数0.1102(As-8.7)As≥50参数β=0.5842(As-21)0.4+0.07886(As-21)21As501f36.1495.7AMs例1：根据下列技术指标，是一个数字FIR低通滤波器：通带截止频率：Ωp=0.2π，通带波动：Rp=0。25dB阻带截止频率：Ωs=0.3π，阻带波动：As=50dBwp=0.2*pi;ws=0.3*pi;%输入滤波器指标tr_width=ws-wp;%计算过渡带宽M=ceil(6.6*pi/tr_width)%计算滤波器阶数Mn=[0:1:M-1];wc=(ws+wp)/2;%近似截止频率wchd=ideal_lp(wc,M);%计算理想hd(n)w_ham=(hamming(M))‘;%计算窗函数(哈明窗)h=hd.*w_ham;%计算h(n)%画图[H,W]=freqz(h,1);subplot(2,2,1);stem(n,hd);title('理想脉冲响应');axis([0,M-1,-0.1,0.3]);gridsubplot(2,2,2);stem(n,w_ham);title('哈明窗');axis([0,M-1,0,1.1]);gridsubplot(2,2,3);stem(n,h);title('实际脉冲响应');axis([0,M-1,-0.1,0.3]);gridsubplot(2,2,4);plot(W/pi,20*log10(abs(H)));title('幅度响应');axis([0,1,-100,10]);grid例2：利用例1给出的设计技术指标，选择凯泽窗,设计出所需的低通滤波器。wp=0.2*pi;ws=0.3*pi;As=50;tr_width=ws-wp;M=ceil((As-7.95)/(14.36*tr_width/(2*pi))+1)n=[0:1:M-1];beta=0.1102*(As-8.7)wc=(ws+wp)/2;hd=ideal_lp(wc,M);w_kai=(kaiser(M,beta))';h=hd.*w_kai;%画图[H,W]=freqz(h,1);subplot(2,2,1);stem(n,hd);title('理想脉冲响应');axis([0,M-1,-0.1,0.3]);gridsubplot(2,2,2);stem(n,w_kai);title('凯泽窗');axis([0,M-1,0,1.1]);gridsubplot(2,2,3);stem(n,h);title('实际脉冲响应');axis([0,M-1,-0.1,0.3]);gridsubplot(2,2,4);plot(W/pi,20*log10(abs(H)));title('幅度响应');axis([0,1,100,10]);grid例3：设计下面的数字带通滤波器.低阻带边缘：ω1s=0.2π，As=60dB低通带边缘：ω1p=0.35π，Rp=1dB高通带边缘：ω2p=0.65π，Rp=1dB高阻带边缘：ω2s=0.8π，As=60dBws1=0.2*pi;wp1=0.35*pi;wp2=0.65*pi;ws2=0.8*pi;As=60;tr_width=min((wp1-ws1),(ws2-wp2))M=ceil(11*pi/tr_width)n=[0:1:M-1];wc1=(ws1+wp1)/2;wc2=(wp2+ws2)/2;hd=ideal_lp(wc2,M)-ideal_lp(wc1,M);w_bla=(blackman(M))';h=hd.*w_bla;%画图[H,W]=freqz(h,1);subplot(2,2,1);stem(n,hd);title('理想脉冲响应');gridsubplot(2,2,2);stem(n,w_bla);title('布莱克曼窗');grid;subplot(2,2,3);stem(n,h);title('实际脉冲响应');grid;subplot(2,2,4);plot(W/pi,20*log10(abs(H)));title('幅度响应（dB）');grid;二、利用fir1函数设计FIR滤波器MATLAB提供函数fir1实现FIR滤波器的窗函数设计。函数格式：b=fir1(n,wn);b=fir1(n,wn,’ftype’);b=fir1(n,wn,window);b=fir1(n,wn,’ftype’,window);fir1函数用经典方法实现加窗线性相位FIR数字滤波器的设计。（1）b=fir1(n,wn);可以设计n阶线性相位的低通FIR数字滤波器，滤波器系数于b中。B(z)=b(1)+b(2)z-1+b(3)z-2+…+b(n+1)z-n注：①此时用默认的哈明窗来设计FIR滤波器。②wn以pi为单位，0≤wn≤1。③当wn=[w1w2]时，则表示设计的滤波器是带通滤波器。④b(n)的下标从1开始。（2）b=fir1(n,wn,’ftype’);可以设计n阶线性相位的高通或带阻FIR数字滤波器，滤波器系数于b中。B(z)=b(1)+b(2)z-1+b(3)z-2+…+b(n+1)z-n当’ftype’为’high’时，设计高通FIR滤波器。当’ftype’为’stop’时，设计带阻FIR滤波器，此时wn=[w1w2]。在设计高通或带阻滤波器时，若n为奇数，则函数fir1会将其自动加1，使其成为偶数。这样就保证了h(n)有奇数个，从而满足线性相位高通和带阻FIR滤波器的设计要求。此时同样默认为哈明窗。（3）b=fir1(n,wn,window);它与(1)的区别在于可以自由地选择不同的窗来完成线性相位FIR滤波器的设计。w=boxcar(M)数组w中返回M点矩形窗函数w=triang(M)数组w中返回M点Bartlett(三角)窗函数w=hanning(M)数组w中返回M点汉宁窗函数w=hamming(M)数组w中返回M点哈明窗函数w=blackman(M)数组w中返回M点布莱克曼窗函数w=kaiser(M,beta)数组w中返回beta值M点凯泽窗函数注意：窗的长度应为M=n+1，它应该与h(n)的点数相同。（4）b=fir1(n,wn,’ftype’,window);它可以用任意的窗函数设计不同的线性相位的高通或带阻FIR滤波器。例4：根据下列技术指标，设计一个数字FIR低通滤波器：wp=0.2π，ws=0.3π,Rp=0.25dB，As=50dB。（1）若选择哈明窗设计：wp=0.2*pi;ws=0.3*pi;tr_width=ws-wp;M=ceil(6.6*pi/tr_width)n=[0:1:M-1];wc=(ws+wp)/2;h=fir1(M,wc/pi);[H,W]=freqz(h,1);plot(W/pi,20*log10(abs(H)));（2）若选择布莱克曼窗设计：wp=0.2*pi;ws=0.3*pi;tr_width=ws-wp;M=ceil(6.6*pi/tr_width)n=[0:1:M-1];wc=(ws+wp)/2;h=fir1(M,wc/pi,blackman(M+1));[H,W]=freqz(h,1);plot(W/pi,20*log10(abs(H)));三、频率采样设计法例5：用频率采样方法设计一个FIR滤波器，其设计指标为wp=0.2π，ws=0.3π,Rp=0.25dB，As=50dB。M=20;alpha=(M-1)/2;l=0:M-1;wl=(2*pi/M)*l;Hrs=[1,1,1,zeros(1,15),1,1];Hdr=[1,1,0,0];wdl=[0,0.25,0.25,1];H=Hrs.*exp(j*angH);h=real(ifft(H,M));[Ha,w]=freqz(h,1);[Hr,ww,a,L]=Hr_Type2(h);%Hr_Type2函数是自定义函数，请参阅附录。subplot(2,2,1);plot(wl(1:11)/pi,Hrs(1:11),'o',wdl,Hdr);axis([0,1,-0.1,1.1]);title('频率样本:M=20')xlabel('频率(单位:Pi)');ylabel('Hr(k)');gridsubplot(2,2,2);stem(l,h);axis([-1,M,-0.1,0.3])title('脉冲响应');xlabel('n');ylabel('h(n)');gridsubplot(2,2,3);plot(ww/pi,Hr,wl(1:11)/pi,Hrs(1:11),'o');axis([0,1,-0.2,1.2]);title('振幅响应')xlabel('频率(单位:Pi)');ylabel('Hr(w)');gridsubplot(2,2,4);plot(w/pi,20*log10(abs(Ha)));axis([0,1,-60,10]);gridtitle('幅度响应');xlabel('频率(单位:Pi)');ylabel('分贝');