數(shù)字信號處理實驗指導(dǎo)書(M)

1、數(shù)字信號處理實驗電子信息科學(xué)與技術(shù)實驗室2007年7月目錄實驗一 離散時間信號的時域表示2實驗二 離散信號的卷積和5實驗三 離散傅立葉變換及其特性驗證7實驗四 信號處理中FFT的應(yīng)用10實驗五 離散系統(tǒng)的Z域分析15實驗六 無限沖激響應(yīng)(IIR)數(shù)字濾波器的三種結(jié)構(gòu)16實驗七 沖激響應(yīng)不變法IIR數(shù)字濾波器設(shè)計20實驗八 雙線性變換法IIR數(shù)字濾波器設(shè)計23實驗一 離散時間信號的時域表示一、實驗?zāi)康?、熟悉Matlab命令，掌握離散時間信號序列的時域表示方法。2、掌握用Matlab描繪二維圖像的方法。3、掌握用Matlab對序列進(jìn)行基本的運算和時域變換的方法。二、實驗原理與計算方法（一）序列的

2、表示方法序列的表示方法有列舉法、解析法和圖形法，相應(yīng)的用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;如果不對向量的位置進(jìn)行定義，則Matlab默認(rèn)該序列的起始位置為n=0。由于內(nèi)存有限，Matlab不能表示一個無限序列。2、解析法對于有解析表達(dá)式的確定信號，首先定義序列的范圍即n的值，然后直接寫出該序列的表

3、達(dá)式，如：例1.2：實現(xiàn)實指數(shù)序列，的Matlab程序為：>>n=0:10;>>x=(0.9).n;例1.3：實現(xiàn)正余弦序列，的Matlab程序為：>>n=5:15;>>x=3*cos(0.1*pi*n+pi/3)+2*sin(0.5*pi*n);3、圖形法在Matlab中用圖形法表示一個序列，是在前兩種表示方法的基礎(chǔ)上將序列的各個量值描繪出來，即首先對序列進(jìn)行定義，然后用相應(yīng)的畫圖語句畫圖，如：例1.4：繪制在例1.1中用列舉法表示的序列的圖形，則在向量定義之后加如下相應(yīng)的繪圖語句：>>stem(n,x);此時得到的圖形的橫坐標(biāo)范圍

4、由向量n的值決定，為-3到4，縱坐標(biāo)的范圍由向量x的值決定，為-1到7。應(yīng)用stem函數(shù)時應(yīng)確保自變量n和函數(shù)值x的個數(shù)相等。此外可用函數(shù)axis(x1,x2,y1,y2)對橫縱坐標(biāo)進(jìn)行限定，以完善圖形，其中x1和x2分別為橫坐標(biāo)的起始和截止位置，y1和y2分別為縱坐標(biāo)的起始和截止位置。也可用xlabel()、ylabel()和title()為該圖添加橫、縱坐標(biāo)說明和標(biāo)題。subplot(m,n,k)函數(shù)可以將當(dāng)前窗口分成m行n列個子窗口，并在第k的子窗口繪圖。窗口的排列順序為從左至右，從上至下分別為1,2,m*n。以上為各個繪圖函數(shù)的基本用法，有關(guān)各函數(shù)的其他參數(shù)可參考Matlab的幫助文件

5、。下面給出產(chǎn)生單位抽樣序列和單位階躍序列的兩個函數(shù)，供參考。例1.5：產(chǎn)生單位抽樣序列的函數(shù)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;該函數(shù)產(chǎn)生一個抽樣位置n0位于n1和n2之間的單位抽樣序列。例1.6：產(chǎn)生單位階躍序列的函數(shù)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;該函數(shù)產(chǎn)生一個起始位置n0位于n1和n2之間的單位階躍序列。注意：由function產(chǎn)生的函數(shù)文件，不能直接運行，并且要放在當(dāng)前路徑下的文件夾里，供其他M文件調(diào)用。（二）序列的基本運算和時域變換1、加法：x1(n)+x2(n)序列的加法運算為對應(yīng)位置處量值的相加，在Matlab中可用運算符“+”實現(xiàn)，但要求參與運算的序列的長度必須相等。如果長度不等或者長度相等但采樣位置不同，則不能直接應(yīng)用該運算符，此時需要先給定參數(shù)使序列具有相同的位置向量和長度。下面給出sigadd函數(shù)

8、實現(xiàn)任意兩序列的加法運算。例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、法運算為對應(yīng)位置處量值的相乘，在Matlab中由數(shù)組運算符“.*”實現(xiàn)，也受到“+”運算符同樣的限制。3、反折：x(n)x(-n)序列的反折指序列的每個量值都對n=0做一個對稱操作，從而得到一個新序列。在Matlab中可由fliplr(x)函數(shù)實現(xiàn)，此時序列位置的反折則由-fliplr(n)實現(xiàn)。4、平移：x(n)x(n-m)平移操作是將序列的每個量值都移動m個位置，在得到的新序列中，量值和原序列相同，只是位置向量n發(fā)生變化，當(dāng)m>0時，表示序列向右平移，此時新序列的位置向量為n+m；當(dāng)m<0時，表示序列向左平移，此時新序列的位置向量為n-m。三、實驗內(nèi)容（1）參考示例程序，產(chǎn)生一

