(成都理工大学)数字信号处理实验

1、本科生实验报告实验课程 数字信号处理 学院名称 信息科学与技术学院 专业名称 物联网工程 学生姓名 曹林鑫 学生学号 201413060301 指导教师 罗耀耀 实验地点 6B607 实验成绩 二一六年 十月 二一六年 十二月实验一 熟悉Matlab环境一、实验目的1. 熟悉MATLAB的主要操作命令。2. 学会简单的矩阵输入和数据读写。3. 掌握简单的绘图命令。4. 用MATLAB编程并学会创建函数。5. 观察离散系统的频率响应。二、实验内容（1）数组的加、减、乘、除和乘方运算。输入A=1 2 3 4，B=3 4 5 6，求C=A+B，D=A-B，E=A.*B，F=A./B，G=A.B并用s

2、tem语句画出A、B、C、D、E、F、G。clear all;a=1 2 3 4; b=3 4 5 6;c=a+b;d=a-b;e=a.*b;f=a./b;g=a.b;n=1:4;subplot(4,2,1);stem(n,a);xlabel('n');xlim(0 5);ylabel('A');subplot(4,2,2);stem(n,b);xlabel('n');xlim(0 5);ylabel('B');subplot(4,2,3);stem(n,c);xlabel('n');xlim(0 5);ylabe

3、l('C');subplot(4,2,4);stem(n,d);xlabel('n');xlim(0 5);ylabel('D');subplot(4,2,5);stem(n,e);xlabel('n');xlim(0 5);ylabel('E');subplot(4,2,6);stem(n,f);xlabel('n');xlim(0 5);ylabel('F');subplot(4,2,7);stem(n,g);xlabel('n');xlim(0 5);ylabe

4、l('G');（2）用MATLAB实现下列序列：1) x(n)=0.8n 0n15 2) x(n)=e(0.2+3j)n 0n15 3) x(n)=3cos(0.125n+0.2)+2sin(0.25n+0.1) 0n15 4) 将3)中的x(n)扩展为以16为周期的函数x16(n)=x(n+16),绘出四个周期。5) 将3)中的x(n)扩展为以10为周期的函数x10(n)=x(n+10),绘出四个周期。解：1) x(n)=0.8n 0n15clear all;xa=0.8.N;figure;subplot(2,1,1);stem(N,xa);xlabel('n'

5、;);xlim(0 16);ylabel('xa'); 2) x(n)=e(0.2+3j)n 0n15xb=exp(0.2+3*j)*N);subplot(2,1,2);stem(N,xb);xlabel('n');xlim(0 16);ylabel('xb');figure; 3) x(n)=3cos(0.125n+0.2)+2sin(0.25n+0.1) 0n15xc=3*cos(0.125*pi*N+0.2*pi)+2*sin(0.25*pi*N+0.1*pi);subplot(3,1,1);stem(N,xc);xlabel('n

6、');xlim(0 16);ylabel('xc');4) 将3)中的x(n)扩展为以16为周期的函数x16(n)=x(n+16),绘出四个周期。 k=0:3;m=0;for i=1:4 for j=1:16 m=m+1; n(m)=N(j)+16*k(i); x16(m)=3*cos(0.125*pi*n(m)+0.2*pi)+2*sin(0.25*pi*n(m)+0.1*pi); endend subplot(3,1,2);stem(n,x16);xlabel('n');ylabel('x16');5) 将3)中的x(n)扩展为以10

7、为周期的函数x10(n)=x(n+10),绘出四个周期。 for j=1:10 x10(j)=x16(j); end for i=1:3 for m=1:10 x10(i*10+m)=x10(m); end endn=1:40;subplot(3,1,3);stem(n,x10);xlabel('n');ylabel('x10');（3）x(n)=1,-1,3,5,产生并绘出下列序列的样本: 1) x1(n)=2x(n+2)-x(n-1)-2x(n)2) clear alln=1:4;T=4;x=1 -1 3 5;x(5:8)=x(1:4);subplot(2,

8、1,1);stem(1:8,x);grid;for i=1:4 if i-1<0 x1(i)=2*x(i+2)-x(i-1)-2*x(i); else x1(i)=2*x(i+2)-x(i-1+T)-2*x(i); endendx1(5:8)=x1(1:4);subplot(2,1,2);stem(1:8,x1);grid; （4）绘出下列时间函数的图形，对x轴、y轴以及图形上方均须加上适当的标注： a) x(t)=sin(2t) 0t10sb) x(t)=cos(100t)sin(t) 0t4sta=0:0.05:10;xa=sin(2*pi*ta);subplot(2,1,1);pl

