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個極點都是單極點，當已經求出