数字信号处理实验报告6

实验报告学院（系）名称：计算机与通信工程学院姓名学号专业班级实验项目实验五FIR数字滤波器的设计课程名称数字信号处理课程代码实验时间2013年05月31日实验地点主校区计算机基础实验室批改意见成绩教师签字：一，实验目的加深理解FIR数字滤波器的时域特性和频率特性，掌握FIR数字滤波器的设计原理与设计方法，以及FIR数字滤波器的应用。二，实验原理FIR数字滤波器总是稳定的系统，且可以设计成具有线性相位的。其在数据通信、图像处理、语言信号处理等实际应用领域中得到广泛的应用。N阶有限冲激响应（FIR）数字滤波器的转移函数为：，系统的单位脉冲响应h（n）是长度为N的有限长因果序列。当满足h(n)=h(N-n-1)的对称条件时，该FIR数字滤波器具有线性相位。FIR数字滤波器的设计方法主要有窗函数法和频率采样法。窗函数法FIR滤波器的冲激响应就是系统函数各次项的系数，所以设计FIR了不起的方法之一就是：从时域出发，截取有限长的一段冲激响应作为H（z）系数，冲激响应长度N就是系统函数H(z)阶数。只要N足够长，并且截取的方法合理，总能够满足频域的要求，这就是FIR滤波器的窗函数设计法。2、频率采样法频率采样法是先对理想频响进行采样，得到样值,再利用插值公式直接求出系统转换函数，以便实现；或者求出频响，以便与理想频响进行比较。在[0,2π]区间上对进行N点采样，等效于时域以N为周期延拓。三，实验内容与方法【例3-5-1】数字滤波器的技术指标如下：、、、，采用窗函数法来设计一个FIR数字滤波器。海明窗和布莱克曼窗均可以提供大于50dB的衰减。如果选用海明窗设计，则它提供较小的过渡带，因此，其具有较小的阶数。尽管在设计中用不到通带波动值，但必须检查设计的实际波动，即验证它是否确实在给定容限内。MATLAB程序如下：wp=0.4*pi;ws=0.6*pi;deltaw=ws-wp;N0=ceil(6.6*pi/deltaw)+1;N=N0+mod(N0+1,2);wdham=(hamming(N))';wc=(ws+wp)/2;tao=(N-1)/2;n=[0:N-1];m=n-tao+eps;hd=sin(wc*m)./(pi*m);h=hd.*wdham[H,w]=freqz(h,[1],1000,'whole');H=(H(1:1:501))';w=(w(1:1:501))';mag=abs(H);db=20*log10((mag+eps)./max(mag));pha=angle(H);grd=grpdelay(h,[1],w);dw=2*pi/1000;Rp=-(min(db(1:wp/dw+1)));As=-round(max(db(ws/dw+1:501)));n=0:N-1;subplot(2,2,1);stem(n,hd,'.');title('理想脉冲响应')axis([0N-1-0.20.3]);ylabel('hd(n)');text(N+1,-0.1,'n')subplot(2,2,2);stem(n,wdham,'.');title('海明窗')axis([0N-101.1]);ylabel('w(n)');text(N+1,0,'n')subplot(2,2,3);stem(n,h,'.');title('实际脉冲响应')axis([0N-1-0.20.3]);xlabel('n');ylabel('h(n)')subplot(2,2,4);plot(w/pi,db);title('幅度响应（单位：dB）');gridaxis([01-10010]);xlabel('频率(单位：pi)');ylabel('分贝')set(gca,'XTickMode','manual','XTick',[0,0.4,0.6,1]);set(gca,'YTickmode','manual','YTick',[-50,0])set(gca,'YTickMode','manual','YTickLabels',['50';'0'])set(gcf,'color','w');计算结果为：N=35，=0.04dB,=52dB,满足要求，并且对通带波动的验证显示为=0.04dB，它也是满足要求的。可以看出用海明窗设计的FIR数字滤波器是满足要求的，且时域和频域的曲线如下图所示。此外，在信号处理工具箱中，MATLAB还直接提供了一个子程序即firl，它利用窗函数法设计FIR滤波器，其标准调用格式为：B=firl(M,WN,‘type’,window)其中：b为待设计的滤波器系数向量，其长度为N=M+1；M为所选的滤波器阶数；为滤波器给定的边缘频率，可以是标量，也可以是一个数组；‘type’为滤波器的类型，如高通、带通、带阻等，缺省时为低通；window为选定的窗函数类型，缺省时为Hamming窗。同上例，MATLAB程序如下：b=firl(34,1/pi,hamming(35));[H,w]=freqz(b,1,512);H_db=20*log10(abs(H));subplot(2,1,1);stem(b);title('hamming');subplot(2,1,2);plot(w,H_db);title('frequency');gridon;得到的结果如下图所示。2、频率采样法【例3-5-2】同例3-5-1数字滤波器的技术指标如下：、、、，利用频率设计法设计一个FIR数字滤波器。MATLAB程序如下：N=35;wp=0.4*pi;ws=0.6*pi;wc=(wp+ws)/2;N=N+mod(N+1,2);N1=fix(wc/(2*pi/N));N2=N-2*N1-1;A=[ones(1,N1+1),zeros(1,N2),ones(1,N1)];theta=-pi*[0:N-1]*(N-1)/N;H=A.*exp(j*theta);h=real(ifft(H));wp1=2*pi/N*fix(wc/(2*pi/N)-1);ws1=wp1+8*pi/N;[H,w]=freqz(h,[1],1000,'whole');H=(H(1:1:501))';w=(w(1:1:501))';mag=abs(H);db=20*log10((mag+eps)./max(mag));pha=angle(H);grd=grpdelay(h,[1],w);N=length(h);L0=(N-1)/2;L=floor(L0);n=1:L+1;ww=[0:511]*pi/512;ifall(abs(h(n)-h(N-n+1))<1e-8)Ar=2*h(n)*cos(((N+1)/2-n)'*ww)-mod(N,2)*h(L+1);type=2-mod(N,2);elseifall(abs(h(n)+h(N-n+1))<1e-8)&(h(L+1)*mod(N-2)<=1e-8)A=2*h(n)*sin(((N+1)/2-n)'*ww);type=4-mod(N,2);elseerror('错误：这不是线性相位滤波器！')enddw=2*pi/1000;Rp=-min(db(1:fix(wp1/dw)+1));As=-round(max(db(fix(ws1/dw)+1:501)));l=0:N-1;wl=(2*pi/N)*l;wdl=[0,wc,wc,2*pi-wc,2*pi-wc,2*pi]/pi;Adl=[1,1,0,0,1,1];subplot(2,2,1);plot(wl(1:N)/pi,A(1:N),'.',wdl,Adl);axis([0,1,-0.1,1.1]);title('频率样本')xlabel('');ylabel('A(k)');set(gca,'XTickMode','manual','XTick',chop([0,0.5,1],2));set(gca,'YTickMode','manual','YTick',[0,1]);gridsubplot(2,2,2);stem(l,h,'.');axis([-1,N,-0.1,0.5]);gridtitle('脉冲响应');ylabel('h(n)');text(N+1,-0.1,'n')subplot(2,2,3);plot(ww/pi,Ar,wl(1:N)/pi,A(1:N),'.');axis([0,1,-0.2,1.2]);title('幅频响应')xlabel('频率（单位：pi）');ylabel('Ar(w)');gridsubplot(2,2,4);plot(w/pi,db);axis([0,1,-50,10]);title('幅度响应');xlabel('频率（单位：pi）');ylabel('分贝数');gridset(gca,'XTickMode','manual','XTick',chop([0,0.5,1],2));set(gca,'YTickMode','manual','YTick',[-As,0]);set(gca,'YTickLabelMode','manual','YTickLabels',['As';'0']);set(gcf,'color','w');经过计算可以得到：=1.6011dB、=24dB,所设计的滤波器的时域与频域特性如下图所示。可以看出，其不满足设计指标。为此，可以在构造的频率响应间断点处增加一个或若干了过度点来改善带外衰减指标。在此选择增加两个过度点、来改善滤波器的设计。MATLAB程序如下：N=35;wp=0.4*pi;ws=0.6*pi;wc=(wp+ws)/2;N=N+mod(N+1,2);N1=fix(wc/(2*pi/N));N2=N-2*N1-1;T1=input('过渡带第一样本T1的值');T2=input('过渡带第二样本T2的值');A=[ones(1,N1),T1,T2,zeros(1,N2-2),T2,T1,ones(1,N1-1)];A=[ones(1,N1+1),zeros(1,N2),ones(1,N1)];theta=-pi*[0:N-1]*(N-1)/N;H=A.*exp(j*theta);h=real(ifft(H));wp1=2*pi/N*fix(wc/(2*pi/N)-1);ws1=wp1+8*pi/N;[H,w]=freqz(h,[1],1000,'whole');H=(H(1:1:501))';w=(w(1:1:501))';mag=abs(H);db=20*log10((mag+eps)./max(mag));pha=angle(H);grd=grpdelay(h,[1],w);N=length(h);L0=(N-1)/2;L=floor(L0);n=1:L+1;ww=[0:511]*pi/512;ifall(abs(h(n)-h(N-n+1))<1e-8)Ar=2*h(n)*cos(((N+1)/2-n)'*ww)-mod(N,2)*h(L+1);type=2-mod(N,2);elseifall(abs(h(n)+h(N-n+1))<1e-8)&(h(L+1)*mod(N-2)<=1e-8)A=2*h(n)*sin(((N+1)/2-n)'*ww);type=4-mod(N,2);elseerror('错误：这不是线性相位滤波器！')enddw=2*pi/1000;Rp=-min(db(1:fix(wp1/dw)+1));As=-round(max(db(fix(ws1/dw)+1:501)));l=0:N-1;wl=(2*pi/N)*l;wdl=[0,wc,wc,2*pi-wc,2*pi-wc,2*pi]/pi;Adl=[1,1,0,0,1,1];subplot(2,2,1);plot(wl(1:N)/pi,A(1:N),'.',wdl,Adl);axis([0,1,-0.1,1.1]);title('频率样本')xlabel('');ylabel('A(k)');set(gca,'XTickMode','manual','XTick',chop([0,0.5,1],2));set(gca,'YTickMode','manual','YTick',[0,1]);gridsubplot(2,2,2);stem(l,h,'.');axis([-1,N,-0.1,0.5]);gridtitle('脉冲响应');ylabel('h(n)');text(N+1,-0.1,'n')subplot(2,2,3);plot(ww/pi,Ar,wl(1:N)/pi,A(1:N),'.');axis([0,1,-0.2,1.2]);title('幅频响应')xlabel('频率（单位：pi）');ylabel('Ar(w)');gridsubplot(2,2,4);plot(w/pi,db);axis([0,1,-50,10]);title('幅度响应');xlabel('频率（单位：pi）');ylabel('分贝数');gridset(gca,'XTickMode','manual','XTick',chop([0,0.5,1],2));set(gca,'YTickMode','manual','YTick',[-As,0]);set(gca,'YTickLabelMode','manual','YTickLabels',['As';'0']);set(gcf,'color','w');运行程序画出幅频特性曲线如图所示。同窗函数一样，MATLAB信号处理工具箱中提供了频率采样法的设计函数fir2。它的典型调用法为：b=fir2(M,f,m)其次，b为待设计的滤波器系数向量，其长度为N=M+1；M为所选的滤波器阶数；f制定归一化的各频率边界频率，从0到1递增，1对应采样频率的一半，及奎奈斯特频率；m为制定各边界处的幅度响应，因此f和m的长度相等。f=[01/pi1/pi1];m=[1100];b=fir2(34,f,m);[H,w]=freqz(b,1,128);H_db=20*log10(abs(H));subplot(2,1,1);stem(b,'.');title('实际脉冲响应');subplot(2,1,2);plot(w/pi,H_db);title('幅度响应');gridonset(gcf,'color','w');得到结果如图所示。四、程序设计实验数字滤波器的技术指标如下：、、、，分别采用窗函数法和频率采样法设计一个FIR数字滤波器。采用窗函数法MATLAB程序如下：wp=0.2*pi;ws=0.3*pi;deltaw=ws-wp;N0=ceil(6.6*pi/deltaw)+1;N=N0+mod(N0+1,2);wdham=(hamming(N))';wc=(ws+wp)/2;tao=(N-1)/2;n=[0:N-1];m=n-tao+eps;hd=sin(wc*m)./(pi*m);h=hd.*wdham[H,w]=freqz(h,[1],1000,'whole');H=(H(1:1:501))';w=(w(1:1:501))';mag=abs(H);db=20*log10((mag+eps)./max(mag));pha=angle(H);grd=grpdelay(h,[1],w);dw=2*pi/1000;Rp=-(min(db(1:wp/dw+1)));As=-round(max(db(ws/dw+1:501)));n=0:N-1;subplot(2,2,1);stem(n,hd,'.');title('理想脉冲响应')axis([0N-1-0.20.3]);ylabel('hd(n)');text(N+1,-0.1,'n')subplot(2,2,2);stem(n,wdham,'.');title('海明窗')axis([0N-101.1]);ylabel('w(n)');text(N+1,0,'n')subplot(2,2,3);stem(n,h,'.');title('实际脉冲响应')axis([0N-1-0.20.3]);xlabel('n');ylabel('h(n)')subplot(2,2,4);plot(w/pi,db);title('幅度响应（单位：dB）');gridaxis([01-10010]);xlabel('频率(单位：pi)');ylabel('分贝')set(gca,'XTickMode','manual','XTick',[0,0.4,0.6,1]);set(gca,'YTickmode','manual','YTick',[-50,0])set(gca,'YTickMode','manual','YTickLabels',['50';'0'])set(gcf,'color','w');程序运行结果如图采用频率采样法MATLAB程序如下：N=35;wp=0.2*pi;ws=0.3*pi;wc=(wp+ws)/2;N=N+mod(N+1,2);N1=fix(wc/(2*pi/N));N2=N-2*N1-1;A=[ones(1,N1+1),zeros(1,N2),ones(1,N1)];theta=-pi*[0:N-1]*(N-1)/N;H=A.*exp(j*theta);h=real(ifft(H));wp1=2*pi/N*fix(wc/(2*pi/N)-1);ws1=wp1+8*pi/N;[H,w]=freqz(h,[1],1000,'whole');H=(H(1:1:501))';w=(w(1:1:501))';mag=abs(H);db=20*log10((mag+eps)./max(mag));pha=angle(H);grd=grpdelay(h,[1],w);N=length(h);L0=(N-1)/2;L=floor(L0);n=1:L+1;ww=[0:511]*pi/512;ifall(abs(h(n)-h(N-n+1))<1e-8)Ar=2*h(n)*cos(((N+1)/2-n)'*ww)-mod(N,2)*h(L+1);type=2-mod(N,2);elseifall(abs(h(n)+h(N-n+1))<1e-8)&(h(L+1)*mod(N-2)<=1e-8)A=2*h(n)*sin(((N+1)/2-n)'*ww);type=4-mod(N,2);elseerror('')enddw=2*pi/1000;Rp=-min(db(1:fix(wp1/dw)+1));As=-round(max(db(fix(ws1/dw)+1:501)));l=0:N-1;wl=(2*pi/N)*l;wdl=[0,wc,wc,2*pi-wc,2*pi-wc,2*pi]/pi;Adl=[1,1,0,0,1,1];subplot(2,2,1);plot(wl(1:N)/pi,A(1:N),'.',wdl,Adl);axis([0,1,-0.1,1.1]);title('频率样本')xlabel('');ylabel('A(k)');set(gca,'XTickMode','manual','XTick',chop([0,0.5,1],2));set(gca,'YTickMode','manual','YTick',[0,1]);gridsubplot(2,2,2);stem(l,h,'.');axis([-1,N,-0.1,0.5]);gridtitle('脉冲响应');ylabel('h(n)');text(N+1,-0.1,'n')subplot(2,2,3);plot(ww/pi,Ar,wl(1:N)/pi,A(1:N),'.');axis([0,1,-0.2,1.2]);title('幅频响应')xlabel('频率（单位：pi）');ylabel('Ar(w)');gridsubplot(2,2,4);plot(w/pi,