11、個有延遲的單位抽樣序列：d (n-11)，繪出序列的圖形。（2）參考示例程序，產(chǎn)生一個向前時移7個時刻的單位階躍序列：u(n+7)，繪出序列的圖形。（3）產(chǎn)生一個指數(shù)為-0.1+(pi/6)*in的復(fù)指數(shù)序列，并繪出序列的實部、虛部、幅度和相位的波形。（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）代數(shù)運算符號和.的區(qū)別是？實驗二 離散信號的卷積和一、實驗?zāi)康?、掌握兩個離散信號卷積和的計算方法和編程技術(shù)。2、進(jìn)一步熟悉用Ma

12、tlab描繪二維圖像的方法。二、實驗原理與計算方法兩個離散序列x(n)與y(n)的卷積和f(n)定義為由于通常信號處理中所碰到的都是有始信號或有限時間信號，因此在實際計算卷積和時，求和是在有限范圍內(nèi)進(jìn)行的。計算過程中上下限的選取和所得結(jié)果的分布區(qū)間取決于參與卷積的兩個序列，下面將分別進(jìn)行討論：1、兩個從n = 0開始的序列和的卷積和 (1)上式右邊因子u(n)表示卷積和的結(jié)果也是一個從n = 0開始的序列。2、從n = n1開始的序列和從n = n2開始的序列的卷積和，其中n1和n2為任意整數(shù)。 (2)上式右邊因子u(n-n1-n2)表示卷積和是一個從n = n1+n2開始的序列。3、從n =

13、 n1開始的長度為N1的加窗序列和從n = n2開始的長度為N2的加窗序列的卷積和，其中則 (3)所得卷積和也是一個加窗序列，從n = n1+ n2開始，長度為N1+ N2-1。Matlab提供了一個內(nèi)部函數(shù)conv(x,h)用來計算兩個有限長度序列的卷積，該函數(shù)得到的卷積結(jié)果默認(rèn)從n=0開始，因此當(dāng)參與卷積的兩個序列的起始位置不是n=0時，則由該函數(shù)得到的計算結(jié)果將出現(xiàn)錯誤，此時需要重新定義結(jié)果的位置向量。由以上卷積運算的原理可知，兩有限長序列卷積后仍為有限長序列，長度為兩序列長度之和減1，結(jié)果的起始位置為兩序列起始位置之和，截止位置為兩序列截至位置之和。據(jù)此，可以得到卷積結(jié)果的位置向量。三

14、、實驗內(nèi)容(1) 根據(jù)（1）式計算兩個從n = 0開始的序列和的卷積和，其中A=40，a = 0.5。取50個樣值點，作出序列、及卷積和f(n)的圖像。(2) 根據(jù)（3）式計算兩個有限長序列和的卷積和，其中作出序列、及卷積和f(n)的圖像。實驗三 離散傅立葉變換及其特性驗證一、實驗?zāi)康?、掌握離散時間傅立葉變換（DTFT）的計算方法和編程技術(shù)。2、掌握離散傅立葉變換（DFT）的計算方法和編程技術(shù)。3、理解離散傅立葉變換（DFT）的性質(zhì)并用Matlab進(jìn)行驗證。二、實驗原理與計算方法1、離散時間傅立葉變換如果序列x(n)滿足絕對可和的條件，即，則其離散時間傅立葉變換定義為： （1）如果x(n)是

15、無限長的，則不能直接用Matlab由x(n)計算X(ejw)，但可以用它來估計X(ejw)表達(dá)式在0,頻率區(qū)間的值并繪制它的幅頻和相頻（或?qū)嵅亢吞摬浚┣€。如果x(n)是有限長的，則可以用Matlab對任意頻率w處的X(ejw)進(jìn)行數(shù)值計算。如果要在0,間按等間隔頻點估計X(ejw)，則（1）式可以用矩陣向量相乘的運算來實現(xiàn)。假設(shè)序列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中，把序列和下標(biāo)排成行向量，對（3）式取轉(zhuǎn)置得：其中nTk是一個N×(M+1)維矩陣。用Matlab實現(xiàn)如下：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函數(shù)dft和idft實現(xiàn)離散傅立葉變換運算。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、離散傅立葉變換的性質(zhì)（1）線性性質(zhì)：注意：若x1(n)和x2(n)分別是N1點和N2點的序列，則選擇N3= max (N1, N2)，將它們作N3點DFT處理。（2） 周期性：離散傅立葉變換(DFT)是周期序列DFS取主值區(qū)間形成的，因此序列及其DFT 具

