《數(shù)字信號處理》上機(jī)實(shí)驗(yàn)指導(dǎo)書

《數(shù)字信號處理》上機(jī)實(shí)驗(yàn)指導(dǎo)書 -PAGE1-《數(shù)字信號處理》上機(jī)實(shí)驗(yàn)指導(dǎo)書一、引言“數(shù)字信號處理”是一門理論和實(shí)驗(yàn)密切結(jié)合的課程，為了深入地掌握課程內(nèi)容，應(yīng)當(dāng)在學(xué)習(xí)理論的同時，做習(xí)題和上機(jī)實(shí)驗(yàn)。上機(jī)實(shí)驗(yàn)不僅可以幫助學(xué)生深入地理解和消化基本理論，而且能鍛煉初學(xué)者的獨(dú)立解決問題的能力。所以，根據(jù)本課程的重點(diǎn)要求編寫了四個實(shí)驗(yàn)。第一章、二章是全書的基礎(chǔ)內(nèi)容，抽樣定理、時域離散系統(tǒng)的時域和頻域分析以及系統(tǒng)對輸入信號的響應(yīng)是重要的基本內(nèi)容。由于第一、二章大部分內(nèi)容已經(jīng)在前期《信號與系統(tǒng)》課程中學(xué)習(xí)完，所以可通過實(shí)驗(yàn)一幫助學(xué)生溫習(xí)以上重要內(nèi)容，加深學(xué)生對“數(shù)字信號處理是通過對輸入信號的一種運(yùn)算達(dá)到處理目的”這一重要概念的理解。這樣便可以使學(xué)生從《信號與系統(tǒng)》課程順利的過渡到本課程的學(xué)習(xí)上來。第三章、四章DFT、FFT是數(shù)字信號處理的重要數(shù)學(xué)工具，它有廣泛的使用內(nèi)容。限于實(shí)驗(yàn)課時，僅采用實(shí)驗(yàn)二“用FFT對信號進(jìn)行譜分析”這一實(shí)驗(yàn)。通過該實(shí)驗(yàn)加深理解DFT的基本概念、基本性質(zhì)。FFT是它的快速算法，必須學(xué)會使用。所以，學(xué)習(xí)完第三、四章后，可安排進(jìn)行實(shí)驗(yàn)二。數(shù)字濾波器的基本理論和設(shè)計(jì)方法是數(shù)字信號處理技術(shù)的重要內(nèi)容。學(xué)習(xí)這一部分時，應(yīng)重點(diǎn)掌握IIR和FIR兩種不同的數(shù)字濾波器的基本設(shè)計(jì)方法。IIR濾波器的單位沖激響應(yīng)是無限長的，設(shè)計(jì)方法是先設(shè)計(jì)模擬濾波器，然后再通過S~Z平面轉(zhuǎn)換，求出相應(yīng)的數(shù)字濾波器的系統(tǒng)函數(shù)。這里的平面轉(zhuǎn)換有兩種方法，即沖激響應(yīng)不變法和雙線性變換法，后者沒有頻率混疊的缺點(diǎn)，且轉(zhuǎn)換簡單，是一種普遍應(yīng)用的方法。學(xué)習(xí)完第六章以后可以進(jìn)行實(shí)驗(yàn)三。FIR濾波器的單位沖激響應(yīng)是有限長的，設(shè)計(jì)濾波器的目的即是求出符合要求的單位沖激響應(yīng)。窗函數(shù)法是一種基本的，也是一種重要的設(shè)計(jì)方法。學(xué)習(xí)完第七章后可以進(jìn)行實(shí)驗(yàn)四。以上所提到的四個實(shí)驗(yàn)，可根據(jù)實(shí)驗(yàn)課時的多少恰當(dāng)安排。例如：實(shí)驗(yàn)一可根據(jù)學(xué)生在學(xué)習(xí)《信號與系統(tǒng)》課程后，掌握的程度來確定是否做此實(shí)驗(yàn)。若時間緊，可以在實(shí)驗(yàn)三、四之中任做一個實(shí)驗(yàn)。建議實(shí)驗(yàn)二必做，因?yàn)镈FT和FFT是學(xué)習(xí)數(shù)字信號處理必須掌握的內(nèi)容。二、關(guān)于使用計(jì)算機(jī)語言由于數(shù)字信號處理實(shí)驗(yàn)的主要目的是驗(yàn)證數(shù)字信號處理的有關(guān)理論，進(jìn)一步理解鞏固所學(xué)理論知識。所以，實(shí)現(xiàn)實(shí)驗(yàn)用的算法語言可以有許多種，但為了提高實(shí)驗(yàn)效率，要求學(xué)生用編程效率比C語言高好幾倍的MATLAB語言。下面介紹MATLAB的主要特點(diǎn)。（有關(guān)MATLAB的啟動、程序運(yùn)行和有關(guān)信號處理工具箱函數(shù)等內(nèi)容將放到最后附錄中介紹。）MATLAB是一種交互式的以矩陣為基本數(shù)據(jù)結(jié)構(gòu)的系統(tǒng)。在生成矩陣對象時，不要求明確的維數(shù)說明。所謂交互式，是指MATLAB的草稿紙編程環(huán)境。即用戶每輸入一條命令并按回車鍵，MATLAB系統(tǒng)便解釋執(zhí)行之，并顯示執(zhí)行結(jié)果。根據(jù)該結(jié)果，用戶立即知道剛輸入的命令的正確性，或利用中間結(jié)果進(jìn)行其他處理等。與C語言或FORTRON語言做科學(xué)數(shù)值計(jì)算的程序設(shè)計(jì)相比較，利用MATLAB可節(jié)省大量的編程時間。將其用于數(shù)字信號處理實(shí)驗(yàn)，則可大大提高實(shí)驗(yàn)效率，在有限的上機(jī)時間內(nèi)，實(shí)驗(yàn)內(nèi)容可增加幾倍。例如，C語言FFT子程序有70多行，而用MATLAB只調(diào)用一個fft函數(shù)即可實(shí)現(xiàn)對序列進(jìn)行FFT計(jì)算。另外，MATLAB的工具箱及圖形顯示（打印）功能，可滿足各層次人員直觀、方便的進(jìn)行分析、計(jì)算和設(shè)計(jì)工作，從而可大大節(jié)省時間。例如，序列的卷積、濾波，系統(tǒng)函數(shù)H(z)的幅頻特性和相頻特性等計(jì)算，均有現(xiàn)成的工具箱函數(shù)。而用其它算法語言完成這些計(jì)算的編程比較麻煩，且程序較長。由于上述特點(diǎn)，在美國一些大學(xué)里，MATLAB已成為輔助教學(xué)的有益工具。MATLAB已成功地用于數(shù)字信號處理課程中的問題分析、實(shí)驗(yàn)、濾波器設(shè)計(jì)及計(jì)算機(jī)模擬。附錄中所介紹的信號處理工具箱函數(shù)及繪圖函數(shù)基本可滿足本教材所要求的上機(jī)實(shí)驗(yàn)需要。對序列進(jìn)行譜分析的MATLAB程序及運(yùn)行結(jié)果見附錄。三、上機(jī)實(shí)驗(yàn)實(shí)驗(yàn)一：信號、系統(tǒng)及系統(tǒng)響應(yīng)1.實(shí)驗(yàn)?zāi)康蘑偈煜みB續(xù)信號經(jīng)過理想抽樣前后的頻譜變化關(guān)系，加深對時域抽樣定理的理解。②熟悉時域離散系統(tǒng)的時域特性。③利用卷積方法觀察分析系統(tǒng)的時域特性。④掌握序列傅里葉變換的計(jì)算機(jī)實(shí)驗(yàn)方法，利用序列的傅里葉變換對連續(xù)信號、離散信號及系統(tǒng)響應(yīng)進(jìn)行頻域分析。2.實(shí)驗(yàn)原理與方法抽樣是連續(xù)信號數(shù)字處理的第一個關(guān)鍵環(huán)節(jié)。對抽樣過程的研究不僅可以了解抽樣前后信號時域和頻域特性發(fā)生的變化以及信號信息不丟失的條件，而且可以加深對傅里葉變換、Z變換和序列傅里葉變換之間關(guān)系式的理解。我們知道，對一個連續(xù)信號xa(t)進(jìn)行理想抽樣的過程可用（1.1）式表示。=δT(t)——（1.1）其中為xa(t)的理想抽樣，δT(t)為周期沖激脈沖，即——（1.2）的傅里葉變換為=——（1.3）（1.3）式表明為的周期延拓，其延拓周期為抽樣角頻率(Ωs=2π/T)。抽樣前后信號的頻譜示意圖見“參考教材數(shù)字信號處理教程（程佩青）圖1-29”。只有滿足抽樣定理時，才不會發(fā)生頻率混疊失真。數(shù)字信號處理教程（程佩青）在計(jì)算機(jī)上用高級語言編程直接按（1.3）式計(jì)算理想抽樣的頻譜很不方便。下面導(dǎo)出用序列的傅里葉變換來計(jì)算的公式。將（1.2）式代入（1.1）式并進(jìn)行傅里葉變換，==

