数字信号处理Matlab课后实验(吴镇扬)

1、数字信号处理实验报告实验一 熟悉MATLAB环境实验二 信号的采样与重建实验三 快速变换及其应用实验四 IIR数字滤波器的设计实验五 FIR数字滤波器的设计实验一 熟悉MATLAB环境一、实验目的(1)熟悉MATLAB的主要操作命令。(2)学会简单的矩阵输入和数据读写。(3)掌握简单的绘图命令。(4)用MATLAB编程并学会创建函数。(5)观察离散系统的频率响应。二、实验内容认真阅读本章附录，在MATLAB环境下重新做一遍附录中的例子，体会各条命令的含义。在熟悉了MATLAB基本命令的基础上，完成以下实验。上机实验内容：(1)数组的加、减、乘、除和乘方运算。输入A=1 2 3 4，B=3 4

2、5 6，求C=A+B，D=A-B，E=A.*B，F=A./B，G=A.B并用stem语句画出A、B、C、D、E、F、G。实验程序：A=1 2 3 4;B=3 4 5 6;n=1:4;C=A+B;D=A-B;E=A.*B;F=A./B;G=A.B;subplot(4,2,1);stem(n,A,'fill');xlabel ('时间序列n');ylabel('A');subplot(4,2,2);stem(n,B,'fill');xlabel ('时间序列n ');ylabel('B');subplo

3、t(4,2,3);stem(n,C,'fill');xlabel ('时间序列n ');ylabel('A+B');subplot(4,2,4);stem(n,D,'fill');xlabel ('时间序列n ');ylabel('A-B');subplot(4,2,5);stem(n,E,'fill');xlabel ('时间序列n ');ylabel('A.*B');subplot(4,2,6);stem(n,F,'fill');x

4、label ('时间序列n ');ylabel('A./B');subplot(4,2,7);stem(n,G,'fill');xlabel ('时间序列n ');ylabel('A.B');运行结果：(2)用MATLAB实现以下序列。a）x(n)=0.8n 0n15实验程序：n=0:15;x=0.8.n;stem(n,x,'fill'); xlabel ('时间序列n ');ylabel('x(n)=0.8n');b）x(n)=e(0.2+3j)n 0n15实验程序

5、：n=0:15;x=exp(0.2+3*j)*n);stem(n,x,'fill'); xlabel ('时间序列n ');ylabel('x(n)=exp(0.2+3*j)*n)');运行结果：a）的时间序列 b）的时间序列c）x(n)=3cos(0.125n+0.2)+2sin(0.25n+0.1) 0n15实验程序：n=0:1:15;x=3*cos(0.125*pi*n+0.2*pi)+2*sin(0.25*pi*n+0.1*pi);stem(n,x,'fill'); xlabel('时间序列n ');yla

6、bel('x(n)=3*cos(0.125*pi*n+0.2*pi)+2*sin(0.25*pi*n+0.1*pi)');运行结果：d）将c)中的x(n)扩展为以16为周期的函数x16(n)=x(n+16),绘出四个周期实验程序：n=0:1:63;x=3*cos(0.125*pi*rem(n,16)+0.2*pi)+2*sin(0.25*pi*rem(n,16)+0.1*pi);stem(n,x,'fill'); xlabel ('时间序列n ');ylabel('x16(n)');e）将c)中的x(n)扩展为以10为周期的函数x

7、10(n)=x(n+10),绘出四个周期实验程序：n=0:1:39;x=3*cos(0.125*pi*rem(n,10)+0.2*pi)+2*sin(0.25*pi*rem(n,10)+0.1*pi);stem(n,x,'fill'); xlabel ('时间序列n ');ylabel('x10(n)');运行结果：d）的时间序列 e）的时间序列(3)x(n)=1,-1,3,5，产生并绘出下列序列的样本。a）x1(n)=2x(n+2)-x(n-1)-2x(n)实验程序：n=0:3;x=1 -1 3 5;x1=circshift(x,0 -2);x