20、有特性和。通常將結(jié)果間的量值表示在k的負(fù)值區(qū)間。（3）對稱性：實序列的離散傅立葉變換可以表示為，其中實部為偶對稱，虛部為奇對稱，幅值為偶對稱，相位為奇對稱。如果序列是實偶對稱序列，則也是實偶對稱，即；如果序列是實奇對稱序列，則是虛奇對稱，即；如果序列是虛偶對稱序列，則也是虛偶對稱，即；如果序列是虛奇對稱序列，則是實奇對稱，即。根據(jù)上述關(guān)系，對于實序列，則有；對于純虛序列，則有。三、實驗內(nèi)容（1）將實指數(shù)函數(shù)抽樣，作DFT，并作出實部、虛部和幅頻、相頻特性曲線。（2）將圖3-2中的兩個連續(xù)函數(shù)經(jīng)抽樣，作DFT，驗證前述的四種奇偶特性，并作出幅頻和相頻特性曲線。0 1 2 te-atu(t)圖3-

21、1連續(xù)時間函數(shù)0 1 2 tx(t)圖3-2 兩個有限時間連續(xù)函數(shù)1 0 1 2 tx(t)1 -1 (a) (b) 實驗四 信號處理中FFT的應(yīng)用一、實驗?zāi)康?、理解用FFT對周期序列進(jìn)行頻譜分析時所面臨的問題并掌握其解決方法。2、掌握用時域窗函數(shù)加權(quán)處理的技術(shù)。3、理解用FFT對非周期信號進(jìn)行頻譜分析所面臨的問題并掌握其解決方法。二、實驗原理與計算方法、對周期序列進(jìn)行頻譜分析應(yīng)注意的問題 k k(a)時域周期整數(shù)倍截斷 (b)時域非周期整數(shù)倍截斷圖 4-1 周期函數(shù)的幅頻曲線對時間序列作FFT時，實際上要作周期延拓（如果取長序列的一段進(jìn)行計算還要先作截斷）。周期序列是無限長時間序列，如果截

22、斷區(qū)間剛好就是該序列周期的整數(shù)倍，那么在進(jìn)行周期延拓后，將還原出原來的周期序列，由此可以較精確地計算出的該周期序列的頻譜。反之，如果截斷區(qū)間并不是該序列周期的整數(shù)倍，那么在進(jìn)行周期延拓后，就不可能還原出原來的周期序列，由此計算出的頻譜與該周期序列的頻譜存在誤差，而且誤差的大小與截斷區(qū)間的選取直接相關(guān)，如圖4-1所示，其中幅度頻譜的量值表示為，以dB(分貝)為單位。這種誤差是由于周期序列與矩形截斷序列相乘在頻域產(chǎn)生二者的頻譜卷積形成的。矩形窗的頻譜是抽樣函數(shù)序列，如圖4-2所示。除了k = 0處主瓣內(nèi)集中了大部分的能量外，兩旁的較小峰值處的旁瓣也分散了一部分能量，它與周期序列頻譜卷積的結(jié)果使原來

23、集中的頻譜展寬，稱為頻率泄漏。 k圖 4-2 矩形窗的頻譜如果對已知周期的信號作頻譜分析，在進(jìn)行時域截斷時，完全可以選取其周期的整數(shù)倍裁取，從而可以避免這種頻率泄漏的發(fā)生。不過，通常需要進(jìn)行頻譜分析的信號是周期未知的信號，或隨機(jī)信號，無法判斷它的周期值，為了盡量避免頻率泄漏對結(jié)果的影響，在作時間截斷時，就應(yīng)選取其頻譜的旁瓣較小的截斷函數(shù)，以減輕泄漏問題。2、時域窗函數(shù)的應(yīng)用作為截斷函數(shù)，矩形窗在作時間截斷時，對所截取區(qū)間內(nèi)的信號不加以任何影響，而其它的窗函數(shù)都將對所截取區(qū)間內(nèi)的信號作加權(quán)處理。除了在實驗二中已經(jīng)介紹過三角窗、Hanning窗和Hamming窗外，常用的窗函數(shù)還有很多，例如Par

24、zen窗、Kaiser窗、Chebyshev窗、Tukey窗、Poisson窗、Caushy窗、Gaussian窗和Blackman窗等等。本次實驗仍是采用實驗二中的幾種窗函數(shù)，但是利用窗函數(shù)作時域加權(quán)截斷。 0 t 0 k(a) 正弦函數(shù)的加權(quán)的非周期時域截斷(b)減小了泄漏的頻譜 圖 4-3 采用Hanning窗加權(quán)后的時域截斷和頻譜圖 4-3 中給出了采用Hanning窗對正弦函數(shù)作非整周期的時域加權(quán)截斷后的波形和頻譜，與圖4-1(b)比較，泄漏已明顯減少。應(yīng)該指出，前面所給出的窗函數(shù)都是定義為以0點為中心、寬度為N +1的加權(quán)函數(shù)，在這里應(yīng)用時，需要將其右移，成為區(qū)間內(nèi)的加權(quán)函數(shù)。3、