9、ot(ta,xa);xlabel('t');ylabel('幅度');tb=0:0.01:4;xb=cos(100*pi*tb).*sin(pi*tb);subplot(2,1,2);plot(tb,xb);xlabel('t');ylabel('幅度'); （5）编写函数stepshift(n0,n1,n2)实现u(n-n0),n1<n0<n2,绘出该函数的图形，起点为n1，终点为n2。n0=5;ns=1;nf=10;%ns为起点;nf为终点；在=n=n0处生成单位阶跃序列n=ns:nf;x=(n-n0)>=0

10、;stem(n,x); （6）给一定因果系统求出并绘制H(z)的幅频响应与相频响应。clear all;b=1,sqrt(2),1;a=1,-0.67,0.9;h,w=freqz(b,a);am=20*log10(abs(h);subplot(2,1,1);plot(w,am);ph=angle(h);subplot(2,1,2);plot(w,ph); （7）计算序列8 -2 -1 2 3和序列2 3 -1 -3的离散卷积，并作图表示卷积结果。clear all;a=8 -2 -1 2 3;b=2 3 -1 -3;c=conv(a,b); %计算卷积M=length(c)-1;n=0:1:M

11、;stem(n,c);xlabel('n');ylabel('幅度'); （8）求以下差分方程所描述系统的单位脉冲响应h(n),0n50 y(n)+0.1y(n-1)-0.06y(n-2)=x(n)-2x(n-1)clear all;N=50;a=1 -2;b=1 0.1 -0.06;x=1 zeros(1,N-1);k=0:1:N-1;y=filter(a,b,x);stem(k,y);xlabel('n');ylabel('幅度 '); 3、 思考题1. 比较实验内容第2题中（4）和（5）两小题的结果，试说明对于周期性信号，应

12、当如何采样，才能保证周期扩展后与原信号保持一致？ 答：对于周期性信号，在进行采样时，其采样周期必须满足采样定理，即采样频率应该大于信号最高频率的两倍，这样才能避免迭混，以便采样后仍能准确的恢复原信号。2. 对于有限长序列，如何用MATLAB计算其DTFT？ 答：用函数Freqz可以计算序列在给定的离散频率点上的DTFT，函数freqz返回的频率响应值为向量H。在H = freqz(num,den,w)中，0到&#960；之间指定的频率集由向量w给出。freqz函数的自变量k就是频率点的总数。3. 对于由两个子系统级联或并联的系统，如何用MATLAB计算它们的幅频响应与相频响应？ 系统的

13、级联或并联实现涉及到了因式分解。在MATLAB中，我们可以用函数roots来实现多项式的因式分解。 实验二 快速Fourier变换(FFT)及其应用一、实验目的1. 在理论学习的基础上，通过本实验，加深对FFT的理解，熟悉MATLAB中的有关函数。2. 熟悉应用FFT对典型信号进行频谱分析的方法。3. 了解应用FFT进行信号频谱分析过程中可能出现的问题，以便在实际中正确应用FFT。4. 熟悉应用FFT实现两个序列的线性卷积和相关的方法。二、实验原理与方法  在各种信号序列中，有限长序列信号处理占有很重要地位，对有限长序列，我们可以使用离散Fouier变换(DFT)。这一变换不但可以很

14、好的反映序列的频谱特性，而且易于用快速算法在计算机上实现，当序列x(n)的长度为N时，它的DFT定义为：反变换为：    有限长序列的DFT是其Z变换在单位圆上的等距采样，或者说是序列Fourier变换的等距采样，因此可以用于序列的谱分析。    FFT并不是与DFT不同的另一种变换，而是为了减少DFT运算次数的一种快速算法。它是对变换式进行一次次分解，使其成为若干小点数的组合，从而减少运算量。常用的FFT是以2为基数的，其长度 。它的效率高，程序简单，使用非常方便，当要变换的序列长度不等于2的整数次方时，为了使用以2为基数的FFT，

15、可以用末位补零的方法，使其长度延长至2的整数次方。     (一)、在运用DFT进行频谱分析的过程中可能产生三种误差：    (1)    混叠    序列的频谱时被采样信号的周期延拓，当采样速率不满足Nyquist定理时，就会发生频谱混叠，使得采样后的信号序列频谱不能真实的反映原信号的频谱。避免混叠现象的唯一方法是保证采样速率足够高，使频谱混叠现象不致出现，即在确定采样频率之前，必须对频谱的性质有所了解，在一般情况下，为了保证高于折叠频率的分量不会出现，在采样前，先用低通模

16、拟滤波器对信号进行滤波。    (2)    泄漏    实际中我们往往用截短的序列来近似很长的甚至是无限长的序列，这样可以使用较短的DFT来对信号进行频谱分析，这种截短等价于给原信号序列乘以一个矩形窗函数，也相当于在频域将信号的频谱和矩形窗函数的频谱卷积，所得的频谱是原序列频谱的扩展。泄漏不能与混叠完全分开，因为泄漏导致频谱的扩展，从而造成混叠。为了减少泄漏的影响，可以选择适当的窗函数使频谱的扩散减至最小。    (3)    栅栏效应

17、0;   DFT是对单位圆上Z变换的均匀采样，所以它不可能将频谱视为一个连续函数，就一定意义上看，用DFT来观察频谱就好像通过一个栅栏来观看一个图景一样，只能在离散点上看到真实的频谱，这样就有可能发生一些频谱的峰点或谷点被“尖桩的栅栏”所拦住，不能别我们观察到。    减小栅栏效应的一个方法就是借助于在原序列的末端填补一些零值，从而变动DFT的点数，这一方法实际上是人为地改变了对真实频谱采样的点数和位置，相当于搬动了每一根“尖桩栅栏”的位置，从而使得频谱的峰点或谷点暴露出来。    (二)、用FFT计算线性卷积&#

18、160;   用FFT可以实现两个序列的圆周卷积。在一定的条件下，可以使圆周卷积等于线性卷积。一般情况，设两个序列的长度分别为N1和N2，要使圆周卷积等于线性卷积的充要条件是FFT的长度NN1N2对于长度不足N的两个序列，分别将他们补零延长到N。    当两个序列中有一个序列比较长的时候，我们可以采用分段卷积的方法。有两种方法：重叠相加法。将长序列分成与短序列相仿的片段，分别用FFT对它们作线性卷积，再将分段卷积各段重叠的部分相加构成总的卷积输出。  重叠保留法。这种方法在长序列分段时，段与段之间保留有互相重叠的部分，在构成总的卷积输

19、出时只需将各段线性卷积部分直接连接起来，省掉了输出段的直接相加。     (三)、用周期图法(平滑周期图的平均法)对随机信号作谱分析    实际中许多信号往往既不具有有限能量，由非周期性的。无限能量信号的基本概念是随机过程，也就是说无限能量信号是一随机信号。周期图法是随机信号作谱分析的一种方法，它特别适用于用FFT直接计算功率谱的估值。    将长度为N的实平稳随机序列的样本x(n)再次分割成K段，每段长度为L,即L=N/K。每段序列仍可表示为：xi(n)=x(n+(i-1)L)，0nL-1，1iK 但是

20、这里在计算周期图之前，先用窗函数w(n)给每段序列xi(n)加权，K个修正的周期图定义为其中U表示窗口序列的能量，它等于在此情况下，功率谱估计量可表示为 三、实验内容及步骤    实验中用到的信号序列：a)       Gaussian序列b)      衰减正弦序列c)      三角波序列d)      反三角波序列  &#