==——（1.4）式中的xa(nT)就是采樣后得到的序列x(n)，即x(n)=xa(nT)x(n)的序列傅里葉變換為X(ejω)=——（1.5）比較（1.5）和（1.4）可知=X(ejω)|ω=ΩT——（1.6）這說明兩者之間只在頻率度量上差一個常數(shù)因子T。實(shí)驗(yàn)過程中應(yīng)注意這一差別。離散信號和系統(tǒng)在時域均可用序列來表示。序列圖形給人以形象直觀的印象，它可加深我們對信號和系統(tǒng)的時域特征的理解。本實(shí)驗(yàn)還將觀察分析幾種信號及系統(tǒng)的時域特性。為了在數(shù)字計(jì)算機(jī)上觀察分析各種序列的頻域特性，通常對X(ejω)在[0，2π]上進(jìn)行M點(diǎn)采樣來觀察分析。對長度為N的有限長序列x(n)，有——（1.7）其中,k=0,1,…,M-1通常M應(yīng)取得大一些，以便觀察譜的細(xì)節(jié)變化。取模||可繪出幅頻待性曲線。一個時域離散線性非移變系統(tǒng)的輸入/輸出關(guān)系為y(n)=x(n)*h(n)=——（1.8）這里，y(n)為系統(tǒng)的輸出序列，x(n)為輸入序列。h(n)、x(n)可以是無限長，也可以是有限長。為了計(jì)算機(jī)繪圖觀察方便，主要討論有限長情況。如果h(n)和x(n)的長度分別為N和M，則y(n)的長度為L=N+M-1。這樣，(1.8)式所描述的卷積運(yùn)算就是序列移位、相乘和累加的過程，所以編程十分簡單。上述卷積運(yùn)算也可以在頻域?qū)崿F(xiàn)(即卷積定理：時域卷積，頻域相乘。)Y(ejω)=X(ejω)H(ejω)——（1.9）(1.9)式右邊的相乘是在各頻點(diǎn){ωk}上的頻譜值相乘。3.實(shí)驗(yàn)內(nèi)容(1)連續(xù)信號分析①連續(xù)時間非周期信號的選擇1單邊指數(shù)脈沖；2雙邊指數(shù)脈沖；3鐘形脈沖；4抽樣脈沖；1.xa(t)=E*exp(-at)(a>0);2.xa(t)=E*exp(-a*abs(t))(a>0)

