成绩：《数字信号处理》作业与上机实验（第七章）班级：电信学号：姓名：任课老师：李宏民完成时间：信息与通信工程学院2015—2016学年第1 学期第7章 有限脉冲响应数字滤波器设计一、教材p238：19,20,21,25,26二、某信号()x t 为：123()0.5cos(2)0.7cos(20.1)0.4cos(2)x t f t f t f t ππππ=+++，其中121100,130,600.f Hz f Hz f Hz ===设计最低阶FIR 数字滤波器，按下图所示对()x t 进行数字滤波处理，实现：(x t ()y t 1）将3f 频率分量以高于50dB 的衰减抑制，同时以低于2dB 的衰减通过1f 和2f 频率分量；2）将1f 和2f 频率分量以高于50dB 的衰减抑制，同时以低于2dB 的衰减通过3f 频率分量；要求：按数字滤波器直接型结构图编写滤波程序，求得()y n ；1)中的FIR 滤波器采用窗函数法设计；2)中的FIR 滤波器采用频率采样法设计。

画出所设计的滤波器频率特性图、信号时域图；给出滤波器设计的MATLAB 代码与滤波器实现的代码；选择合适的信号采样周期T 。

3）与第6章作业2的IIR 滤波方法进行比较研究。

一、19、 Fs=80000; fp=15000;fs=20000;rs=40;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs; rp=-20*log10(1-0.02);rs=40; [N1,wpo]=ellipord(wp/pi,ws/pi,rp,rs); [B,A]=ellip(N1,rp,rs,wpo); [Hk,wk]=freqz(B,A,500);Bt=ws-wp;alph=0.5842*(rs-21)^0.4+0.07886*(rs-21); M=ceil((rs-8)/2.285/Bt) wc=(wp+ws)/2/pi;hn=fir1(M,wc,kaiser(M+1,alph)); [Hk1,wk1]=freqz(hn,1,500); figure(1);plot(wk1/pi,20*log10(abs(Hk1)),'k'); hold on plot(wk/pi,20*log10(abs(Hk)),'r--'); hold off legend('FIR 滤波器,'IIR 滤波器');axis([0,1,-80,5]);xlabel('w/\pi');ylabel('幅度/dB'); title('损耗函数'); figure(2)plot(wk1/pi,angle(Hk1)/pi,'k'); hold on plot(wk/pi,angle(Hk)/pi,'r--'); hold off legend('FIR 滤波器','IIR 滤波器');xlabel('w/\pi');ylabel('相位/\pi'); title('相频特性曲线');0.20.40.60.810w/幅度/d B损耗函数0.20.40.60.81-1-0.8-0.6-0.4-0.200.20.40.60.81w/π相位/π相频特性曲线20、N=21;n=0:20;wc=pi/4;hn1=fir1(N-1,wc,'s',hanning(N)); hn2=fir1(N-1,wc,'s',hamming(N)); hn3=fir1(N-1,wc,'s',boxcar(N)); hn4=fir1(N-1,wc,'s',blackman(N)); figure(1)plot(n,hn1,'*b');hold on ;plot(n,hn2,'--','linewidth',2); plot(n,hn3,'r:','linewidth',3); plot(n,hn4);hold off ;xlabel('n');ylabel('h(n)');legend('汉宁窗','哈明窗','矩形窗','布莱克曼窗'); title('单位冲击响应'); figure(2)[Hk1,wk1]=freqz(hn1,1,500);plot(wk1/pi,20*log10(abs(Hk1)),'*b');hold on [Hk2,wk2]=freqz(hn2,1,500);plot(wk2/pi,20*log10(abs(Hk2)),'--','linewidth',2); [Hk3,wk3]=freqz(hn3,1,500);plot(wk3/pi,20*log10(abs(Hk3)),'r:','linewidth',3);[Hk4,wk4]=freqz(hn4,1,500);plot(wk4/pi,20*log10(abs(Hk4)));hold offlegend('汉宁窗','哈明窗','矩形窗','布莱克曼窗'); axis([0,1,-80,30]);xlabel('w/\pi');ylabel('幅度'); title('四种低通滤波器的损耗函数');5101520nh (n )单位冲击响应0.20.40.60.81w/幅度四种低通滤波器的损耗函数21、N=21;n=0:20;wc=pi/4;hn1=fir1(N-1,wc,'high',hanning(N)); hn2=fir1(N-1,wc,'high',hamming(N)); hn3=fir1(N-1,wc,'high',boxcar(N)); hn4=fir1(N-1,wc,'high',blackman(N)); figure(1)plot(n,hn1,'*b');hold on ;plot(n,hn2,'--','linewidth',2); plot(n,hn3,'r:','linewidth',3); plot(n,hn4);hold off ;xlabel('n');ylabel('h(n)');legend('汉宁窗','哈明窗','矩形窗','布莱克曼窗'); title('单位冲击响应'); figure(2)[Hk1,wk1]=freqz(hn1,1,500);plot(wk1/pi,20*log10(abs(Hk1)),'*b');hold on [Hk2,wk2]=freqz(hn2,1,500);plot(wk2/pi,20*log10(abs(Hk2)),'--','linewidth',2); [Hk3,wk3]=freqz(hn3,1,500);plot(wk3/pi,20*log10(abs(Hk3)),'r:','linewidth',3); [Hk4,wk4]=freqz(hn4,1,500);plot(wk4/pi,20*log10(abs(Hk4)));hold offlegend('汉宁窗','哈明窗','矩形窗','布莱克曼窗'); axis([0,1,-80,30]);xlabel('w/\pi'); ylabel('幅度'); title('四种窗的损耗函数');5101520nh (n )单位冲击响应0102030w/幅度四种窗的损耗函数25、wp=0.6*pi;ws=0.45*pi;rp=0.2;rs=45;Bt=wp-ws;N1=ceil(6.2*pi/Bt); N11=N1+mod(N1+1,2); disp('N11='),disp(N11) N2=ceil(11*pi/Bt); N22=N2+mod(N2+1,2); disp('N22='),disp(N22),alph=0.5842*(rs-21)^0.4+0.07886*(rs-21); M=ceil((rs-8)/2.285/Bt); disp('M='),disp(M) wc=(wp+ws)/2/pi;hn1=fir1(N11-1,wc,'high',hamming(N11)); [Hk1,w]=freqz(hn1,1);hn2=fir1(N22-1,wc,'high',blackman(N22)); [Hk2,w]=freqz(hn2,1);hn3=fir1(M-1,wc,'high',kaiser(M,alph)); [Hk3,w]=freqz(hn3,1); figure(1)plot(w/pi,20*log10(Hk1),w/pi,20*log10(Hk2),'-.',w/pi,20*log10(Hk3),'--');grid onxlabel('w/\pi');ylabel('幅度');legend('哈明窗', '布莱克曼窗', '凯塞窗') title('三种窗FIR 高通损耗函数曲线') figure(2)stem(hn1)title('哈明窗单位脉冲响应h(n)') figure(3) stem(hn2)title('布莱克曼窗单位脉冲响应h(n)') figure(4) stem(hn3)title('凯塞窗单位脉冲响应h(n)')0.20.40.60.81-160-140-120-100020w/幅度三种窗FIR 高通损耗函数曲线51015202530354045-0.4-0.3-0.2-0.100.10.20.30.40.5哈明窗单位脉冲响应h(n)1020304050607080布莱克曼窗单位脉冲响应h(n)5101520253035-0.4-0.3-0.2-0.100.10.20.30.40.5凯塞窗单位脉冲响应h(n)26、wp1=0.55*pi;wp2=0.7*pi;ws1=0.45*pi;ws2=0.8*pi;Bt=wp2-wp1; N=ceil(6.2*pi/Bt);wc=[(wp1+ws1)/2/pi,(ws2+wp2)/2/pi]; hn=fir1(N-1,wc,hanning(N)); [Hk1,wk1]=freqz(hn,1,500); n=0:N-1; figure(1); stem(n,hn);xlabel('n');ylabel('h(n)'); title('单位冲击响应'); figure(2);plot(wk1/pi,20*log10(abs(Hk1))); axis([0,1,-100,5]);xlabel('w/\pi');ylabel('·幅度'); title('损耗函数');grid on051015202530354045nh (n )单位冲击响应-1000w/?幅度损耗函数二、(1) Fs=3500; fp=130; fs=600; rs=50;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs; Bt=ws-wp;alph=0.5842*(rs-21)^0.4+0.07886*(rs-21); N=ceil((rs-8)/2.285/Bt) wc=(wp+ws)/2/pi;hn=fir1(N,wc,kaiser(N+1,alph)) figure(1) freqz(hn,1); N=500;n=0:N-1; T=1/Fs;t=n*T;x=0.5*cos(2*pi*100*t)+0.7*cos(2*pi*130*t+0.1*pi)+0.4*cos(2*pi*600*t); figure(2),plot(n,x);title('信号x(n)');ylabel('·幅值');xlabel('n');m1=0;m2=0;m3=0;m4=0;m5=0;m6=0;m7=0;m8=0;m9=0;m10=0;m11=0;m12=0;m13=0;m14=0; m15=0;m16=0;m17=0;m18=0;m19=0;m20=0;m21=0;m22=0; for m=1:500y(m)=hn(1)*x(m)+m1*hn(2)+m2*hn(3)+m3*hn(4)+m4*hn(5)+m5*hn(6)+...m6*hn(7)+m7*hn(8)+m8*hn(9)+m9*hn(10)+m10*hn(11)+m11*hn(12)+m12*hn(13)+...m13*hn(14)+m14*hn(15)+m15*hn(16)+m16*hn(17)+m17*hn(18)+m18*hn(19)+... m19*hn(20)+m20*hn(21)+m21*hn(22)+m22*hn(23);m22=m21;m21=m20;m20=m19;m19=m18;m18=m17;m17=m16;m16=m15;m15=m14;m14=m13;m13=m12;m12=m11;m11=m10;m10=m9;m9=m8;m8=m7;m7=m6;m6=m5; m5=m4;m4=m3;m3=m2;m2=m1;m1=x(m);end figure(3) plot(n,y); xlabel('n'); ylabel('y(n)');title('直接网络型y(n)');-800-600-400-2000Normalized Frequency (⨯π rad/sample)P h a s e (d e g r e e s )-150-100050Normalized Frequency (⨯π rad/sample)M a g n i t u d e (d B )0100200300400500012信号x(n)?幅值n01002003004005001ny (n )直接网络型y(n)(2) T=0.3; Fs=3500;fp=600;fs=100;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs; Bt=wp-ws;m=1;N=ceil((m+1)*2*pi/Bt+1)N=N+mod(N+1,2);Np=fix(wp/(2*pi/N));Ns=N-2*Np-1;Hw=[zeros(1,Np+1),ones(1,Ns),zeros(1,Np)]; Hw(Np+2)=T;Ak(N-Np)=T; thetak=-pi*(N-1)*(0:N-1)/N; Hdk=Hw.*exp(1j*thetak); hn=real(ifft(Hdk)) figure(1) freqz(hn,1); N=500;n=0:N-1; T=1/Fs;t=n*T;x=0.5*cos(2*pi*100*t)+0.7*cos(2*pi*130*t+0.1*pi)+0.4*cos(2*pi*600*t); figure(2),plot(n,x);title('x(n)');ylabel('·幅度');xlabel('n');m1=0;m2=0;m3=0;m4=0;m5=0;m6=0;m7=0;m8=0;m9=0;m10=0;m11=0;m12=0;m13=0;m14=0; for m=1:500y(m)=hn(1)*x(m)+m1*hn(2)+m2*hn(3)+m3*hn(4)+m4*hn(5)+m5*hn(6)+...m6*hn(7)+m7*hn(8)+m8*hn(9)+m9*hn(10)+m10*hn(11)+m11*hn(12)+m12*hn(13)+... m13*hn(14)+m14*hn(15);m14=m13;m13=m12;m12=m11;m11=m10;m10=m9;m9=m8;m8=m7;m7=m6;m6=m5; m5=m4;m4=m3;m3=m2;m2=m1;m1=x(m);end figure(3) plot(n,y); xlabel('n'); ylabel('y(n)');title('直接网络型y(n)');00.20.40.60.81-1000-5000500Normalized Frequency (⨯π rad/sample)P h a s e (d e g r e e s )0.20.40.60.81-400-300-200-1000100Normalized Frequency (⨯π rad/sample)M a g n i t u d e (d B )0100200300400500012x(n)?幅度n01002003004005000ny (n )直接网络型y(n)(3) T=1;Fs=3500;fp=130;fs=600;wp1=(2*pi*fp)/Fs/T;ws1=(2*pi*fs)/Fs/T; rp=2;rs=50;[N,wc]=ellipord(wp1,ws1,rp,rs,'s') [B,A]=ellip(N,rp,rs,wc,'s'); [Bz,Az]=impinvar(B,A,1/T) [Hk,wk]=freqz(Bz,Az,500);wp=(2*pi*fp)/Fs;ws=(2*pi*fs)/Fs; Bt=ws-wp;alph=0.5842*(rs-21)^0.4+0.07886*(rs-21); N=ceil((rs-8)/2.285/Bt) wc=(wp+ws)/2/pi;hn=fir1(N,wc,kaiser(N+1,alph)); [Hk1,wk1]=freqz(hn,1,500); figure(1);plot(wk1/pi,20*log10(abs(Hk1)),'k'); hold onplot(wk/pi,20*log10(abs(Hk)),'r--','linewidth',2); hold off legend('FIR 滤波器','IIR 滤波器');axis([0,1,-80,5]);xlabel('w/\pi');ylabel('·幅度'); title('损耗函数'); figure(2)plot(wk1/pi,angle(Hk1)/pi,'k'); hold onplot(wk/pi,angle(Hk)/pi,'--','linewidth',2); hold off legend('FIR 滤波器','IIR 滤波器');xlabel('w/\pi');ylabel('相位/\pi'); title('相频特性曲线');w/π?幅度损耗函数01w/π相位/π相频特性曲线。