21、160; 上机实验内容：(1)、观察高斯序列的时域和幅频特性，固定信号xa(n)中参数p=8，改变q的值，使q分别等于2，4，8，观察它们的时域和幅频特性，了解当q取不同值时，对信号序列的时域幅频特性的影响；固定q=8，改变p，使p分别等于8，13，14，观察参数p变化对信号序列的时域及幅频特性的影响，观察p等于多少时，会发生明显的泄漏现象，混叠是否也随之出现？记录实验中观察到的现象，绘出相应的时域序列和幅频特性曲线。% %定义高斯序列% function Xa,Fa =gauss(p,q)% n=0:15;% Xa(n+1)=exp(-(n+1-p).2./q);% F=fft(Xa);%

22、Fa=abs(F);clear all;% p=8,q=2 %Xa1,Fa1= gauss(8,2);k=0:15;subplot(5,2,1);plot(k,Xa1);xlabel('n');ylabel('时域特性');text(10,0.5,'p=8,q=2');subplot(5,2,2);plot(k,Fa1);xlabel('n');ylabel('幅频特性');text(8,3,'p=8,q=2');% p=8,q=4 %Xa2,Fa2= gauss(8,4);subplot(5,2,

23、3);plot(k,Xa2);xlabel('n');ylabel('时域特性');text(10,0.5,'p=8,q=4');subplot(5,2,4);plot(k,Fa2);xlabel('n');ylabel('幅频特性');text(8,3,'p=8,q=4');% p=8,q=8 %Xa3,Fa3= gauss(8,8);subplot(5,2,5);plot(k,Xa3);xlabel('n');ylabel('时域特性');text(10,0.5,

24、'p=8,q=8');subplot(5,2,6);plot(k,Fa3);xlabel('n');ylabel('幅频特性');text(8,3,'p=8,q=8');% p=13,q=8 %Xa4,Fa4= gauss(13,8);subplot(5,2,7);plot(k,Xa4);xlabel('n');ylabel('时域特性');text(10,0.5,'p=13,q=8');subplot(5,2,8);plot(k,Fa4);xlabel('n');y

