数字信号处理实验指导书(M)

《数字信号处理实验指导书(M)》由会员分享，可在线阅读，更多相关《数字信号处理实验指导书(M)（28页珍藏版）》请在装配图网上搜索。

1、数字信号处理实验电子信息科学与技术实验室2007年7月目录实验一 离散时间信号的时域表示2实验二 离散信号的卷积和5实验三 离散傅立叶变换及其特性验证7实验四 信号处理中FFT的应用10实验五 离散系统的Z域分析15实验六 无限冲激响应(IIR)数字滤波器的三种结构16实验七 冲激响应不变法IIR数字滤波器设计20实验八 双线性变换法IIR数字滤波器设计23 实验一 离散时间信号的时域表示一、实验目的1、熟悉Matlab命令，掌握离散时间信号序列的时域表示方法。2、掌握用Matlab描绘二维图像的方法。3、掌握用Matlab对序列进行基本的运算和时域变换的方法。二、实验原理与计算方法（一）序列

2、的表示方法序列的表示方法有列举法、解析法和图形法，相应的用Matlab也可以有这样几种表示方法，分别介绍如下：1、列举法在Matlab中，用一个列向量来表示一个有限长序列，由于一个列向量并不包含位置信息，因此需要用表示位置的n和表示量值的x两个向量来表示任意一个序列，如：例1.1：n=-3,-2,-1,0,1,2,3,4;x=2,1,-1,0,1,4,3,7;如果不对向量的位置进行定义，则Matlab默认该序列的起始位置为n=1。由于内存有限，Matlab不能表示一个无限序列。2、解析法对于有解析表达式的确定信号，首先定义序列的范围即n的值，然后直接写出该序列的表达式，如：例1.2：实现实指数

3、序列，的Matlab程序为：n=0:10;x=(0.9).n;例1.3：实现正余弦序列，的Matlab程序为：n=5:15;x=3*cos(0.1*pi*n+pi/3)+2*sin(0.5*pi*n);3、图形法在Matlab中用图形法表示一个序列，是在前两种表示方法的基础上将序列的各个量值描绘出来，即首先对序列进行定义，然后用相应的画图语句画图，如：例1.4：绘制在例1.1中用列举法表示的序列的图形，则在向量定义之后加如下相应的绘图语句：stem(n,x);此时得到的图形的横坐标范围由向量n的值决定，为-3到4，纵坐标的范围由向量x的值决定，为-1到7。应用stem函数时应确保自变量n和函数

