数字信号处理实验

1、数字信号处理实验姓名：学号：班级： 专业：通信工程实验 1 基 2-FFT 算法实现一、实验目的掌握基 2-FFT 的原理及具体实现方法。编程实现基 2-FFT 算法。加深理解 FFT 算法的特点。二、实验设备与环境计算机、MATLAB 软件环境。三、实验基础理论FFT 是 DFT 的一种快速算法，能使 DFT 的计算大大简化，运算时间缩短。FFT 利用了的三个固有特性，即对称性、周期性和可约性，将长序列的 DFT 分解为短序列的 DFT，合并了DFT 运算中的某些项，从而减少了 DFT 的运算量。FFT 算法基本上可以分为两大类，即按时间抽取法和按频率抽取法。在实现 FFT 算法时，要重点考

2、虑两个问题，注意数据的读取和存储：（1）输入输出的排序；（2）碟形运算的实现。按时间抽取算法中输入反序输出顺序，按频率抽取算法中输入顺序输出反序；运算过程中的每一级都由 N/2 个碟形运算构成，每一个碟形运算单元中，两个节点变量运算后得到的结果为下一列相同位置的节点变量，而和其他节点变量无关，可以采用原位运算，节省存储单元。另外，碟形运算中的复系数可以存储为能及时查阅的系数表，这样可以节约计算量，但是需要 N/2 个复数存储器。MATLAB 中提供了用于计算 FFT 的函数 fft，可将实验中所得到的结果与利用 MATLAB 中fft 函数计算的结果相比较，以此验证结果的正确性。四、实验内容编