25、label('幅频特性');text(8,3,'p=13,q=8');% p=14,q=8 %Xa5,Fa5= gauss(14,8);subplot(5,2,9);plot(k,Xa5);xlabel('n');ylabel('时域特性');text(10,0.5,'p=14,q=8');subplot(5,2,10);plot(k,Fa5);xlabel('n');ylabel('幅频特性');text(8,3,'p=14,q=8');(2)、观察衰减正弦序列xb

26、(n)的时域和幅频特性，a=0.1，f=0.0625，检查谱峰出现位置是否正确，注意频谱的形状，绘出幅频特性曲线，改变f，使f分别等于0.4375和0.5625，观察这两种情况下，频谱的形状和谱峰出现位置，有无混叠和泄漏现象？说明产生现象的原因。% 定义衰减正弦序列% function Xb,Fb = downsin(a,f)% n=0:15;% Xb(n+1)=exp(-a.*n).*sin(2*pi*f.*n); %自然对数的底:e=:2.71828 18284 59045 23536 % F = fft(Xb);% Fb=abs(F);clear all;k=0:15;% a=0.1,f

27、=0.0.0625 %Xb,Fb=downsin(0.1,0.0625);subplot(3,2,1); plot(k,Xb);xlabel('n');ylabel('时域特性');text(8,0.5,'a=0.1,f=0.0625');subplot(3,2,2); plot(k,Fb);xlabel('n');ylabel('幅值特性');text(10,3,'a=0.1,f=0.0625');% a=0.1,f=0.4375 %Xb1,Fb1=downsin(0.1,0.4375);subp