25、對非周期序列進(jìn)行頻譜分析應(yīng)注意的問題（）混疊一般非周期信號作FFT之前要進(jìn)行時域采樣和周期延拓(無限長時間信號還應(yīng)先截斷再延拓)。根據(jù)Fourier變換理論，經(jīng)等周期的沖激采樣后，離散序列的頻譜是原信號頻譜以為周期的周期延拓。按照Nyquist采樣定理，由采樣引起周期延拓后頻譜之間不發(fā)生混疊的條件是：(1)原信號應(yīng)該是有限帶寬信號，設(shè)其頻帶寬度為fm；(2)頻譜的周期，即采樣周期應(yīng)滿足Nyquist 條件。 0 n 0 N/2 N (a)時域按周期Ts采樣 (b)頻域一個周期內(nèi)在N/2附近出現(xiàn)混疊 圖 4-4 非周期函數(shù)采樣后的幅頻曲線由于實際上有限長時間信號的FT是頻域的無限函數(shù)，因此采樣所

26、得的離散序列的頻譜必定產(chǎn)生混疊，減小采樣周期只能減小而不能消除混疊。對于時間有限函數(shù)，當(dāng)采樣周期較大時，也會在FFT得到的頻域出現(xiàn)混疊，形成頻譜失真，造成頻譜分析結(jié)果與原信號的實際頻譜的差異，也無法恢復(fù)出原信號。當(dāng)然，實際工作中只要采樣和截斷產(chǎn)生的誤差在許可的范圍內(nèi)就行了，但應(yīng)該認(rèn)識到混疊是引起頻譜分析誤差的一個主要原因。還應(yīng)該注意的是，離散Fourier變換的頻域也是周期化的，區(qū)間內(nèi)的樣點實際上是負(fù)頻率區(qū)的量值，因此如果出現(xiàn)混疊，就將在一個周期內(nèi)出現(xiàn)，并發(fā)生在附近的區(qū)域，如圖4-4所示。要減少混疊，就要盡量減小采樣周期。（）泄漏周期函數(shù)截斷引起的頻率泄漏問題，在非周期函數(shù)截斷處理后同樣存在，

27、這種誤差是由于采樣后的離散序列與矩形截斷序列相乘在頻域造成二者的頻譜卷積形成的。矩形窗的頻譜是抽樣函數(shù)序列，它與離散序列頻譜卷積的結(jié)果使原來集中于每一個樣點處的頻譜展寬，其影響在高頻區(qū)(接近N/2的樣點)特別明顯，如圖4-5所示。同樣，為了盡量避免頻率泄漏對結(jié)果的影響，在對非周期函數(shù)作時間截斷時，除盡量增加截斷序列的寬度外，也應(yīng)選取其頻譜的旁瓣較小的截斷函數(shù)，以減輕泄漏問題。 0 n 0 N/2 N (a)時域截斷 (b)頻域一個周期內(nèi)在N/2附近出現(xiàn)泄漏 圖 4-5 函數(shù)采樣后作截斷的幅頻曲線在選取了適當(dāng)?shù)拇昂瘮?shù)后，應(yīng)當(dāng)使窗函數(shù)的寬度與被處理的序列長度相同，如果作變換前還需要補(bǔ)零（例如為了作

28、卷積運算或避免柵欄效應(yīng)），則應(yīng)將原序列與窗函數(shù)相乘后再補(bǔ)零，即補(bǔ)零的樣點不用窗函數(shù)加權(quán)處理。（）柵欄效應(yīng)非周期信號應(yīng)具有連續(xù)的頻譜，在對作抽樣后進(jìn)行DFT，得到的是離散的頻譜。如果排除混疊和泄漏等誤差的影響，所得的結(jié)果也只是的連續(xù)頻譜上的個樣值。這就象通過柵欄上的等間距縫隙觀看到的另一邊的景象，故此稱柵欄效應(yīng)。被柵欄遮住的景象中有可能存在與顯現(xiàn)出的頻譜差異較大的變化，即顯示信號特征的頻譜分量。為了使被柵欄遮住的部分能盡可能地顯現(xiàn)出來，可以采用增加頻域樣點密度的方法，即在不增加信號采樣點的情況下，用時域補(bǔ)零加寬變換尺度N來實現(xiàn)，稱為補(bǔ)零重構(gòu)。例如原來信號采樣得到12個樣點，在其后面再加上4個零，