3、程实现序列长度为N=8 的按时间抽取的基 2-FFT 算法。给定一个 8 点序列，采用编写的程序计算其 DFT，并与 MATLAB 中 fft 函数计算的结果相比较，以验证结果的正确性。编程实现序列长度为 N=8 的按频率抽取的基 2-FFT 算法。给定一个 8 点序列，采用编写的程序计算其 DFT，并与 MATLAB 中 fft 函数计算的结果相比较，以验证结果的正确性。将上述 FFT 程序推广到序列长度为N=2v 的情况，要求利用原位运算。五、实验代码及分析实验代码： x=1+2j 0.5+3j 5+4j -2+3j 6-3j 5+1j 9 -1-2j;n=bin2dec(fliplr(d

4、ec2bin(1:8-1,3)+1;x0=x;for l=1:3for m=1:2(l-1) for k=1:2(3-l)t=x0(k+(m-1)*2(4-l)+2(3-l);x0(k+(m-1)*2(4-l)+2(3-l)=(x0(k+(m-1)*2(4-l)-t)*exp(-2*1j*pi*(k-1)/(8/2(l-1); x0(k+(m-1)*2(4-l)=x0(k+(m-1)*2(4-l)+t;endendendy1=fft(x)for i=1:8y0(i)=x0(n(i);endy0y1 =Columns 1 through 523.5000 + 8.0000i1.4749 +10.

5、7678i-4.0000 -13.5000i-1.5754 + 7.0104i18.5000- 2.0000iColumns 6 through 8-3.4749 + 7.2322i -10.0000 + 3.5000i -16.4246 - 5.0104iy0 =Columns 1 through 523.5000 + 8.0000i1.4749 +10.7678i-4.0000 -13.5000i-1.5754 + 7.0104i18.5000- 2.0000iColumns 6 through 8-3.4749 + 7.2322i -10.0000 + 3.5000i -16.4246

6、- 5.0104i分析：y0 和y1 两者结果完全相同，说明结果正确。实验代码：x=1+2j 0.5+3j 5+4j -2+3j 6-3j 5+1j 9 -1-2j;n=bin2dec(fliplr(dec2bin(1:8-1,3)+1;x0=x;for l=1:3for m=1:2(l-1) for k=1:2(3-l)t=x0(k+(m-1)*2(4-l)+2(3-l);x0(k+(m-1)*2(4-l)+2(3-l)=(x0(k+(m-1)*2(4-l)-t)*exp(-2*1j*pi*(k-1)/(8/2(l-1); x0(k+(m-1)*2(4-l)=x0(k+(m-1)*2(4-l

7、)+t;end endendy1=fft(x)for i=1:8y0(i)=x0(n(i);endy0y1 =Columns 1 through 523.5000 + 8.0000i1.4749 +10.7678i-4.0000 -13.5000i-1.5754 + 7.0104i18.5000- 2.0000iColumns 6 through 8-3.4749 + 7.2322i -10.0000 + 3.5000i -16.4246 - 5.0104iy0 =Columns 1 through 523.5000 + 8.0000i1.4749 +10.7678i-4.0000 -13.5

8、000i-1.5754 + 7.0104i18.5000- 2.0000iColumns 6 through 8-3.4749 + 7.2322i -10.0000 + 3.5000i -16.4246 - 5.0104i分析：两种算法的结果相同，说明结果相同。实验代码：按时间抽取 xn= 1,4,7,5,9,2,3,8,6,5,2,7,6,4,13,2; M=nextpow2(length(xn); N=2M; for m=0:N/2-1WN(m+1)=exp(-j*2*pi/N)m; endDF1=xn,zeros(1,N-length(xn);H=0;for I=0:N-1;if IK=

9、N/2;while H=K;H=H-K;K=K/2;endH=H+K;endfor G=1:M; F=2(G-1); for S=0:F-1;P=2(M-G)*S; for K=S:2G:N-2;T=DF1(K+1)+DF1(K+F+1)*WN(P+1); DF1(K+F+1)=DF1(K+1)-DF1(K+F+1)*WN(P+1); DF1(K+1)=T;end endend DF1 DF1 =Columns 1 through 1084.0000-5.0057-7.8833-5.0141-3.0131-5.2100-7.9867-5.006010.0000-4.5526Columns 11

10、 through 16-8.1429-4.9962-2.9869-5.2093-7.9872-5.0060F1=fft(xn) F1 =Columns 1 through 1084.0000-5.0057-7.8833-5.0141-3.0131-5.2100-7.9867-5.006010.0000-4.5526Columns 11 through 16-8.1429-4.9962-2.9869-5.2093-7.9872-5.0060按频率抽取 xn= 1,4,7,5,9,2,3,8,6,5,2,7,6,4,13,2;M=nextpow2(length(xn);N=2M; M=log2(N

11、); for k1=0:M-1 D=2k1;E=N/2k1; F=N/2(k1+1); G=N/2(k1+1)-1;Wn=exp(-j*2*pi/E); for g=1:DH1=(g-1)*E;H2=(g-1)*E+F;for r=0:G;k=r+1; xn(k+H1)=xn(k+H1)+xn(k+H2);xn(k+H2)=xn(k+H1)-xn(k+H2)-xn(k+H2)*Wnr;end endendn1=fliplr(dec2bin(0:N-1);n2=bin2dec(n1);for i=1:NXk(i)=xn(n2(i)+1); end XkXk =Columns 1 through

12、1084.0000-5.0057-7.8833-5.0141-3.0131-5.2100-7.9867-5.006010.0000-4.5526Columns 11 through 16-8.1429-4.9962-2.9869-5.2093-7.9872-5.0060F2=fft(xn) F2 =Columns 1 through 1084.0000-5.0057-7.8833-5.0141-3.0131-5.2100-7.9867-5.006010.0000 -4.5526Columns 11 through 16-8.1429-4.9962-2.9869-5.2093-7.9872-5.

13、0060六、 收获与感悟通过本次数字信号处理实验，学习了基 2-FFT 算法的实现，通过对 MATLAB 中 fft 函数的运行结果与所编写程序结果的对比，了解了按时间抽取算法和按频率抽取算法的异同，按时间抽取算法，两种方式的结果相同，而按频率抽取算法，两种方式的结果不同。并且知道了碟形图在 MATLAB 中的实现方法实验 2 利用 DFT 分析信号频谱一、实验目的加深对 DFT 原理的理解。应用 DFT 分析信号的频谱。深刻理解利用 DFT 分析信号频谱的原理，分析实现过程中出现的现象及解决方法。二、实验设备与环境计算机、MATLAB 软件环境。三、实验基础理论DFT 和 DTFT 的关系有

14、限长序列x(n)(0 n N 1)的离散时间傅里叶变换X()在频率区间（0 2） 的 N 个等间隔分布的点k = 2k/N(0 k N 1)上的N 个取样值可以由下式表示：由上式可知，序列 x(n)的 N 点 DFTX(k),实际上就是 x(n)序列的 DTFT 在 N 个等间隔频率点 k = 2k/N(0kN 1)上样本 X(k)。利用 DFT 求 DTFT方法 1：其中(x)为内插函数方法 2：然而在实际 MATLAB 计算中，上述插值运算不见得是最好的办法。由于 DFT是 DTFT 的取样值，其相邻两个频率样本点的间距为 2/N，所以如果我们增加数据的长度N，使得到的 DFT 谱线就更加

15、精细，其包络就越接近 DTFT 的结果，这样就可以利用 DFT来近似计算DTFT，如果没有更多的数据，可以通过补零来增加数据长度。利用 DFT 分析连续时间信号的频谱采用计算机分析连续时间信号的频谱，第一部就是把连续时间信号离散化。这里需要两个操作：一是采样，二是截断。对于连续时间非周期信号()，按采样间隔T 进行采样，截取长度为 M，那么对()进行 N 点频域采样，得到因此，可以将利用 DFT 分析连续非周期信号频谱的步骤归纳如下：（1）确定时域采样间隔 T，得到离散序列 x(n)；（2）确定截取长度 M，得到 M 点离散序列() = ()(),这里()为窗函数。确定频域采样点数 N，要求

16、NM。利用 FFT 计算离散序列的 N 点 DFT，得到()。（5）根据式由()计算()采样点的近似值。采用上述方法计算()的频谱，需要注意如下三个问题：频谱混叠。如果不满足采样定理的条件，频谱会出现混叠误差。对于频谱无限宽的信号，应考虑覆盖大部分主要频率分量的范围。栅栏效应和频谱分辨率。使用 DFT 计算频谱，得到的结果只是 N 个频谱样本值， 样本值之间的频谱是未知的，像通过一个栅栏观察频谱，称为“栅栏效应”。频谱分辨率与记录长度成反比，要提高频谱分辨率，就要增加记录时间。频谱泄露。对信号阶段会把窗函数的频谱引入信号频谱，造成频谱泄露，解决这个问题的主要办法是采用旁瓣小的窗函数，频谱泄露和

17、窗函数均会引起误差。因此，要合理选取采样间隔和截取长度，必要时还考虑加适当的窗。对于连续时间周期信号，我们在采用计算机进行计算时，也总是要进行截断，序列总是有限长的，仍然可以采用上述方法近似计算。可能用到的 MATLAB 函数与代码实验中 DFT 运算可采用 MATLAB 中提供的函数 fft 来实现。DTFT 可以利用 MATLAB 矩阵运算的方法进行计算X() = = , , , 1 2 四、实验内容=1121.已知 x(n)=2,-1,1,1,完成如下要求：计算器 DTFT，并画出*, +区间的波形。计算 4 点 DFT，并把结果显示在(1)所画的图形中。对 x(n)补零，计算 64 点

18、 DFT，并显示结果。根据实验结果，分析是否可以由 DFT 计算 DTFT，如果可以，如何实现？考察序列0 n10 时，用 DFT 估计 x(n)的频谱；将 x(n)补零加长到长度为 100 点序列用 DFT估计 x(n)的频谱，要求画出相应波形。0 n 100 时，用 DFT 估计 x(n)的频谱，并画出波形。根据实验结果，分析怎样提高频谱分辨率。已知信号 ，其中1 =1, 2 = 2, 3 = 3,从 x(n)的表达式可以看出，它包含三个频率的正弦波，但是，从其时域波形来看，似乎是一个正弦信号，利用 DFT 做频谱分析，确定适合的参数，使得到的的频率分辨率复合需要。利用 DFT 近似分析连

19、续时间信号 的频谱（幅度谱）。分析采用不同的采样间隔和截取长度进行计算的结果，并最终确定适合的参数。五、实验结果及分析1.（1）实验代码：x=2 -1 1 1;n=0:3;w=-pi:0.01*pi:pi;X=x*exp(-j*n*w);subplot(211);plot(w,abs(X);xlabel(Omega/pi);title(Magnitude);axis tightsubplot(212);plot(w,angle(X)/pi);xlabel(Omega/pi);title(Phase);axis tight波形图：Magnitude321-3-2-10123/ Phase0.50

20、-0.5-3-2-10123/实验代码：fft(x);k=0:3;subplot(211);hold onstem(k,abs(ans);subplot(212);hold onstem(k,angle(ans)/pi);计算结果图示如下：Magnitude321-3-2-10123/ Phase0.50-0.5-3-2-10123/实验代码：fft(x,64);k=0:63;subplot(211);stem(k,abs(ans);xlabel(k);title(Magnitude);axis tightsubplot(212);stem(k,angle(ans);xlabel(k);tit

21、le(Phase);axis tight波形图如下：Magnitude32100102030405060k Phase10-10102030405060k分析：可以用 DFT 计算 DTFT。实现方法：利用对 x(n)补零，提高采样密度，当补零的个数足够多时，即可用 DFT 的结果近似 DTFT。2.（1）当 0 n 10时，实验代码:n=0:10;x=cos(0.48*pi*n)+cos(0.52*pi*n);X=fft(x);subplot(211);k=0:10;stem(k,abs(X);xlabel(k);title(Magnitude);axis tightsubplot(212)

22、;stem(k,angle(X);xlabel(k);title(Phase);axis tightx(n)的频谱为：Magnitude86420012345678910k Phase10-1012345678910k当 x(n)的长度变为 100 时，实验代码：X1=fft(x,100);subplot(211);k=0:99;stem(k,abs(X1);xlabel(k);title(Magnitude);axis tightsubplot(212);stem(k,angle(X1);xlabel(k);title(Phase);axis tightx(n)的频谱为：Magnitude1

23、0500102030405060708090k Phase10-10102030405060708090k（2）当0 n 100时，实验代码:n=0:100;x=cos(0.48*pi*n)+cos(0.52*pi*n);X=fft(x);subplot(211);k=0:100;stem(k,abs(X);xlabel(k);title(Magnitude);axis tightsubplot(212);stem(k,angle(X);xlabel(k);title(Phase);axis tight频谱波形：Magnitude402000102030405060708090100k Pha

24、se210-1-20102030405060708090100k（3）分析：由频谱波形可以看出，通过增加补零的个数可以提高频谱分辨率，但是不能提高分辨力，提高分辨力需要通过延长所选信号的长度来实现。实验代码： t=0:0.01:1f1=1;f2=2;f3=3;x=0.15*sin(2*pi*f1.*t)+sin(2*pi*f2.*t)-0.1*sin(2*pi*f3.*t);plot(t,x);grid;时域波形图为：1.510.50-0.5-1-1.500.10.20.30.40.50.60.70.80.91 n=0:2000; x=0.15*sin(2*pi*f1.*n)+sin(2*pi

25、*f2.*n)-0.1*sin(2*pi*f3.*n); X=fft(x); figure(2); stem(n,X,filled);-10 x 10420-2-4-6-8-1002004006008001000 1200 1400 1600 1800 2000分析：由图中频谱可知，存在 3 个频谱分量。实验代码：（1）采样区间：0，100，采样间隔为 1 时： n=0:100; x=exp(-0.1*n); X1=fft(x); k=0:100; stem(k,abs(X1); xlabel(k);波形图：1210864200102030405060708090100k（2）采样区间：0，1

26、00，采样间隔为 2 时： n=0:2:100; x=exp(-0.1*n); X2=fft(x); k=0:50; stem(k,abs(X2); xlabel(k);波形图：654321005101520253035404550k采样区间：0，50，采样间隔为 1 时： n=0:50; x=exp(-0.1*n); X3=fft(x); k=0:50; stem(k,abs(X3); xlabel(k);波形图：12108642005101520253035404550k采样区间：0，50，采样间隔为 2 时： n=0:2:50; x=exp(-0.1*n); X4=fft(x); k=0

27、:25; stem(k,abs(X4); xlabel(k);波形图：65432100510152025k（5）采样区间：0，150，采样间隔为 1 时： n=0:150; x=exp(-0.1*n); X5=fft(x); k=0:150; stem(k,abs(X5); xlabel(k);波形图：121086420050100150k（6）采样区间：0，150，采样间隔为 2 时： n=0:2:150; x=exp(-0.1*n); X5=fft(x); k=0:75; stem(k,abs(X5); xlabel(k);波形图：654321001020304050607080k分析：由

28、上述对比发现，采样区间取0，100，采样间隔取 2 时，结果最为接近 x(t)的频谱。六、收获及体会通过本次实验，学习到了如何用 DFT 分析频谱，知道了如何提高频谱的分辨力和分辨率，理解力二者的区别，并且能够利用 DFT 来近似DTFT 的结果。通过观察频谱波形，能够看到频谱泄漏现象，使得在时域无法观察到的波形，能够在频域有所体现。不同的采样间隔和信号截取长度，所得到的 DFT 结果也有所不同，需要选择合适的参数，才能得到完整清晰的频谱。实验中加深了对 Matlab 进行数字信号分析的理解，能够通过更加直观的方式，了解到课本上的知识，使得对课本知识的学习有了更深刻的体会。实验 5 脉冲响应不

29、变法设计 IIR 数字滤波器一、实验目的掌握利用脉冲响应不变法设计 IIR 数字滤波器的原理及具体方法。加深理解数字滤波器和模拟滤波器之间的技术指标转化。掌握脉冲响应不变法设计 IIR 数字滤波器的优缺点及适用范围。二、实验上机环境计算机、MATLAB 软件环境三、实验原理1、基本原理从时域响应出发，使数字滤波器的单位脉冲响应 h(n)模仿模拟滤波器的单位冲激响应(),h(n)等于()的取样值。2、变换方法思路：将H a(s)进行部分分式展开对 H a(s)进行拉式变换对ha (t )时域采样得到h(n)对h(n)进行 z 变换3.设计步骤确定数字滤波器的性能指标p ,st , Rp , As

30、 。将数字滤波器频率指标转换成响应的模拟滤波器频率指标 ppT ststT根据指标 p ， st ， Rp 和 As 设计模拟滤波器 H a (s) 。将 Ha (s) 展成部分分式形式NH a (s)=k 1Ak。s pkk把模拟极点 p 转换成数字极点e pkT ，得到数字滤波器NH (z) Akk 1 1 e pkT z1可见 H a (s) 至 H(z)间的变换关系为1s sk11 eskT z1zz eskT在 MATLAB 中有两种方法可以实现实现上述变换。方法 1：利用 residue 函数和 residuez 函数实现脉冲响应不变变换法，实用方法如下：实现多项式形式和部分分式形

31、式之间的转换。实现多项式形式和部分分式形式之间的转换。方法 2：matlab 中提供了 impinvar 函数采用脉冲响应不变法实现模拟滤波器到数字滤波器的变换，其使用如下：bz,az=impinvar(b,a,fs)采用脉冲响应不变法将模拟滤波器系统函数的系数向量 b 和 a 变换成为数字滤波器系统函数的系数向量 bz 和 az，fs 为采样频率（默认为 1）。bz,az=impinvar(b,a) 采样频率默认为 1 的情况下，采用脉冲响应不变法将模拟滤波器变换为数字滤波器。四、实验内容设采样频率为 fs =4kHz，采用脉冲响应不变法设计一个三阶巴特沃斯数字低通滤波器，其 3dB 截止频

32、率为 fc=1kHz。设采样频率为 fs=10kHz，设计数字低通滤波器，满足如下指标通带截止频率：fp=1kHz，通带波动：Rp=1dB阻带截止频率：fst=1.5kHz，阻带衰减：As=15dB要求分别采用巴特沃斯、切比雪夫 I 型、切比雪夫 II 型和椭圆模拟原型滤波器及脉冲响应不变法进行设计。结合实验结果，分别讨论采用上述设计的数字滤波器是否都能满足给定指标要求，分析脉冲响应不变法设计 IIR 数字滤波器的优缺点及适用范围。五、实验结果及分析1.实验代码：fs=4000;fc=1000;Wc=fc*2*pi; b=Wc3;a=1 2*Wc 2*Wc2 Wc3;bz,az=impinva

33、r(b,a,fs) bz =0.0000az =0.58130.21141.0000-0.39840.2475-0.0432w=0:500*pi/500; H,w=freqz(bz,az);subplot(221);plot(w/pi,abs(H);grid on;xlabel(omega(pi);ylabel(|H(ejomega)|);subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H);grid on;xlabel(omega(pi);ylabel(|H(ejomega)|);subplot(223);plot(w/pi,angle(H)/p

34、i);grid on;xlabel(omega(pi),ylabel(Phase of H(ejomega);grd=grpdelay(bz,az,w);subplot(224);plot(w/pi,grd);grid on;xlabel(omega(pi),ylabel(Group delay); G=tf(bz,az,0.1,variable,z-1) G =0.5813 z-1 + 0.2114 z-2- 1 - 0.3984 z-1 + 0.2475 z-2 - 0.04321 z-3波形图如下：1|H(ej)|0.80.60.40.200.51()0|H(ej)|-5-10-1500

35、.51()1Phase of H(ej)0.50-0.5-100.51()2Group delay10-100.51()2.（1）巴特沃斯滤波器实验代码： fs=10000; fp=1000; fst=1500; Rp=1; As=15; Ap=Rp; Wp=fp*2*pi; Wst=fst*2*pi; N=ceil(log10(10(Ap/10)-1)/(10(As/10)-1)/(2*log10(Wp/Wst); Wc=Wp/(10(Ap/10)-1)(1/2/N); b,a=butter(N,Wc,s); bz,az=impinvar(b,a,fs) bz =-0.00000.00060

36、.01010.01610.00410.00010az =1.0000-3.36355.0684-4.27592.1066-0.57060.0661 w=0:500*pi/500; H,w=freqz(bz,az); subplot(2,2,1);plot(w/pi,abs(H); grid on;xlabel(omega(pi);ylabel(|H(ejomega)|); subplot(2,2,2);plot(w/pi,20*log10(abs(H)/max(abs(H); grid on;xlabel(omega(pi);ylabel(|H(ejomega)|); subplot(2,2,

37、3);plot(w/pi,angle(H)/pi); grid on;xlabel(omega(pi),ylabel(Phase of H(ejomega); grd=grpdelay(bz,az,w); subplot(2,2,4);plot(w/pi,grd); grid on;xlabel(omega(pi),ylabel(Group delay); G=tf(bz,az,0.1,variable,z-1) G =z-5-5.457e-15 + 0.000631 z-1 + 0.0101 z-2 + 0.01614 z-3 + 0.004101 z-4 + 0.0001033-1 - 3

38、.364 z-1 + 5.068 z-2 - 4.276 z-3 + 2.107 z-4 - 0.5706 z-5 + 0.06607 z-6Sample time: 0.1 seconds Discrete-time transfer function. 波形图：1|H(ej)|0.5000.51()0|H(ej)|-20-40-60-8000.51()1Phase of H(ej)0.50-0.5-100.51()10Group delay5000.51()由图线可知，结果满足设计指标要求。切比雪夫 I 型实验代码：fs=10000; fp=1000; fst=1500; Rp=1; As

39、=15; Ap=Rp; epc=sqrt(10(Ap/10)-1); Wp=fp*2*pi; Wc=Wp; Wst=fst*2*pi; N=ceil(acosh(sqrt(10(0.1*As)-1)/epc)/acosh(Wst/Wc); b,a=cheby1(N,Rp,Wp,s); bz,az=impinvar(b,a,fs) bz =0.00000.00540.01810.00400az =1.0000-3.05913.8323-2.29190.5495 w=0:500*pi/500; H,w=freqz(bz,az); subplot(2,2,1);plot(w/pi,abs(H); g

40、rid on;xlabel(omega(pi);ylabel(|H(ejomega)|); subplot(2,2,2);plot(w/pi,20*log10(abs(H)/max(abs(H); grid on;xlabel(omega(pi);ylabel(|H(ejomega)|); subplot(2,2,3);plot(w/pi,angle(H)/pi); grid on;xlabel(omega(pi),ylabel(Phase of H(ejomega); grd=grpdelay(bz,az,w); subplot(2,2,4);plot(w/pi,grd); grid on;

41、xlabel(omega(pi),ylabel(Group delay); G=tf(bz,az,0.1,variable,z-1) G =3.411e-17 + 0.005373 z-1 + 0.0181 z-2 + 0.003985 z-3-1 - 3.059 z-1 + 3.832 z-2 - 2.292 z-3 + 0.5495 z-4Sample time: 0.1 seconds Discrete-time transfer function. 波形图：1.5|H(ej)|10.5000.51()0|H(ej)|-20-40-60-8000.51()1Phase of H(ej)0

42、.50-0.5-100.51()15Group delay105000.51()由图线可知，满足设计指标要求。切比雪夫 II 型实验代码： fs=10000; fp=1000; fst=1500; Rp=1; As=15; epc=sqrt(10(Rp/10)-1); Wp=fp*2*pi; Wst=fst*2*pi; Wc=Wst; N=ceil(acosh(sqrt(10(0.1*As)-1)/epc)/acosh(Wc/Wp); b,a=cheby2(N,As,Wst,s); bz,az=impinvar(b,a,fs) bz =-0.42140.9724-0.81190.42050.0

43、000az =1.0000-1.66581.4289-0.51930.0935 w=0:500*pi/500; H,w=freqz(bz,az);subplot(2,2,1);plot(w/pi,abs(H);grid on;xlabel(omega(pi);ylabel(|H(ejomega)|);subplot(2,2,2);plot(w/pi,20*log10(abs(H)/max(abs(H);grid on;xlabel(omega(pi);ylabel(|H(ejomega)|);subplot(2,2,3);plot(w/pi,angle(H)/pi);grid on;xlabe

44、l(omega(pi),ylabel(Phase of H(ejomega);grd=grpdelay(bz,az,w);subplot(2,2,4);plot(w/pi,grd);grid on;xlabel(omega(pi),ylabel(Group delay);G=tf(bz,az,0.1,variable,z-1) G =-0.4214 + 0.9724 z-1 - 0.8119 z-2 + 0.4205 z-3 + 1.663e-06 z-4-1 - 1.666 z-1 + 1.429 z-2 - 0.5193 z-3 + 0.09352 z-4Sample time: 0.1

45、seconds Discrete-time transfer function.波形图:1.50|H(ej)|H(ej)|1-50.5000.51()-1000.51()1Phase of H(ej)0.50-0.5-100.51()5Group delay0-500.51()由图线可知，阻带衰减小于 15dB，不满足设计要求。椭圆模拟原型滤波器实验代码： fs=10000; fp=1000; fst=1500; Rp=1; As=15; epc=sqrt(10(Rp/10)-1); Wp=fp*2*pi; Wst=fst*2*pi; Wc=Wp; A=10(As/20); k1=epc/(s

46、qrt(A2-1); k=Wp/Wst; N=ceil(ellipke(k)*ellipke(sqrt(1-k12)/ellipke(k1)/ellipke(sqrt(1-k2); bz,az=ellip(N,Rp,As,fp/fs*pi) bz =0.18550.04960.04960.1855az =1.0000-1.41071.2357-0.3549 w=0:500*pi/500; H,w=freqz(bz,az);subplot(2,2,1);plot(w/pi,abs(H);grid on;xlabel(omega(pi);ylabel(|H(ejomega)|);subplot(2

47、,2,2);plot(w/pi,20*log10(abs(H)/max(abs(H);grid on;xlabel(omega(pi);ylabel(|H(ejomega)|);subplot(2,2,3);plot(w/pi,angle(H)/pi);grid on;xlabel(omega(pi),ylabel(Phase of H(ejomega);grd=grpdelay(bz,az,w);subplot(2,2,4);plot(w/pi,grd);grid on;xlabel(omega(pi),ylabel(Group delay);G=tf(bz,az,0.1,variable,

48、z-1) G =0.1855 + 0.04956 z-1 + 0.04956 z-2 + 0.1855 z-3-1 - 1.411 z-1 + 1.236 z-2 - 0.3549 z-3Sample time: 0.1 seconds Discrete-time transfer function. 波形图：1|H(ej)|0.5000.51()0|H(ej)|-20-40-60-8000.51()1Phase of H(ej)0.50-0.5-100.51()15Group delay105000.51()由图线可知，满足设计指标要求。脉冲相应不变法： 实验代码： fs=10000; fp

49、=1000; fst=1500; Rp=1; As=15; Ap=Rp; Wp=fp*2*pi; Wst=fst*2*pi; N=ceil(log10(10(Ap/10)-1)/(10(As/10)-1)/(2*log10(Wp/Wst); Wc=Wp/(10(Ap/10)-1)(1/2/N); b=Wc3; a=1 2*Wc 2*Wc2 Wc3; bz,az=impinvar(b,a,fs) bz =00.10570.0663az =1.0000-1.64911.0663-0.2450 w=0:500*pi/500; H,w=freqz(bz,az);subplot(2,2,1);plot(

50、w/pi,abs(H);grid on;xlabel(omega(pi);ylabel(|H(ejomega)|);subplot(2,2,2);plot(w/pi,20*log10(abs(H)/max(abs(H);grid on;xlabel(omega(pi);ylabel(|H(ejomega)|);subplot(2,2,3);plot(w/pi,angle(H)/pi);grid on;xlabel(omega(pi),ylabel(Phase of H(ejomega);grd=grpdelay(bz,az,w);subplot(2,2,4);plot(w/pi,grd);gr

51、id on;xlabel(omega(pi),ylabel(Group delay);G=tf(bz,az,0.1,variable,z-1) G =0.1057 z-1 + 0.06633 z-2-1 - 1.649 z-1 + 1.066 z-2 - 0.245 z-3Sample time: 0.1 secondsDiscrete-time transfer function.波形图：1|H(ej)|0.5000.51()0|H(ej)|-20-40-6000.51()1Phase of H(ej)0.50-0.5-100.51()4Group delay20-200.51()由图线可知

52、，满足设计指标要求。分析：脉冲响应不变法设计 IIR 数字滤波器的优缺点：优点：脉冲响应不变法的频率坐标的变换是线性的。因此，如果模拟滤波器的的频率响应是限带于折叠频率以内的话，则通过变换后所得的数字滤波器的频率响应可以不失真的反映原响应与频率的关系。缺点：脉冲响应不变法最大的缺点是有频谱的周期延拓效应。其次，脉冲响应不变法所得的数字滤波器会出现频率的混叠失真。第三，数字滤波器在 T 很小（或者说取样频率很高）时，可能具有不希望的高增益。适用范围：脉冲响应不变法设计 IIR 数字滤波器只能用于限带的频率响应，如衰减特性很好的低通或带通。六、收获与体会通过本次实验，学会了利用用脉冲响应不变法设计

53、 IIR 数字滤波器的原理及具体方法， 通过比较不同类型的数字滤波器，了解了它们各自的优缺点以及适用范围，在相同的技术指标下，不同滤波器的效果也有较大差异，不是所有类型都可以满足要求。了解了如何将所给的技术指标，转化为设计过程中的参数，包括模拟域和数字域的转化，加深了对课本内容的理解与认识，进一步掌握了 IIR 数字滤波器的设计要点。实验 7窗函数法设计 FIR 数字滤波器一、实验目的掌握窗函数法设计 FIR 数字滤波器的原理和具体方法二、实验设备与环境计算机、Matlab 软件环境三、实验基础理论1、基本原理窗函数设计法的基本思想为，首先选择一个适当的理想的滤波器H()，然后用窗函数截取它的

54、单位脉冲响应h()，得到线性相位和因果的 FIR 滤波器，这种方法的重点是选择一个合适的窗函数和理想滤波器，使设计的滤波器的单位脉冲响应逼近理想滤波器的单位脉冲响应。2、设计步骤给定理想滤波器的频率响应H()，在通带上具有单位增益和线性相位，在阻带上具有零响应。一个带宽为( Wp=0.2*pi; Wst=0.3*pi; tr_width=Wst-Wp; Wn=(Wp+Wst)/2; Wn1=Wn/pi; Rp=0.25; As=50; N=ceil(1.8*pi/tr_width)+1; h=fir1(N,Wn1,boxcar(N+1); H,w=freqz(h,1); figure(1);

55、plot(w/pi,20*log10(abs(H)/max(abs(H); grid on; xlabel(omegapi);ylabel(dB); figure(2);plot(w/pi,angle(H); xlabel(omegapi);ylabel(rad); grid on;频率特性曲线：0-10-20-30dB-40-50-60-70-80-9000.10.20.30.40.50.60.70.80.914321rad0-1-2-3-400.10.20.30.40.50.60.70.80.91分析：阻带衰减小于 50dB，故不满足设计指标要求。汉宁窗：实验代码： Wp=0.2*pi;W

56、st=0.3*pi;tr_width=Wst-Wp;Wn=(Wp+Wst)/2;Wn1=Wn/pi;Rp=0.25;As=50;N=ceil(6.2*pi/tr_width)+1;h=fir1(N,Wn1,hanning(N+1); H,w=freqz(h,1);figure(1);plot(w/pi,20*log10(abs(H)/max(abs(H);grid on;xlabel(omegapi);ylabel(dB);figure(2);plot(w/pi,angle(H);xlabel(omegapi);ylabel(rad);grid on;频率特性曲线：0-50dB-100-15000.10.20.30.40.50.60.70.80.914321rad0-1-2-3-400.10.20.30.40.50.60.70.80.91分析：由图可知，满足设计指标要求。海明窗：实验代码： Wp=0.2*pi;Wst=0.3*pi;tr_width=Wst-Wp;Wn=(Wp+Ws