8、2=circshift(x,0 1);x3=2*x1-x2-2*x;stem(x3,'fill'); xlabel ('时间序列n ');ylabel('x1(n)=2x(n+2)-x(n-1)-2x(n)');b）实验程序：n=0:3;x=1 -1 3 5;x1=circshift(x,0 1);x2=circshift(x,0 2);x3=circshift(x,0 3);x4=circshift(x,0 4);x5=circshift(x,0 5);xn=1*x1+2*x2+3*x3+4*x4+5*x5;stem(xn,'fill&

9、#39;); xlabel('时间序列n ');ylabel('x2(n)=x(n-1)+2x(n-2)+3x(n-3)+4x(n-4)+5x(n-5)');运行结果： a）的时间序列b）的时间序列(4)绘出时间函数的图形，对x轴、y轴图形上方均须加上适当的标注。 a) x(t)=sin(2t) 0t10s b) x(t)=cos(100t)sin(t) 0t4s实验程序：clc;t1=0:0.001:10;t2=0:0.01:4;xa=sin(2*pi*t1);xb=cos(100*pi*t2).*sin(pi*t2);subplot(2,1,1);plot(

10、t1,xa);xlabel ('t');ylabel('x(t)');title('x(t)=sin(2*pi*t)');subplot(2,1,2);plot(t2,xb);xlabel ('t');ylabel('x(t)');title('x(t)=cos(100*pi*t2).*sin(pi*t2)');运行结果：(5)编写函数stepshift(n0,n1,n2)实现u(n-n0),n1<n0<n2,绘出该函数的图形，起点为n1，终点为n2。实验程序：clc;n1=input(

11、'请输入起点：');n2=input('请输入终点:');n0=input('请输入阶跃位置:');n=n1:n2;x=n-n0>=0;stem(n,x,'fill');xlabel('时间序列n');ylabel('u(n-n0)');请输入起点：2请输入终点：8请输入阶跃位置：6运行结果： （5）运行结果 （6）运行结果(6)给一定因果系统求出并绘制H(z)的幅频响应与相频响应。实验程序：a=1 -0.67 0.9;b=1 sqrt(2) 1;h w=freqz(b,a);fp=20*lo