②用MATLAB編制程序開始結(jié)束開始結(jié)束產(chǎn)生信號xa(t)求連續(xù)時間非周期信號的傅氏變換Xa(jΩ)分別繪制xa(t)，|Xa(jf)|圖形圖1.1連續(xù)信號分析程序框圖③時域觀察，頻域分析調(diào)整信號xa(t)參數(shù)，觀察時域波形變化。分析頻域波形，選取近似帶限頻率值fmax，確定對xa(t)→x(n)=xa(nT)的抽樣周期T(T=1/fs)，不失真的條件應(yīng)該滿足抽樣頻率fs≥2fmax（即：時域抽樣定理）。例如圖1.2給出了單邊減幅余弦信號xa(t)=exp(-10*t).*cos(2*pi*25*t)的時域和幅頻特性曲線，由此圖可以近似選取xa(t)的抽樣頻率fs≥2fmax=2×150=300Hz。圖1.2xa(t)的時域和幅頻特性曲線(2)離散信號分析①離散時間非周期信號x(n)的生成對前一個實(shí)驗(yàn)中的連續(xù)信號xa(t)抽樣，用于產(chǎn)生實(shí)驗(yàn)中要用到的信號序列x(n)=xa(nT)。根據(jù)上一個實(shí)驗(yàn)分析結(jié)果選取近似帶限頻率值fmax，按照抽樣頻率fs=2fmax、fs>2fmax和fs<2fmax三種情況編制MATLAB源程序(T=1/fs)。

②用MATLAB編制程序開始結(jié)束對xa(t)抽樣產(chǎn)生序列x(n)求序列x(n)的傅氏變換X(e開始結(jié)束對xa(t)抽樣產(chǎn)生序列x(n)求序列x(n)的傅氏變換X(ejω)分別繪制x(n)，|X(ej2πfT)|圖形確定抽樣頻率fs,產(chǎn)生離散時間變量t=nT圖1.3離散信號分析程序框圖③時域觀察，頻域分析調(diào)整抽樣頻率fs=2fmax、fs>2fmax和fs<2fmax，觀察時域波形變化，分析頻域波形。根據(jù)時域、頻域的變化驗(yàn)證時域抽樣定理X(ej2πfT)=。改變程序中頻率軸變量為：f(Hz)模擬頻率、Ω(rad/s)模擬角頻率和ω(rad)數(shù)字頻率，觀察頻率特性的周期變化，并用所學(xué)理論對以上各種情況加以解釋。例如圖1.4給出了單邊減幅余弦抽樣信號x(n)的時域和幅頻特性曲線，由此圖可以觀察分析選取抽樣頻率fs=2fmax=2×150=300Hz時的時域和幅頻特性。圖1.4xa(t)的時域和幅頻特性曲線

(3)系統(tǒng)響應(yīng)分析①生成實(shí)驗(yàn)用的輸入序列x(n)和系統(tǒng)單位沖激響應(yīng)序列h(n)

輸入序列：x(n)=R10(n)

單位沖激響應(yīng)序列：h(n)=δ(n)+2.5δ(n-1)+2.5δ(n-2)+δ(n-3)②時域離散信號、系統(tǒng)和系統(tǒng)響應(yīng)分析

觀察系統(tǒng)h(n)對輸入信號x(n)=R10(n)的響應(yīng)特性y(n)。即利用線性卷積求系統(tǒng)響應(yīng)y(n)，調(diào)用有限長序列線性卷積子程序，用于完成兩個給定長度的序列的卷積?？梢灾苯诱{(diào)用MATLAB語言中的卷積函數(shù)conv。conv用于兩個有限長度序列的卷積，它假定兩個序列都從n=0開始。

調(diào)用格式如下：y=conv(x,h)

改變輸入信號x(n)長度，打印輸出y(n)圖形并判斷y(n)圖形及其非零值序列長度是否與理論計(jì)算結(jié)果一致。③卷積定理的驗(yàn)證

a.調(diào)用序列博里葉變換數(shù)值計(jì)算子程序，求輸出y(n)的頻率響應(yīng)Y(ejωk),觀察|Y(ejωk)|特性曲線。

b.調(diào)用序列博里葉變換數(shù)值計(jì)算子程序，分別求輸入序列x(n)的頻率響應(yīng)X(ejωk)和離散系統(tǒng)的沖激響應(yīng)序列h(n)的頻率響應(yīng)H(ejωk)，計(jì)算Y(ejωk)=X(ejωk)×H(ejωk)，并繪出|Y(ejωk)|曲線。

c.將a.與b.分別計(jì)算的幅頻特性曲線|Y(ejωk)|進(jìn)行比較，驗(yàn)證時域卷積定理。

圖1.5給出系統(tǒng)響應(yīng)分析實(shí)驗(yàn)程序流程框圖，供實(shí)驗(yàn)者參考。開始結(jié)束產(chǎn)生輸入信號序列x(n)和單位沖激響應(yīng)序列h(n)用序列傅氏變換數(shù)值計(jì)算，求X(ejω開始結(jié)束產(chǎn)生輸入信號序列x(n)和單位沖激響應(yīng)序列h(n)用序列傅氏變換數(shù)值計(jì)算，求X(ejωk)并繪制其圖形分別繪制y(n)，|Y(ejωk)|圖形用序列傅氏變換數(shù)值計(jì)算，求H(ejωk)并繪制其圖形求Y(ejωk)=X(ejωk)×H(ejωk)計(jì)算離散系統(tǒng)輸出信號序列y(n)=x(n)*h(n)用序列傅氏變換數(shù)值計(jì)算，求Y(ejωk)繪制y(n)的Y(ejωk)圖形與Y(ejωk)=X(ejωk)×H(ejωk)的圖形圖1.5系統(tǒng)響應(yīng)分析實(shí)驗(yàn)的程序框圖4.思考題在分析理想抽樣序列特性的實(shí)驗(yàn)中，采樣頻率不同時，相應(yīng)理想抽樣序列的傅里葉變換頻譜的數(shù)字頻率度量是否都相同？它們所對應(yīng)的模擬頻率是否相同？為什么？5.實(shí)驗(yàn)報(bào)告要求(1)簡述實(shí)驗(yàn)?zāi)康募皩?shí)驗(yàn)原理。(2)按實(shí)驗(yàn)要求附上實(shí)驗(yàn)過程中編制的MATLAB源程序及實(shí)驗(yàn)中涉及的所有信號序列、系統(tǒng)單位沖激響應(yīng)及系統(tǒng)響應(yīng)序列的時域和幅頻特性曲線，并對所得結(jié)果進(jìn)行分析和解釋。(3)結(jié)合所學(xué)理論知識總結(jié)實(shí)驗(yàn)中的主要結(jié)論。(4)簡要回答思考題。實(shí)驗(yàn)二：用FFT做譜分析1.實(shí)驗(yàn)?zāi)康?1)進(jìn)一步加深DFT算法原理和基本性質(zhì)的理解（因?yàn)镕FT只是DFT的一種快速算法，所以FFT的運(yùn)算結(jié)果必然滿足DFT的基本性質(zhì)）。(2)熟悉FFT算法原理和FFT子程序的應(yīng)用。(3)學(xué)習(xí)用FFT對連續(xù)信號和時域離散信號進(jìn)行譜分析的方法，了解可能出現(xiàn)的分析誤差及其原因，以便在實(shí)際中正確應(yīng)用FFT。2.實(shí)驗(yàn)步驟(1)復(fù)習(xí)DFT的定義、性質(zhì)和用DFT作譜分析的有關(guān)內(nèi)容。(2)復(fù)習(xí)按時間抽選法FFT算法原理及相應(yīng)的運(yùn)算流圖(3)編制信號產(chǎn)生子程序，產(chǎn)生以下典型信號供譜分析用：x1(n)=R4(n)x1=ones(1,4)=[1,1,1,1]x2(n)=x2=[1,2,3,4,4,3,2,1…0]x3(n)=x4(n)=cos(π/4)nx5(n)=sin(π/8)nx6(t)=cos8πt+cos16πt+cos20πt應(yīng)當(dāng)注意，如果給出的是連續(xù)信號xa(t)，則首先要根據(jù)其最高頻率確定抽樣頻率fs以及由頻率分辨率選擇抽樣點(diǎn)數(shù)N，然后對其進(jìn)行軟件抽樣(即計(jì)算x(n)=xa(nT)，

0≤n≤N-1)，產(chǎn)生對應(yīng)序列x(n)。對信x6(t)，頻率分辨率的選擇要以能分辨開其中的三個頻率對應(yīng)的譜線為準(zhǔn)則。對周期序列，最好截取周期的整數(shù)倍進(jìn)行譜分析，否則有可能產(chǎn)生較大的分析誤差。請實(shí)驗(yàn)者根據(jù)DFT的隱含周期性思考這個問題。

