北理工大学数字信号处理实验报告.docx

实验1 利用DFT分析信号频谱一、实验目的1.加深对DFT原理的理解。2.应用DFT分析信号的频谱。3.深刻理解利用DFT分析信号频谱的原理，分析实现过程中出现的现象及解决方法。二、实验设备与环境计算机、MATLAB软件环境。三、实验原理1. DFT与 DTFT的关系有限长序列x(n)(0 nN-1) 的离散时间傅里叶变换X (ej) 在频率区间(02)的N个等间隔分布的点k=2k/N(0kN-1)上的N个取样值可以由下式表示：Xej=2k/N=k=0N-1xne-j2Nkn=X(k) 0kN-1由上式可知，序列x(n)的N 点DFT X(k)，实际上就是x(n)序列的DTFT在N 个等间隔频率点k=2k/N(0kN-1)上样本X(k)。2. 利用DFT 求DTFT法一：由 X (k)恢复出 X (ej )的方法：Xej=n=-xne-jnxej=k=1NXk(-2N)=sinN2Nsin(/2)法二：然而在实际 MATLAB计算中，上述插值运算不见得是最好的办法。由于DFT是 DTFT的取样值，其相邻两个频率样本点的间距为2/N，所以如果我们增加数据的长度N，使得到的DFT谱线就更加精细，其包络就越接近DTFT的结果，这样就可以利用DFT来近似计算DTFT。如果没有更多的数据，可以通过补零来增加数据长度。3.利用DFT分析连续时间信号的频谱采用计算机分析连续时间信号的频谱，第一步就是把连续时间信号离散化，这里需要进行两个操作：一是采样，二是截断。将利用DFT分析连续非周期信号频谱的步骤归纳如下：（1）确定时域采样间隔T，得到离散序列x(n)；（2）确定截取长度M，得到M点离散序列xMn=xnw(n)，这里w(n)为窗函数。（3）确定频域采样点数 N，要求NM。（4）利用FFT计算离散序列的 N点DFT，得到xMn。（5）根据式（2-6）由xMn计算 Xa(j)采样点的近似值。采用上述方法计算xa(t)的频谱，需要注意如下三个问题：（1）频谱混叠。如果不满足采样定理的条件，频谱会出现混叠误差。对于频谱无限宽的信号，应考虑覆盖大部分主要频率分量的范围。（2）栅栏效应和频谱分辨率。使用DFT计算频谱，得到的结果只是 N个频谱本值，样本值之间的频谱是未知的，像通过一个栅栏观察频谱，称为“栅栏效应”。频谱分辨率与记录长度成反比，要提高频谱分表率，就要增加记录时间。（3）频谱泄露。对信号截断会把窗函数的频谱引入信号频谱，照成频谱泄露 。解决这个问题的主要办法是采用旁瓣小的窗函数，频谱泄露和窗函数均会引起误差 。因此，要合理选取采样间隔和截取长度，必要时还需考虑加适当的窗。对于连续时间周期信号，我们在采用计算机进行计算时，也总是要进行截断，序列总是有限长的，仍然可以采用上述方法近似计算。四、实验内容1. 已知 x(n)=1 , 1 , 1 , 2，完成如下要求：（1）计算其 DTFT，并画出-pi,pi区间的波形。源代码：n=0:3;x=2 -1 1 1;w=-pi:0.01*pi:pi;X=x*exp(-1j*n*w);subplot(211);plot(w,abs(X);xlabel(Omega/pi);title(Magnitude);axis(-pi pi 0 4);subplot(212);plot(w,angle(X)/pi);xlabel(Omega/pi);axis(-pi pi -1 1);title(Phase); 实验结果：（2）计算 4 点 DFT，并把结果显示在（1）所画的图形中。源代码：n=0:3;x=2 -1 1 1;w=-pi:0.01*pi:pi;X=x*exp(-1j*n*w);subplot(211);plot(w,abs(X);xlabel(Omega/pi);title(Magnitude);axis(-pi pi 0 4);hold on;w1=-pi:0.5*pi:0.5*pi;X1=x*exp(-1j*n*w1);stem(w1,abs(X1),filled);subplot(212);plot(w,angle(X)/pi);hold on;stem(w1,angle(X1)/pi,filled);xlabel(Omega/pi);axis(-pi pi -1 4);title(Phase);实验结果：（3）对 x(n)补零，计算64点DFT，并显示结果。源代码：n=0:3;n1=0:63;x=2 -1 1 1;x1=2 -1 1 1 zeros(1,60);w=-pi:0.01*pi:pi;w1=-pi:pi/32:pi-pi/32;X=x*exp(-1j*n*w);X1=x1*exp(-1j*n1*w1);subplot(211);plot(w,abs(X);xlabel(Omega/pi);title(Magnitude);axis(-pi pi 0 4);hold on;stem(w1,abs(X1),filled);subplot(212);plot(w,angle(X)/pi);hold on;stem(w1,angle(X1)/pi,filled);xlabel(Omega/pi);axis(-pi pi -1 4);title(Phase);（4）根据实验结果，分析是否可以由 DFT计算 DTFT，如果可以，如何实现。 可以用DFT计算DTFT，由上图看出，N值越大，DFT的包络就越像DTFT的曲线图。所以，我们可以很直接的得到这个答案，那就是，当我们取得N足够大的时候，我们便可以用DFT来表示DTFT了。2. 考察序列xn=cos0.48n+cos(0.52n)（1）0n10 时，用DFT 估计x(n)的频谱；将x(n)补零加长到长度为100 点序列用DFT 估计x(n)的频谱。要求画出相应波形。 源代码：n1=0:10;x1=(cos(0.48*pi.*n1)+cos(0.52*pi.*n1);X1=fft(x1);subplot(211);xlabel(Omega/pi);title(Magnitude);stem(n1,abs(X1),filled);n=0:100;ni = 0:9;x=(cos(0.48*pi.*ni)+cos(0.52*pi.*ni), zeros(1,91);X=fft(x);subplot(212);xlabel(Omega/pi);title(Magnitude);stem(n/10,abs(X),filled);（2）0n100 时，用DFT 估计x(n)的频谱，并画出波形。 源代码：n1=0:100;x1=cos(0.48*pi.*n1)+cos(0.52*pi.*n1);X1=fft(x1);xlabel(Omega/pi);title(Magnitude);stem(n1/100*2*pi,abs(X1),filled);（3）根据实验结果，分析怎样提高频谱分辨率。加大n，n加的足够大的时候，信号的频谱看起来更像是原来实际信号的频谱。3. 已知信号xt=0.15sin2f1t+sin(2f2t) 0.1sin(2f3t) ，其中f1=1Hz，f2=2Hz ， f3=3Hz 。从x(t)的表达式可以看出，它包含三个频率的正弦波，但是，从其时域波形（图E2-1）来看，似乎是一个正弦信号，利用DFT 做频谱分析，确定适合的参数，使得到的频谱的频率分辨率符合需要。源代码：n=0:0.1:4.9;x=0.15*sin(2*pi*n)+sin(4*pi*n)-0.1*sin(6*pi*n);X=fft(x);stem(n,abs(X),filled);4.利用DFT 近似分析连续时间信号 xt=e0.1t u(t) 的频谱（幅度谱）。分析采用不同的采样间隔和截取长度进行计算的结果，并最终确定适合的参数。源代码：t=0:0.1:99.9;w=-pi:0.01:pi;x=exp(0.1*t);X=x*exp(-1j*t*w);plot(w,abs(X);五、 实验中的主要结论、遇到的问题及解决方法，收获和体会。主要结论：1）可以用DFT来替代DTFT，但是有个条件，就是，在DFT里面去的值要足够的大，DFT的包络才可以代表DTFT的曲线。2）当时域信号离散的时候，我们可以加大时间的取样，并且时间取得越长，离散信号的频谱就更接近真实连续信号的频谱。遇到的问题：在用MATLAB做DFT的时候，一度时间做不出来，信号出来的时候相位有延时，后来发现是时间取样不对。收获和体会：用MATLAB实实在在地看到了DFT和DTFT的区别和联系，对课本上的知识不再仅仅只是只能想象的阶段了，我们还可以，用计算机软件来仿真，这给了我们一种思路，一种从其他方面解决问题的方法。实验三 IIR数字滤波器设计一、实验目的1. 掌握利用脉冲响应不变法和双线性变换法设计IIR数字滤波器的原理及具体方法。2. 加深理解数字滤波器和模拟滤波器之间的技术指标转化。3. 掌握脉冲响应不变法和双线性变换法设计IIR数字滤波器的优缺点及使用范围。二、实验设备与环境计算机、MATLAB软件环境三、实验原理 IIR滤波器具有无限长持续时间脉冲响应，而模拟滤波器一般都具有无限长的脉冲响应，因此它与模拟滤波器相匹配。IIR滤波器设计的基本方法就是先设计一个合适的模拟滤波器，然后利用复值映射把模拟滤波器变换成数字滤波器。 1.数字滤波器和模拟滤波器的一些指标通带、过渡带、通带响应中的容限、阻带的容限、通带波动、阻带衰减等等。2.模拟原型滤波器IIR滤波器设计方法由已有的模拟滤波器得到数字滤波器，我们将这些模拟滤波器称为原型滤波器。常用的模拟原型滤波器有巴特沃斯（Butterworth）滤波器、切比雪夫(Chebyshev，型和型)滤波器和椭圆（Ellipse）滤波器等。(1) 巴特沃斯滤波器(2) 切比雪夫低通滤波器a、切比雪夫I低通滤波器b、切比雪夫II低通滤波器(3) 椭圆滤波器(4)三种滤波器的比较如果给定相同的设计指标，选用椭圆滤波器所要求的阶数 N最低，切比雪夫滤波器次之，巴特沃斯滤波器最高；如果要求的阶数相同，切比雪夫滤波器的过渡带比巴特沃斯滤波器陡，椭圆滤波器的过渡带又比切比雪夫滤波器陡。然而，从通带的相位响应来看，椭圆滤波器虽然提供了最优的幅度平方响应，但通带上的相位响应非线性较大，而巴特沃斯滤波器在通带上具有相当的线性相位，切比雪夫滤波器的相位特征介于两者之间。所以在实际设计中，选用何种滤波器应视实际用途和指标要求而定。3模拟滤波器到数字滤波器的变换从模拟滤波器到数字滤波器的变换就是要由 H a (s)进一步求得 H(z)，也就是由 s平面到 z平面的变换，这种变换要求数字滤波器能模仿模拟滤波器的特性。我们最常用的从模拟滤波器到数字滤波器的变换方法有两种：脉冲响应不变法和双线性变换法。3.1脉冲响应不变法基本原理:从时域响应出发，使数字滤波器的单位脉冲响应 h(n)模仿模拟滤波器的单位冲激响应ha(t)，h(n)等于ha(t)的取样值。3.2双线性变换法基本原理:从频率响应出发，直接使数字滤波器的频率响应逼近模拟滤波器的频率响应，进而求得数字滤波器的系统函数 H(z)。4. 直接利用 MATLAB函数设计 IIR数字滤波器前面我们介绍了 butter函数、cheby1函数、cheby2函数、ellip函数设计模拟滤波器的方法，这些函数也可直接用于设计 IIR数字滤波器，由于双线性变换法的优点，这些函数都采用双线性变换。这些函数的具体使用方法如下：b,a = butter(N,Wn)设计低通巴特沃斯数字滤波器b,a = cheby1(N,Rp,Wn)设计低通切比雪夫型数字滤波器b,a = cheby2(N,As,Wn,)设计低通切比雪夫型数字滤波器b,a = ellip(N,Rp,As,Wn)设计低通椭圆数字滤波器b,a = butter(N,Wn,ftype)设计巴特沃斯数字滤波器b,a = cheby1(N,Rp,Wn,ftype)设计切比雪夫型数字滤波器b,a = cheby2(N,As,Wn,ftype)设计切比雪夫型数字滤波器b,a = ellip(N,Rp,As,Wn,ftype)设计椭圆数字滤波器ftype表示数字滤波器的类型，可选选项有high、low、stop，分别表示高通、低通、带阻。上述函数中的参数N和Wn可以采用函数buttord、函数cheb1ord、函数cheb2ord、函数ellipord求得，这些函数能在已知设计指标时给出滤波器的阶数N和截止频率Wn。这些函数的具体使用方法如下：N,Wn = buttord(Wp,Wst,Rp,As)根据数字滤波器设计指标给出数字巴特沃斯滤波器的阶数 N和截止频率Wn；N,Wn = cheb1ord(Wp,Wst,Rp,As) 根据数字滤波器设计指标给出数字切比雪夫型滤波器的阶数N和截止频率Wn；N,Wn = cheb2ord(Wp,Wst,Rp,As) 根据数字滤波器设计指标给出数字切比雪夫型滤波器的阶数N和截止频率Wn；N,Wn = ellipord(Wp,Wst,Rp,As)根据数字滤波器设计指标给出数字椭圆滤波器的阶数 N和截止频率Wn。其中，Wp和 Wst分别表示通带和阻带截止频率，Rp表示通带波动，As表示阻带衰减。四、实验内容设频率为f=10kHz，设计数字低通滤波器，要求：通带截止频率：fp=1kHz，通带波动：Rp=1dB阻带截止频率：fst=1.5kHz，阻带衰减：As=15dB要求分别设计巴特沃斯、切比雪夫I型、切比雪夫II型和椭圆模拟原型滤波器，并分别结合脉冲响应不变法和双线性变换法进行设计。结合实验结果，分别讨论采用上述方法设计的数字滤波器是否都能满足给定指标要求，分析脉冲响应不变法和双线性变换法设计IIR数字滤波器的优缺点及使用范围。（1） 巴特沃斯wp=2*pi*1000; ws=2*pi*1500; Rp=1; As=15; Fs=10000; OmegaP=wp;OmegaS=ws; ep=sqrt(10(Rp/10)-1); Ripple=sqrt(1/(1+ep*ep); Attn=1/(10(As/20); N,OmegaC=buttord(OmegaP,OmegaS,Rp,As,s);z0,p0,k0=buttap(N); p=p0*OmegaC;z=z0*OmegaC; k=k0*OmegaCN; ba=k*real(poly(z); aa=real(poly(p); bd,ad=impinvar(ba,aa,Fs);H,w=freqz(bd,ad,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(bd,ad,w); subplot(221);plot(w/pi,mag); title(A); ylabel(|H|);grid on; subplot(2,2,3);plot(w/pi,db); title(A(dB); xlabel(f(pi);ylabel(dB);axis(0,1,-100,50); grid on; subplot(222);plot(w/pi,pha/pi); title(Phase); xlabel();ylabel(Phase(pi);axis(0,1,-1,1); grid on; subplot(224);plot(w/pi,grd); title(t); xlabel(pi);ylabel(x);axis(0,1,0,20); grid on; （2） 切比雪夫I型wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;Fs=10000;OmegaP=wp*Fs; OmegaS=ws*Fs; ep=sqrt(10(Rp/10)-1); Ripple=sqrt(1/(1+ep*ep); Attn=1/(10(As/20); N,OmegaC=cheb1ord(OmegaP,OmegaS,Rp,As,s); z0,p0,k0=cheb1ap(N,Rp); p=p0*OmegaC;z=z0*OmegaC; k=k0*OmegaCN; ba=k*real(poly(z); aa=real(poly(p); bd,ad=impinvar(ba,aa,Fs); H,w=freqz(bd,ad,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(bd,ad,w); subplot(2,2,1);plot(w/pi,mag); title(|A|); xlabel();ylabel(|H|);axis(0,1,0,1.1); set(gca,XTickMode,manual,XTick,0,0.4,0.6,1);grid on; set(gca,XTickMode,manual,YTick,0,Attn,Ripple,1); subplot(2,2,3);plot(w/pi,db);title(|A|(dB); xlabel(Frep(pi);ylabel(dB);axis(0,1,-100,50); grid on; subplot(2,2,2);plot(w/pi,pha/pi); title(Phase); xlabel();ylabel(Phase(pi);axis(0,1,-1,1); grid on; subplot(2,2,4);plot(w/pi,grd);title(lay off); xlabel(Frep(pi);ylabel(x);axis(0,1,0,20); grid on; （3） 切比雪夫II型wp=0.2*pi; ws=0.3*pi; Rp=1; As=15;Fs=10000;T=1/Fs; OmegaP=(2/T)*tan(wp/2);OmegaS=(2/T)*tan(ws/2); ep=sqrt(10(Rp/10)-1); Ripple=sqrt(1/(1+ep*ep);Attn=1/(10(As/20); N,OmegaC=cheb2ord(OmegaP,OmegaS,Rp,As,s);z0,p0,k0=cheb2ap(N,As);p=p0*OmegaC;z=z0*OmegaC; k=k0*OmegaCN;ba0=real(poly(z0);ba0=k0*ba0;aa0=real(poly(p0);ba=real(poly(z);ba=k*ba;aa=real(poly(p); bd,ad=bilinear(ba,aa,Fs);bd1,ad1=bilinear(ba0,aa0,Fs/OmegaC);subplot(1,1,1) H,w=freqz(bd,ad,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(bd,ad,w);subplot(2,2,1);plot(w/pi,mag); title(A); xlabel();ylabel(|A|);set(gca,XTickMode,manual,XTick,0,0.4,0.6,1);grid on; set(gca,XTickMode,manual,YTick,0,Attn,Ripple,1); subplot(2,2,3);plot(w/pi,db);title(|A|(dB); xlabel(freq(rad/s);ylabel(dB);grid on; subplot(2,2,2);plot(w/pi,pha/pi); title(Phase); xlabel();ylabel(pi);grid on; subplot(2,2,4);plot(w/pi,grd); title(lay off); xlabel(freq(rad/s);ylabel();grid on;（4） 椭圆模拟原型wp=0.2*pi; ws=0.3*pi; Rp=1; As=15;Fs=10000;T=1/Fs; OmegaP=(2/T)*tan(wp/2);OmegaS=(2/T)*tan(ws/2); ep=sqrt(10(Rp/10)-1); Ripple=sqrt(1/(1+ep*ep);Attn=1/(10(As/20); N,OmegaC=ellipord(OmegaP,OmegaS,Rp,As,s);z0,p0,k0=ellipap(N,Rp,As);p=p0*OmegaC;z=z0*OmegaC; k=k0*OmegaCN;ba0=real(poly(z0);ba0=k0*ba0;aa0=real(poly(p0);ba=real(poly(z);ba=k*ba;aa=real(poly(p); bd,ad=bilinear(ba,aa,Fs);bd1,ad1=bilinear(ba0,aa0,Fs/OmegaC);subplot(1,1,1) H,w=freqz(bd,ad,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(bd,ad,w);subplot(2,2,1);plot(w/pi,mag); title(A); xlabel();ylabel(|A|);set(gca,XTickMode,manual,XTick,0,0.4,0.6,1);grid on; set(gca,XTickMode,manual,YTick,0,Attn,Ripple,1); subplot(2,2,3);plot(w/pi,db);title(|A|(dB); xlabel(freq(rad/s);ylabel(dB);grid on; subplot(2,2,2);plot(w/pi,pha/pi); title(Phase); xlabel();ylabel(pi);grid on; subplot(2,2,4);plot(w/pi,grd); title(lay off); xlabel(freq(rad/s);ylabel();grid on;（5） 脉冲响应不变法和双线性变换法设计IIR数字滤波器的优缺点及使用范围脉冲响应不变法：优点：频率坐标的变换是线性的，若模拟滤波器的频响是限带于折叠频率以内的话，则通过变换后所得的数字滤波器的频响可以不失真的反应源响应与频率的关系。缺点：频谱的周期延拓效应，因此当模拟滤波器的频响不是限带于折叠频率以内的话，就不可避免的会产生频谱混叠失真。使用范围：当模拟滤波器频响在折叠频率以上衰减很大时，采用脉冲响应不变法设计的滤波器才满足要求。双线性变换法：优点：消除了脉冲响应不变法产生的频谱混叠失真。缺点：频率和之间有严重的非线性。使用范围：被变换的连续时间系统的响应必须是分段常数型的，不然数字域上的频率与连续时间域上的频谱会有很大的失真存在。五、主要结论、遇到的问题及解决方法，收获和体会。主要结论：同样的系统用巴特沃斯型滤波器，所需的阶数最少，并且也是最好看的。切比雪夫I型在通带会有选择性的滤波，切比雪夫II型在阻带会有一定的滤波效应，而椭圆滤波器的通带和阻带都会有一定的选择滤波效应。遇到的问题及解决办法：刚刚入手切比雪夫型的时候，不知道他的函数都是什么，而且不知道怎么用，后来看了实验教程里的巴特沃斯的滤波设计，就慢慢学会了。收获和体会：MATLAB功能好强大，又学会了几种用MATLAB分析方法。实验4 窗函数法设计FIR数字滤波器一、 实验目的(1) 掌握用窗函数法设计FIR数字滤波器的原理和方法。(2) 熟悉线性相位FIR数字滤波器特性。(3) 了解各种窗函数对滤波特性的影响。二、实验设备与环境计算机、MATLAB软件环境。三、 实验原理 滤波器的理想频率响应函数为Hd(ej)，则其对应的单位脉冲响应为hdn=12-piHdejejd窗函数设计法的基本原理是用有限长单位脉冲响应序列h(n)逼hd(n)。由于hd(n)往往是无限长序列，且是非因果的，所以用窗函数。w(n)将hd(n)截断，并进行加权处理： h(n) = hd(n) w(n) h(n)就作为实际设计的FIR数字滤波器的单位脉冲响应序列，其频率响应函数H(ej)为Hej=n=0N-1hnej用窗函数法设计的滤波器性能取决于窗函数w(n)的类型及窗口长度N的取值。设计过程中，要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度N。用窗函数法设计的滤波器性能取决于窗函数(n)的类型及窗口长度N的取值。设计过程中，要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度N。各种类型的窗函数可达到的阻带最小衰减和过渡带宽度可参见课本。 这样选定窗函数类型和长度N后，求出单位脉冲响应hn=hdnn，并求出H(ej)。H(ej)是否满足要求，要进行验算。一般在h(n)尾部加零使长度满足2的整数次幂，以便用FFT计算Hej。如果要观察细节，补零点数增多即可。如果不满足要求，则要重新选择窗函数类型和长度N，再次验算，直至满足要求。如果要求线性相位特性，则h(n)还必须满足：hn=h(N-1-n)，根据式中的正负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。要根据所设计的滤波特性正确选择其中一类。四、实验内容及步骤1、设计一个数字低通FIR滤波器，其技术指标如下：p=0.2，Rp=0.25dBst=0.3，As=50dB分别采用矩形窗、汉宁窗、海明窗、布莱克曼窗、凯瑟窗设计该滤波器。结合实验结果，分别讨论采用上述方法设计的数字滤波器是否都能满足给定指标要求。矩形窗：Wp=0.2*pi;Wst=0.3*pi;Rp=0.25;As=50;tr_width=Wst-Wp;N=ceil(1.8*pi/tr_width)+1;n=0:(N-1);Wc=(Wp+Wst)/2;alpha=(N-1)/2;hd=(Wc/pi)*sinc(Wc/pi)*(n-alpha);w_boxcar=boxcar(N);h=hd.*w_boxcar;subplot(221);stem(n,hd,filled);axis tight;xlabel(n);ylabel(hd(n);Hr,w1=zerophase(h);subplot(222);plot(w1/pi,Hr);axis tight;xlabel(omega/pi);ylabel(H(omega);subplot(223);stem(n,h,filled);axis tight;xlabel(n);ylabel(h(n);H,w=freqz(h,1);subplot(224);plot(w/pi,20*log10(abs(H)/max(abs(H);xlabel(omega/pi);ylabel(dB);grid on;由最后一张图看出，当0.3之后的区域内还有大于50dB的存在，所以，不满足设计要求。汉宁窗：Wp=0.2*pi;Wst=0.3*pi;Rp=0.25;As=50;tr_width=Wst-Wp;N=ceil(6.2*pi/tr_width)+1;n=0:(N-1);Wc=(Wp+Wst)/2;alpha=(N-1)/2;hd=(Wc/pi)*sinc(Wc/pi)*(n-alpha);w_hanning=hanning(N);h=hd.*w_hanning;subplot(221);stem(n,hd,filled);axis tight;xlabel(n);ylabel(hd(n);Hr,w1=zerophase(h);subplot(222);plot(w1/pi,Hr);axis tight;xlabel(omega/pi);ylabel(H(omega);subplot(223);stem(n,h,filled);axis tight;xlabel(n);ylabel(h(n);H,w=freqz(h,1);subplot(224);plot(w/pi,20*log10(abs(H)/max(abs(H);xlabel(omega/pi);ylabel(dB);grid on; 满足设计要求。海明窗：Wp=0.2*pi;Wst=0.3*pi;Rp=0.25;As=50;tr_width=Wst-Wp;N=ceil(6.6*pi/tr_width)+1;n=0:(N-1);Wc=(Wp+Wst)/2;alpha=(N-1)/2;hd=(Wc/pi)*sinc(Wc/pi)*(n-alpha);w_hamming=hamming(N);h=hd.*w_hamming;subplot(221);stem(n,hd,filled);axis tight;xlabel(n);ylabel(hd(n);Hr,w1=zerophase(h);subplot(222);plot(w1/pi,Hr);axis tight;xlabel(omega/pi);ylabel(H(omega);subplot(223);stem(n,h,filled);axis tight;xlabel(n);ylabel(h(n);H,w=freqz(h,1);subplot(224);plot(w/pi,20*log10(abs(H)/max(abs(H);xlabel(omega/pi);ylabel(dB);grid on; 满足设计要求。布莱克曼窗：Wp=0.2*pi;Wst=0.3*pi;Rp=0.25;As=50;tr_width=Wst-Wp;N=ceil(11*pi/tr_width)+1;n=0:(N-1);Wc=(Wp+Wst)/2;alpha=(N-1)/2;hd=(Wc/pi)*sinc(Wc/pi)*(n-alpha);w_blackman=blackman(N);h=hd.*w_blackman;subplot(221);stem(n,hd,filled);axis tight;xlabel(n);ylabel(hd(n);Hr,w1=zerophase(h);subplot(222);plot(w1/pi,Hr);axis tight;xlabel(omega/pi);ylabel(H(omega);subplot(223);stem(n,h,filled);axis tight;xlabel(n);ylabel(h(n);H,w=freqz(h,1);subplot(224);plot(w/pi,20*log10(abs(H)/max(abs(H);xlabel(omega/pi);ylabel(dB);grid on; 满足设计要求。凯瑟窗：Wp=0.2*pi;Wst=0.3*pi;Rp=0.25;As=50;tr_width=Wst-Wp;N=ceil(As-7.95)/(2.285*tr_width)+1;n=0:(N-1);Wc=(Wp+Wst)/2;alpha=(N-1)/2;hd=(Wc/pi)*sinc(Wc/pi)*(n-alpha);w_kaiser=kaiser(N);h=hd.*w_kaiser;subplot(221);stem(n,hd,filled);axis tight;xlabel(n);ylabel(hd(n);Hr,w1=zerophase(h);subplot(222);plot(w1/pi,Hr);axis tight;xlabel(omega/pi);ylabel(H(omega);subplot(223);stem(n,h,filled);axis tight;xlabel(n);ylabel(h(n);H,w=freqz(h,1);subplot(224);plot(w/pi,20*log10(abs(H)/max(abs(H);xlabel(omega/pi);ylabel(dB);grid on; 不满足设计要求。2.设计一个数字带通FIR滤波器，其技术指标如下下阻带边缘：st1=0.2，As=60dB，下通带边缘：p1=0.35，RP=1dB， 上通带边缘：p2=0.65，RP=1dB， 上阻带边缘：st2=0.8，As=60dB，Wp=0.15*pi;Wst=0.3*pi;Rp=1;As=60;tr_width=Wst-Wp;N=ceil(6.6*pi/tr_width)+1;n=0:(N-1);Wc=(Wp+Wst)/2;alpha=(N-1)/2;hd=(Wc/pi)*sinc(Wc/pi)*(n-alpha);w_hamming=hamming(N);h1=hd.*w_hamming;h=h1.*cos(0.5*pi*n);subplot(221);stem(n,hd,filled);axis tight;xlabel(n);ylabel(hd(n);Hr,w1=zerophase(h);subplot(222);plot(w1/pi,Hr);axis tight;xlabel(omega/pi);ylabel(H(omega);subplot(223);stem(n,h,filled);axis tight;xlabel(n);ylabel(h(n);H,w=freqz(h,1);subplot(224);plot(w/pi,20*log10(abs(H)/max(abs(H);xlabel(omega/pi);ylabel(dB);grid on;3.采用频率取样设计法设计 FIR数字低通滤波器，满足以下指标p=0.2，Rp=0.25dBst=0.3，As=50dB（1） 取 N = 20，过渡带没有样本。源代码：N=20;alpha=(N-1)/2;l=0:N-1;w1=(2*pi/N)*l;Hrs=1,1,1,zeros(1,15),1,1;Hdr=1,1,0,0;wd1=0,0.25,0.25,1;k1=0:floor(N-1)/2);k2=floor(N-1)/2)+1:N-1;angH=-alpha*(2*pi)/N*k1,alpha*(2*pi)/N*(N-k2);H=Hrs.*exp(1j*angH);h=ifft(H,N);w=0:500*pi/500;H=freqz(h,1,w); Hr,wr=zerophase(h);subplot(221);plot(wd1,Hdr,w1(1:11)/pi,Hrs(1:11),o);axis(0,1,-0.1,1.1);xlabel(omega(pi);ylabel(H(k);subplot(222);stem(l,h,filled);axis(0,N-1,-0.1,0.3);xlabel(n);ylabel(h(n);subplot(223);plot(wr/pi,Hr,w1(1:11)/pi,Hrs(1:11),o);axis(0,1,-0.2,1.2);xlabel(omega(pi);ylabel(Hr(w);subplot(224);plot(w/pi,20*log10(abs(H)/max(abs(H);axis(0,1,-50,5);grid;xlabel(omega(pi);ylabel(dB);（2） 取 N = 40，过渡带有一个样本，T = 0.39。源代码：N=40;alpha=(N-1)/2;l=0:N-1;w1=(2*pi/N)*l;Hrs=1,1,1,1,1,0.39,zeros(1,29),0.39,1,1,1,1;Hdr=1,1,0.39,0,0;wd1=0,0.2,0.25,0.3,1;k1=0:floor(N-1)/2);k2=floor(N-1)/2)+1:N-1;angH=-alpha*(2*pi)/N*k1,alpha*(2*pi)/N*(N-k2);H=Hrs.*exp(1j*angH);h=ifft(H,N);w=0:500*pi/500;H=freqz(h,1,w); Hr,wr=zerophase(h);subplot(221);plot(wd1,Hdr,w1(1:21)/pi,Hrs(1:21),o);axis(0,1,-0.1,1.1);xlabel(omega(pi);ylabel(H(k);subplot(222);stem(l,h,filled);axis(0,N-1,-0.1,0.3);xl