28、lot(3,2,3); plot(k,Xb1);xlabel('n');ylabel('时域特性');text(8,0.5,'a=0.1,f=0.4375');subplot(3,2,4); plot(k,Fb1);xlabel('n');ylabel('幅值特性');text(10,3,'a=0.1,f=0.4375');% a=0.1,f= 0.5625 %Xb2,Fb2=downsin(0.1,0.5625);subplot(3,2,5); plot(k,Xb2);xlabel('n&

29、#39;);ylabel('时域特性');text(8,0.5,'a=0.1,f=0.5625');subplot(3,2,6); plot(k,Fb2);xlabel('n');ylabel('幅值特性');text(10,3,'a=0.1,f=0.5625');    (3)、观察三角波和反三角波序列的时域和幅频特性，用N=8点FFT分析信号序列xc(n)和xd(n)的幅频特性，观察两者的序列形状和频谱曲线有什么异同？绘出两序列及其幅频特性曲线。在xc(n)和xd(n)末尾补零，用

30、N=16点FFT分析这两个信号的幅频特性，观察幅频特性发生了什么变化？两情况的FFT频谱还有相同之处吗？这些变化说明了什么？clear all;n=0:3;k=1:8;%定义三角波序列Xc(n+1) = n;Xc(n+5) =4-n;%定义反三角波序列Xd(n+1) = 4-n;Xd(n+5) =n;% 三角波特性 %subplot(2,2,1);plot(k-1,Xc);xlabel('n');ylabel('时域特性');text(1,3,'三角波');subplot(2,2,2);plot(k-1,abs(fft(Xc);xlabel(&#

31、39;k');ylabel('幅频特性');text(4,10,'三角波');% 反三角波特性 %subplot(2,2,3);plot(k-1,Xd);xlabel('n');ylabel('时域特性');text(3,3,'反三角波');subplot(2,2,4);plot(k-1,abs(fft(Xd);xlabel('k');ylabel('幅频特性');text(4,10,'反三角波');%末尾补0，计算32点FFTXc(9:32)=0;Xd(9:

32、32)=0;k=1:32;figure;% 三角波特性 %subplot(2,2,1);plot(k-1,Xc);xlabel('n');ylabel('时域特性');text(1,3,'三角波');subplot(2,2,2);plot(k-1,abs(fft(Xc);xlabel('k');ylabel('幅频特性');text(4,10,'三角波');% fan三角波特性 %subplot(2,2,3);plot(k-1,Xd);xlabel('n');ylabel('

33、时域特性');text(3,3,'反三角波');subplot(2,2,4);plot(k-1,abs(fft(Xd);xlabel('k');ylabel('幅频特性');text(4,10,'反三角波');    (4)、一个连续信号含两个频率分量，经采样得x(n)=sin2*0.125n+cos2*(0.125+f)n n=0,1,N-1已知N=16，f分别为1/16和1/64，观察其频谱；当N=128时，f不变，其结果有何不同，为什么？clear all;% N = 16 %N=16;

34、detf=1/16;n=0:N-1;x1(n+1)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n);detf = 1/64;x2(n+1)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n);% N = 16,detf = 1/16 %subplot(2,2,1);stem(n,x1);hold;plot(n,x1);xlabel('n');ylabel('时域特性');text(6,1,'N=16,detf=1/16');subplot(2,2,2);stem(n,

35、abs(fft(x1);xlabel('n');ylabel('幅值特性');text(6,4,'N=16,detf=1/16');% N = 16,detf = 1/64 %subplot(2,2,3);stem(n,x2);xlabel('n');ylabel('时域特性');text(6,1,'N=16,detf=1/64');subplot(2,2,4);stem(n,abs(fft(x2);xlabel('n');ylabel('幅值特性');text(6,

36、4,'N=16,detf=1/64');% N = 128 %N=128;detf=1/16;n=0:N-1;x3(n+1)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n);detf = 1/64;x4(n+1)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n);% N = 128,detf = 1/16 %figure;subplot(2,2,1);stem(n,x3);xlabel('n');ylabel('时域特性');axis(0 128 -2 2);tex

37、t(6,1.5,'N=128,detf=1/16');subplot(2,2,2);stem(n,abs(fft(x3);xlabel('n');ylabel('幅值特性');axis(0 128 -10 70);text(40,60,'N=128,detf=1/16');% N = 128,detf = 1/64 %subplot(2,2,3);stem(n,x3);xlabel('n');ylabel('时域特性');axis(0 128 -2 2);text(6,1.5,'N=128,

38、detf=1/16');subplot(2,2,4);stem(n,abs(fft(x4);xlabel('n');ylabel('幅值特性');axis(0 128 -10 70);text(40,60,'N=128,detf=1/16');(5)、用FFT分别实现xa(n)（p8，q2）和 xb(n)（a0.1，f0.0625）的16点圆周卷积和线性卷积。clear all;N=16;n=0:N-1;p=8;q=2;Xa(n+1)=exp(-(n-p).2./q);a=0.1;f=0.0625;Xb(n+1)=exp(-a.*n).*

39、sin(2*pi*f.*n);%16点循环卷积Fa=fft(Xa); Fb=fft(Xb);Fx=Fa.*Fb;X51=ifft(Fx);stem(n,X51);%16点线性卷积Xa(N+1:2*N-1)=0;Xb(N+1:2*N-1)=0;Fa=fft(Xa); Fb=fft(Xb);Fc=Fa.*Fb;X52=ifft(Fc);figure;stem(1:2*N-1,X52);（6）用FFT分别计算xa(n)(p=8,q=2)和xb(n)(a=0.1,f=0.0625)的16点循环相关和线性相关，问一共有多少种结果，他们之间有何异同点。clear all;N=16;n=0:N-1;p=8;

40、q=2;Xa(n+1)=exp(-(n-p).2./q);a=0.1;f=0.0625;Xb(n+1)=exp(-a.*n).*sin(2*pi*f.*n); N=length(Xa); %16点循环相关Fa=fft(Xa,2*N); Fb=fft(Xb,2*N);Fx=conj(Fa).*Fb;X71=real(ifft(Fx);X71=X71(N+2:2*N) X71(1:N);n=(-N+1):(N-1); stem(n,X71);%16点线性相关Xa(N+1:2*N-1)=0;Xb(N+1:2*N-1)=0;Fa=fft(Xa); Fb=fft(Xb);Fc=conj(Fa).*Fb;

41、X72=real(ifft(Fc);figure;stem(1:2*N-1,X72);（7）用FFT分别计算xa(n)(p=8,q=2)和xb(n)(a=0.1,f=0.0625)的自相关函数。clear all;N=16;n=0:N-1;p=8;q=2;Xa(n+1)=exp(-(n-p).2./q);a=0.1;f=0.0625;Xb(n+1)=exp(-a.*n).*sin(2*pi*f.*n); %自然对数的底:e=:2.71828 18284 59045 23536 N=length(Xa); % Xa(n) 16点自相关Fa=fft(Xa,2*N); Fb=fft(Xb,2*N);

42、F1=conj(Fa).*Fa;X81=real(ifft(F1);X81=X81(N+2:2*N) X81(1:N);n=(-N+1):(N-1); subplot(2,1,1);stem(n,X81);xlabel('n'); ylabel('幅度');% Xb(n) 16点自相关Fb=fft(Xb,2*N);F2=conj(Fb).*Fb;X82=real(ifft(F2);X82=X82(N+2:2*N) X82(1:N);% n=(-N+1):(N-1); subplot(2,1,2);stem(n,X82);xlabel('n');

43、ylabel('幅度');四、思考题 实验中的信号序列和，在单位圆上的Z变换频谱和会相同吗？如果不同，说出哪一个低频分量更多一些，为什么？答：xc和xd的幅频和相频完全一样。两者为同一个信号序列。只是采样周期的起始位置不同，相差相位180度而已。但由FFT对其所得的图像去不同。 因此，两者在单位圆上的z变换频谱完全一样；并且直流分量大小完全一样。实验三 IIR数字滤波器的设计一、实验目的 1. 掌握双线性变换法及脉冲响应不变法设计IIR数字滤波器的具体设计方法及其原理，熟悉用双线性变换法及脉冲响应不变法设计低通、高通和带通IIR数字滤波器的MATLAB编程。 2. 观

44、察双线性变换及脉冲响应不变法设计的滤波器的频域特性，了解双线性变换法及脉冲响应不变法的特点。 3. 熟悉Butterworth滤波器、Chebyshev滤波器和椭圆滤波器的频率特性。 2、 实验原理 1 脉冲响应不变法用数字滤波器的单位脉冲响应序列模仿模拟滤波器的冲激响应,让正好等于的采样值，即，其中为采样间隔，如果以及分别表示的拉式变换及的Z变换，则2双线性变换法    S平面与z平面之间满足以下映射关系：s平面的虚轴单值地映射于z平面的单位圆上，s平面的左半平面完全映射到z平面的单位圆内。双线性变换不存在混叠问题。    双线性变

45、换是一种非线性变换 ，这种非线性引起的幅频特性畸变可通过预畸而得到校正。    IIR低通、高通、带通数字滤波器设计采用双线性原型变换公式：变换类型           变换关系式     备   注 低通 高通   带通 ：带通的上下边带临界频率 以低通数字滤波器为例，将设计步骤归纳如下：1. 确定数字滤波器的性能指标：通带临界频率、阻带临界频率、通带波动、阻带内的最小衰减、采样周期、采样频率 ；2.  

46、确定相应的数字角频率，； 3.  计算经过预畸的相应模拟低通原型的频率，；4. 根据c和r计算模拟低通原型滤波器的阶数N，并求得低通原型的传递函数； 5. 用上面的双线性变换公式代入，求出所设计的传递函数； 6.  分析滤波器特性，检查其指标是否满足要求。 三、实验内容：(1)，,;设计一切比雪夫高通滤波器，观察其通带损耗和阻带衰减是否满足要求。wc=2*1000*tan(2*pi*300/(2*1000);wt=2*1000*tan(2*pi*200/(2*1000); %对临界频率进行预畸变；N,wn=cheb1ord(wc,wt,0.8,20,'s')

47、; %给定通带边界频率wc和阻带边界频率wr、通带波动0.8dB和最小阻带衰减20dB，求满足指标的模拟滤波器的最低阶数N和通带边界频率wn；B,A=cheby1(N,0.8,wn,'high','s'); %给定模拟滤波器的最低阶数N、通带边界频率wn和通带波动0.8dB，设计N阶模拟巴特沃斯高通滤波器，B和A分别表示系统函数的分子和分母多项式的系数；num,den=bilinear(B,A,1000); %给定模拟滤波器系统函数H(s)=B(s)/A(s)和采样频率1000Hz，根据双线性变换法求出数字滤波器的系统函数H(z)=B(z)/A(z)，num和d

48、en分别为数字滤波器的分子分母多项式系数；h,w=freqz(num,den); %计算以(num,den)为参数的滤波器的频率向量w和复频率响应向量h；f=w/pi*500;plot(f,20*log10(abs(h);axis(0,500,-80,10);grid;xlabel('频率/Hz');ylabel('幅度/dB');分析：f=200Hz时阻带衰减大于30dB，通过修改axis(0,fs/2,-80,10)为axis(200,fs/2,-1,1)发现通带波动rs满足<0.8。（2），,;分别用脉冲响应不变法及双线性变换法设计一巴特沃思数字低通

49、滤波器，观察所设计数字滤波器的幅频特性曲线，记录带宽和衰减量，检查是否满足要求。比较这两种方法的优缺点。Wc=2*pi*200;Wr=2*pi*300;Rp=1;Rs=25;N,Wn=buttord(Wc,Wr,Rp,Rs,s); B,A=butter(N,Wn,s); num1,den1=impinvar(B,A,1000); h1,w=freqz(num1,den1);B,A=butter(N,Wn,s);num2,den2=bilinear(B,A,1000); h2,w=freqz(num2,den2);f=w/pi*500;plot(f,abs(h1),-.,f,abs(h2),-)

50、;grid;xlabel(频率/Hz);ylabel(幅度);bz1 =0.0000 0.0002 0.0153 0.0995 0.1444 0.0611 0.0075 0.0002 0.0000 0az1 =1.0000 -1.9199 2.5324 -2.2053 1.3868 -0.6309 0.2045 -0.0450 0.0060 -0.0004因此脉冲响应不变法的系统函数为:bz2 =0.0179 0.1072 0.2681 0.3575 0.2681 0.1072 0.0179az2 =1.0000 -0.6019 0.9130 -0.2989 0.1501 -0.0208 0.

51、0025因此双线性变换法的系统函数为:分析：脉冲响应不变法的N=9，双线性变换法的N=6，由图知它们都满足要求，但脉冲响应的衰减较快，双线性变换的过渡带窄一些，且阶数比脉冲小，容易实现。（3）利用双线性变换法分别设计满足下列指标的巴特沃思型、切比雪夫型数字低通滤波器，并作图验证设计结果：fc=1.2kHz，0.5dB，fr2kHz, At40dB,fs=8kHz。比较这三种滤波器的阶数。wc=2*1000*tan(2*pi*1200/(2*8000);wt=2*1000*tan(2*pi*2000/(2*8000); %对临界频率进行预畸变；N,wn= buttord (wc,wt,0.5,4

52、0,'s'); %给定通带边界频率wc和阻带边界频率wr、通带波动0.5dB和最小阻带衰减40dB，求满足指标的模拟滤波器的最低阶数N和通带边界频率wn；B,A= butter (N, wn, 's'); %给定模拟滤波器的最低阶数N、通带边界频率wn和通带波动0.5dB，设计N阶模拟巴特沃斯低通滤波器，B和A分别表示系统函数的分子和分母多项式的系数；num,den=bilinear(B,A,8000); %给定模拟滤波器系统函数H(s)=B(s)/A(s)和采样频率8000Hz，根据双线性变换法求出数字滤波器的系统函数H(z)=B(z)/A(z)，num和de

53、n分别为数字滤波器的分子分母多项式系数；h,w=freqz(num,den); %计算以(num,den)为参数的滤波器的频率向量w和复频率响应向量h；f=w/pi*500;plot(f,20*log10(abs(h);axis(0,400,-200,10);grid;xlabel('频率/Hz');ylabel('幅度/dB');title（巴特沃思滤波器）;wc=2*1000*tan(2*pi*1200/(2*8000);wt=2*1000*tan(2*pi*2000/(2*8000); %对临界频率进行预畸变；N,wn=cheb1ord(wc,wt,0.5

54、,40,'s'); %给定通带边界频率wc和阻带边界频率wr、通带波动0.5dB和最小阻带衰减40dB，求满足指标的模拟滤波器的最低阶数N和通带边界频率wn；B,A=cheby1(N,0.5,wn, 's'); %给定模拟滤波器的最低阶数N、通带边界频率wn和通带波动0.5dB，设计N阶模拟切比雪夫低通滤波器，B和A分别表示系统函数的分子和分母多项式的系数；num,den=bilinear(B,A,8000); %给定模拟滤波器系统函数H(s)=B(s)/A(s)和采样频率8000Hz，根据双线性变换法求出数字滤波器的系统函数H(z)=B(z)/A(z)，num

55、和den分别为数字滤波器的分子分母多项式系数；h,w=freqz(num,den); %计算以(num,den)为参数的滤波器的频率向量w和复频率响应向量h；f=w/pi*500;plot(f,20*log10(abs(h);axis(0,400,-200,10);grid;xlabel('频率/Hz');ylabel('幅度/dB');title（切比雪夫滤波器）;分析：切比雪夫滤波器的阶数为5；巴特沃思滤波器的阶数为9。（4）分别用脉冲响应不变法及双线性不变法设计一巴特沃思型数字带通滤波器，已知fs=30kHz，其等效的模拟滤波器指标为3dB，2kHzf3k

56、Hz；At5dB，f6kHz；At20dB，f1.5kHz。w1=2*30000*tan(2*pi*2000/(2*30000);w2=2*30000*tan(2*pi*3000/(2*30000);wr=2*30000*tan(2*pi*6000/(2*30000); %对临界频率进行预畸变；N,wn=buttord(w1 w2,1 wr,3,5,'s'); %给定通带边界频率带w1 w2和阻带边界频率带1 wr、通带波动3dB和最小阻带衰减5dB，求满足指标的模拟滤波器的最低阶数N和3dB边界频率wn；B,A=butter(N,wn,'s'); %给定模拟滤

57、波器的最低阶数N和3dB边界频率wn，设计N阶模拟巴特沃斯低通滤波器，B和A分别表示系统函数的分子和分母多项式的系数；num,den=bilinear(B,A,30000); %给定模拟滤波器系统函数H(s)=B(s)/A(s)和采样频率30000Hz，根据双线性变换法求出数字滤波器的系统函数H(z)=B(z)/A(z)，num和den分别为数字滤波器的分子分母多项式系数；h,w=freqz(num,den); %计算以(num,den)为参数的滤波器的频率向量w和复频率响应向量h；f=w/pi*200;plot(f,20*log10(abs(h);axis(-40,40,-30,10); %给定幅频响应的频率和幅度的范围；grid;xlabel('频率/kHz');ylabel('幅度/dB');4、 思考题1. 双线性变换法中和之间的关系是非线性的，在实验中你注意到这种非线性关系了吗？从那几种数字滤波器的幅频特性曲线中可以观察到这种非线性关系？ 答：在双线性变化法中，模拟频率与数字频率不再是线性关系，所以一个线性相位模拟滤波器经双线性变换后，得到的数字滤波器不再保持原有的线性相位了，在每一幅使用了双线性变换的图中，可以看到在采样频率一半处，幅度为零，这