開始結(jié)束開始結(jié)束讀入長度N調(diào)用信號產(chǎn)生子程序產(chǎn)生實(shí)驗(yàn)信號調(diào)用繪圖子程序(函數(shù))繪制時間序列波形圖調(diào)用FFT子程序(函數(shù))計(jì)算信號的DFT調(diào)用繪圖子程序(函數(shù))繪制|X(k)|曲線

圖2.1主程序框圖(5)按實(shí)驗(yàn)內(nèi)容要求，上機(jī)實(shí)驗(yàn)，并寫出實(shí)驗(yàn)報(bào)告。3．上機(jī)實(shí)驗(yàn)內(nèi)容(1)對2中所給出的信號逐個進(jìn)行譜分析。下面給出針對各信號的FFT變換區(qū)間N以及對連續(xù)信號x6(t)的抽樣頻率fs，供實(shí)驗(yàn)時參考。x1(n),x2(n),x3(n),x4(n),x5(n)：N=8,16x6(t)：fs=64(Hz),N=16,32,64(2)令x(n)=x4(n)+x5(n)，用FFT計(jì)算8點(diǎn)和16點(diǎn)離散傅里葉變換，X(k)=DFT[x(n)]并根據(jù)DFT的對稱性，由X(k)求出X4(k)=DFT[x4(n)]和X5(k)=DFT[x5(n)]，并與(1)中所得結(jié)果比較。［提示：取N=16時，x4(n)=x4(N-n),x5(n)=-x5(N-n)。］(3)令x(n)=x4(n)+jx5(n)，重復(fù)(2)。4.思考題(1)在N=8時，x2(n)和x3(n)的幅頻特性會相同嗎？為什么？N=16呢？(2)如果周期信號的周期預(yù)先不知道，如何用FFT進(jìn)行譜分析？5.實(shí)驗(yàn)報(bào)告要求(1)簡述實(shí)驗(yàn)原理及目的。(2)結(jié)合實(shí)驗(yàn)中所得給定典型序列幅頻特性曲線，與理論結(jié)果比較，并分析說明誤差產(chǎn)生的原因以及用FFT作譜分析時有關(guān)參數(shù)的選擇方法。(3)總結(jié)實(shí)驗(yàn)所得主要結(jié)論。(4)簡要回答思考題。實(shí)驗(yàn)三：用雙線性變換法設(shè)計(jì)IIR數(shù)字濾波器1.實(shí)驗(yàn)?zāi)康?1)熟悉用雙線性變換法設(shè)計(jì)IIR數(shù)字濾波器的原理與方法。(2)掌握數(shù)字濾波器的計(jì)算機(jī)仿真方法。(3)通過觀察對實(shí)際心電圖信號的濾波作用，獲得數(shù)字濾波的感性知識。2.實(shí)驗(yàn)內(nèi)容(1)用雙線性變換法設(shè)計(jì)一個巴特沃斯低通IIR數(shù)字濾波器。設(shè)計(jì)指標(biāo)參數(shù)為：在通帶內(nèi)頻率低于0.2π時，最大衰減小于1dB；在阻帶內(nèi)[0.3π,π]頻率區(qū)間上，最小衰減大于15dB。(2)以0.02π為抽樣間隔，打印出數(shù)字濾波器在頻率區(qū)間[0,π/2]上的幅頻響應(yīng)特性曲線。(3)用所設(shè)計(jì)的濾波器對實(shí)際心電圖信號抽樣序列（在本實(shí)驗(yàn)后面給出）進(jìn)行仿真濾波處理，并分別打印出濾波前后的心電圖信號波形圖，觀察總結(jié)濾波作用與效果。3.實(shí)驗(yàn)步驟(1)復(fù)習(xí)有關(guān)巴特沃斯模擬濾波器設(shè)計(jì)和用雙線性變換法設(shè)計(jì)IIR數(shù)字濾波器的內(nèi)容，按照“2.實(shí)驗(yàn)內(nèi)容(1)”的要求設(shè)計(jì)滿足指標(biāo)的數(shù)字濾波器函數(shù)H(z)?，F(xiàn)給出滿足本實(shí)驗(yàn)要求的數(shù)字濾波器系統(tǒng)函數(shù)：H(z)=