12、g(abs(h);subplot(2,1,1);plot(w,fp);xlabel('时间序列t');ylabel('幅频特性');xp=angle(h);subplot(2,1,2);plot(w,xp);xlabel('时间序列t');ylabel('相频特性');运行结果：(右上图)实验二 信号的采样与重建一、 实验目的1 学习本章内容的基础上，通过实验加强有关信号采样与重建的基本概念，熟悉相关的Matlab函数。通过观察采样信号的混叠现象，进一步理解奈奎斯特采样频率的意义。通过实验，了解数字信号采样率转换过程中的频率特征。

13、对实际的音频文件作内插和抽取操作，体会低通滤波器在内插和抽取中的作用。二、 实验内容1.一个信号是三个正弦信号的和，正弦信号的频率为50Hz、500Hz、1000Hz，该信号以8kHz采样。用适当数量的样本画出该信号。2.一个信号是三个正弦信号的和，正弦信号的频率为50Hz、500Hz、1000Hz，该信号以800Hz采样。用适当数量的样本画出该信号，并讨论信号的混叠状况。3.令，其中，即每个周期内有16个点。试利用MATLAB编程实现：作M=4倍的抽取，使每个周期变成4点。作L=3倍的差值，使每个周期变成48点。4.输入信号x（n）为归一化频率分别是，的正弦信号相加而成，N=50，内插因子为

14、5，抽取因子为3，给出按有理因子5/3作采样率转换的输入输出波形。5.常见的音频文件采样率为44.1kHz。请找一个wav格式、采样率为44.1kHz的音频文件，用MATLAB编写程序，把它转换成采样率为48kHz、32kHz、22.05Khz、16kHz和8KHz的音频文件，用播放器分别进行播放，比较音质的变化，并解释原因。三、 实验程序和结果分析1.一个信号是三个正弦信号的和，正弦信号的频率为50Hz、500Hz、1000Hz，该信号以8kHz采样。用适当数量的样本画出该信号。答：n=0:99;x=sin(0.00625*2*pi*n)+sin(0.0625*2*pi*n)+sin(0.1

15、25*2*pi*n);stem(n,x);title('抽样图像');xlabel('n');ylabel('幅值');h,w=freqz(x);figure(2);plot(w(1:256),abs(h(1:256);title('幅频响应图像');xlabel('/omega');ylabel('幅度');分析：模拟信号 (A/D) 变换转换为数字信号，每隔一个采样频率 fs=1/8k，为保证采样后信号的频谱形状不失真，采样频率必须大于信号中最高频率成分的两倍，该采样频率符合奈奎斯特抽样定理，抽

16、样前后的频谱不会失真。由频谱图可以看出，当w=0，0.4，0.8处幅度最大值。n取值0：99范围，x（n）本身是一个带限信号，在w=0.8后幅度随着w的增大幅度逐渐衰减。2.一个信号是三个正弦信号的和，正弦信号的频率为50Hz、500Hz、1000Hz，该信号以800Hz采样。用适当数量的样本画出该信号，并讨论信号的混叠状况。Matlab编程：n=0:99;x=sin(0.0625*2*pi*n)+sin(0.625*2*pi*n)+sin(1.25*2*pi*n);stem(n,y);title('抽样图像');xlabel('n');ylabel('

17、幅值');h,w=freqz(x);figure(2);plot(w(1:512),abs(h(1:512);title('幅频响应图像');xlabel('/omega');ylabel('幅度');分析:对连续信号进行等间隔采样时，如果不能满足采样定理，采样后信号的频率就会重叠，即高于采样频率一半的频率成分将被重建成低于采样频率一半的信号。这种频谱的重叠导致的失真称为混叠，而重建出来的信号称为原信号的混叠替身，因为这两个信号有同样的样本值。在本实验中因为采样频率为800hz，所以会发生部分混叠，主要由500Hz、1000Hz的信号产生

18、。因为采样频率小于奈奎斯特采样频率，引起混叠，采样后信号的频谱形状失真，如若要采样后的频谱不发生失真，采样频率必须大于信号中最高频率成分的两倍2 （或 2）。（对取样频率的要求，即取样频率要足够大，采得的样值要足够多，才能恢复原信号。）即采样频率最小要为2000Hz3. 令，其中，即每个周期内有16个点。试利用MATLAB编程实现：作M=4倍的抽取，使每个周期变成4点。作L=3倍的插值，使每个周期变成48点。clear;N=52; M=4;f=1;fs=16;n=0:N-1; x=cos(2*pi*(f/fs)*n);stem(n,x(1:N);grid;title('原图像'

19、);xlabel('n');ylabel('幅度');y=decimate(x,M);figure(2);m=0:N/M-1;stem(m,y(1:N/M); title('抽样图像');xlabel('n');ylabel('幅度');grid;clear;N=50;L=3;f=1;fs=16;n=0:N-1;x=cos(2*pi*(f/fs)*n);figure(1);stem(n,x(1:N);title('inputsequence');xlabel('n');ylabel(

20、'幅度');grid;y=interp(x,L);figure(2);m=0:N*L-1;stem(m,y(1:N*L);title('outputsequence');xlabel('n');ylabel('幅度');grid;分析：减少抽样率以去掉过多数据的过程称为信号的抽取，增加抽样率以增加数据的过程称为信号的插值。抽取、插值及其二者相结合的使用便可实现信号抽样率的转换。抽取：如果要把采样速率减小M倍(M是整数)，可以把原始的采样序列每隔M-1个点取一个点，形成新的采样序列，该过程称为M倍抽取,此题中M=4，即从原信号中每4

21、个点抽取一个值组成新的序列。从图可以看出原信号有16个抽样值，而抽取后的信号一个周期只有4个抽样值。要求抽样频率fs必须满足fs>2Mfc³。但是，如果fs>2Mfc³的条件不能得到满足，那么频谱将发生混迭，因此也就无法重建出x (t)。插值：如果要把采样速率增大L倍，可以把原始信号相邻两个点之间取L-1个零点，增大抽样频率，它可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。插值会让信号频谱周期拓展，但由于输出信号点数增加，采样率增加，所以频率分量只会周期存在但不会发生混叠。如题图，信号进行3倍插值，频率增加3倍，可以从图抽样点数看出。每个周期的

22、抽样点数为16的3倍，48。4.输入信号x（n）为归一化频率分别是，的正弦信号相加而成，N=50，内插因子为5，抽取因子为3，给出按有理因子5/3作采样率转换的输入输出波形。clear;N=50;M=3;L=5;f1=0.04;f2=0.3;n=0:N-1;x=sin(2*pi*f1*n)+sin(2*pi*f2*n);figure(1);stem(n,x(1:N);title('input sequence');xlabel('n');ylabel('幅度')y=resample(x,L,M);figure(2);stem(n,y(1:N);t

23、itle('output sequence');xlabel('n');ylabel('幅度')figure(3);h,w=freqz(x);plot(w,abs(h);title('输入幅频响应');xlabel('w');ylabel('幅度');figure(4);h,w=freqz(y);plot(w,abs(h);title('输出幅频响应');xlabel('w');ylabel('幅度');分析：对给定的信号x (n)，若希望将抽样率转变

24、为数L / M倍，先将x (n)作M倍的抽取，再作L倍的插值来实现，或是先作L倍的插值，再作M倍的抽取。一般来说，抽取使x (n)的数据点减少，会产生信息的丢失，因此，合理的方法是先对信号作插值，然后再抽取。从幅频响应的图做对比可以看出，抽样速率转变后，频谱压缩了，相应的最大幅度值发生叠加变大。5.常见的音频文件采样率为44.1kHz。请找一个wav格式、采样率为44.1kHz的音频文件，用MATLAB编写程序，把它转换成采样率为48kHz、32kHz、22.05Khz、16kHz和8KHz的音频文件，用播放器分别进行播放，比较音质的变化，并解释原因。y,fs,bits=wavread(

25、9;1.wav'); c=resample(y,12,11); wavwrite(y,48000,bits,'2.wav'); sound(c,48000); figure; plot(c); title('采样频率提高12/11倍的时域图');y,fs,bits=wavread('1.wav'); c=resample(y,8,11); wavwrite(y,32000,bits,'3.wav'); sound(c,32000); figure; plot(c); title('采样频率提高8/11倍的时域图&#

26、39;);y,fs,bits=wavread('1.wav'); c=resample(y,1,2); wavwrite(y,22050,bits,'4.wav'); sound(c,22050); figure; plot(c); title('采样频率提高1/2倍的时域图');y,fs,bits=wavread('1.wav'); c=resample(y,4,11); wavwrite(y,16000,bits,'5.wav'); sound(c,16000); figure; plot(c); title(

27、'采样频率提高4/11倍的时域图') ;y,fs,bits=wavread('1.wav'); c=resample(y,2,11); wavwrite(y,8000,bits,'6.wav'); sound(c,8000); figure; plot(c); title('采样频率提高2/11倍的时域图');分析：采样频率越高所能描述的声波频率就越高。采样率决定声音频率的范围（相当于音调）当频率越低时，抽样点数越少，声音音调降低，音质越差。当频率越高，抽样点数越多，声音音调随之上升，音质也随之上升变好。音频在转换时候会产生失真，

28、采样率只能在采样的时候设置，之后被禁止改动，声音本身是连续的声波，是一组连续的模拟信号，量化编码后变成波形文件，此题中本省的采样率为44.1KHz，采样率就是每一秒对这个波形进行采样的次数，次数越多即采样率越高，则采样后的波形与模拟信号的波形越接近。如果采样点增加，声音就会被拉长，也就声音变慢，减少采样点，音速则增加。从采样后的波形中看或者听音频文件。可以得出，当抽样率变小时候，声音被拉长，音质变差。时域图变窄。四、 思考题1. 试说明对于周期信号，应当如何采样，才能保证周期扩展后与原信号保持一 致。答：整周期采样，即每周期完整地采集若干（通常是8、16及其倍数）点，就可以保证周期扩展后与原来

29、信号保持一致。2. 模拟抗混叠滤波器的指标是如何确定的，欠采样的情况下是否需要模拟抗混叠滤波器？答：带内波纹度：通带的幅值精度指标，例如：带内波纹度为±0.1dB时，对幅值精度的影响 约为±1%。 阻带下降斜率：滤波器在截至频率开始下降，下降斜率越大越好。一般采用每个倍频 程 的下降分贝数衡量。 滤波器落差：带通到带阻差值的分贝数dB。欠采样的情况下，没有混叠也就不需要模拟抗混叠滤波器。3. 抽取是否会造成信号频谱成分的丢失？为什么还要这么做？答：会造成信号频谱成分的丢失，这么做可以减少信号所占空间和信号本身的频率。实验三 快速变换及其应用 一、实验目的(1)在理

30、论学习的基础上，通过本实验，加深对FFT的理解，熟悉FFT子程序。(2)熟悉应用FFT对典型信号进行频谱分析的方法。(3)了解应用FFT进行信号频谱分析过程中可能出现的问题以便在实际中正确应用FFT。(4)熟悉应用FFT实现两个序列的线性卷积的方法。(5)对DCT变换用作数据压缩有初步的认识。二、实验内容(1)观察高斯序列的时域和幅频特性。实验程序：clc;n=0:1:15;% p=8,q=2p=8;q=2;x1=exp(-(n-p).2./q);fp1=fft(x1);fp1=abs(fp1);subplot(5,2,1);plot(x1);xlabel('n');ylabe

31、l('时域特性');title('p=8,q=2');subplot(5,2,2);plot(fp1);xlabel('n');ylabel('幅频特性');title('p=8,q=2');% p=8,q=4p=8;q=4;x2=exp(-(n-p).2./q);fp2=fft(x2);fp2=abs(fp2);subplot(5,2,3);plot(x2);xlabel('n');ylabel('时域特性');title('p=8,q=4');subplot(5,

32、2,4);plot(fp2);xlabel('n');ylabel('幅频特性');title('p=8,q=4');% p=8,q=8p=8;q=8;x3=exp(-(n-p).2./q);fp3=fft(x3);fp3=abs(fp3);subplot(5,2,5);plot(x3);xlabel('n');ylabel('时域特性');title('p=8,q=8');subplot(5,2,6);plot(fp3);xlabel('n');ylabel('幅频特性&#

33、39;);title('p=8,q=8');% p=13,q=8p=13,q=8;x4=exp(-(n-p).2./q);fp4=fft(x4);fp4=abs(fp4);subplot(5,2,7);plot(x4);xlabel('n');ylabel('时域特性');title('p=13,q=8');subplot(5,2,8);plot(fp4);xlabel('n');ylabel('幅频特性');title('p=13,q=8');% p=14,q=8p=14,q=8;

34、x5=exp(-(n-p).2./q);fp5=fft(x5);fp5=abs(fp5);subplot(5,2,9);plot(x5);xlabel('n');ylabel('时域特性');title('p=14,q=8');subplot(5,2,10);plot(fp5);xlabel('n');ylabel('幅频特性');title('p=14,q=8');(2) 观察衰减正弦序列xb(n)的时域和幅频特性。a=0.1,f=0.0625,检查谱峰出现位置是否正确，注意频谱的形状，绘出幅频特

35、性曲线，改变f，使f 分别等于0.4375 和0.5625,观察这两种情况下频谱的形状和谱峰出现为止，又无混叠和泄露现象？说明产生现象的原因。实验程序：clc;n=0:1:15;a=0.1;f1=0.0625;f2=0.4375;f3=0.5625;x1=exp(-a*n).*sin(2*pi*f1.*n);x2=exp(-a*n).*sin(2*pi*f2.*n);x3=exp(-a*n).*sin(2*pi*f3.*n);fp1=fft(x1);fp1=abs(fp1);fp2=fft(x2);fp2=abs(fp2);fp3=fft(x3);fp3=abs(fp3);subplot(3,

36、2,1);plot(x1);xlabel('n');ylabel('时域特性');title('a=0.1,f=0.0625');subplot(3,2,2);plot(fp1);xlabel('n');ylabel('幅频特性');title('a=0.1,f=0.0625');subplot(3,2,3);plot(x2);xlabel('n');ylabel('时域特性');title('a=0.1,f=0.4375');subplot(3,2,

37、4);plot(fp2);xlabel('n');ylabel('幅频特性');title('a=0.1,f=0.4375');subplot(3,2,5);plot(x3);xlabel('n');ylabel('时域特性');title('a=0.1,f=0.5625');subplot(3,2,6);plot(fp3);xlabel('n');ylabel('幅频特性');title('a=0.1,f=0.5625');运行结果：(3)观察三角波

38、和反三角波序列的时域和幅频特性。实验程序：N=8时：clc;n=0:1:7;% x1为三角波序列时域特性，x2为反三角波时域特性x1=zeros(size(n);x2=zeros(size(n);xa=(n>=0)&(n<=3);x1(xa)=n(xa);x2(xa)=4-n(xa);xa=(n>=4)&(n<=7);x1(xa)=8-n(xa);x2(xa)=n(xa)-4; fp1=fft(x1);fp1=abs(fp1);fp2=fft(x2);fp2=abs(fp2);% fp1为三角波序列幅频特性，fp2为反三角波幅频特性subplot(2,2

39、,1);stem(n,x1,'.');xlabel('n');ylabel('时域特性');title('三角波序列');subplot(2,2,2);stem(n,fp1,'.');xlabel('n');ylabel('幅频特性');title('三角波序列');subplot(2,2,3);stem(n,x2,'.');xlabel('n');ylabel('时域特性');title('反三角波序列'

40、;);subplot(2,2,4);stem(n,fp2,'.');xlabel('n');ylabel('幅频特性');title('反三角波序列');实验程序：N=16时：clc;n=0:1:16;% x1为三角波序列时域特性，x2为反三角波时域特性x1=zeros(size(n);x2=zeros(size(n);xa=(n>=0)&(n<=3);x1(xa)=n(xa);x2(xa)=4-n(xa);xa=(n>=4)&(n<=7);x1(xa)=8-n(xa);x2(xa)=n(x

41、a)-4; xa=(n>=8)&(n<=15);x1(xa)=0;x2(xa)=0; fp1=fft(x1);fp1=abs(fp1);fp2=fft(x2);fp2=abs(fp2);% fp1为三角波序列幅频特性，fp2为反三角波幅频特性subplot(2,2,1);stem(n,x1,'.');xlabel('n');ylabel('时域特性');title('三角波序列');subplot(2,2,2);stem(n,fp1,'.');xlabel('n');ylabel

42、('幅频特性');title('三角波序列');subplot(2,2,3);stem(n,x2,'.');xlabel('n');ylabel('时域特性');title('反三角波序列');subplot(2,2,4);stem(n,fp2,'.');xlabel('n');ylabel('幅频特性');title('反三角波序列');运行结果：(4)一个连续信号含两个频率分量，经采样得x(n)=sin2*0.125n+cos2*(

43、0.125+f)n n=0,1,N-1已知N=16，f分别为1/16和1/64，观察其频谱；当N=128时，f不变，其结果有何不同。实验程序：N=16时：clc;n1=0:1:15;% N=16时，vf1=1/16,vf2=1/64两种情况下vf1=1/16;vf2=1/64;x1=sin(2*pi*0.125.*n1)+cos(2*pi*(0.125+vf1).*n1); %时域特性x2=sin(2*pi*0.125.*n1)+cos(2*pi*(0.125+vf2).*n1);fp1=fft(x1);fp1=abs(fp1); fp2=fft(x2);fp2=abs(fp2);subplo

44、t(2,2,1);stem(n1,x1,'.');xlabel('n');ylabel('时域特性');title('N=16 vf=1/16');subplot(2,2,2);stem(n1,fp1,'.');xlabel('n');ylabel('幅频特性');title('N=16 vf=1/16');subplot(2,2,3);stem(n1,x2,'.');xlabel('n');ylabel('时域特性')

45、;title('N=16 vf=1/64');subplot(2,2,4);stem(n1,fp2,'.');xlabel('n');ylabel('幅频特性');title('N=16 vf=1/64');运行结果：N=128时：clc;n2=0:1:127;% N=32时，vf1=1/16,vf2=1/64两种情况下vf1=1/16;vf2=1/64;x3=sin(2*pi*0.125.*n2)+cos(2*pi*(0.125+vf1).*n2); %时域特性x4=sin(2*pi*0.125.*n2)+cos

46、(2*pi*(0.125+vf2).*n2);fp3=fft(x3);fp3=abs(fp3); %幅频特性fp4=fft(x4);fp4=abs(fp4);subplot(2,2,1);stem(n2,x3,'.');xlabel('n');ylabel('时域特性');title('N=128 vf=1/16');subplot(2,2,2);stem(n2,fp3,'.');xlabel('n');ylabel('幅频特性');title('N=128 vf=1/16&

47、#39;);subplot(2,2,3);stem(n2,x4,'.');xlabel('n');ylabel('幅频特性');title('N=128 vf=1/64');subplot(2,2,4);stem(n2,fp4,'.');xlabel('n');ylabel('时域特性');title('N=128 vf=1/64');运行结果：(5)用FFT分别计算xa(n) (p=8,q=2)和xb(n) (a=0.1,f=0.0625)的16点循环卷积和线性卷积

48、。实验程序：clc;n1=0:1:15;n2=0:1:30;p=8;q=2;a=0.1;f=0.0625;xa=exp(-(n1-p).2)/q);xb=exp(-a.*n1).*sin(2*pi*f.*n1);fa=fft(xa);fb=fft(xb);circle=fa.*fb; %圆周卷积line=conv(xa,xb); %线性卷积subplot(2,2,1);stem(n1,xa,'.');ylabel('时域特性');title('高斯序列xa');subplot(2,2,2);stem(n1,xb,'.');ylab

49、el('时域特性');title('衰减正弦序列xb');subplot(2,2,3);stem(n1,circle,'.');ylabel('幅频特性');title('xa与xb的16点循环卷积');subplot(2,2,4);stem(n2,line,'.');ylabel('幅频特性');title(' xa与xb的线性卷积');运行结果：(6)产生一512点的随机序列xe(n)，并用xc(n)和xe(n)做线性卷积，观察卷积前后xe(n)频谱的变化。要求将

50、xe(n)分成8段，分别采用重叠相加法和重叠保留法。实验程序：clc;n1=0:1:3;n2=4:7;n=0:518;k=1:7;xe=rand(1,512);xc1=n1;xc2=8-n2;xc=xc1,xc2;xe1=k-k;xe2=xe1,xe;%重叠相加法yn=zeros(1,519);yn3=zeros(1,519);for j=0:7 %重叠相加法%重叠保留法 xj=xe(64*j+1:64*(j+1); xj1=xe2(64*j+1:64*j+71); xak=fft(xj,71); xak1=fft(xj1); xck=fft(xc,71); xck1=fft(xc,71);

51、yn1=ifft(xak.*xck); yn2=ifft(xak1.*xck1); %每段卷积的结果 temp=zeros(1,519); temp1=zeros(1,519); temp(64*j+1:64*j+71)=yn1; temp1(64*j+1:64*j+64)=yn2(8:71); yn=yn+temp; yn3=yn3+temp1; %将每段的卷积结果加到yn中end;figure(1);subplot(2,1,1);plot(n,yn);xlabel('n');ylabel('y(n) ');title('xc(n)与xe(n)的线性卷

52、积(重叠相加法)');subplot(2,1,2);plot(n,abs(fft(yn);xlabel('k');ylabel('Y(k)');axis(0,600,0,300);figure(2);subplot(2,1,1);plot(n,yn3);xlabel('n');ylabel('y(n)');title(' xc(n)与xe(n)的线性卷积(重叠保留法)');subplot(2,1,2);plot(n,abs(fft(yn3);xlabel('k');ylabel('Y

53、(k)的幅频特性');axis(0,600,0,300);运行结果：重叠相加法：重叠保留法：(7) 对做16点DCT变换，取前m个系数做16点IDCT，得到不同m下的，求归一化的均方误差，绘出随m变化的关系并观察不同的m下的。实验结果：实验四 IIR数字滤波器的设计一 、实验目的 (1) 掌握双线性变换法及脉冲相应不变法设计IIR 数字滤波器的具体设计方法及其原理，熟 悉用双线性变换法及脉冲响应不变法设计低通、高通和带通IIR 数字滤波器的计算机编程。 (2) 观察双线性变换及脉冲响应不变法设计的滤波器的频域特性，了解双线性变换法及脉冲 响应不变法的特点。 (3) 熟悉巴特沃思滤波器、

54、切比雪夫滤波器和椭圆滤波器的频率特性。 二 、实验内容(1)，,，设计一切比雪夫高通滤波器，观察其通带损耗和阻带衰减是否满足要求。实验程序：clc;fc=300;Ap=0.8;fr=200;At=20;T=1000;wc=2*T*tan(2*pi*fc/(2*T);wt=2*T*tan(2*pi*fr/(2*T);N,wn=cheb1ord(wc,wt,Ap,At,'s');B,A=cheby1(N,0.5,wn,'high','s');num,den=bilinear(B,A,T);h,w=freqz(num,den);f=w/(2*pi)*T

55、;plot(f,20*log10(abs(h);axis(0,500,-80,10);grid;xlabel('频率/Hz');ylabel('幅度/dB');title('切比雪夫高通滤波器');运行结果：分析：f=200Hz时阻带衰减大于30dB，通过修改axis(0,fs/2,-80,10)为axis(200,fs/2,-1,1)发现通带波动rs满足<0.8。(2)，,，分别用脉冲响应不变法及双线性变换法设计一巴特沃思数字低通滤波器，观察所设计数字滤波器的幅频特性曲线，记录带宽和衰减量，检查是否满足要求。比较这两种方法的优缺点。实验程

56、序：clc;fs=1000;fc=200;fr=300;T=0.001;wp1=2*pi*fc;wr1=2*pi*fr;N1,wn1=buttord(wp1,wr1,1,25,'s');B1,A1=butter(N1,wn1,'s');num1,den1=impinvar(B1,A1,fs);%脉冲响应不变法h1,w=freqz(num1,den1);wp2=2*fs*tan(2*pi*fc/(2*fs);wr2=2*fs*tan(2*pi*fr/(2*fs);N2,wn2=buttord(wp2,wr2,1,25,'s');B2,A2=butt

57、er(N2,wn2,'s');num2,den2=bilinear(B2,A2,fs);%双线性变换法h2,w=freqz(num2,den2);f=w/(2*pi)*fs;plot(f,20*log10(abs(h1),'-.',f,20*log10(abs(h2),'-');axis(0,500,-80,10);grid;xlabel('频率/Hz');ylabel('幅度/dB');title('巴特沃思数字低通滤波器');legend('脉冲响应不变法','双线性变换

58、法');运行结果：(3)利用双线性变换法分别设计满足下列指标的巴特沃思型、切比雪夫型和椭圆型数字低通滤波器，并作图验证设计结果: f =1.2kHz，0.5dB，f =2kHz,At40dB，f =8kHz。 比较这三种滤波器的阶数。实验程序：clc;wc=2*pi*1200;wr=2*pi*2000;rp=0.5;rs=40;fs=8000; w1=2*fs*tan(wc/(2*fs); w2=2*fs*tan(wr/(2*fs); Nb,wn=buttord(w1,w2,rp,rs,'s') %巴特沃思B,A=butter(Nb,wn,'s'); n

59、um1,den1=bilinear(B,A,fs); h1,w=freqz(num1,den1); Nc,wn=cheb1ord(w1,w2,rp,rs,'s') %切比雪夫B,A=cheby1(Nc,rp,wn,'s'); num2,den2=bilinear(B,A,fs); h2,w=freqz(num2,den2); Ne,wn=ellipord(w1,w2,rp,rs,'s') %椭圆型B,A=ellip(Ne,rp,rs,wn,'low','s'); num3,den3=bilinear(B,A,fs)

60、; h3,w=freqz(num3,den3); f=w/pi*4000; plot(f,20*log10(abs(h1),'-',f,20*log10(abs(h2),'-',f,20*log10(abs(h3),':'); axis(0,3000,-100,10);grid; xlabel('频率/Hz'); ylabel('幅度/dB'); title('三种数字低通滤波器'); legend('巴特沃思数字低通滤波器','切比雪夫数字低通滤波器','椭圆数字低通滤波器',3);运行结果：(4)分别用脉冲响应不变法和双线性变换法设计一巴特沃思型数字带通滤波器，已知f=30kHz，其等效的模拟滤波器指标为 3dB，2kHzf3kHz；At5dB，f 6kHz；At20dB，f1.5kHz。实验程序：clc;wc=2*pi*2000,2*pi*3000;wr=2*pi*1500,2*pi*6000;rp=3;rs=20;fs=30000; N,wn=buttord(wc,wr,rp,rs,'s'); B,A=butter(N,wn,