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

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=0。由于内存有限，Matlab不能表示一个无限序列。2、解析法对于有解析表达式的确定信号，首先定义序列的范围即n的值，然后直接写出该序列的表

3、达式，如：例1.2：实现实指数序列，的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);此时得到的图形的横坐标范围

4、由向量n的值决定，为-3到4，纵坐标的范围由向量x的值决定，为-1到7。应用stem函数时应确保自变量n和函数值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的帮助文件

5、。下面给出产生单位抽样序列和单位阶跃序列的两个函数，供参考。例1.5：产生单位抽样序列的函数impseq(n0,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 < n1) | (n0 > n2) | (n1 > n2)error('arguments must satisfy n1 <= n0 <= n2')endn = n1:n2;x = (n

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

7、:n2;%x = zeros(1,(n0-n1), ones(1,(n2-n0+1);x = (n-n0) >= 0;该函数产生一个起始位置n0位于n1和n2之间的单位阶跃序列。注意：由function产生的函数文件，不能直接运行，并且要放在当前路径下的文件夹里，供其他M文件调用。（二）序列的基本运算和时域变换1、加法：x1(n)+x2(n)序列的加法运算为对应位置处量值的相加，在Matlab中可用运算符“+”实现，但要求参与运算的序列的长度必须相等。如果长度不等或者长度相等但采样位置不同，则不能直接应用该运算符，此时需要先给定参数使序列具有相同的位置向量和长度。下面给出sigadd函数

8、实现任意两序列的加法运算。例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 can be different from n1)%n = min(min(n1),min(n2):max(max(n1),max(n2); %

9、duration of y(n)y1 = zeros(1,length(n); y2 = y1; % initializationy1(find(n>=min(n1)&(n<=max(n1)=1)=x1; % x1 with duration of yy2(find(n>=min(n2)&(n<=max(n2)=1)=x2; % x2 with duration of yy = y1+y2; % sequence addition其中x1和x2为参与加法运算的两序列，n1和n2分别为x1和x2的位置向量。2、乘法：x1(n)·x2(n)序列的乘

10、法运算为对应位置处量值的相乘，在Matlab中由数组运算符“.*”实现，也受到“+”运算符同样的限制。3、反折：x(n)x(-n)序列的反折指序列的每个量值都对n=0做一个对称操作，从而得到一个新序列。在Matlab中可由fliplr(x)函数实现，此时序列位置的反折则由-fliplr(n)实现。4、平移：x(n)x(n-m)平移操作是将序列的每个量值都移动m个位置，在得到的新序列中，量值和原序列相同，只是位置向量n发生变化，当m>0时，表示序列向右平移，此时新序列的位置向量为n+m；当m<0时，表示序列向左平移，此时新序列的位置向量为n-m。三、实验内容（1）参考示例程序，产生一

11、个有延迟的单位抽样序列：d (n-11)，绘出序列的图形。（2）参考示例程序，产生一个向前时移7个时刻的单位阶跃序列：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、进一步熟悉用Ma

12、tlab描绘二维图像的方法。二、实验原理与计算方法两个离散序列x(n)与y(n)的卷积和f(n)定义为由于通常信号处理中所碰到的都是有始信号或有限时间信号，因此在实际计算卷积和时，求和是在有限范围内进行的。计算过程中上下限的选取和所得结果的分布区间取决于参与卷积的两个序列，下面将分别进行讨论：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 =

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

14、、实验内容(1) 根据（1）式计算两个从n = 0开始的序列和的卷积和，其中A=40，a = 0.5。取50个样值点，作出序列、及卷积和f(n)的图像。(2) 根据（3）式计算两个有限长序列和的卷积和，其中作出序列、及卷积和f(n)的图像。实验三 离散傅立叶变换及其特性验证一、实验目的1、掌握离散时间傅立叶变换（DTFT）的计算方法和编程技术。2、掌握离散傅立叶变换（DFT）的计算方法和编程技术。3、理解离散傅立叶变换（DFT）的性质并用Matlab进行验证。二、实验原理与计算方法1、离散时间傅立叶变换如果序列x(n)满足绝对可和的条件，即，则其离散时间傅立叶变换定义为： （1）如果x(n)是

15、无限长的，则不能直接用Matlab由x(n)计算X(ejw)，但可以用它来估计X(ejw)表达式在0,频率区间的值并绘制它的幅频和相频（或实部和虚部）曲线。如果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维矩阵：

16、将k和n排成列向量，则在Matlab中，把序列和下标排成行向量，对（3）式取转置得：其中nTk是一个N×(M+1)维矩阵。用Matlab实现如下：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%

17、-% Xk = dft(xn,N)% Xk = DFT coeff. array over 0 <= k <= N-1% xn = N-point 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

18、 = xn * WNnk; % row vector for DFT coefficientsfunction xn = idft(Xk,N)% Computes Inverse 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

19、for kWN = exp(-j*2*pi/N); % Wn factornk = n'*k; % creates a N by N matrix of nk valuesWNnk = WN . (-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 具

20、有特性和。通常将结果间的量值表示在k的负值区间。（3）对称性：实序列的离散傅立叶变换可以表示为，其中实部为偶对称，虚部为奇对称，幅值为偶对称，相位为奇对称。如果序列是实偶对称序列，则也是实偶对称，即；如果序列是实奇对称序列，则是虚奇对称，即；如果序列是虚偶对称序列，则也是虚偶对称，即；如果序列是虚奇对称序列，则是实奇对称，即。根据上述关系，对于实序列，则有；对于纯虚序列，则有。三、实验内容（1）将实指数函数抽样，作DFT，并作出实部、虚部和幅频、相频特性曲线。（2）将图3-2中的两个连续函数经抽样，作DFT，验证前述的四种奇偶特性，并作出幅频和相频特性曲线。0 1 2 te-atu(t)图3-

21、1连续时间函数0 1 2 tx(t)图3-2 两个有限时间连续函数1 0 1 2 tx(t)1 -1 (a) (b) 实验四 信号处理中FFT的应用一、实验目的1、理解用FFT对周期序列进行频谱分析时所面临的问题并掌握其解决方法。2、掌握用时域窗函数加权处理的技术。3、理解用FFT对非周期信号进行频谱分析所面临的问题并掌握其解决方法。二、实验原理与计算方法、对周期序列进行频谱分析应注意的问题 k k(a)时域周期整数倍截断 (b)时域非周期整数倍截断图 4-1 周期函数的幅频曲线对时间序列作FFT时，实际上要作周期延拓（如果取长序列的一段进行计算还要先作截断）。周期序列是无限长时间序列，如果截

22、断区间刚好就是该序列周期的整数倍，那么在进行周期延拓后，将还原出原来的周期序列，由此可以较精确地计算出的该周期序列的频谱。反之，如果截断区间并不是该序列周期的整数倍，那么在进行周期延拓后，就不可能还原出原来的周期序列，由此计算出的频谱与该周期序列的频谱存在误差，而且误差的大小与截断区间的选取直接相关，如图4-1所示，其中幅度频谱的量值表示为，以dB(分贝)为单位。这种误差是由于周期序列与矩形截断序列相乘在频域产生二者的频谱卷积形成的。矩形窗的频谱是抽样函数序列，如图4-2所示。除了k = 0处主瓣内集中了大部分的能量外，两旁的较小峰值处的旁瓣也分散了一部分能量，它与周期序列频谱卷积的结果使原来

23、集中的频谱展宽，称为频率泄漏。 k图 4-2 矩形窗的频谱如果对已知周期的信号作频谱分析，在进行时域截断时，完全可以选取其周期的整数倍裁取，从而可以避免这种频率泄漏的发生。不过，通常需要进行频谱分析的信号是周期未知的信号，或随机信号，无法判断它的周期值，为了尽量避免频率泄漏对结果的影响，在作时间截断时，就应选取其频谱的旁瓣较小的截断函数，以减轻泄漏问题。2、时域窗函数的应用作为截断函数，矩形窗在作时间截断时，对所截取区间内的信号不加以任何影响，而其它的窗函数都将对所截取区间内的信号作加权处理。除了在实验二中已经介绍过三角窗、Hanning窗和Hamming窗外，常用的窗函数还有很多，例如Par

24、zen窗、Kaiser窗、Chebyshev窗、Tukey窗、Poisson窗、Caushy窗、Gaussian窗和Blackman窗等等。本次实验仍是采用实验二中的几种窗函数，但是利用窗函数作时域加权截断。 0 t 0 k(a) 正弦函数的加权的非周期时域截断(b)减小了泄漏的频谱 图 4-3 采用Hanning窗加权后的时域截断和频谱图 4-3 中给出了采用Hanning窗对正弦函数作非整周期的时域加权截断后的波形和频谱，与图4-1(b)比较，泄漏已明显减少。应该指出，前面所给出的窗函数都是定义为以0点为中心、宽度为N +1的加权函数，在这里应用时，需要将其右移，成为区间内的加权函数。3、

25、对非周期序列进行频谱分析应注意的问题（）混叠一般非周期信号作FFT之前要进行时域采样和周期延拓(无限长时间信号还应先截断再延拓)。根据Fourier变换理论，经等周期的冲激采样后，离散序列的频谱是原信号频谱以为周期的周期延拓。按照Nyquist采样定理，由采样引起周期延拓后频谱之间不发生混叠的条件是：(1)原信号应该是有限带宽信号，设其频带宽度为fm；(2)频谱的周期，即采样周期应满足Nyquist 条件。 0 n 0 N/2 N (a)时域按周期Ts采样 (b)频域一个周期内在N/2附近出现混叠 图 4-4 非周期函数采样后的幅频曲线由于实际上有限长时间信号的FT是频域的无限函数，因此采样所

26、得的离散序列的频谱必定产生混叠，减小采样周期只能减小而不能消除混叠。对于时间有限函数，当采样周期较大时，也会在FFT得到的频域出现混叠，形成频谱失真，造成频谱分析结果与原信号的实际频谱的差异，也无法恢复出原信号。当然，实际工作中只要采样和截断产生的误差在许可的范围内就行了，但应该认识到混叠是引起频谱分析误差的一个主要原因。还应该注意的是，离散Fourier变换的频域也是周期化的，区间内的样点实际上是负频率区的量值，因此如果出现混叠，就将在一个周期内出现，并发生在附近的区域，如图4-4所示。要减少混叠，就要尽量减小采样周期。（）泄漏周期函数截断引起的频率泄漏问题，在非周期函数截断处理后同样存在，

27、这种误差是由于采样后的离散序列与矩形截断序列相乘在频域造成二者的频谱卷积形成的。矩形窗的频谱是抽样函数序列，它与离散序列频谱卷积的结果使原来集中于每一个样点处的频谱展宽，其影响在高频区(接近N/2的样点)特别明显，如图4-5所示。同样，为了尽量避免频率泄漏对结果的影响，在对非周期函数作时间截断时，除尽量增加截断序列的宽度外，也应选取其频谱的旁瓣较小的截断函数，以减轻泄漏问题。 0 n 0 N/2 N (a)时域截断 (b)频域一个周期内在N/2附近出现泄漏 图 4-5 函数采样后作截断的幅频曲线在选取了适当的窗函数后，应当使窗函数的宽度与被处理的序列长度相同，如果作变换前还需要补零（例如为了作

28、卷积运算或避免栅栏效应），则应将原序列与窗函数相乘后再补零，即补零的样点不用窗函数加权处理。（）栅栏效应非周期信号应具有连续的频谱，在对作抽样后进行DFT，得到的是离散的频谱。如果排除混叠和泄漏等误差的影响，所得的结果也只是的连续频谱上的个样值。这就象通过栅栏上的等间距缝隙观看到的另一边的景象，故此称栅栏效应。被栅栏遮住的景象中有可能存在与显现出的频谱差异较大的变化，即显示信号特征的频谱分量。为了使被栅栏遮住的部分能尽可能地显现出来，可以采用增加频域样点密度的方法，即在不增加信号采样点的情况下，用时域补零加宽变换尺度N来实现，称为补零重构。例如原来信号采样得到12个样点，在其后面再加上4个零，

29、使序列的总长度为16个样点。这样处理的结果原来信号的采样间隔和频率都没有改变，设采样频率为，经补零重构之后，采样频率仍然为，但是原来频域样点间宽度为/12，经补零重构之后频域样点间宽度为/16。这就使补零重构之后频域样点密度增加，而且显示出原来没有显露的一些频率位置的频谱。三、实验内容（1）将数字频率f=1/53，样点数为256的余弦序列作FFT，画出频谱曲线，观察并记录频率泄漏现象，然后用Hamming窗和三角窗作加权截断，观察并记录泄漏的衰减。（2）将幅度为1，周期为2的方波信号，按Ts=1/37 s的间距抽样，取样点数N=256作FFT，画出频谱曲线，然后用Hamming窗和三角窗作加权

30、截断，观察并记录作不同的加权截断引起的频谱差异。（3）将单边指数函数抽样截断后作FFT，首先选取不同的抽样周期s，并取N = 256，观察频谱混叠。然后作不同宽度的截断，选取矩形窗宽为4，8，16，32等，并保持N = 128，观察频谱泄漏。说明在什么情况下混叠和泄漏基本消失。（4）计算下面三个正弦函数的组合的频谱其中频率，令，抽样周期。分别取N=32,64,128将其截断后作FFT，观察和记录混叠和泄漏形态。分别采取补零加宽和增加截取时间宽度的方法作出频谱图，了解栅栏效应和频率分辨力的意义。实验五 离散系统的Z域分析一、实验目的1、掌握离散序列z变换的计算方法。2、掌握离散系统系统函数零极点

31、的计算方法和零极点图的绘制方法，并能根据零极点图分析系统的因果性和稳定性。3、掌握利用Matlab进行z反变换的计算方法。二、实验原理与计算方法1、z变换离散序列x(n)的z变换定义为：。在Matlab中可以利用符号表达式计算一个因果序列的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)定义为系统函数，它是单位抽样

32、响应h(n)的z变换，即对于线性移不变系统，若n<0时，h(n)=0，则系统为因果系统；若，则系统稳定。由于h(n)为因果序列，所以H(z)的收敛域为收敛圆外部区域，因此H(z)的收敛域为收敛圆外部区域时，系统为因果系统。因为，若z=1时H(z)收敛，即，则系统稳定，即H(z)的收敛域包括单位圆时，系统稳定。因此因果稳定系统应满足的条件为：，即系统函数H(z)的所有极点全部落在z平面的单位圆之内。3、Matlab中系统函数零极点的求法及零极点图的绘制方法Matlab中系统函数的零点和极点可以用多项式求根函数roots ()来实现，调用该函数的命令格式为：p=roots(A)。其中A为待求

33、根多项式的系数构成的行向量，返回向量p是包含该多项式所有根位置的列向量。如：求多项式的根的Matlab命令为：A=1 3/4 1/8;p=roots(A)运行结果为：p=-0.5000-0.2500需要注意的是，离散系统的系统函数可能有两种形式，一种是分子和分母多项式均按z的降幂次序排列，另一种是分子分母多项式均按z-1的升幂次序排列，两种方式在构造多项式系数向量时稍有不同。若H(z)是按z的降幂次序排列，则系数向量一定要由多项式的最高幂次开始，一直到常数项，缺项用0补齐，如，其分子多项式的系数向量应为：B=1 0 2 0；分母多项式的系数向量应为：A=1 3 2 2 1。若H(z)是按z-1

34、的升幂次序排列，则分子和分母多项式系数向量的维数一定要相同，缺项用0补齐，否则零点和极点就可能被漏掉。如，其分子多项式的系数向量应为：B=1 1 0；分母多项式的系数向量应为：A=1 1/2 1/4。可利用Matlab中的zplane函数实现系统函数的零极点图的绘制。该函数的调用方法为：zplane(B,A);其中B、A为系统函数分子分母多项式的系数向量。4、z反变换的计算方法z反变换可由部分分式展开法求得。由于指数序列anu(n)的z变换为，因此求反变换时，通常对进行展开：其中称为有理函数的留数。分两种情况进行讨论：（1）X(z)的所有极点均为单实极点此时，则X(z)的z反变换为：（2）X(

35、z)有共轭极点设X(z)有一对共轭极点，则，其中留数的计算方法与单极点相同，即，r2=r1 *因此，只要求出部分分式展开的系数（留数），就可以直接求出X(z)的z反变换x(n)。在Matlab中可利用函数residue()求解。令B和A分别是的分子和分母多项式构成的系数向量，则函数r,p,k=residue (B,A)将产生三个向量r、p、k，其中r为包含部分分式展开系数ri(i=1,2,N)的列向量，p为包含所有极点的行向量，k为包含部分分式展开的多项式项的系数cj(j=1,2,M-N)的列向量，若MN，则k为空阵。用residue()函数求出部分分式展开的系数后，便可根据其极点位置分布情况

36、直接求出X(z)的反变换x(n)。如：已知，求其z反变换x(n)。首先利用residue()函数求出的部分分式展开的系数和极点，相应的Matlab命令为：B=0 1 0;A=1 3 2;r,p,k=residue (B,A)运行结果为：r = 2 -1p = -2 -1k = 由以上结果可得：；即X(z)只有两个单极点，其z反变换为：。已知，求其z反变换x(n)。利用residue()函数求出的部分分式展开的系数和极点，可得：B=0 0 1 1;A=1 -2 2 -1;r,p,k=residue (B,A)r = 2.0000 -1.0000 + 0.0000i -1.0000 - 0.000

37、0ip = 1.0000 0.5000 + 0.8660i 0.5000 - 0.8660ik = 可见，包含一对共轭极点，用abs()和angle()函数即可求出共轭极点的模和相位，相应命令为：p1=abs(p')p1 = 1.0000 1.0000 1.0000a1=angle(p')/pia1 = 0 -0.3333 0.3333即共轭极点为：，则，其z反变换为：三、实验内容（1）求下列序列的z变换：2-nu(n)；-(1/2)n；(1/2)n+(1/3)n（2）已知两个离散系统的系统函数分别为： ；分别求出各系统的零极点，绘制零极点图，分析系统的稳定性；求出各系统单位抽

38、样响应。实验六 无限冲激响应(IIR)数字滤波器的三种结构一、实验目的1、掌握IIR滤波器的直接II型、级联型和并联型三种结构的基本原理和特点。2、掌握利用Matlab实现IIR滤波器的三种结构的程序设计方法，并能够进行三者之间的相互转换。3、掌握滤波器频响特性的绘制方法。二、实验原理与计算方法按照结构划分数字滤波器，有递归式和非递归式两种。递归式数字滤波器的差分方程为 (1)其中至少有一个.非递归式数字滤波器的差分方程为 (2)可以看出递归式数字滤波器的响应不仅与激励有关，而且与以前的输出信号有关；而非递归式数字滤波器的响应仅只与激励有关。按照单位样值响应划分数字滤波器，则有无限冲激响应(I

39、IR)和有限冲激响应(FIR)之分。IIR滤波器是递归式的，差分方程如(1)式所示，FIR滤波器一般是非递归式的，差分方程如(2)式所示。IIR滤波器常用的典型结构有直接II型、级联型和并联型，分别介绍如下：、直接II型（也称为正准型结构）根据(1)式，IIR滤波器的传输函数为 (3)其中已假设(1)式中的，对于其它情况，则可令相应的某些系数为零。令则有(4)由此可以得到相应的时域中激励与响应之间的关系为 (5)其中是与(4)式中的相应的中间函数序列。x(n) y2(n) b0 y(n) z-1 a1 b1 z-1 a2 b2 z-1 aN-1 bN-1 z-1 aN bN图 6-1 直接II

40、型结构信号流由(5)式确定的直接II型的信号流图如图6-1所示，其中将中间的两条延时链合并为一条，实际的信号流将按（5）式分成两个延时链独立运行。编程时，用三个数组分别存放系数和，。由图可以看出，沿中间的延时链自上向下传播过程中将逐级向右移位，每一级向左边与的乘积按累加，再和相加，得到的沿中间的延时链又自上向下传播逐级向右移位，每一级向右边与的乘积按累加，其结果就是对应于激励的响应。然后向右移位一个单位时间，输入激励计算下一个时刻的响应。直接II型结构具有简单直观的典型网络结构形式，在计算机上很容易实现。但是它对系数的量值变化比较敏感，直接确定了系统零、极点的位置，从而影响到系统的性能。尤其当

41、阶数N较高时，系统对系数的字长效应很敏感，产生的误差也较大。、级联型结构由于当直接II型结构传输函数阶数增加时，系数量化引起的误差影响到滤波器的性能，因此要采用其它形式的结构。将(3)式的传输函数分子和分母进行因式分解，即用它的零、极点表示为 (6)其中，分子和分母中的实系数二阶因子分别对应于共轭零、极点。可以将上式分子和分母中单根一阶因子作为二阶因子的一种特例，那么(6)式就可以表示为M个实系数二阶基本节级联的形式 (7) (8) x(n) h1(n) b01 hM(n) b0M y(n) z-1 z-1 a11 b11 a1M b1M z-1 z-1 a21 b21 a2M b2M 图 6

42、-2 级联型结构信号流其中为滤波器的二阶基本节，为取整结果。则分别是第k个基本节的分母、分子系数，它们只关系到滤波器的某一对零、极点。级联型结构如图6-2所示，它是一系列二阶基本节的级联，每一个可以用直接II型结构实现。级联型结构的特点是对滤波器性能的调整比较方便，调整系数，只单独涉及到第k级零、极点，而不会影响到其它任一级的零、极点，因而可以独立地控制滤波器的各零、极点的分布。、并联型结构这种结构将传输函数展开为部分分式，即表示为若干一阶和二阶基本节网络与一个常数之和 (9)其中，同样也可以统一表示为二阶基本节的形式(10) B0x(n) y1(n) b01 y(n) z-1 a11 b11

43、 z-1 a21 yM(n) z-1 a1M b1M z-1 a2M 图 6-3 并联型结构信号流并联型结构信号流如图6-3所示，其中二阶基本节网络可以用直接II型结构实现，程序设计也可参考直接型II结构的方法。并联型结构也可以单独调整极点位置，但却不能象级联型结构那样直接控制零点的分布。因为并联型结构各二阶基本节网络的零点并不是整个系统函数的零点。因此，当要准确传输零点时，以采用级联型结构为宜。不过，由于并联型基本节之间互不影响，所以运算误差比级联型的要小一些。通常，IIR滤波器用系统函数的有理式（直接结构形式）描述，Matalb中提供相应的函数可以把直接型结构转换成级联和并联型结构。在Ma

44、tlab中，直接型结构由两个行向量描述，b包含系数bk，a包含系数ak，可由filter函数实现。要由已知的直接型结构实现级联型结构，可利用函数dir2cas根据直接型的系数bk和ak得到系数b0，Bk,i和Ak,i，该函数把矢量b和a转换成K×3维矩阵B和A。首先计算b0，它等于b0/ a0，（a01），然后通过给较短的矢量添零，使矢量b和a一样长。这就保证每一个双二阶环节的分子和分母均不为零。接着计算多项式B(z)和A(z)的根，用cplxpair函数把这些根以共轭复根对的次序排列，最后用poly函数把每一对根再转换成二阶分子或分母多项式。级联型由casfiltr函数实现，它把每

45、个二阶环节的系数存放在矩阵B和A中，放在一个循环里，而在此循环中采用filter函数。输入被乘以b0，每个滤波器的输出作为下一级滤波器的输入，最后一个滤波器的输出即为总的输出。函数cas2dir可以把级联形式转换成直接形式。这是一个包含几次多项式乘法的简单运算，在一个K次的循环中实用conv函数。类似的，利用dir2par函数可以把直接型系数bk和ak转换成并联型系数Bk,i和Ak,i，并在parfiltr函数中使用这些系数以实现并联形式。而par2dir函数则可用来实现并联形式转换成直接形式。三种结构的IIR滤波器均可通过以下方式计算单位抽样响应。delta=impseq(0,0,7);hd

46、ir=filter(b,a,delta)hcas=casfiltr(b0,B,A,delta)hpar=parfiltr(C,B,A,delta)传输函数为的滤波器的频率特性为。由于我们所讨论的滤波器都是稳定系统，其所有极点都在z平面的单位圆内，单位样值响应是一个无限长的衰减序列，所以要截取有限长的一段来作频率特性分析。截取的长度越长，逼近的程度越高。如果截取长度为N，滤波器的频率特性由Fourier变换得 (11)则滤波器的幅频特性和相频特性可表示为 (12) (13)三、实验内容设三阶滤波器的传输函数为激励信号为根据所给定的滤波器系统函数和参数值，用三种结构实现该IIR数字滤波器结构，分别

47、使激励和通过该滤波器，求出相应于不同激励的响应，作出输出响应信号的时间曲线和幅频、相频特性曲线，并判断为何种滤波器(低通、高通、带通)。实验七 冲激响应不变法IIR数字滤波器设计一、实验目的1、掌握构成一个频率响应与给定的滤波特性相接近的模拟滤波器的设计原理。2、掌握用冲激响应不变法设计IIR数字滤波器的基本原理和算法。3、了解数字滤波器和模拟滤波器的频率响应特性，掌握相应的计算方法，分析用冲激响应不变法获得的数字滤波器频率响应特性中出现的混叠现象。二、实验原理与计算方法、冲激响应不变法设计IIR数字滤波器的基本原理和算法采用冲激响应不变法设计数字滤波器，就是使其单位样值响应与相应的模拟滤波器的冲激响应在抽样点处的量值相等，即 (1)其中T为抽样周期。因此用冲激响应不变法设计IIR数字滤波器的基本步骤，就是首先根据设计要求确定相应的模拟滤波器的传递函数，经Laplace反变换求出冲激响应后，对它进行抽样得到的等于数字滤波器的单位样值响应，再经z变换所得就是数字滤波器的传递函数。如果模拟滤波器的传递函数的N个极点都是单极点，则可以将写成部分分式展开的形式 (2)那么，经Laplace反变换求出的模拟滤波器的冲激响应为相对应的数字滤波器的单位样值响应为对上式作z变换，得 (3)由上面的推导可见，只要模拟滤波器的传递函数的N个极点都是单极点，当已经求出