=——（3.1）

式中Hk(z)=,k=1,2,3——（3.2）A=0.09036

B1=1.2686,C1=－0.7051

B2=1.0106,C2=－0.3583

B3=0.9044,C3=－0.2155由(3.1)式和(3.2)式可見，濾波器H(z)由三個二階濾波器H1(z)、H2(z)和H3(z)級聯(lián)組成，如圖3.1所示。H1H1(z)H2(z)H3(z)x(n)y1(n)y2(n)y3(n)=y(n)H(z)圖3.1濾波器H(z)的組成(2)編寫濾波器仿真程序，計(jì)算H(z)對心電圖信號抽樣序列x(n)的響應(yīng)序列y(n)。設(shè)yk(n)為第k級二階濾波器Hk(z)的輸出序列，yk-1(n)為輸入序列，如圖3.1所示。由(3.2)式可得到差分方程：yk(n)=Ayk-1(n)+2Ayk-1(n-1)+Ayk-1(n-2)+Bkyk(n-1)+Ckyk(n-2)——（3.3）當(dāng)k=1時，yk-1(n)=x(n)。所以H(z)對x(n)的總響應(yīng)序列y(n)可以用順序迭代算法得到。即依次對k=1，2，3，求解差分方程(3.3)，最后得到y(tǒng)3(n)=y(n)。仿真程序就是實(shí)現(xiàn)上述求解差分方程和順序迭代算法的通用程序。也可以直接調(diào)用MATLABfilter函數(shù)實(shí)現(xiàn)仿真。(3)在通用計(jì)算機(jī)上運(yùn)行仿真濾波程序，并調(diào)用通用繪圖子程序，完成“2.實(shí)驗(yàn)內(nèi)容(2)和(3)”。本實(shí)驗(yàn)要用的MATLAB繪圖函數(shù)見附錄。4.思考題用雙線性變換法設(shè)計(jì)數(shù)字濾波器過程中，變換公式s=中T的取值，對設(shè)計(jì)結(jié)果有無影響？為什么？5.實(shí)驗(yàn)報(bào)告要求(1)簡述實(shí)驗(yàn)?zāi)康募霸怼?2)由所打印的|H(ejω)|特性曲線及設(shè)計(jì)過程簡述雙線性變換法的特點(diǎn)。(3)對比濾波前后的心電圖信號波形，說明數(shù)字濾波器的濾波過程與濾波作用。(4)簡要回答思考題。6.心電圖信號抽樣序列x(n)人體心電圖信號在測量過程中往往受到工業(yè)高頻干擾，所以必須經(jīng)過低通濾波處理后，才能作為判斷心臟功能的有用信息。下面給出一實(shí)際心電圖信號抽樣序列樣本x(n)，其中存在高頻干擾。在實(shí)驗(yàn)中，以x(n)作為輸入序列，濾除其中的干擾成分。｛x(n)｝=｛-4，-2，0，-4，-6，-4，-2，-4，-6，-6，

-4，-4，-6，-6，-2，6，12，8，0，-16，

-38，-60，-84，-90，-66，-32，-4，-2，-4，8，

12，12，10，6，6，6，4，0，0，0，

0，0，-2，-4，0，0，0，-2，-2，0，

0，-2，-2，-2，-2，0，｝實(shí)驗(yàn)四：用窗函數(shù)法設(shè)計(jì)FIR數(shù)字濾波器1.實(shí)驗(yàn)?zāi)康?1)掌握用窗函數(shù)法設(shè)計(jì)FIR數(shù)字濾波器的原理和方法。(2)熟悉線性相位FIR數(shù)字濾波器特性。(3)了解各種窗函數(shù)對濾波特性的影響。2.實(shí)驗(yàn)原理與方法如果所希望的濾波器的理想頻率響應(yīng)函數(shù)為Hd(ejω)，則其對應(yīng)的單位脈沖響應(yīng)為hd(n)=——（4.1）窗函數(shù)設(shè)計(jì)法的基本原理是用有限長單位脈沖響應(yīng)序列h(n)逼hd(n)。由于hd(n)往往是無限長序列，且是非因果的，所以用窗函數(shù)。w(n)將hd(n)截?cái)啵⑦M(jìn)行加權(quán)處理，得到：h(n)=hd(n)w(n)——（4.2）h(n)就作為實(shí)際設(shè)計(jì)的FIR數(shù)字濾波器的單位脈沖響應(yīng)序列，其頻率響應(yīng)函數(shù)H(ejω)為H(ejω)=——（4.3）式中，N為所選窗函數(shù)w(n)的長度。由第七章可知，用窗函數(shù)法設(shè)計(jì)的濾波器性能取決于窗函數(shù)w(n)的類型及窗口長度N的取值。設(shè)計(jì)過程中，要根據(jù)對阻帶最小衰減和過渡帶寬度的要求選擇合適的窗函數(shù)類型和窗口長度N。各種類型的窗函數(shù)可達(dá)到的阻帶最小衰減和過渡帶寬度見第七章。這樣選定窗函數(shù)類型和長度N后，求出單位脈沖響應(yīng)h(n)=hd(n)·w(n)，并按式(4.3)求出H(ejω)。H(ejω)是否滿足要求，要進(jìn)行驗(yàn)算。一般在h(n)尾部加零使長度滿足2的整數(shù)次冪，以便用FFT計(jì)算H(ejω)。如果要觀察細(xì)節(jié)，補(bǔ)零點(diǎn)數(shù)增多即可。如果H(ejω)不滿足要求，則要重新選擇窗函數(shù)類型和長度N，再次驗(yàn)算，直至滿足要求。如果要求線性相位特性，則h(n)還必須滿足：h(n)=±h(N-1–n)根據(jù)上式中的正、負(fù)號和長度N的奇偶性又將線性相位FIR濾波器分成四類。要根據(jù)所設(shè)計(jì)的濾波特性正確選擇其中一類。

例如，要設(shè)計(jì)線性相位低通特性，可選擇h(n)=h(N－1－n)一類，