29、使序列的總長度為16個樣點。這樣處理的結(jié)果原來信號的采樣間隔和頻率都沒有改變，設(shè)采樣頻率為，經(jīng)補(bǔ)零重構(gòu)之后，采樣頻率仍然為，但是原來頻域樣點間寬度為/12，經(jīng)補(bǔ)零重構(gòu)之后頻域樣點間寬度為/16。這就使補(bǔ)零重構(gòu)之后頻域樣點密度增加，而且顯示出原來沒有顯露的一些頻率位置的頻譜。三、實驗內(nèi)容（1）將數(shù)字頻率f=1/53，樣點數(shù)為256的余弦序列作FFT，畫出頻譜曲線，觀察并記錄頻率泄漏現(xiàn)象，然后用Hamming窗和三角窗作加權(quán)截斷，觀察并記錄泄漏的衰減。（2）將幅度為1，周期為2的方波信號，按Ts=1/37 s的間距抽樣，取樣點數(shù)N=256作FFT，畫出頻譜曲線，然后用Hamming窗和三角窗作加權(quán)

30、截斷，觀察并記錄作不同的加權(quán)截斷引起的頻譜差異。（3）將單邊指數(shù)函數(shù)抽樣截斷后作FFT，首先選取不同的抽樣周期s，并取N = 256，觀察頻譜混疊。然后作不同寬度的截斷，選取矩形窗寬為4，8，16，32等，并保持N = 128，觀察頻譜泄漏。說明在什么情況下混疊和泄漏基本消失。（4）計算下面三個正弦函數(shù)的組合的頻譜其中頻率，令，抽樣周期。分別取N=32,64,128將其截斷后作FFT，觀察和記錄混疊和泄漏形態(tài)。分別采取補(bǔ)零加寬和增加截取時間寬度的方法作出頻譜圖，了解柵欄效應(yīng)和頻率分辨力的意義。實驗五 離散系統(tǒng)的Z域分析一、實驗?zāi)康?、掌握離散序列z變換的計算方法。2、掌握離散系統(tǒng)系統(tǒng)函數(shù)零極點

31、的計算方法和零極點圖的繪制方法，并能根據(jù)零極點圖分析系統(tǒng)的因果性和穩(wěn)定性。3、掌握利用Matlab進(jìn)行z反變換的計算方法。二、實驗原理與計算方法1、z變換離散序列x(n)的z變換定義為：。在Matlab中可以利用符號表達(dá)式計算一個因果序列的z變換。其命令格式為：syms n; f=(1/2)n+(1/3)n;ztrans(f)2、離散系統(tǒng)的系統(tǒng)函數(shù)及因果穩(wěn)定的系統(tǒng)應(yīng)滿足的條件一個線性移不變離散系統(tǒng)可以用它的單位抽樣響應(yīng)h(n)來表示其輸入與輸出關(guān)系，即y(n)= x(n)* h(n)對該式兩邊取z變換，得： Y(z)= X(z)· H(z)則： 將H(z)定義為系統(tǒng)函數(shù)，它是單位抽樣

32、響應(yīng)h(n)的z變換，即對于線性移不變系統(tǒng)，若n<0時，h(n)=0，則系統(tǒng)為因果系統(tǒng)；若，則系統(tǒng)穩(wěn)定。由于h(n)為因果序列，所以H(z)的收斂域為收斂圓外部區(qū)域，因此H(z)的收斂域為收斂圓外部區(qū)域時，系統(tǒng)為因果系統(tǒng)。因為，若z=1時H(z)收斂，即，則系統(tǒng)穩(wěn)定，即H(z)的收斂域包括單位圓時，系統(tǒng)穩(wěn)定。因此因果穩(wěn)定系統(tǒng)應(yīng)滿足的條件為：，即系統(tǒng)函數(shù)H(z)的所有極點全部落在z平面的單位圓之內(nèi)。3、Matlab中系統(tǒng)函數(shù)零極點的求法及零極點圖的繪制方法Matlab中系統(tǒng)函數(shù)的零點和極點可以用多項式求根函數(shù)roots ()來實現(xiàn)，調(diào)用該函數(shù)的命令格式為：p=roots(A)。其中A為待求

33、根多項式的系數(shù)構(gòu)成的行向量，返回向量p是包含該多項式所有根位置的列向量。如：求多項式的根的Matlab命令為：A=1 3/4 1/8;p=roots(A)運行結(jié)果為：p=-0.5000-0.2500需要注意的是，離散系統(tǒng)的系統(tǒng)函數(shù)可能有兩種形式，一種是分子和分母多項式均按z的降冪次序排列，另一種是分子分母多項式均按z-1的升冪次序排列，兩種方式在構(gòu)造多項式系數(shù)向量時稍有不同。若H(z)是按z的降冪次序排列，則系數(shù)向量一定要由多項式的最高冪次開始，一直到常數(shù)項，缺項用0補(bǔ)齊，如，其分子多項式的系數(shù)向量應(yīng)為：B=1 0 2 0；分母多項式的系數(shù)向量應(yīng)為：A=1 3 2 2 1。若H(z)是按z-1

