电气与计算机工程学院

ADDINCNKISM.UserStyle中山大学南方学院中山大学南方学院电气与计算机工程学院《数字信号处理与实践》实验讲义2018年8月1日目录实验一 3实验二 8实验三 15实验四 22实验五 32实验六 39实验一用MATLAB产生时域离散信号一、实验目的了解常用时域离散信号及其特点掌握用MATLAB产生时域离散信号的方法二、实验原理1、时域离散信号的概念在时间轴的离散点上取值的信号，称为离散时间信号。通常，离散时间信号用x(n)表示，其幅度可以在某一范围内连续取值。由于信号处理设备或装置（如计算机、专用的信号处理芯片等）均以有限位的二进制数来表示信号的幅度，因此，信号的幅度也必须离散化。我们把时间和幅度均取离散值的信号称为时域离散信号或数字信号。在MATLAB语言中，时域离散信号可以通过编写程序直接产生。2、常用时域离散信号的生成1)单位抽样序列单位抽样序列的表示式为或以下三段程序分别用不同的方法来产生单位抽样序列。例1-1用MATLAB的关系运算式来产生单位抽样序列。n1=-5;n2=5;n0=0;n=n1:n2;x=[n==n0];stem(n,x,'filled');axis([n1,n2,0,1.1*max(x)]);xlabel('时间(n)');ylabel('幅度x(n)');title('单位脉冲序列');运行结果如图1-1所示：图1-1单位脉冲序列 2)单位阶跃序列单位阶跃序列表示式为 或以下三段程序分别用不同的方法来产生单位阶跃序列。例1-2用MATLAB的关系运算式来产生单位阶跃序列。n1=-2;n2=8;n0=0;n=n1:n2;x=[n>=n0];stem(n,x,'filled');axis([n1,n2,0,1.1*max(x)]);xlabel('时间(n)');ylabel('幅度x(n)');title('单位阶跃序列');box运行结果如图1-2所示：图1-2单位阶跃序列 3)实指数序列实指数序列的表示式为x(n)=an其中a为实数例1-3编写产生a=1/2和a=2的实指数连续信号和离散序列的程序n1=-10;n2=10;a1=0.5;a2=2;na1=n1:0;x1=a1.^na1;na2=0:n2;x2=a2.^na2;subplot(2,2,1);plot(na1,x1);(图形窗口)title('实指数信号（a<1）');subplot(2,2,3);stem(na1,x1,'filled');title('实指数序列（a<1）');subplot(2,2,2);plot(na2,x2);title('实指数信号（a>1）');subplot(2,2,4);stem(na2,x2,'filled');title('实指数序列（a<1）');box程序运行结果如图1-3所示：图1-3实指数序列4）正（余）弦序列正（余）弦序列的表示式为x(n)=Umsin(ω0n+Θ)例1-4已知一时域周期性正弦信号的频率为1Hz，振幅值为1V。编写程序在图形窗口上显示两个周期的信号波形，并对该信号的一个周期进行32点采样获得离散信号。f=1;Um=1;nt=2;N=32;T=1/f;dt=T/N;n=0:nt*N-1;tn=n*dt;x=Um*sin(2*f*pi*tn);subplot(2,1,1);plot(tn,x);axis([0,nt*T,1.1*min(x),1.1*max(x)]);ylabel('x(t)');subplot(2,1,2);stem(tn,x);axis([0,nt*T,1.1*min(x),1.1*max(x)]);ylabel('x(n)');box程序运行结果如图1-4所示图1-4正弦信号5）矩形波序列MATLAB提供有专门函数square用于产生矩形波。其调用格式如下：x=square(t)：类似于sin(t),产生周期为2π，幅值为±1的方波。x=square(t，duty)：产生指定周期的矩形波，其中duty用于指定占空比。将square的参数t换成n，且n取整数，则可以获得矩形序列。例1-5一个周期性矩形信号频率为5kHz，信号幅度在0～2V之间，占空比为0.25，。编写程序生成该信号，要求在图形窗口上显示2个周期的信号波形；对信号的一个周期进行16点采样获得离散信号。f=5000;nt=2;N=16;T=1/f;dt=T/N;n=0:nt*N-1;tn=n*dt;x=square(2*f*pi*tn,25)+1;subplot(2,1,1);plot(tn,x);axis([0,nt*T,1.1*min(x),1.1*max(x)]);ylabel('x(t)');subplot(2,1,2);stem(tn,x);axis([0,nt*T,1.1*min(x),1.1*max(x)]);ylabel('x(n)');box程序运行结果如图1-5所示;图1-5矩形信号的采样三、操作步骤及注意事项任务一：阅读并上机验证实验原理部分的例题程序，理解每一条语句的含义。改变例题中的有关参数（如信号的频率、周期、幅度、显示时间的取值范围、采样点数等），观察对信号波形的影响。任务二：编写程序，产生以下离散序列：（1）f(n)=δ(n)(-3<n<4)（2）f(n)=u(n)(-5<n<5)（3）f(n)=e(0.1+j1.6∏)n(0<n<16)（4）f(n)=3sin(nП/4)(0<n<20)任务三：一个连续的周期性方波信号频率为200Hz，信号幅度在-1～+1V之间，要求在图形窗口上显示其两个周期的波形。以4kHz的频率对连续信号进行采样，编写程序生成连续信号和其采样获得的离散信号波形。实验报告1、列写调试通过的实验程序，打印实验程序产生的曲线图形。2、思考题：通过例题程序，你发现采样频率Fs、采样点数N、采样时间间隔dt在程序编写中有怎样的联系？使用时需注意什么问题？实验二离散LSI系统的时域分析一、实验目的（1）加深对离散系统的差分方程、单位脉冲响应、单位阶跃响应和卷积分析方法的理解。（2）初步了解用MATLAB语言进行离散时间系统时域分析的基本方法。（3）掌握求解离散时间系统的单位脉冲响应、单位阶跃响应、线性卷积以及差分方程的程序的编写方法，了解常用子函数的调用格式。二、实验原理1、离散LSI系统的响应与激励由离散时间系统的时域分析方法可知，一个离散LSI系统的响应与激励可以用如下框图表示：其输入、输出关系可用以下差分方程描述：2、用函数impz和dstep求解离散系统的单位脉冲响应和单位阶跃响应。例2-1已知描述某因果系统的差分方程为6y(n)+2y(n-2)=x(n)+3x(n-1)+3x(n-2)+x(n-3)满足初始条件y(-1)=0，x(-1)=0，求系统的单位脉冲响应和单位阶跃响应。解：将y(n)项的系数a0进行归一化，得到y(n)+1/3y(n-2)=1/6x(n)+1/2x(n-1)+1/2x(n-2)+1/6x(n-3)分析上式可知，这是一个3阶系统，列出其bk和ak系数：a0=1，a,1=0，a,2=1/3，a,3=0b0=1/6，b,1=1/2，b,2=1/2，b,3=1/6程序清单如下：a=[1,0,1/3,0];b=[1/6,1/2,1/2,1/6];N=32;n=0:N-1;hn=impz(b,a,n);gn=dstep(b,a,n);subplot(1,2,1);stem(n,hn,'k');title('系统的单位序列响应');ylabel('h(n)');xlabel('n');axis([0,N,1.1*min(hn),1.1*max(hn)]);subplot(1,2,2);stem(n,gn,'k');title('系统的单位阶跃响应');ylabel('g(n)');xlabel('n');axis([0,N,1.1*min(gn),1.1*max(gn)]);程序运行结果如图2-1所示：图2-1单位序列响应和单位阶跃响应3、用函数filtic和filter求解离散系统的单位序列响应和单位阶跃响应。例2-2已知描述某因果系统的差分方程为6y(n)-2y(n-4)=x(n)-3x(n-2)+3x(n-4)-x(n-6)，满足初始条件y(-1)=0，x(-1)=0，求系统的单位脉冲响应和单位阶跃响应。时间轴上N取32点。注意：原式非标准形式，必须化为标准形式后再列出系数b,a。程序清单如下：x01=0;y01=0;a=[1,0,0,0,-1/3,0,0];b=[1/6,0,-1/2,0,1/2,0,-1/6];N=32;n=0:N-1;xi=filtic(b,a,0);x1=[n==0];hn=filter(b,a,x1,xi);x2=[n>=0];gn=filter(b,a,x2,xi);subplot(1,2,1);stem(n,hn,'k');title('系统的单位序列响应');ylabel('h(n)');xlabel('n');axis([0,N,1.1*min(hn),1.1*max(hn)]);subplot(1,2,2);stem(n,gn,'k');title('系统的单位阶跃响应');ylabel('g(n)');xlabel('n');axis([0,N,1.1*min(gn),1.1*max(gn)]);程序运行结果如图2-2所示：图2-2单位序列响应和单位阶跃响应4、用MATLAB实现线性卷积1）用函数conv进行卷积运算：求解两个序列的卷积和，关键在于如何确定卷积结果的时宽区间。MATLAB提供的求卷积函数conv默认两个序列的序号均从n=0开始，卷积结果y对应的序列的序号也从n=0开始。例2-3已知两个序列f1=0.8n（0<n<20），f2=u(n)（0<n<10），求两个序列的卷积和。n1=0:20;f1=0.8.^n1;subplot(2,2,1);stem(n1,f1,'filled');title('f1(n)');n2=0:10;N2=length(n2);f2=ones(1,N2);subplot(2,2,2);stem(n2,f2,'filled');title('f2(n)');y=conv(f1,f2);subplot(2,1,2);stem(y,'filled');程序运行结果如图2-3所示：图2-3卷积示意图5、离散LSI系统时域响应的求解：MATLAB提供了多种方法求解离散LSI系统的响应：1）用conv函数进行卷积积分，求任意输入的系统零状态响应；例2-4已知描述某因果系统的差分方程为6y(n)+2y(n-2)=x(n)+3x(n-1)+3x(n-2)+x(n-3)，满足初始条件y(-1)=0，x(-1)=0。在该系统的输入端加一个矩形脉冲序列，其占空比为0.25，一个周期取16个采样点，求该系统的响应。程序清单如下：N=16;n=0:N-1;x=[ones(1,N/4),zeros(1,3*N/4)];subplot(2,2,1);stem(n,x);a=[1,0,1/3,0,];b=[1/6,1/2,1/2,1/6];hn=impz(b,a,n);subplot(2,2,2);stem(n,hn);y=conv(x,hn);subplot(2,1,2);stem(y);程序运行结果如图2-4所示：图2-4零状态响应3）用filtic和filter函数求任意输入的系统完全响应。例2-5已知描述某系统的差分方程为y(n)-1.5y(n-1)+0.5y(n-2)=x(n)n≥0，满足初始条件y(-1)=4，y(-2)=10，求系统输入为x(n)=(0.25)nu(n)时的零输入、零状态及全响应。解为了更深入地理解filtic和filter函数的用法，先用经典法求得系统完全响应表达式：，并编写程序绘出其图形，以便与用MATLAB函数求解的结果进行对比。程序清单如下：a=[1,-1.5,0.5];b=1;N=20;n=0:N-1;x=0.25.^n;x0=zeros(1,N);y01=[4,10];xi=filtic(b,a,y01);y0=filter(b,a,x0,xi);xi0=filtic(b,a,0);y1=filter(b,a,x,xi0);y=filter(b,a,x,xi);y2=((1/3)*(1/4).^n+(1/2).^n+(2/3)).*ones(1,N);subplot(2,3,1);stem(n,x);title('输入信号x(n)');subplot(2,3,2);stem(n,y0);title('系统的零输入响应');subplot(2,3,3);stem(n,y1);title('系统的零状态响应');subplot(2,2,3);stem(n,y);title('用filter求得的完全响应');subplot(2,2,4);stem(n,y2);title('经典法求得的完全响应');程序运行结果如图2-5所示：图2-5系统响应对比图三、操作步骤及注意事项任务一：已知描述某离散LSI系统的差分方程为2y(n)-3y(n-1)+y(n-2)=x(n-1)，分别用impz和dstep函数、filtic和filter函数两种方法求解系统的单位序列响应和单位阶跃响应。任务二：编写程序描绘下列序列的卷积波形：（1）f1(n)=u(n),f2(n)=u(n-2),(0≤n<10)（2）x(n)=sin(n/2),h(n)=(0.5)n(-3≤n≤4П)实验报告1、列写调试通过的实验程序，打印实验程序产生的曲线图形。2、列出本实验提出的有关MATLAB函数在调用时应注意的问题。实验三离散LSI系统的频域分析一、实验目的加深对离散系统变换域分析——z变换的理解，掌握使用MATLAB进行z变换和逆z变换的常用函数的用法。（2）了解离散系统的零极点与系统因果性和稳定性的关系，熟悉使用MATLAB进行离散系统的零极点分析的常用函数的用法。二、实验原理1、z变换和逆z变换（1）用ztrans函数求无限长序列的z变换。该函数只给出z变换的表达式，而没有给出收敛域。另外，由于这一函数还不尽完善，有的序列的z变换还不能求出，逆z变换也存在同样的问题。例3-1求以下各序列的z变换x1(n)=anx2(n)=nx3(n)=n(n-1)/2x4(n)=ejωonx5(n)=1/[n(n-1)]程序清单如下：symsw0nza;x1=n*a^n;X1=ztrans(x1)x2=sin(w0*n);X2=ztrans(x2)x3=exp(-a*n)*sin(w0*n);X3=ztrans(x3)程序运行结果如下：X1=z/a/(z/a-1)X2=z/(z-1)^2X3=1/2*z*(z+1)/(z-1)^3-1/2*z/(z-1)^2X4=z/exp(i*w0)/(z/exp(i*w0)-1)X5=z/(z-1)-ztrans(1/n,n,z)（2）用iztrans函数求无限长序列的逆z变换。2、离散系统的零极点分析（系统极点位置对系统响应的影响）例3-2研究z右半平面的实数极点对系统的影响。已知系统的零极点增益模型分别为：求这些系统的零极点分布图以及系统的单位序列响应，判断系统的稳定性。程序清单如下：z1=[0]';p1=[0.85]';k=1;[b1,a1]=zp2tf(z1,p1,k);subplot(3,2,1);zplane(z1,p1);title('极点在单位圆内');subplot(3,2,2);impz(b1,a1,20);z2=[0]';p2=[1]';[b2,a2]=zp2tf(z2,p2,k);subplot(3,2,3);zplane(z2,p2);title('极点在单位圆上');subplot(3,2,4);impz(b2,a2,20);z3=[0]';p3=[1.5]';[b3,a3]=zp2tf(z3,p3,k);subplot(3,2,5);zplane(z3,p3);title('极点在单位圆外');subplot(3,2,6);impz(b3,a3,20);程序运行结果如图3-1所示。由图可见，这三个系统的极点均为实数且处于z平面的右半平面。由图可知，当极点位于单位圆内，系统的单位序列响应随着频率的增大而收敛；当极点位于单位圆上，系统的单位序列响应为等幅振荡；当极点位于单位圆外，系统的单位序列响应随着频率的增大而发散。由此可知系统1、2为稳定系统。图3-1零极点分布例3-3已知某离散时间系统的系统函数为求该系统的零极点及零极点分布图，并判断系统的因果稳定性。程序清单如下：b=[0.2,0.1,0.3,0.1,0.2];a=[1,-1.1,1.5,-0.7,0.3];rz=roots(b)rp=roots(a)subplot(2,1,1);zplane(b,a);title('系统的零极点分布图');subplot(2,1,2);impz(b,a,20);title('系统的单位序列响应');xlabel('n');ylabel('h(n)');程序运行结果如下：rz=-0.5000+0.8660i-0.5000-0.8660i0.2500+0.9682i0.2500-0.9682irp=0.2367+0.8915i0.2367-0.8915i0.3133+0.5045i0.3133-0.5045i图3-2零极点分布由零极点分布图可见，该系统的所有极点均在单位圆内，因此该系统是一个因果稳定系统。3、离散系统的频率响应（1）离散系统的频率响应的基本概念已知稳定系统传递函数的零极点增益模型为则系统的频响函数为其中，系统的幅频特性为系统的相频特性为由以上各式可见，系统函数与频率响应有着密切的联系。适当地控制系统函数的零极点分布，可以改变离散系统的频响特性：①在原点（z=0）处的零点或极点至单位圆的距离始终保持不变，其值|ejω|=1，所以，对幅度响应不起作用；②单位圆附近的零点对系统幅度响应的谷值位置及深度有明显影响；③单位圆内且靠近单位圆附近的极点对系统幅度的峰值位置及大小有明显的影响。（2）系统的频响特性分析例3-4已知某离散时间系统的系统函数为求该系统在0～П频率范围内的相对幅频响应与相频响应。程序清单如下：b=[0.1321,0,-0.3963,0,0.3963,0,-0.1321];a=[1,0,0.34319,0,0.60439,0,0.20407];freqz(b,a);程序运行结果如图3-5所示。该系统是一个IIR数字带通滤波器。其中幅频特性采用归一化的相对幅度值，以分贝（dB）为单位。图3-3幅频响应与相频响应（3）一个求解频率响应的实用函数。在实际使用freqz进行离散系统频响特性分析时。通常需要求解幅频响应、相频响应、群时延，幅频响应又分为绝对幅频和相对幅频两种表示方法。下面定义函数freqz_m，利用该函数，可方便求出上述各项。freqz_m函数定义如下：function[db,mag,pha,grd，w]=freqz_m(b,a);[H,w]=freqz(b,a,1000,'whole');H=(H(1:501))';w=(w(1:501))';mag=abs(H);db=20*log10((mag+eps)/max(mag));pha=angle(H);grd=grpdelay(b,a,w);例3-5已知某离散时间系统的系统函数为求该系统在0～П频率范围内的绝对幅频响应与相频响应、相对幅频响应与相频响应及群时延。程序清单如下：b=[0.1321,0,0.3963,0,0.3963,0,0.1321];a=[1,0,-0.34319,0,0.60439,0,-0.20407];[db,mag,pha,grd,w]=freqz_m(b,a);subplot(2,2,1);plot(w/pi,mag);gridaxis([0,1,1.1*min(mag),1.1*max(mag)]);title('幅频特性（V）');xlabel('\omega/\pi');ylabel('幅度(V)');subplot(2,2,2);plot(w/pi,pha);grid;axis([0,1,1.1*min(pha),1.1*max(pha)]);xlabel('\omega/\pi');ylabel('相位');title('相频特性');subplot(2,2,3);plot(w/pi,db);gridaxis([0,1,-100,5]);title('幅频特性（dB）');subplot(2,2,4);plot(w/pi,grd);gridaxis([0,1,0,10])title('群时延');程序运行结果如图3-8所示：图3-4频率响应三、操作步骤及注意事项任务一：求一下系统函数所描述的离散系统的零极点分布图，并判断系统的稳定性（1）（2）任务二：已知某离散时间系统的系统函数为求该系统在0～П频率范围内的绝对幅频响应与相频响应、相对幅频响应与相频响应及群时延。实验报告1、列写调试通过的实验程序，打印实验程序产生的曲线图形。2、思考题：系统函数零极点的位置与系统单位序列响应有何关系？实验四DFS、DFT与FFT一、实验目的加深对周期序列DFS、有限长序列DFT和FFT的基本概念及理论的理解。（2）MATLAB语言求解DFS、DFT、FFT的以及相应反变换的方法。二、实验原理1、周期序列的离散傅里叶级数（DFS）（1）DFS的基本概念：离散时间序列x(n)满足x(n)=x(n+rN)，称为离散周期序列，用表示。其中，N为信号的周期，x(n)称为离散周期序列的主值。周期序列可以用离散傅里叶级数（DFS）表示：其中，是周期序列DFS第k次谐波分量的系数，也称为周期序列的频谱，可表示为以上两式也是周期序列的一对傅里叶级数变换对。令，以上DFS变换对又可以写成：与连续周期信号的傅里叶级数相比，周期序列的离散傅里叶级数有以下特点：①连续周期信号的傅里叶级数由无穷多个与基波频率成整数倍的谐波分量叠加而成，而周期为N的周期序列的傅里叶级数仅有N个独立的谐波分量。②周期序列的频谱也是一个以N为周期的周期序列。（2）周期序列的DFS和IDFS例4-1已知一个周期性矩形序列的脉冲宽度占整个周期的1/4，一个周期的采样点数为16点，编程显示3个周期的序列波形，并：①用傅里叶级数求信号的幅度谱和相位谱。②求傅里叶级数逆变换的图形，并与原序列进行比较。程序清单如下：N=16;xn=[ones(1,N/4),zeros(1,3*N/4)];xn=[xn,xn,xn];n=0:3*N-1;k=0:3*N-1;Xk=xn*exp(-j*2*pi/N).^(n'*k);x=(Xk*exp(j*2*pi/N).^(n'*k))/N;subplot(2,2,1);stem(n,xn);title('x(n)');axis([-1,3*N,1.1*min(xn),1.1*max(xn)]);subplot(2,2,2);stem(n,abs(x));title('IDFS|X(k)|');axis([-1,3*N,1.1*min(x),1.1*max(x)]);subplot(2,2,3),stem(k,abs(Xk));title('|X(k)|');axis([-1,3*N,1.1*min(abs(Xk)),1.1*max(abs(Xk))]);subplot(2,2,4),stem(k,angle(Xk));title('arg|X(k)|');axis([-1,3*N,1.1*min(angle(Xk)),1.1*max(angle(Xk))]);程序运行结果如图4-1所示：图4-1幅度谱与相位谱由离散傅里叶级数逆变换图形可见，与原序列相比，幅度扩大了32倍。这是因为周期序列为原主值序列周期的3倍，做逆变换时未做处理。可将逆变换程序改为：x=(Xk*exp(j*2*pi/N).^(n'*k))/3*3*N;由上例可见，周期序列的DFS和IDFS是依据变换公式编程的，无论信号序列如何变化，求解的公式总是一样的。因此，可将其编写成通用子程序：①离散傅里叶级数变换通用子程序dfs.mfunction[Xk]=dfs(xn,N)n=0:N-1;k=0:N-1;WN=exp(-j*2*pi/N);nk=n'*k;Xk=xn*WN.^nk;②离散傅里叶级数逆变换通用子程序idfs.mfunction[xn]=idfs(Xk,N)n=0:N-1;k=0:N-1;WN=exp(j*2*pi/N);nk=n'*k;xn=(Xk*WN.^nk)/N;例4-2利用上述两个子程序，重做例4-1程序清单如下：N=16;xn=[ones(1,N/4),zeros(1,3*N/4)];Xk=dfs(xn,N);x=idfs(Xk,N);subplot(2,2,1);stem(n,xn);title('x(n)');axis([-1,3*N,1.1*min(xn),1.1*max(xn)]);subplot(2,2,2);stem(n,abs(x));title('IDFS|X(k)|');axis([-1,3*N,1.1*min(x),1.1*max(x)]);subplot(2,2,3),stem(k,abs(Xk));title('|X(k)|');axis([-1,3*N,1.1*min(abs(Xk)),1.1*max(abs(Xk))]);subplot(2,2,4),stem(k,angle(Xk));title('arg|X(k)|');axis([-1,3*N,1.1*min(angle(Xk)),1.1*max(angle(Xk))]);程序运行结果如图4-2所示。由于子程序仅适用于对主值区间进行变换，周期次数无法传递给子程序，因此程序执行结果仅显示一个周期的变换情况。图4-2幅度谱与相位谱（2）周期重复次数对序列频谱的影响理论上讲，周期序列不满足绝对可积条件，因此不能用傅里叶级数来表示。实际处理时可先取K个周期进行处理，然后令K趋于无穷大，分析其极限情况。根据这一分析思路，可以观察序列由非周期到周期变化时，频谱由连续谱逐渐向离散谱过渡的过程。2、离散傅里叶变换（DFT）（1）DFT与IDFT在实际中常常使用有限长序列。如果有限长序列为x(n)，则该序列的离散傅里叶变换对可表示为从离散傅里叶变换定义式可以看出，有限长序列在时域上是离散的，在频域上也是离散的，式中，即仅在单位圆上N个等间距的点上取值，这为使用计算机进行处理带来了方便。由有限长序列的傅里叶变换和逆变换定义可知，DFT和DFS的变换公式非常相似，因此，在程序编写上也基本一致。例4-3已知x(n)=[0,1,2,3,4,5,6,7]，求其DFT和IDFT。要求：①画出序列傅里叶变换对应的|X(k)|和arg[X(k)]图形。②画出x(n)图形，并与IDFT[X(k)]图形进行比较。程序清单如下：xn=[0,1,2,3,4,5,6,7];N=length(xn);n=0:N-1;k=0:N-1;Xk=xn*exp(-j*2*pi/N).^(n'*k);x=(Xk*exp(j*2*pi/N).^(n'*k))/N;subplot(2,2,1);stem(n,xn);title('x(n)');axis([-1,N,1.1*min(xn),1.1*max(xn)]);subplot(2,2,2);stem(n,abs(x));title('IDFT|X(k)|');axis([-1,N,1.1*min(x),1.1*max(x)]);subplot(2,2,3),stem(k,abs(Xk));title('|X(k)|');axis([-1,N,1.1*min(abs(Xk)),1.1*max(abs(Xk))]);subplot(2,2,4),stem(k,angle(Xk));title('arg|X(k)|');axis([-1,N,1.1*min(angle(Xk)),1.1*max(angle(Xk))]);程序运行结果如图4-3所示。由图可见，与周期序列不同，有限长序列本身是仅有N点的离散序列，相当于周期序列的主值部分。因此，其频谱也对应序列的主值部分，是长度为N的离散序列。图4-3DFT和IDFT（2）DFT与DFS的联系将周期序列的傅里叶级数变换对和有限长序列的离散傅里叶变换对进行比较可见，两者的区别仅仅是将周期序列换成了有限长序列x(n)，同时，由于式中的周期性，因而有限长序列的离散傅里叶变换实际上隐含着周期性。（3）DFT与DTFT的联系若离散时间非周期序列为x(n)，则它的离散傅里叶变换（DTFT）对定义为其中称为序列的频谱。可表示为，称为序列的幅度谱，称为序列的相位谱。由DTFT的定义可见，序列在时域是离散的、非周期的，在频域是连续的、周期的。与有限长序列相比，仅在单位圆上取值，X(k)是在单位圆上N个等间距的点上取值。因此，连续谱可由离散谱X(k)经插值后得到。3、FFT（1）MATLAB提供的FFT函数由理论学习可知，DFT是唯一在时域和频域均离散的变换方法，它适用于有限长序列。尽管这种变换方法是可以用于数值计算的，但如果只是简单地按照定义进行数据处理，当序列长度很大时，将占用很大的内存空间，且运算时间很长。快速傅里叶变换是用于DFT运算的高效快速算法的统称，FFT只是其中的一种。FFT主要有时域抽取算法和频域抽取算法，基本思想是将一个长度为N的序列分解成多个段序列，如基2算法、基4算法等，大大地缩短了DFT的时间。有关详细理论可参考教材。MATLAB提供了进行FFT的函数fft和ifft分别用于计算DFT和IDFT。（2）用FFT进行频谱分析①对有限长序列进行谱分析一个序号从n1到n2的时域有限长序列x(n)，它的频谱定义为它的离散傅里叶变换，且在Nyquist频率范围内有界并连续。序列的长度为N，则N=n2-n1+1。计算x(n)的DFT得到的是的N个样本点。其中数字频率为式中：为数字频率的分辨率；k取对应-（N-1）/2到（N-1）/2区间的整数。在实际使用中，往往要求计算出信号以模拟频率为横坐标的频谱，此时对应的模拟频率为式中：D为模拟频率的分辨率或频率间隔；Ts为采样信号的周期，Ts=1/Fs；定义信号的长度L=NTs。在使用FFT进行DFT的高校运算时，一般不直接用n从n1到n2的x(n)，而是取的主值区间（n=0,1,…,N-1）的数据，经FFT将产生N个数据，定位在k=0,1,…,N-1的数字频率点上，即对应[0,2]。如果要显示[-，]范围的频谱，则可以使用fftshift(X)进行位移。例4-4已知有限长序列x(n)=[1，2，3，2，1]，其采样频率Fs=10Hz。请使用FFT计算其频谱。程序清单如下Fs=10;xn=[1,2,3,2,1];N=length(xn);D=2*pi*Fs/N;k=floor(-(N-1)/2:(N-1)/2);X=fftshift(fft(xn,N));subplot(1,2,1);plot(k*D,abs(X),'o:');title('幅度频谱');xlabel('rad/s');subplot(1,2,2);plot(k*D,angle(X),'o:');title('相位频谱');xlabel('rad/s');程序运行结果如图4-4所示：图4-4FFT由图4-4可知，当有限长序列的长度N=5时，频谱的样本点数也为5，频率点之间的间距非常大，即分辨率很低。及时使用了plot命令的插值功能，显示出的曲线仍是断续的，与真实曲线有较大误差。改变分辨率的基本方法是给输入序列补零，即增加频谱的密度。这种方法只是改善了图形的是在分辨率，并不增加频谱的细节信息。将上述有限长序列x(n)[1,2,3,2,1]末尾补零到N=1000点，将程序改为：Fs=10;N=1000;xn=[1,2,3,2,1];Nx=length(xn);xn=[1,2,3,2,1,zeros(1,N-Nx-1)];D=2*pi*Fs/N;k=floor(-(N-1)/2:(N-1)/2);X=fftshift(fft(xn,N));subplot(1,2,1);plot(k*D,abs(X));title('幅度频谱');xlabel('rad/s');subplot(1,2,2);plot(k*D,angle(X));title('相位频谱');xlabel('rad/s');程序的运行结果如图4-5所示，由图可见，图形的分辨率提高，曲线几乎是连续的频谱了。图4-5补零FFT②对无限长序列进行谱分析用FFT进行无限长序列的频谱分析，首先要将无限长序列截断成一个有限长序列。序列长度的取值对频谱有较大的影响，带来的问题是引起频谱的泄漏和波动。三、操作步骤及注意事项任务一：已知有限长序列x(n)=[1,0.5,0,0.5,1,1,0.5,0]，要求：①求该序列的DFT、IDFT的图形；②用FFT算法求该序列的DFT、IDFT的图形；③假定采用频率Fs=20Hz，序列长度N分别取8、32和64，用FFT计算其幅度谱和相位谱。实验报告1、列写调试通过的实验程序，打印实验程序产生的曲线图形。2、思考题：DFS、DFT、FFT有何联系？实验五IIR数字滤波器的设计一、实验目的掌握双线性变换法及脉冲相应不变法设计IIR数字滤波器的具体设计方法及其原理，熟悉用双线性变换法及脉冲响应不变法设计低通、高通和带通IIR数字滤波器的计算机编程。掌握用MATLAB产生时域离散信号的方法（2）观察双线性变换及脉冲响应不变法设计的滤波器的频域特性，了解双线性变换法及脉冲响应不变法的特点。（3）熟悉巴特沃思滤波器、切比雪夫滤波器和椭圆滤波器的频率特性。二、实验原理1．脉冲响应不变法所谓脉冲响应不变法就是使数字滤波器的单位脉冲响应序列h(n)等于模拟滤波器的单位冲激响应和ha(t)的采样值，即：h(n)=ha(nT)，其中，T为采样周期。在脉冲响应不变法中，模拟角频率和数字角频率的变换关系为：，可见，Ω和ω之间的变换关系为线性的。在MATLAB中，可用函数impinvar实现从模拟滤波器到数字滤波器的脉冲响应不变映射，调用格式为：[B,A]=impinvar(b,a,fs1)[B,A]=impinvar(b,a)其中，b、a分别为模拟滤波器的分子和分母多项式系数向量；fs1为采样频率（Hz），缺省值fs=1Hz；B、A分别为数字滤波器分子和分母多项式系数向量。2．双线性变换法：由于s平面和z平面的单值双线性映射关系为，其中T为采样周期。因此，若已知模拟滤波器的传递函数，将上式代入即可得到数字滤波器的系统函数H(z)。在双线性变换中，模拟角频率和数字角频率的变换关系为：，可见，Ω和ω之间的变换关系为非线性的。在MATLAB中，可用函数bilinear实现从模拟滤波器到数字滤波器的双线性变换映射，调用格式为：[B,A]=bilinear(b,a,fs1)3．滤波器设计（1）定技术指标转换为模拟滤波器设计性能指标。（2）估计满足性能指标的模拟相应滤波器性能阶数和截止频率。利用MATLAB中buttord、cheb1ord、cheb2ord、ellipord等函数，对于模拟滤波器：[n,Wn]=buttord(Ωp,Ωs,αp,αs,’s’)其中，Ωp为通带边界频率，rad/s；Ωs为阻带边界频率，rad/s;αp为带通波动，dB;αs为阻带衰减，dB；‘s’表示为模拟滤波器；函数返回值N为模拟滤波器的最小阶数；Ωc为模拟滤波器的截止频率(-3dB频率)，rad/s。函数适用低通、高通、带通、带阻滤波器。对于数字滤波器，[n,Wn]=buttord(Wp,Ws,Rp,Rs),其中Wp，Ws是归一化数字频率，.例如，设计一个低通数字滤波器，采样频率为1000Hz，通带截止频率为40Hz，衰减为3dB，阻带截止频率为150Hz，衰减为60dB,则命令为：Wp=40/500,Ws=150/500;[n,Wn]=buttord(Wp,Ws,3,60);利用buttord计算阶数N和通带截止频率Wn。[B,A]=BUTTER(N,Wn,'high')用来设计高通滤波器

[B,A]=BUTTER(N,Wn,'low')designsalowpassfilter.--低通滤波器

[B,A]=BUTTER(N,Wn)--带通滤波器y=filter(B,A,x)得到滤波器系数之后可以直接用。例5-1，，，,，设计一切比雪夫高通滤波器，观察其通带损耗和阻带衰减是否满足要求。实验程序：clc;fc=300;Ap=0.8;fr=200;At=20;T=1000;%这里的T表示的是采样频率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;plot(f,20*log10(abs(h)));axis([0,500,-80,10]);grid;xlabel('频率/Hz');ylabel('幅度/dB');title('切比雪夫高通滤波器');运行结果：图5-1切比雪夫滤波器分析：f=200Hz时阻带衰减大于30dB，通过修改axis([0,fs/2,-80,10])为axis([200,fs/2,-1,1])发现通带波动rs满足<0.8。例5-2，，，,，分别用脉冲响应不变法及双线性变换法设计一巴特沃思数字低通滤波器，观察所设计数字滤波器的幅频特性曲线，记录带宽和衰减量，检查是否满足要求。比较这两种方法的优缺点。实验程序：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]=butter(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('脉冲响应不变法','双线性变换法');运行结果：图5-2巴特沃斯滤波器例5-3利用双线性变换法分别设计满足下列指标的巴特沃思型、切比雪夫型和椭圆型数字低通滤波器，并作图验证设计结果:f=1.2kHz，δ≤0.5dB，f=2kHz,At≥40dB，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');[num1,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);[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);运行结果：图5-3三种滤波器对比也可以直接用butter数字滤波器来做：w1=1200/4000;w2=2000/4000;[Nb,wn]=buttord(w1,w2,rp,rs);[B,A]=butter(Nb,wn,'low');[h1,w]=freqz(B,A);例5-4分别用脉冲响应不变法和双线性变换法设计一巴特沃思型数字带通滤波器，已知f=30kHz，其等效的模拟滤波器指标为δ＜3dB，2kHz＜f≤3kHz；At≥5dB，f≥6kHz；At≥20dB，f≤1.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,'s');[num1,den1]=impinvar(B,A,fs);%脉冲响应不变法[h1,w]=freqz(num1,den1);w1=2*fs*tan(2*pi*2000/(2*fs));w2=2*fs*tan(2*pi*3000/(2*fs));wr1=2*fs*tan(2*pi*1500/(2*fs));wr2=2*fs*tan(2*pi*6000/(2*fs));[N,wn]=buttord([w1,w2],[wr1,wr2],rp,rs,'s');[B,A]=butter(N,wn,'s');[num2,den2]=bilinear(B,A,fs);%双线性变换法[h2,w]=freqz(num2,den2);f=w/pi*15000;plot(f,20*log10(abs(h1)),'-.',f,20*log10(abs(h2)),'-');axis([500,7000,-30,10]);grid;xlabel('频率/Hz');ylabel('幅度、dB');title('巴特沃思数字低通滤波器');legend('脉冲响应不变法','双线性变换法',1);运行结果：图5-4两种方法对比图三、操作步骤及注意事项任务一：设计模拟巴特沃斯低通滤波器，fp=300Hz，αp=1dB，fs=800Hz，αs=20dB。任务二：fp=0.１kHZ,αp=1dB，fs=0.3kHZ，αs=25dB，T=1ms；分别用脉冲响应不变法和双线性变换法设计一个Butterworth数字低通滤波器。实验报告列写调试通过的实验程序，打印实验程序产生的曲线图形。实验六FIR数字滤波器的设计一、实验目的掌握用窗函数法，频率采样法及优化设计法设计FIR滤波器的原理及方法，熟悉响应的计算机编程；（2）熟悉线性相位FIR滤波器的幅频特性和相频特性；二、实验原理窗函数设计法是一种把一个长序列变成有限长的短序列的设计方法，是在时域进行的。用窗函数法设计FIR数字滤波器时，先根据Wc和N求出相应的的理想滤波器的单位脉冲响应hd(n)。因为hd(n)一般是非因果的，且无限长，物理上是不可实现的。为此可选择适当的窗函数w(n)截取有限长的hd(n)，即h(n)=hd(n)w(n)，只要阶数足够长，截取的方法合理，总能够满足频域的要求。实际中常用的窗函数有矩形（Boxcar）窗、三角（Bartlett）窗、汉宁(Hanning)窗、汉明(Hamming)窗和布莱克曼（Blackman）窗。这些窗函数各有优缺点,所以要根据实际情况合理选择窗函数类型。1．窗函数法设计线性相位FIR滤波器的一般步骤为：（1）确定理想滤波器的特性；（2）由Hd（ejw）求出；（3）根据过渡带宽度和阻带最小衰减，借助窗函数确定窗的形式及N的大小，即选择适当的窗函数，并根据线性相位条件确定窗函数的长度N；在MATLAB中，可由w=boxcar(N)（矩形窗）、w=hanning(N)（汉宁窗）、w=hamming(N)（汉明窗）、w=Blackman(N)（布莱克曼窗）、w=Kaiser(N,beta)（凯塞窗）等函数来实现窗函数设计法中所需的窗函数。（4）由h(n)=.w(n),0≤n≤N-1，得出单位脉冲响应h(n);（5）对h(n)作离散时间傅立叶变换，得到H()。2．在MATLAB中，可以用b=fir1(N,Wn,’ftype’,taper)等函数辅助设计FIR数字滤波器。N代表滤波器阶数；Wn代表滤波器的截止频率(归一化频率)，当设计带通和带阻滤波器时，Wn为双元素相量；ftype代表滤波器类型，如’high’高通，’stop’带阻等；taper为窗函数，默认为海明窗，窗函数实现需要用窗函数blackman,hamming,hanningchebwin,kaiser产生。例6-1.N=45，计算并画出矩形窗、汉明窗、布莱克曼窗的归化幅度谱，并比较各自的特点。实验程序：N=45;%矩形窗window1=boxcar(N);wvtool(window1);%汉明窗window2=hamming(N);wvtool(window2);%布莱克曼窗window3=black