而不能選h(n)=－h(huán)(N－1－n)一類。3.實(shí)驗(yàn)內(nèi)容及步驟(1)復(fù)習(xí)用窗函數(shù)法設(shè)計(jì)FIR數(shù)字濾波器一節(jié)內(nèi)容，閱讀本實(shí)驗(yàn)原理，掌握設(shè)計(jì)步驟。(2)編寫程序。①編寫能產(chǎn)生矩型窗、升余弦窗、改進(jìn)升余弦窗和二階升余弦窗的窗函數(shù)子程序。②編寫主程序。主程序框圖如圖4.1所示，僅供參考。其中幅度特性要求用dB表示。

開始讀入窗口長度N開始讀入窗口長度N調(diào)用子程序(函數(shù))計(jì)算H(k)=DFT[h(n)]結(jié)束計(jì)算hd(n)調(diào)用繪圖子程序(函數(shù))繪制|H(k)|幅度相位曲線調(diào)用窗函數(shù)子程序求w(n)計(jì)算h(n)=hd(n)w(n)設(shè)H(k)=DFT[h(n)]H(k)=HR(k)+jHI(k)|H(k)|=畫圖時，用20lg|H(k)|打印幅度特性。第k點(diǎn)對應(yīng)的頻率ωk=。為使曲線包絡(luò)更接近H(ejω)的幅度特性曲線，DFT變換區(qū)間要選大些。例如窗口長度N=33時，可通過在h(n)末尾補(bǔ)零的方法，使長度變?yōu)?4，再進(jìn)行64點(diǎn)DFT，則可得到更精確的幅度衰減特性曲線。FFT子程序可調(diào)用MALTAB函數(shù)fft實(shí)現(xiàn)。(3)上機(jī)實(shí)驗(yàn)內(nèi)容。①用升余弦窗設(shè)計(jì)一線性相位低通FIR數(shù)字濾波器，截止頻率ωc=rad。窗口長度N=15,33。要求在兩種窗口長度情況下，分別求出h(n)，打印出相應(yīng)的幅頻特性和相頻特性曲線，觀察3dB帶寬20dB帶寬?？偨Y(jié)窗口長度N對濾波特性的影響。設(shè)計(jì)低通FIR數(shù)字濾波器時，一般以理想低通濾波特性為逼近函數(shù)Hd(ejω)，即其中=(N-1)/2==②n=33，ωc=π/4，用四種窗函數(shù)設(shè)計(jì)線性相位低通濾波器。繪制相應(yīng)的幅頻特性曲線，觀察3dB和20dB帶寬以及阻帶最小衰減，比較四種窗函數(shù)對濾波器特性的影響。4.思考題(1)如果給定通帶截止頻率和阻帶截止頻率以及阻帶最小衰減，如何用窗函數(shù)法設(shè)計(jì)線性相位低通濾波器？寫出設(shè)計(jì)步驟。(2)如果要求用窗函數(shù)法設(shè)計(jì)帶通濾波器，且給定上、下邊帶截止頻率為ω1和ω2，試求理想帶通的單位脈沖響應(yīng)hd(n)。5.實(shí)驗(yàn)報(bào)告要求(1)簡述實(shí)驗(yàn)?zāi)康募霸怼?2)按照實(shí)驗(yàn)步驟及要求，比較各種情況下的濾波性能，說明窗口長度N和窗函數(shù)類型對濾波特性的影響。(3)總結(jié)用窗函數(shù)法設(shè)計(jì)FIR濾波器的主要特點(diǎn)。(4)簡要回答思考題。附錄：實(shí)驗(yàn)用MATLAB語言工具箱函數(shù)簡介MATLAB(ANTrixLABoratory)是Mathwords公司于1982年推出的一套高性能的數(shù)值計(jì)算和可視化軟件。它集數(shù)值分析、矩陣運(yùn)算、信號處理和圖形顯示于一體，構(gòu)成了一個方便且界面友好的用戶環(huán)境。MATLAB的推出得到了各個領(lǐng)域?qū)＜覍W(xué)者的廣泛關(guān)注。其強(qiáng)大的擴(kuò)展功能為各個領(lǐng)域的應(yīng)用提供了基礎(chǔ)。由各個領(lǐng)域的專家學(xué)者相繼推出了十幾種MATLAB工具箱，信號處理工具箱就是其中之一，借助這一工具箱進(jìn)行信號處理方面的研究，可大大節(jié)省時間。這些工具箱函數(shù)的使用可通過Help命令得到。為了便于使用，下面簡要介紹MATLAB的啟動、編程和運(yùn)行程序的過程，最后給出實(shí)驗(yàn)中可能用到的工具箱函數(shù)。(考慮到上機(jī)方便，附錄中外文均采用正白體。)1、啟動MATLAB系統(tǒng)及編程在Windows桌面上雙擊MATLAB快捷方式圖標(biāo)，啟動MATLAB系統(tǒng)，Windows工作平臺上自動彈出一個MALAB窗口，其中有CommandWindow窗口，稱其為MATLAB命令窗口。該窗口是用戶與MATLAB解釋器進(jìn)行通信的工作環(huán)境，所有的MATLAB函數(shù)、命令及程序都要在該命令窗口下運(yùn)行。在命令窗口中，用戶可以發(fā)出MATLAB命令。每條命令輸入并接回車鍵后，MATLAB系統(tǒng)便解釋并執(zhí)行之，顯示命令執(zhí)行結(jié)果。如果不想顯示中間結(jié)果，可在相應(yīng)命令后加分號即可。應(yīng)當(dāng)注意，在命令窗口輸人命令只適合一些簡單的運(yùn)算程序，不便于編輯修改較復(fù)雜的程序。這是因?yàn)橹苯釉诿畲翱谳斎氲拿钜唤?jīng)執(zhí)行后，只記錄其運(yùn)行結(jié)果變量，如果有錯誤，就得再輸入一遍（當(dāng)然，有較簡便的重新輸入方法）。這樣一來調(diào)試程序太麻煩。M文件為這種編輯運(yùn)行程序過程提供了方便。在MATLAB系統(tǒng)中，根據(jù)功能可將MATLAB系統(tǒng)所使用的外部文件分成幾類，并用不同的擴(kuò)展名作為標(biāo)識。(1)M文件以字母m為擴(kuò)展名，一般為ASCII碼文本文件，可以用任何文本編輯器進(jìn)行編輯。(2)MAT文件是MATLAB系統(tǒng)的二進(jìn)制數(shù)據(jù)文件，用于保存MATLAB系統(tǒng)所使用的數(shù)據(jù)。(3)MEX文件是經(jīng)過MATLAB編譯系統(tǒng)編譯的函數(shù)二進(jìn)制可執(zhí)行文件，可直接調(diào)入MATLAB系統(tǒng)中運(yùn)行，運(yùn)行速度比M文件快得多，所以常把已經(jīng)調(diào)試好且比較大的M文件編譯成MEX文件，供以后使用。M文件又分為程序M文件和函數(shù)M文件。下面簡要介紹實(shí)驗(yàn)過程中要用到的程序M文件的編輯與運(yùn)行方法。如前述，M文件為ASCII文本文件，可以用任何文本編輯器進(jìn)行編輯。本節(jié)只介紹在MATLAB命令窗口狀態(tài)下M文件的編輯運(yùn)行與修改。1)新建M文件在MATLAB命令窗口狀態(tài)，選File-New-Mfile，MATLAB系統(tǒng)彈出文件編輯窗口，類似于Windows的記事本（低版本MATLAB就是調(diào)用Windows記事本），這時一切編輯就如同記事本一樣。2)運(yùn)行與修改M文件MATLAB采用的是行命令模式，用戶每輸入一條命令，回車后MATLAB就解釋并執(zhí)行這條命令，再根據(jù)要求顯示運(yùn)算結(jié)果。此外，MATLAB也可執(zhí)行某個文件中的MATLAB語句序列（類似于DOS的批處理文件）。這兩種模式一起構(gòu)成了MATLAB的解釋環(huán)境。當(dāng)一個M文件編輯存盤后，在MATLAB命令窗口輸入M文件名并回車，則MATLAB系統(tǒng)就逐行解釋并執(zhí)行該M文件中的命令序列（即程序）。若發(fā)現(xiàn)錯誤，則打開M文件修改，再運(yùn)行。應(yīng)當(dāng)注意，如果已編輯好的M文件不在MATLAB系統(tǒng)搜索路徑中，用戶應(yīng)當(dāng)用path命令在MATLAB的搜索路徑中添加新的搜索路徑。例如，M文件存在C：\MYFILES路徑，下面的命令：path(‘C:MYFILES’，path)；將搜索路徑改為在搜索完當(dāng)前目錄之后，先搜索C:\MYFILES，再在當(dāng)前的搜索路徑目錄中搜索。2、實(shí)驗(yàn)用MATLAB函數(shù)簡介為了便于實(shí)驗(yàn)時查閱，下面簡要介紹幾種實(shí)驗(yàn)中可能要用到的MATLAB函數(shù)。1）abs功能：求絕對值（模值）。格式：y=abs(x)說明：y=abs(x)用于計(jì)算x的絕對值，當(dāng)x為復(fù)數(shù)時，得到的是復(fù)數(shù)的模值。當(dāng)x為字符串時，abs(x)得到字符串的各個字符的ASCII碼，例如，x=’123’，則abs(x)得到：4950512）angle功能：求相角。格式：φ=angle(h)說明：φ=angle(h)用于求復(fù)矢量或復(fù)矩陣的相角（以弧度為單位），相角介于-π和+π之間。例如，對復(fù)數(shù)h可用兩種方法表示：h=x+jy=mejφ則m和φ可由x+jy表示格式的復(fù)數(shù)求得m=abs(h)φ=angle(h)當(dāng)然，也可由m和φ求取h=x+jy格式中的x和y：h=m.*exp(i*φ)x=real(h)y=imag(h)3）conv功能：求卷積。格式：c=conv(a,b)說明：conv(a,b)用于求矢量a和b的卷積，即c(n)=a(k+1)b(n-k),n=1,2,…式中N為矢量a和b的最大長度。例如，當(dāng)a=[123],b=[156]時，則c=conv(a,b)c=413282718此函數(shù)可直接用于求兩個有限長序列的卷積。設(shè)x(n)和h(n)的長度分別為M和N，則計(jì)算二者卷積的MATLAB語句如下：y=conv(x,h)y的長度為N+M-1。4)filter功能：利用IIR濾波器或FIR濾波器對數(shù)據(jù)進(jìn)行濾波。格式：y=filter(b,a,x)[y,zf]=filter(b,a,x)y=filter(b,a,x,zi)說明：filter利用數(shù)字濾波器對數(shù)據(jù)進(jìn)行濾波，其實(shí)現(xiàn)采用直接五型結(jié)構(gòu)，因而適用于IIR和FIR兩種濾波器。濾波器的系統(tǒng)函數(shù)為即濾波器系數(shù)a=[a0a1a2…aN],b=[b0b1…bM]，輸入序列矢量為X。這里的標(biāo)準(zhǔn)形式為a0=1，如果輸入矢量a時，a0≠1，則MATLAB將自動進(jìn)行歸一化系數(shù)的操作；如果a0=0，則給出出錯信息。y=filter(b,a,x)利用給定系數(shù)矢量a和b對x中的數(shù)據(jù)進(jìn)行濾波，結(jié)果放入y矢量中，y的長度取max(N,M)。y=filter(b,a,x,zi)可在zi中指定x的初始狀態(tài)。[y,zf]=filter(b,a,x,zi)除了得到矢量y外，還得到x的最終狀態(tài)矢量zf。例C1設(shè)時域離散系統(tǒng)可由如下差分方程描述：y(n)=y(n-1)-0.9y(n-2)+x(n),n(1)計(jì)算n=0,…，120上的單位脈沖響應(yīng)；(2)計(jì)算n=0,…，120上的單位階躍響應(yīng)。解可利用filter函數(shù)直接求解。MATLAB程序?yàn)閑xc1.m。%Examplecla=[1,-1,0.9];b=1;x=zeros(1,120);n=[0:119];x(1)=1;%x(n)=δ(n)h=filter(b,a,x);subplot(2,1,1);stem(n,h,'.');title('ImpulseResponse');xlabel('n');ylabel('h(n)');x=ones(1,120);%x(n)=u(n)s=filter(b,a,x);subplot(2,1,2);stem(n,s,'.');title('StepResponse');xlabel('n');ylabel('s(n)');系統(tǒng)的單位脈沖響應(yīng)h(n)和單位階躍響應(yīng)s(n)如圖C．1所示。仿照該例題，可調(diào)用filter函數(shù)進(jìn)行實(shí)驗(yàn)三中的IIR濾波仿真。