34、的升冪次序排列，則分子和分母多項式系數(shù)向量的維數(shù)一定要相同，缺項用0補(bǔ)齊，否則零點和極點就可能被漏掉。如，其分子多項式的系數(shù)向量應(yīng)為：B=1 1 0；分母多項式的系數(shù)向量應(yīng)為：A=1 1/2 1/4?？衫肕atlab中的zplane函數(shù)實現(xiàn)系統(tǒng)函數(shù)的零極點圖的繪制。該函數(shù)的調(diào)用方法為：zplane(B,A);其中B、A為系統(tǒng)函數(shù)分子分母多項式的系數(shù)向量。4、z反變換的計算方法z反變換可由部分分式展開法求得。由于指數(shù)序列anu(n)的z變換為，因此求反變換時，通常對進(jìn)行展開：其中稱為有理函數(shù)的留數(shù)。分兩種情況進(jìn)行討論：（1）X(z)的所有極點均為單實極點此時，則X(z)的z反變換為：（2）X(

35、z)有共軛極點設(shè)X(z)有一對共軛極點，則，其中留數(shù)的計算方法與單極點相同，即，r2=r1 *因此，只要求出部分分式展開的系數(shù)（留數(shù)），就可以直接求出X(z)的z反變換x(n)。在Matlab中可利用函數(shù)residue()求解。令B和A分別是的分子和分母多項式構(gòu)成的系數(shù)向量，則函數(shù)r,p,k=residue (B,A)將產(chǎn)生三個向量r、p、k，其中r為包含部分分式展開系數(shù)ri(i=1,2,N)的列向量，p為包含所有極點的行向量，k為包含部分分式展開的多項式項的系數(shù)cj(j=1,2,M-N)的列向量，若MN，則k為空陣。用residue()函數(shù)求出部分分式展開的系數(shù)后，便可根據(jù)其極點位置分布情況

36、直接求出X(z)的反變換x(n)。如：已知，求其z反變換x(n)。首先利用residue()函數(shù)求出的部分分式展開的系數(shù)和極點，相應(yīng)的Matlab命令為：B=0 1 0;A=1 3 2;r,p,k=residue (B,A)運行結(jié)果為：r = 2 -1p = -2 -1k = 由以上結(jié)果可得：；即X(z)只有兩個單極點，其z反變換為：。已知，求其z反變換x(n)。利用residue()函數(shù)求出的部分分式展開的系數(shù)和極點，可得：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()函數(shù)即可求出共軛極點的模和相位，相應(yīng)命令為：p1=abs(p')p1 = 1.0000 1.0000 1.0000a1=angle(p')/pia1 = 0 -0.3333 0.3333即共軛極點為：，則，其z反變換為：三、實驗內(nèi)容（1）求下列序列的z變換：2-nu(n)；-(1/2)n；(1/2)n+(1/3)n（2）已知兩個離散系統(tǒng)的系統(tǒng)函數(shù)分別為： ；分別求出各系統(tǒng)的零極點，繪制零極點圖，分析系統(tǒng)的穩(wěn)定性；求出各系統(tǒng)單位抽

38、樣響應(yīng)。實驗六 無限沖激響應(yīng)(IIR)數(shù)字濾波器的三種結(jié)構(gòu)一、實驗?zāi)康?、掌握IIR濾波器的直接II型、級聯(lián)型和并聯(lián)型三種結(jié)構(gòu)的基本原理和特點。2、掌握利用Matlab實現(xiàn)IIR濾波器的三種結(jié)構(gòu)的程序設(shè)計方法，并能夠進(jìn)行三者之間的相互轉(zhuǎn)換。3、掌握濾波器頻響特性的繪制方法。二、實驗原理與計算方法按照結(jié)構(gòu)劃分?jǐn)?shù)字濾波器，有遞歸式和非遞歸式兩種。遞歸式數(shù)字濾波器的差分方程為 (1)其中至少有一個.非遞歸式數(shù)字濾波器的差分方程為 (2)可以看出遞歸式數(shù)字濾波器的響應(yīng)不僅與激勵有關(guān)，而且與以前的輸出信號有關(guān)；而非遞歸式數(shù)字濾波器的響應(yīng)僅只與激勵有關(guān)。按照單位樣值響應(yīng)劃分?jǐn)?shù)字濾波器，則有無限沖激響應(yīng)(I

39、IR)和有限沖激響應(yīng)(FIR)之分。IIR濾波器是遞歸式的，差分方程如(1)式所示，F(xiàn)IR濾波器一般是非遞歸式的，差分方程如(2)式所示。IIR濾波器常用的典型結(jié)構(gòu)有直接II型、級聯(lián)型和并聯(lián)型，分別介紹如下：、直接II型（也稱為正準(zhǔn)型結(jié)構(gòu)）根據(jù)(1)式，IIR濾波器的傳輸函數(shù)為 (3)其中已假設(shè)(1)式中的，對于其它情況，則可令相應(yīng)的某些系數(shù)為零。令則有(4)由此可以得到相應(yīng)的時域中激勵與響應(yīng)之間的關(guān)系為 (5)其中是與(4)式中的相應(yīng)的中間函數(shù)序列。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、型結(jié)構(gòu)信號流由(5)式確定的直接II型的信號流圖如圖6-1所示，其中將中間的兩條延時鏈合并為一條，實際的信號流將按（5）式分成兩個延時鏈獨立運行。編程時，用三個數(shù)組分別存放系數(shù)和，。由圖可以看出，沿中間的延時鏈自上向下傳播過程中將逐級向右移位，每一級向左邊與的乘積按累加，再和相加，得到的沿中間的延時鏈又自上向下傳播逐級向右移位，每一級向右邊與的乘積按累加，其結(jié)果就是對應(yīng)于激勵的響應(yīng)。然后向右移位一個單位時間，輸入激勵計算下一個時刻的響應(yīng)。直接II型結(jié)構(gòu)具有簡單直觀的典型網(wǎng)絡(luò)結(jié)構(gòu)形式，在計算機(jī)上很容易實現(xiàn)。但是它對系數(shù)的量值變化比較敏感，直接確定了系統(tǒng)零、極點的位置，從而影響到系統(tǒng)的性能。尤其當(dāng)

