數(shù)字信號處理MATLAB實(shí)現(xiàn)及綜合應(yīng)用

第7章數(shù)字信號處理MATLAB實(shí)現(xiàn)及綜合應(yīng)用7.1數(shù)字信號處理MATLAB實(shí)現(xiàn)7.2在雙音撥號系統(tǒng)中的應(yīng)用7.3在wav信號分析方面的應(yīng)用7.4在自適應(yīng)濾波器方面的應(yīng)用隨著信息學(xué)科和計(jì)算機(jī)學(xué)科的飛速發(fā)展，數(shù)字信號處理的重要性日益顯著，數(shù)字信號的處理已成為電子信息學(xué)科中的工程技術(shù)人員必不可少的知識。由于數(shù)字信號處理的概念比較抽象，再加上其數(shù)值計(jì)算又比較繁瑣，因此，國外很早就開始把MATLAB用于數(shù)字信號處理的教學(xué)過程，并取得了很好的經(jīng)驗(yàn)，現(xiàn)在MATLAB已經(jīng)成為解決數(shù)字信號處理問題的公認(rèn)的標(biāo)準(zhǔn)軟件。為了使讀者更好地掌握數(shù)字信號處理的基本概念和基本理論，本章首先在7.1節(jié)給出前面各章的MATLAB典型例題。通過這些例題，幫助讀者弄清一些重要的概念，加深對課程基本概念和基本原理的理解，學(xué)會利用計(jì)算機(jī)軟件解決在理論學(xué)習(xí)中不易解決的問題。然后在7.2節(jié)~7.4節(jié)分別給出了數(shù)字信號處理技術(shù)在通信、信號分析、自適應(yīng)濾波等方面的實(shí)際應(yīng)用，并給出了MATLAB仿真程序和結(jié)果分析，這些內(nèi)容可以擴(kuò)寬讀者數(shù)字信號處理技術(shù)應(yīng)用的知識面，促進(jìn)理論與實(shí)際相結(jié)合。本章和下一章上機(jī)實(shí)驗(yàn)都是課程學(xué)習(xí)的重要實(shí)踐環(huán)節(jié)，它不僅能幫助同學(xué)們理解和掌握課程的重要概念和基本內(nèi)容，而且也有益于讀者深入學(xué)習(xí)和應(yīng)用信號處理知識來解決實(shí)際問題。7.1數(shù)字信號處理MATLAB實(shí)現(xiàn)7.1.1MATLAB簡介7.1.2離散時(shí)間信號與系統(tǒng)分析應(yīng)用舉例7.1.3離散傅里葉變換應(yīng)用舉例7.1.4IIR濾波器設(shè)計(jì)常用函數(shù)及舉例7.1.5FIR濾波器設(shè)計(jì)常用函數(shù)及舉例7.1.1MATLAB簡介MATLAB是由美國MathWorks公司推出的軟件產(chǎn)品，“MATLAB”是MatrixLaboratory（矩陣實(shí)驗(yàn)室）的縮寫。MATLAB是一個(gè)完整的并可擴(kuò)展的計(jì)算機(jī)環(huán)境，是一種進(jìn)行科學(xué)和工程計(jì)算的交互式程序語言。它的基本數(shù)據(jù)單元是不需要指定維數(shù)的矩陣，解決同樣的數(shù)值計(jì)算問題，使用MATLAB效率要提高許多倍。MATLAB采用開放式的環(huán)境，你可以讀到它的原碼、了解它的算法，并能改變當(dāng)前的函數(shù)或增添自己編寫的函數(shù)。目前，在國內(nèi)外大學(xué)和研究機(jī)構(gòu)中，MATLAB是一種非常流行的計(jì)算機(jī)語言，許多重要的學(xué)術(shù)刊物上發(fā)表的論文均是用MATLAB來分析計(jì)算和繪制各種圖形。它還是一種有力的教學(xué)工具，在大學(xué)的線性代數(shù)、自動控制理論、數(shù)字信號處理、動態(tài)系統(tǒng)仿真等課程的教學(xué)中，已成為標(biāo)準(zhǔn)的教學(xué)工具。1．變量2．基本操作符3．輸入、輸出語句4．庫函數(shù)5．繪圖操作表7-1 MATLAB的常用繪圖命令及功能實(shí)際的圖形顯示有多種情況，下面分別加以說明。（1）需要用多個(gè)圖形窗顯示多個(gè)圖形時(shí)，應(yīng)在圖形顯示語句之前加上圖形窗指令figure；（2）在一個(gè)圖形窗中顯示多個(gè)圖形時(shí)，應(yīng)在圖形顯示語句之前調(diào)用圖形窗分塊函數(shù)subplot()；（3）在一個(gè)圖形坐標(biāo)系中顯示多條曲線時(shí)，要在plot()語句前先執(zhí)行指令holdon用以保持當(dāng)前曲線，顯示完成之后，再執(zhí)行指令holdoff；（4）需要選擇坐標(biāo)軸尺寸時(shí)，在plot()語句之后調(diào)用axis()函數(shù)。信號處理工具箱覆蓋了經(jīng)典信號處理的大多數(shù)內(nèi)容，是一個(gè)非常優(yōu)秀的算法研究、輔助設(shè)計(jì)工具。它提供了大量的工具函數(shù)、濾波器分析和設(shè)計(jì)工具，常用的函數(shù)有以下幾種。（1）設(shè)計(jì)IIR濾波器相關(guān)函數(shù)buttord、butter、cheblord、cheb2ord、ellipord、cheby1、cheby2等。6．用于數(shù)字信號處理的工具箱——SIGNAL（2）設(shè)計(jì)FIR濾波器相關(guān)函數(shù)fir1、fir2、firls、remezord、remez等；窗函數(shù)：boxcar、hanning、hamming、blackman、kaiser。（3）傅里葉變換相關(guān)函數(shù)FFT（一維快速傅里葉變換）、IFFT（一維快速傅里葉逆變換）、FFT2（二維快速傅里葉變換）、abs（絕對值或模）、angle（相角）等。（4）頻率變換相關(guān)函數(shù)freqz、freqs、lp2lp、lp2hp、lp2bp、lp2bs。（5）頻率響應(yīng)相關(guān)函數(shù)freqz、freqs。7.1.2離散時(shí)間信號與系統(tǒng)分析應(yīng)用舉例1．zplane(b,a)函數(shù)或pzplotz(b,a)函數(shù)2．freqz函數(shù)3．filter函數(shù)圖7-1零極點(diǎn)圖圖7-2頻率響應(yīng)圖圖7-3仿真波形7.1.3離散傅里葉變換應(yīng)用舉例在MATLAB中，可用函數(shù)dft來實(shí)現(xiàn)離散傅里葉變換運(yùn)算。而且，在實(shí)用中，MATLAB還提供了內(nèi)部函數(shù)來快速進(jìn)行離散傅里葉變換運(yùn)算。它采用了優(yōu)化的算法，并且程序是用目的碼編寫的，使得它有極高的運(yùn)算速度，命名為fft，反變換是ifft。調(diào)用方式為y=fft(x)或y=fft(x,N)。圖7-44點(diǎn)DTFT和DFT圖圖7-5DFT圖（N=8，N=16）圖7-6x(n)的前10點(diǎn)數(shù)據(jù)對應(yīng)的x(n)、X(ejw)、X(k)