4、值x的个数相等。此外可用函数axis(x1,x2,y1,y2)对横纵坐标进行限定，以完善图形，其中x1和x2分别为横坐标的起始和截止位置，y1和y2分别为纵坐标的起始和截止位置。也可用xlabel()、ylabel()和title()为该图添加横、纵坐标说明和标题。subplot(m,n,k)函数可以将当前窗口分成m行n列个子窗口，并在第k的子窗口绘图。窗口的排列顺序为从左至右，从上至下分别为1,2,m*n。以上为各个绘图函数的基本用法，有关各函数的其他参数可参考Matlab的帮助文件。下面给出产生单位抽样序列和单位阶跃序列的两个函数，供参考。例1.5：产生单位抽样序列的函数impseq(n0

5、,n1,n2)。function x,n = impseq(n0,n1,n2)% Generates x(n) = delta(n-n0); n1 = n,n0 = n2% -% x,n = impseq(n0,n1,n2)%if (n0 n2) | (n1 n2)error(arguments must satisfy n1 = n0 = n2)endn = n1:n2;x = (n-n0) = = 0;该函数产生一个抽样位置n0位于n1和n2之间的单位抽样序列。例1.6：产生单位阶跃序列的函数stepseq(n0,n1,n2)。function x,n = stepseq(n0,n1,n2

6、)% Generates x(n) = u(n-n0); n1 = n,n0 = n2% -% x,n = stepseq(n0,n1,n2)%if (n0 n2) | (n1 n2)error(arguments must satisfy n1 = n0 = 0;该函数产生一个起始位置n0位于n1和n2之间的单位阶跃序列。注意：由function产生的函数文件，不能直接运行，并且要放在当前路径下的文件夹里，供其他M文件调用。（二）序列的基本运算和时域变换1、加法：x1(n)+x2(n)序列的加法运算为对应位置处量值的相加，在Matlab中可用运算符“+”实现，但要求参与运算的序列的长度必须相

7、等。如果长度不等或者长度相等但采样位置不同，则不能直接应用该运算符，此时需要先给定参数使序列具有相同的位置向量和长度。下面给出sigadd函数实现任意两序列的加法运算。例1.7：function y,n = sigadd(x1,n1,x2,n2)% implements y(n) = x1(n)+x2(n)% % y,n = sigadd(x1,n1,x2,n2)% y = sum sequence over n, which includes n1 and n2% x1 = first sequence over n1% x2 = second sequence over n2 (n2 ca

8、n be different from n1)%n = min(min(n1),min(n2):max(max(n1),max(n2); % duration of y(n)y1 = zeros(1,length(n); y2 = y1; % initializationy1(find(n=min(n1)&(n=min(n2)&(n0时，表示序列向右平移，此时新序列的位置向量为n+m；当m0时，表示序列向左平移，此时新序列的位置向量为n-m。三、实验内容（1）参考示例程序，产生一个有延迟的单位抽样序列：d (n-11)，绘出序列的图形。（2）参考示例程序，产生一个向前时移7个时刻的单位阶跃序列

9、：u(n+7)，绘出序列的图形。（3）产生一个指数为-0.1+(pi/6)*in的复指数序列，并绘出序列的实部、虚部、幅度和相位的波形。（4）已知x(n)=1,2,3,4,5,6,7, 6,5,4,3,2,1，参考示例程序，绘出下列序列的波形。a.x1(n)=2x(n-5)-3x(n+4)b.x2(n)=x(3-n)+x(n)x(n-2)四、思考（1）代数运算符号和.的区别是？实验二 离散信号的卷积和一、实验目的1、掌握两个离散信号卷积和的计算方法和编程技术。2、进一步熟悉用Matlab描绘二维图像的方法。二、实验原理与计算方法两个离散序列x(n)与y(n)的卷积和f(n)定义为由于通常信号处

10、理中所碰到的都是有始信号或有限时间信号，因此在实际计算卷积和时，求和是在有限范围内进行的。计算过程中上下限的选取和所得结果的分布区间取决于参与卷积的两个序列，下面将分别进行讨论：1、两个从n = 0开始的序列和的卷积和 (1)上式右边因子u(n)表示卷积和的结果也是一个从n = 0开始的序列。2、从n = n1开始的序列和从n = n2开始的序列的卷积和，其中n1和n2为任意整数。 (2)上式右边因子u(n-n1-n2)表示卷积和是一个从n = n1+n2开始的序列。3、从n = n1开始的长度为N1的加窗序列和从n = n2开始的长度为N2的加窗序列的卷积和，其中则 (3)所得卷积和也是一个

11、加窗序列，从n = n1+ n2开始，长度为N1+ N2-1。Matlab提供了一个内部函数conv(x,h)用来计算两个有限长度序列的卷积，该函数得到的卷积结果默认从n=0开始，因此当参与卷积的两个序列的起始位置不是n=0时，则由该函数得到的计算结果将出现错误，此时需要重新定义结果的位置向量。由以上卷积运算的原理可知，两有限长序列卷积后仍为有限长序列，长度为两序列长度之和减1，结果的起始位置为两序列起始位置之和，截止位置为两序列截至位置之和。据此，可以得到卷积结果的位置向量。三、实验内容(1) 根据（1）式计算两个从n = 0开始的序列和的卷积和，其中A=40，a = 0.5。取50个样值点

12、，作出序列、及卷积和f(n)的图像。(2) 根据（3）式计算两个有限长序列和的卷积和，其中作出序列、及卷积和f(n)的图像。实验三 离散傅立叶变换及其特性验证一、实验目的1、掌握离散时间傅立叶变换（DTFT）的计算方法和编程技术。2、掌握离散傅立叶变换（DFT）的计算方法和编程技术。3、理解离散傅立叶变换（DFT）的性质并用Matlab进行验证。二、实验原理与计算方法1、离散时间傅立叶变换如果序列x(n)满足绝对可和的条件，即，则其离散时间傅立叶变换定义为： （1）如果x(n)是无限长的，则不能直接用Matlab由x(n)计算X(ejw)，但可以用它来估计X(ejw)表达式在0,频率区间的值并

13、绘制它的幅频和相频（或实部和虚部）曲线。如果x(n)是有限长的，则可以用Matlab对任意频率w处的X(ejw)进行数值计算。如果要在0,间按等间隔频点估计X(ejw)，则（1）式可以用矩阵向量相乘的运算来实现。假设序列x(n)在(即不一定在0, N-1)有N个样本，要估计下列各点上的X(ejw)：它们是0,之间的(M+1)个等间隔频点，则（1）式可写成： （2）将x(nl)和X(ejwk)分别排列成向量x和X，则有： X=Wx （3）其中W是一个(M+1)N维矩阵：将k和n排成列向量，则在Matlab中，把序列和下标排成行向量，对（3）式取转置得：其中nTk是一个N(M+1)维矩阵。用Mat

14、lab实现如下：k=0:M; n=n1:n2;X=x*(exp(-j*pi/M).(n*k);2、离散傅立叶变换一个有限长序列的离散傅立叶变换对定义为： （4） （5）以列向量x和X形式排列x(n)和X(k)，则式（4）、（5）可写成：X=WNx其中矩阵WN由下式给出：可由下面的Matlab函数dft和idft实现离散傅立叶变换运算。function Xk = dft(xn,N)% Computes Discrete Fourier Transform% -% Xk = dft(xn,N)% Xk = DFT coeff. array over 0 = k = N-1% xn = N-poin

15、t finite-duration sequence% N = Length of DFT%n = 0:1:N-1; % row vector for nk = 0:1:N-1; % row vecor for kWN = exp(-j*2*pi/N); % Wn factornk = n*k; % creates a N by N matrix of nk valuesWNnk = WN . nk; % DFT matrixXk = xn * WNnk; % row vector for DFT coefficientsfunction xn = idft(Xk,N)% Computes I

16、nverse Discrete Transform% -% xn = idft(Xk,N)% xn = N-point sequence over 0 = n = N-1% Xk = DFT coeff. array over 0 = k = N-1% N = length of DFT%n = 0:1:N-1; % row vector for nk = 0:1:N-1; % row vecor for kWN = exp(-j*2*pi/N); % Wn factornk = n*k; % creates a N by N matrix of nk valuesWNnk = WN . (-

17、nk); % IDFT matrixxn = (Xk * WNnk)/N; % row vector for IDFT values3、离散傅立叶变换的性质（1）线性性质：注意：若x1(n)和x2(n)分别是N1点和N2点的序列，则选择N3= max (N1, N2)，将它们作N3点DFT处理。（2） 周期性：离散傅立叶变换(DFT)是周期序列DFS取主值区间形成的，因此序列及其DFT 具有特性和。通常将结果间的量值表示在k的负值区间。（3）对称性：实序列的离散傅立叶变换可以表示为，其中实部为偶对称，虚部为奇对称，幅值为偶对称，相位为奇对称。如果序列是实偶对称序列，则也是实偶对称，即；如果序列

18、是实奇对称序列，则是虚奇对称，即；如果序列是虚偶对称序列，则也是虚偶对称，即；如果序列是虚奇对称序列，则是实奇对称，即。根据上述关系，对于实序列，则有；对于纯虚序列，则有。三、实验内容（1）将实指数函数抽样，作DFT，并作出实部、虚部和幅频、相频特性曲线。（2）将图3-2中的两个连续函数经抽样，作DFT，验证前述的四种奇偶特性，并作出幅频和相频特性曲线。0 1 2 te-atu(t)图3-1连续时间函数0 1 2 tx(t)图3-2 两个有限时间连续函数1 0 1 2 tx(t)1 -1 (a) (b) 实验四 信号处理中FFT的应用一、实验目的1、理解用FFT对周期序列进行频谱分析时所面临的

19、问题并掌握其解决方法。2、掌握用时域窗函数加权处理的技术。3、理解用FFT对非周期信号进行频谱分析所面临的问题并掌握其解决方法。二、实验原理与计算方法、对周期序列进行频谱分析应注意的问题 k k(a)时域周期整数倍截断 (b)时域非周期整数倍截断图 4-1 周期函数的幅频曲线对时间序列作FFT时，实际上要作周期延拓（如果取长序列的一段进行计算还要先作截断）。周期序列是无限长时间序列，如果截断区间刚好就是该序列周期的整数倍，那么在进行周期延拓后，将还原出原来的周期序列，由此可以较精确地计算出的该周期序列的频谱。反之，如果截断区间并不是该序列周期的整数倍，那么在进行周期延拓后，就不可能还原出原来的

20、周期序列，由此计算出的频谱与该周期序列的频谱存在误差，而且误差的大小与截断区间的选取直接相关，如图4-1所示，其中幅度频谱的量值表示为，以dB(分贝)为单位。这种误差是由于周期序列与矩形截断序列相乘在频域产生二者的频谱卷积形成的。矩形窗的频谱是抽样函数序列，如图4-2所示。除了k = 0处主瓣内集中了大部分的能量外，两旁的较小峰值处的旁瓣也分散了一部分能量，它与周期序列频谱卷积的结果使原来集中的频谱展宽，称为频率泄漏。 k图 4-2 矩形窗的频谱如果对已知周期的信号作频谱分析，在进行时域截断时，完全可以选取其周期的整数倍裁取，从而可以避免这种频率泄漏的发生。不过，通常需要进行频谱分析的信号是周

21、期未知的信号，或随机信号，无法判断它的周期值，为了尽量避免频率泄漏对结果的影响，在作时间截断时，就应选取其频谱的旁瓣较小的截断函数，以减轻泄漏问题。2、时域窗函数的应用作为截断函数，矩形窗在作时间截断时，对所截取区间内的信号不加以任何影响，而其它的窗函数都将对所截取区间内的信号作加权处理。除了在实验二中已经介绍过三角窗、Hanning窗和Hamming窗外，常用的窗函数还有很多，例如Parzen窗、Kaiser窗、Chebyshev窗、Tukey窗、Poisson窗、Caushy窗、Gaussian窗和Blackman窗等等。本次实验仍是采用实验二中的几种窗函数，但是利用窗函数作时域加权截断。

22、 0 t 0 k(a) 正弦函数的加权的非周期时域截断(b)减小了泄漏的频谱 图 4-3 采用Hanning窗加权后的时域截断和频谱图 4-3 中给出了采用Hanning窗对正弦函数作非整周期的时域加权截断后的波形和频谱，与图4-1(b)比较，泄漏已明显减少。应该指出，前面所给出的窗函数都是定义为以0点为中心、宽度为N +1的加权函数，在这里应用时，需要将其右移，成为区间内的加权函数。3、对非周期序列进行频谱分析应注意的问题（）混叠一般非周期信号作FFT之前要进行时域采样和周期延拓(无限长时间信号还应先截断再延拓)。根据Fourier变换理论，经等周期的冲激采样后，离散序列的频谱是原信号频谱以

23、为周期的周期延拓。按照Nyquist采样定理，由采样引起周期延拓后频谱之间不发生混叠的条件是：(1)原信号应该是有限带宽信号，设其频带宽度为fm；(2)频谱的周期，即采样周期应满足Nyquist 条件。 0 n 0 N/2 N (a)时域按周期Ts采样 (b)频域一个周期内在N/2附近出现混叠 图 4-4 非周期函数采样后的幅频曲线由于实际上有限长时间信号的FT是频域的无限函数，因此采样所得的离散序列的频谱必定产生混叠，减小采样周期只能减小而不能消除混叠。对于时间有限函数，当采样周期较大时，也会在FFT得到的频域出现混叠，形成频谱失真，造成频谱分析结果与原信号的实际频谱的差异，也无法恢复出原信

24、号。当然，实际工作中只要采样和截断产生的误差在许可的范围内就行了，但应该认识到混叠是引起频谱分析误差的一个主要原因。还应该注意的是，离散Fourier变换的频域也是周期化的，区间内的样点实际上是负频率区的量值，因此如果出现混叠，就将在一个周期内出现，并发生在附近的区域，如图4-4所示。要减少混叠，就要尽量减小采样周期。（）泄漏周期函数截断引起的频率泄漏问题，在非周期函数截断处理后同样存在，这种误差是由于采样后的离散序列与矩形截断序列相乘在频域造成二者的频谱卷积形成的。矩形窗的频谱是抽样函数序列，它与离散序列频谱卷积的结果使原来集中于每一个样点处的频谱展宽，其影响在高频区(接近N/2的样点)特别

25、明显，如图4-5所示。同样，为了尽量避免频率泄漏对结果的影响，在对非周期函数作时间截断时，除尽量增加截断序列的宽度外，也应选取其频谱的旁瓣较小的截断函数，以减轻泄漏问题。 0 n 0 N/2 N (a)时域截断 (b)频域一个周期内在N/2附近出现泄漏 图 4-5 函数采样后作截断的幅频曲线在选取了适当的窗函数后，应当使窗函数的宽度与被处理的序列长度相同，如果作变换前还需要补零（例如为了作卷积运算或避免栅栏效应），则应将原序列与窗函数相乘后再补零，即补零的样点不用窗函数加权处理。（）栅栏效应非周期信号应具有连续的频谱，在对作抽样后进行DFT，得到的是离散的频谱。如果排除混叠和泄漏等误差的影响，

26、所得的结果也只是的连续频谱上的个样值。这就象通过栅栏上的等间距缝隙观看到的另一边的景象，故此称栅栏效应。被栅栏遮住的景象中有可能存在与显现出的频谱差异较大的变化，即显示信号特征的频谱分量。为了使被栅栏遮住的部分能尽可能地显现出来，可以采用增加频域样点密度的方法，即在不增加信号采样点的情况下，用时域补零加宽变换尺度N来实现，称为补零重构。例如原来信号采样得到12个样点，在其后面再加上4个零，使序列的总长度为16个样点。这样处理的结果原来信号的采样间隔和频率都没有改变，设采样频率为，经补零重构之后，采样频率仍然为，但是原来频域样点间宽度为/12，经补零重构之后频域样点间宽度为/16。这就使补零重构

27、之后频域样点密度增加，而且显示出原来没有显露的一些频率位置的频谱。三、实验内容（1）将数字周期Ts=1/53，样点数为256的余弦序列作FFT，画出频谱曲线，观察并记录频率泄漏现象，然后用Hamming窗和三角窗作加权截断，观察并记录泄漏的衰减。（2）将幅度为1，周期为2的方波信号，按Ts=1/37 s的间距抽样，取样点数N=256作FFT，画出频谱曲线，然后用Hamming窗和三角窗作加权截断，观察并记录作不同的加权截断引起的频谱差异。（3）将单边指数函数抽样截断后作FFT，首先选取不同的抽样周期s，并取N = 256，观察频谱混叠。然后作不同宽度的截断，选取矩形窗宽为4，8，16，32等，

28、并保持N = 128，观察频谱泄漏。说明在什么情况下混叠和泄漏基本消失。（4）计算下面三个正弦函数的组合的频谱其中频率，令，抽样周期。分别取N=32,64,128将其截断后作FFT，观察和记录混叠和泄漏形态。分别采取补零加宽和增加截取时间宽度的方法作出频谱图，了解栅栏效应和频率分辨力的意义。实验五 离散系统的Z域分析一、实验目的1、掌握离散序列z变换的计算方法。2、掌握离散系统系统函数零极点的计算方法和零极点图的绘制方法，并能根据零极点图分析系统的因果性和稳定性。3、掌握利用Matlab进行z反变换的计算方法。二、实验原理与计算方法1、z变换离散序列x(n)的z变换定义为：。在Matlab中可

29、以利用符号表达式计算一个因果序列的z变换。其命令格式为：syms n; f=(1/2)n+(1/3)n;ztrans(f)2、离散系统的系统函数及因果稳定的系统应满足的条件一个线性移不变离散系统可以用它的单位抽样响应h(n)来表示其输入与输出关系，即y(n)= x(n)* h(n)对该式两边取z变换，得： Y(z)= X(z) H(z)则： 将H(z)定义为系统函数，它是单位抽样响应h(n)的z变换，即对于线性移不变系统，若nN=3;OmegaC=0.5; b,a=u_buttap(N,Omegac); C,B,A=sdir2cas(b,a)（2）按给定技术指标设计Butterworth模拟低

30、通滤波器函数b,a=afd_butt(Wp,Ws,Rp,As)用来实现按给定技术指标设计Butterworth模拟低通滤波器；函数db,mag,pha,w=freqs_m(b,a,0.5*pi)绘制模拟滤波器频响特性，得出衰减值,幅频特性,相频特性和自变量w；函数H,w=freqs(b,a,0.5*pi)绘制模拟滤波器频响特性,给出复函数H,再由abs(H)和angle(H)求出幅频和相频特性；ha,x,t=impulse(b,a)求出模拟滤波器冲激响应。Wp=0.2*pi;Ws=0.3*pi;Rp=7;As=16; b,a=afd_butt(Wp,Ws,Rp,As) db,mag,pha,w

31、=freqs_m(b,a,0.5*pi) ha,x,t=impulse(b,a)（3）数字滤波器设计Butterworth函数b,a=imp_incr(c,d,T)用来实现冲激响应不变的映射，c、d分别为模拟滤波器系统函数Ha(s)的分子、分母的系数；T为抽样周期；%数字滤波器技术指标wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;%对应的模拟滤波器技术指标T=1;Wp=wp/T;Ws=ws/T;cs,cd=afd_butt(Wp,Ws,Rp,As); b,a=imp_invr(cs,ds,T)%并联形式的系统函数系数C,B,A=dir2par(b,a)函数H,w=freqz(b

32、,a)绘制数字滤波器频响特性,给出复函数H,再由abs(H)和angle(H)求出幅频和相频特性。（4）Chebyshev原型Matlab提供了函数z,p,k=cheblap(N,Rp)用来设计N阶通带波动为Rp的归一化Chebyshev模拟低通滤波器；函数b,a=u_chb1ap(N,Omegac)给出未归一化的Chebyshev模拟低通滤波器原型。（5）按给定技术指标设计Chebysheb模拟低通滤波器函数b,a=afd_chb1(Wp,Ws,Rp,As)用来实现按给定技术指标设计Chebyshev模拟低通滤波器；函数db,mag,pha,w=freqs_m(b,a,0.5*pi)绘制频响特性；ha,x,t=impulse(b,a)求出模拟滤波器冲激响应。Wp=0.2*pi;Ws=0.3*pi;Rp=1;As=16; b,a=afd_chb1(Wp,Ws,Rp,As) db,mag,pha,w=freqs_m(b,a,0.5*pi) ha,x,t=impulse(b,a)（6）数字滤波器设计Chebyshev函数b,a