41、階數(shù)N較高時，系統(tǒng)對系數(shù)的字長效應(yīng)很敏感，產(chǎn)生的誤差也較大。、級聯(lián)型結(jié)構(gòu)由于當(dāng)直接II型結(jié)構(gòu)傳輸函數(shù)階數(shù)增加時，系數(shù)量化引起的誤差影響到濾波器的性能，因此要采用其它形式的結(jié)構(gòu)。將(3)式的傳輸函數(shù)分子和分母進(jìn)行因式分解，即用它的零、極點表示為 (6)其中，分子和分母中的實系數(shù)二階因子分別對應(yīng)于共軛零、極點?？梢詫⑸鲜椒肿雍头帜钢袉胃浑A因子作為二階因子的一種特例，那么(6)式就可以表示為M個實系數(shù)二階基本節(jié)級聯(lián)的形式 (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 級聯(lián)型結(jié)構(gòu)信號流其中為濾波器的二階基本節(jié)，為取整結(jié)果。則分別是第k個基本節(jié)的分母、分子系數(shù)，它們只關(guān)系到濾波器的某一對零、極點。級聯(lián)型結(jié)構(gòu)如圖6-2所示，它是一系列二階基本節(jié)的級聯(lián)，每一個可以用直接II型結(jié)構(gòu)實現(xiàn)。級聯(lián)型結(jié)構(gòu)的特點是對濾波器性能的調(diào)整比較方便，調(diào)整系數(shù)，只單獨涉及到第k級零、極點，而不會影響到其它任一級的零、極點，因而可以獨立地控制濾波器的各零、極點的分布。、并聯(lián)型結(jié)構(gòu)這種結(jié)構(gòu)將傳輸函數(shù)展開為部分分式，即表示為若干一階和二階基本節(jié)網(wǎng)絡(luò)與一個常數(shù)之和 (9)其中，同樣也可以統(tǒng)一表示為二階基本節(jié)的形式(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 并聯(lián)型結(jié)構(gòu)信號流并聯(lián)型結(jié)構(gòu)信號流如圖6-3所示，其中二階基本節(jié)網(wǎng)絡(luò)可以用直接II型結(jié)構(gòu)實現(xiàn)，程序設(shè)計也可參考直接型II結(jié)構(gòu)的方法。并聯(lián)型結(jié)構(gòu)也可以單獨調(diào)整極點位置，但卻不能象級聯(lián)型結(jié)構(gòu)那樣直接控制零點的分布。因為并聯(lián)型結(jié)構(gòu)各二階基本節(jié)網(wǎng)絡(luò)的零點并不是整個系統(tǒng)函數(shù)的零點。因此，當(dāng)要準(zhǔn)確傳輸零點時，以采用級聯(lián)型結(jié)構(gòu)為宜。不過，由于并聯(lián)型基本節(jié)之間互不影響，所以運算誤差比級聯(lián)型的要小一些。通常，IIR濾波器用系統(tǒng)函數(shù)的有理式（直接結(jié)構(gòu)形式）描述，Matalb中提供相應(yīng)的函數(shù)可以把直接型結(jié)構(gòu)轉(zhuǎn)換成級聯(lián)和并聯(lián)型結(jié)構(gòu)。在Ma

44、tlab中，直接型結(jié)構(gòu)由兩個行向量描述，b包含系數(shù)bk，a包含系數(shù)ak，可由filter函數(shù)實現(xiàn)。要由已知的直接型結(jié)構(gòu)實現(xiàn)級聯(lián)型結(jié)構(gòu)，可利用函數(shù)dir2cas根據(jù)直接型的系數(shù)bk和ak得到系數(shù)b0，Bk,i和Ak,i，該函數(shù)把矢量b和a轉(zhuǎn)換成K×3維矩陣B和A。首先計算b0，它等于b0/ a0，（a01），然后通過給較短的矢量添零，使矢量b和a一樣長。這就保證每一個雙二階環(huán)節(jié)的分子和分母均不為零。接著計算多項式B(z)和A(z)的根，用cplxpair函數(shù)把這些根以共軛復(fù)根對的次序排列，最后用poly函數(shù)把每一對根再轉(zhuǎn)換成二階分子或分母多項式。級聯(lián)型由casfiltr函數(shù)實現(xiàn)，它把每

45、個二階環(huán)節(jié)的系數(shù)存放在矩陣B和A中，放在一個循環(huán)里，而在此循環(huán)中采用filter函數(shù)。輸入被乘以b0，每個濾波器的輸出作為下一級濾波器的輸入，最后一個濾波器的輸出即為總的輸出。函數(shù)cas2dir可以把級聯(lián)形式轉(zhuǎn)換成直接形式。這是一個包含幾次多項式乘法的簡單運算，在一個K次的循環(huán)中實用conv函數(shù)。類似的，利用dir2par函數(shù)可以把直接型系數(shù)bk和ak轉(zhuǎn)換成并聯(lián)型系數(shù)Bk,i和Ak,i，并在parfiltr函數(shù)中使用這些系數(shù)以實現(xiàn)并聯(lián)形式。而par2dir函數(shù)則可用來實現(xiàn)并聯(lián)形式轉(zhuǎn)換成直接形式。三種結(jié)構(gòu)的IIR濾波器均可通過以下方式計算單位抽樣響應(yīng)。delta=impseq(0,0,7);hd

46、ir=filter(b,a,delta)hcas=casfiltr(b0,B,A,delta)hpar=parfiltr(C,B,A,delta)傳輸函數(shù)為的濾波器的頻率特性為。由于我們所討論的濾波器都是穩(wěn)定系統(tǒng)，其所有極點都在z平面的單位圓內(nèi)，單位樣值響應(yīng)是一個無限長的衰減序列，所以要截取有限長的一段來作頻率特性分析。截取的長度越長，逼近的程度越高。如果截取長度為N，濾波器的頻率特性由Fourier變換得 (11)則濾波器的幅頻特性和相頻特性可表示為 (12) (13)三、實驗內(nèi)容設(shè)三階濾波器的傳輸函數(shù)為激勵信號為根據(jù)所給定的濾波器系統(tǒng)函數(shù)和參數(shù)值，用三種結(jié)構(gòu)實現(xiàn)該IIR數(shù)字濾波器結(jié)構(gòu)，分別

47、使激勵和通過該濾波器，求出相應(yīng)于不同激勵的響應(yīng)，作出輸出響應(yīng)信號的時間曲線和幅頻、相頻特性曲線，并判斷為何種濾波器(低通、高通、帶通)。實驗七 沖激響應(yīng)不變法IIR數(shù)字濾波器設(shè)計一、實驗?zāi)康?、掌握構(gòu)成一個頻率響應(yīng)與給定的濾波特性相接近的模擬濾波器的設(shè)計原理。2、掌握用沖激響應(yīng)不變法設(shè)計IIR數(shù)字濾波器的基本原理和算法。3、了解數(shù)字濾波器和模擬濾波器的頻率響應(yīng)特性，掌握相應(yīng)的計算方法，分析用沖激響應(yīng)不變法獲得的數(shù)字濾波器頻率響應(yīng)特性中出現(xiàn)的混疊現(xiàn)象。二、實驗原理與計算方法、沖激響應(yīng)不變法設(shè)計IIR數(shù)字濾波器的基本原理和算法采用沖激響應(yīng)不變法設(shè)計數(shù)字濾波器，就是使其單位樣值響應(yīng)與相應(yīng)的模擬濾波器的沖激響應(yīng)在抽樣點處的量值相等，即 (1)其中T為抽樣周期。因此用沖激響應(yīng)不變法設(shè)計IIR數(shù)字濾波器的基本步驟，就是首先根據(jù)設(shè)計要求確定相應(yīng)的模擬濾波器的傳遞函數(shù)，經(jīng)Laplace反變換求出沖激響應(yīng)后，對它進(jìn)行抽樣得到的等于數(shù)字濾波器的單位樣值響應(yīng)，再經(jīng)z變換所得就是數(shù)字濾波器的傳遞函數(shù)。如果模擬濾波器的傳遞函數(shù)的N個極點都是單極點，則可以將寫成部分分式展開的形式 (2)那么，經(jīng)Laplace反變換求出的模擬濾波器的沖激響應(yīng)為相對應(yīng)的數(shù)字濾波器的單位樣值響應(yīng)為對上式作z變換，得 (3)由上面的推導(dǎo)可見，只要模擬濾波器的傳遞函數(shù)的N個極點都是單極點，當(dāng)已經(jīng)求出