圖7-7x(n)補(bǔ)零至100點(diǎn)對應(yīng)的x(n)、X(ejw)、X(k)

圖7-9128點(diǎn)x(n)對應(yīng)的x(n)、X(ejw)、X(k)7.1.4IIR濾波器設(shè)計(jì)常用函數(shù)及舉例1．最小階數(shù)選擇函數(shù)[n,wn]=buttord/cheb1ord/cheb2ord/ellipord(wp,ws,Rp,Rs,'s')1．最小階數(shù)選擇函數(shù)[z,p,k]=buttap(n)/cheb1ap(n,Rp)/cheb2ap(n,Rs)/ellipap(n,Rp,Rs)2．模擬低通原型函數(shù)（1）低通到低通的頻率變換[b1,a1]=lp2lp(b,a,w0)，其中，w0為低通濾波器的截止頻率（rad/s）。（2）低通到高通的頻率變換[b1,a1]=lp2hp(b,a,w0)，其中，w0為高通濾波器的截止頻率（rad/s）。3．頻率變換函數(shù)（3）低通到帶通的頻率變換[b1,a1]=lp2bp(b,a,w0,Bw)，其中，w0為帶通濾波器的中心頻率，Bw為帶通濾波器的帶寬。當(dāng)濾波器通帶的下截止頻率為w1，上截止頻率為w2時(shí)，w0=sqrt(w1*w2)，Bw=w2-w1。（4）低通到帶阻的頻率變換[b1,a1]=lp2bs(b,a,w0,Bw)，其中，w0為帶阻濾波器的中心頻率，Bw為帶阻濾波器的帶寬。當(dāng)濾波器通帶的下截止頻率為w1，上截止頻率為w2時(shí)，w0=sqrt(w1*w2)，Bw=w2-w1。MATLAB信號處理工具箱還提供了模擬濾波器設(shè)計(jì)的完全工具函數(shù)：butter、cheby1、cheby2、ellip、besself。用戶只需一次調(diào)用就可自動完成全部設(shè)計(jì)過程，編程十分簡單。這些工具函數(shù)既適用于模擬濾波器設(shè)計(jì)，也適用于數(shù)字濾波器。4．完全設(shè)計(jì)函數(shù)（1）巴特沃思濾波器：[b,a]=butter(n,wn,'ftype','s')其中，n為濾波器階數(shù)；wn為濾波器截止頻率；'s'為模擬濾波器，缺省時(shí)為數(shù)字濾波器。（2）切比雪夫?yàn)V波器：[b,a]=cheby1(n,Rp,wn,'ftype','s')或[b,a]=cheby2(n,Rs,wn,'ftype','s')（3）橢圓濾波器：[b,a]=ellip(n,Rp,Rs,wn,'ftype','s')[bz,az]=impinvar/bilinear(b,a,fs)5．脈沖響應(yīng)和雙線性變換函數(shù)其中，b、a分別為模擬濾波器的分子和分母多項(xiàng)式系數(shù)向量；fs為采樣頻率，單位為Hz，缺省時(shí)fs默認(rèn)為1Hz；bz，az分別為數(shù)字濾波器分子和分母多項(xiàng)式系數(shù)向量。圖7-10不同采樣頻率對所設(shè)計(jì)數(shù)字濾波器頻率響應(yīng)的影響圖7-11ChebyshevⅠ型數(shù)字高通濾波器幅頻特性圖7-127階ChebyshevⅡ型數(shù)字低通濾波器幅度和相位圖7.1.5FIR濾波器設(shè)計(jì)常用函數(shù)及舉例1．窗函數(shù)法設(shè)計(jì)FIR濾波器（2）計(jì)算截止頻率為wc的理想低通濾波器的單位脈沖響應(yīng)hd(n)，語句如下hd=ideal_lp(wc,N)（3）h(n)為h=hd.*win。（4）已知h(n)，為了得到頻率響應(yīng)的幅度、相位及群時(shí)延，需用語句如下[db,mag,pha,grd,w]=freqz_m(h,1)圖7-13設(shè)計(jì)的實(shí)際低通濾波器結(jié)果（1）h(n)為h=real(ifft(H,N))。（2）MATLAB程序中，可以用freqz計(jì)算頻率響應(yīng)，但不能從中確定振幅響應(yīng)。這是因?yàn)樵贛ATLAB函數(shù)中，不存在與abs函數(shù)功能相當(dāng)而能找到振幅的函數(shù)。但在MATLAB中，提供了4種擴(kuò)展函數(shù)，可以計(jì)算FIR4種類型濾波器的振幅響應(yīng)。●Hr_type1●Hr_type2●Hr_type3●Hr_type42．頻率采樣法設(shè)計(jì)FIR濾波器前面介紹的窗函數(shù)法設(shè)計(jì)FIR濾波器，完全是根據(jù)窗函數(shù)法的設(shè)計(jì)思路和設(shè)計(jì)流圖一步一步進(jìn)行的，這有助于掌握方法的原理。在實(shí)際應(yīng)用中，還可以采用另外一種方法——直接利用MATLAB中的內(nèi)部函數(shù)來進(jìn)行FIR濾波器的設(shè)計(jì)。MATLAB信號處理工具箱提供了兩個(gè)用窗函數(shù)法設(shè)計(jì)FIR濾波器的函數(shù)，見表7-2。值得一提的是，F(xiàn)IR濾波器的窗函數(shù)設(shè)計(jì)法不僅可以設(shè)計(jì)標(biāo)準(zhǔn)頻率響應(yīng)（低通、高通、帶通及帶阻）濾波器，還可設(shè)計(jì)任意頻率響應(yīng)的多帶濾波器。3．直接采用MATLAB內(nèi)部函數(shù)設(shè)計(jì)FIR濾波器表7-2 用窗函數(shù)法設(shè)計(jì)FIR濾波器的內(nèi)部函數(shù)函數(shù)名函