圖C．1系統(tǒng)的單位脈沖響應(yīng)h(n)和階躍響應(yīng)s(n)5)freqz功能：數(shù)字濾波器的頻率響應(yīng)。格式：[h,w]=freqz(b,a,n)[h,f]=freqz(b,a,n,Fs)h=freqz(b,aw)h=freqz(b,a,f,Fs)freqz(b,a)說明：freqz用于計(jì)算數(shù)字濾波器的頻率響應(yīng)函數(shù)H(ejω)。H(z)具有如下形式:H(z)=矢量a=[1a1a2…aN]和b=[b0b1…bM]分別表示分母多項(xiàng)式A(z)和分子多項(xiàng)式B(z)的系數(shù)。[h,w]=freqz(b,a,n)可得到數(shù)字濾波器的n點(diǎn)頻率響應(yīng)值，這n個點(diǎn)均勻地分布在[0，π]上，并將這n個頻點(diǎn)的頻率記錄在w中，相應(yīng)的頻響值記錄在h中。要求n為大于零的整數(shù)，最好為2的整數(shù)次冪，以便采用FFT計(jì)算，以提高速度。[h,f]=freqz(b,a,n,Fs)用于對H(ejω)在［0，F(xiàn)s/2］上等間隔采樣n點(diǎn)，采樣點(diǎn)頻率及相應(yīng)頻響值分別記錄在f和h中。由用戶指定Fs（以Hz為單位）值。h=freqz(b，a，w)用于對H(ejω)在[0，2π]上進(jìn)行采樣，采樣頻率點(diǎn)由矢量w指定。h=freqz(b,a,f,Fs)用于對H(ejω)在[0，F(xiàn)s]上采樣，采樣頻率點(diǎn)由矢量f給定。freqz(b,a,n)用于在當(dāng)前圖形窗口中繪制出幅頻和相頻特性曲線。6)impz功能：計(jì)算H(z)相應(yīng)的單位脈沖響應(yīng)h(n)。格式：[h,t]=impz(b,a)[h,t]=impz(b,a,z)[h,t]=impz(b,a,n,Fs)impz(b,a)其中a和b分別為H(z)的分母和分子多項(xiàng)式系數(shù)。該函數(shù)的詳細(xì)說明請看MATLABHelp命令。7）fft功能：一維快速傅里葉變換(FFT)。格式：y=fft(x)y=fft(x,n)說明：fft函數(shù)用于計(jì)算矢量或矩陣的離散傅里葉變換。y=fft(x)利用FFT算法計(jì)算矢量x的離散傅里葉變換，當(dāng)x為矩陣時，y為矩陣x每一列的FFT。當(dāng)x長度為2的整數(shù)次冪時，fft采用基2FFT算法，否則采用稍慢的混合基算法。y=fft(x,n)采用n點(diǎn)FFT。當(dāng)x長度小于n時，fft函數(shù)自動在x尾部補(bǔ)零，以構(gòu)成n點(diǎn)數(shù)據(jù)；當(dāng)x的長度大于n時，fft截取x的前面n點(diǎn)數(shù)據(jù)進(jìn)行FFT。8）ifft功能：一維逆快速傅里葉變換（IFFT）。格式：y=ifft(x)y=ifft(x,n)9)plot功能：線型繪圖函數(shù)。格式：plot(v)plot(x,y)說明：plot(v)是一種最簡單的調(diào)用方式，v是長度為n的數(shù)值向量。Plot(v)的作用是在坐標(biāo)系中順序地用直線連接頂點(diǎn){i,v(i),i=1,2,…,n}生成一條折（曲）線。當(dāng)向量元素充分多時，既可生成一條光滑的曲線。在實(shí)驗(yàn)中，若FFT點(diǎn)數(shù)足夠多時，用plot打印的幅頻特性就很接近x(ejω)連續(xù)曲線。Plot(x,y)稱為參數(shù)式，參數(shù)x和y都是長度為n的向量。plot(x,y)在坐標(biāo)系中生成順序連接頂點(diǎn){x(i),y(i),i=1,2,…,n}的折（曲）線，這種調(diào)用可被用來生成參數(shù)方程的圖形。10）stem功能：繪制離散序列圖。格式：stem(y)stem(x,y)stem(…,’線端符號’)stem(…,’線型’)stem(…,’線型’,’線端符號’)說明：stem(y)和stem(x,y)分別與plot(x)和plot(x,y)的繪圖規(guī)則相同，只是stem繪制的是離散序列圖（或稱為“桿狀”圖），序列線端為圓圈。實(shí)驗(yàn)中用于繪制時域序列x(n)的波形圖和序列的離散傅里葉變換x(k)的幅度圖。后面三種調(diào)用方式中的‘線型’用于指定圖中序列線的線型和顏色，“線端符號”用于指定序列線端點(diǎn)表示序列值的繪圖形狀。線型、顏色和線端符號如表C.1所示。表C.1線型、顏色與線端符號線型顏色線端符號類型符號類型符號類型符號實(shí)線(缺省)-紅r實(shí)點(diǎn).點(diǎn)線:綠g星號*點(diǎn)虛線-.藍(lán)b圓圈o虛線--黑k三角形^表中只列出幾種共實(shí)驗(yàn)用，其他更多的符號和顏色清查閱MATLAB書。11）subplot功能：多坐標(biāo)設(shè)置與定位當(dāng)前坐標(biāo)系。格式：subplot(m,n,k)說明：s