數(shù)

功

能fir1用窗函數(shù)法設(shè)計(jì)具有標(biāo)準(zhǔn)頻率響應(yīng)的FIR濾波器fir2用窗函數(shù)法設(shè)計(jì)具有任意頻率響應(yīng)的多帶FIR濾波器圖7-14例7-10運(yùn)行結(jié)果圖（1）函數(shù)fir1例7-11設(shè)計(jì)一個(gè)低通濾波器，通帶0～200Hz，阻帶300～1000Hz，通帶波動7%，阻帶衰減60dB，采樣頻率為2000Hz。圖7-15例7-11濾波器幅度和相位圖（2）函數(shù)fir2圖7-16被污染信號的頻譜特性圖7-17設(shè)計(jì)的多帶濾波器幅頻特性圖7-18原始信號、污染信號及濾波后信號的比較7.2在雙音撥號系統(tǒng)中的應(yīng)用7.2.1電話的雙音撥號系統(tǒng)7.2.2Goertzel算法7.2.3檢測DTMF信號的DFT參數(shù)選擇雙音多頻（DualToneMultiFrequency，DTMF）是用按鍵進(jìn)行電話撥號的制式，它不單單用在電話中，還可以用于傳輸十進(jìn)制數(shù)據(jù)的其他通信系統(tǒng)中。DTMF也廣泛應(yīng)用于電子郵件和銀行系統(tǒng)。在這些系統(tǒng)中，用戶可以從電話發(fā)送DTMF信號來選擇菜單語音進(jìn)行操作。圖7-19雙頻撥號頻率分配7.2.1電話的雙音撥號系統(tǒng)7.2.2Goertzel算法圖7-20計(jì)算X(k)的系統(tǒng)結(jié)構(gòu)7.2.3檢測DTMF信號的DFT參數(shù)選擇DTMF信號的參數(shù)選擇要考慮多方面的因素。首先電話數(shù)字化的采樣頻率已知為8kHz，實(shí)際上為了抑止語音干擾，檢測系統(tǒng)除了規(guī)定的8個(gè)頻率之外，還要檢測它們的2次倍頻處的DFT幅度。這些2倍頻處幅度一般很小，如果基頻和倍頻分量同時(shí)都大，那就可能是外來聲音的干擾。表7-4 各個(gè)撥號頻率及N＝205時(shí)其對應(yīng)的DFT序號8個(gè)基頻Hz準(zhǔn)確k值最近整數(shù)k值絕對誤差2次諧波頻率Hz準(zhǔn)確k值最近整數(shù)k值絕對誤差69717.861180.139139435.024350.02477019.531200.269154038.692390.30885221.833220.167170442.813430.18794124.113240.113188247.285470.285120930.981310.019241860.752610.248133634.235340.235267267.134670.134147737.848380.152295474.219740.219163341.846420.154326682.058820.0587.2.4DTMF信號檢測的MATLAB仿真圖7-21發(fā)送端產(chǎn)生的雙頻信號時(shí)域波形圖7-22接收端雙頻信號在8個(gè)近似基頻上的DFT幅度7.3在wav信號分析方面的應(yīng)用7.3.1wav文件的一次性傅里葉 變換7.3.2wav文件的分段傅里葉分析本節(jié)將通過分析計(jì)算機(jī)中的wav文件來進(jìn)一步討論數(shù)字信號處理中的信號分析方法。首先介紹幾個(gè)常用的函數(shù)：[x，fs，bits]＝waveread（‘filename’）——這是一個(gè)MATLAB中讀取wav文件的數(shù)據(jù)函數(shù)。其中的x表示一長串的數(shù)據(jù)，一般是兩列（立體聲）；fs是該wav文件在采集時(shí)用的采樣頻率；bits是指在進(jìn)行A/D轉(zhuǎn)換時(shí)用的量化位長（一般是8bits或16bits）。[d]＝fft(w,N)——這是MATLAB中快速傅里葉變換（FFT）函數(shù)的一種輸入輸出形式。w是一列波形數(shù)據(jù)；N是指用多少點(diǎn)的FFT，一般選擇2的整次方數(shù)（如16，128，1024等），因?yàn)檫@樣就可以使用優(yōu)化的蝶形算法；d是頻域的輸出。由于FFT（DFT）的對稱性，且輸入的是實(shí)數(shù)，F(xiàn)FT結(jié)果的復(fù)數(shù)序列是共軛反對稱的，所以它們模是大小對稱的，一般來說只用其一半的數(shù)據(jù)就可以了。sound（w，fs，bits）——和前面的waveread一樣的參數(shù)表示，它將序列的數(shù)據(jù)通過聲卡轉(zhuǎn)化為聲音。下面先介紹wav文件的一次性傅里葉變換。7.3.1wav文件的一次性傅里葉變換1．聲波主要頻率的分析圖7-23信號的時(shí)域波形圖7-24信號頻域的幅值2．聲音信號的重構(gòu)前面介紹了對wav文件直接進(jìn)行一次性傅里葉變換的分析，但存在兩個(gè)問題：第一，對分析的wav文件直接一次性做FFT，點(diǎn)數(shù)太大，而在實(shí)際運(yùn)用中DSP的FFT的點(diǎn)數(shù)是有限的，一般只能達(dá)到千點(diǎn)；第二，語音是分音節(jié)的，應(yīng)該對它分段分析，所以下面引入分段處理數(shù)據(jù)，并連接每個(gè)段的頻域結(jié)果。7.3.2