1、 - 29 -DTMF電話中的撥號信號發(fā)生與檢測系統(tǒng)仿真摘 要雙音多頻（DTMF）信號是音頻電話中的撥號信號，由美國AT&T貝爾公司實驗室研制，并用于電話網(wǎng)絡(luò)中。這種信號制式具有很高的撥號速度，且容易自動檢測識別，很快就代替了原有的用脈沖計數(shù)方式的撥號制式。這種雙音多頻信號制式不僅用在電話網(wǎng)絡(luò)中，還可以用于傳輸十進制數(shù)據(jù)的其他通信系統(tǒng)中，用于電子郵件和銀行系統(tǒng)中。這些系統(tǒng)中用戶可以用電話發(fā)送DTMF信號選擇語音菜單進行操作。 DTMF信號系統(tǒng)是一個典型的小型信號處理系統(tǒng)，它要用數(shù)字方法產(chǎn)生模擬信號并進行傳輸，其中還用到了D/A變換器；在接收端用A/D變換器將其轉(zhuǎn)換成數(shù)字信號，并進行數(shù)字信號處理

2、。本文著重介紹了一種快速的頻率檢測算法GEORTZEL算法的原理，以及該算法在雙音多頻撥號檢測系統(tǒng)中的應(yīng)用。并根據(jù)此算法設(shè)計出了相應(yīng)的濾波器，利用MATLAB在計算機上對雙音多頻的信令產(chǎn)生及檢測接收系統(tǒng)進行了仿真測試。關(guān)鍵詞：雙音多頻，Goertzel算法，數(shù)字信號處理，傅里葉變換；目 錄摘要 11.緒論31.1 引言31.2 數(shù)字信號簡介31.3 數(shù)字濾波器41.4 頻率分析4 1.4.1 傅里葉級數(shù)及傅里葉變換 5 1.4.2 Z變換及系統(tǒng)函數(shù)5 1.4.3 離散系統(tǒng)傅里葉變換（DFT）62.雙音多頻（DTMF）通信系統(tǒng)6 2.1 雙音多頻（DTMF）信號簡介6 2.2 雙音多頻（DTMF

3、）信號的產(chǎn)生7 2.3 雙音多頻（DTMF）信號的檢測83.戈澤爾(Goertzel)算法 103.1 戈澤爾(Goertzel)算法原理 103.2 戈澤爾(Goertzel)算法程序 125.雙音多頻（DTMF）系統(tǒng)的仿真136.結(jié)束語 15參考文獻 16 1 緒論1.1 引言雙音多頻（DTMF）信號，由美國貝爾實驗室研制，用于音頻電話網(wǎng)絡(luò)中的撥號信號。一方面這種信號具有非常高的撥號速度，另一方面它便于自動檢測識別及電話業(yè)務(wù)的拓展，所以很快代替了原有的脈沖計數(shù)方式的撥號形式。自1876年美國貝爾發(fā)明電話以來，電話交換技術(shù)發(fā)展主要經(jīng)歷了三個階段：人工交換、機電交換、電子交換。早在1878年就

4、出現(xiàn)了人工供電制交換機，它借助話務(wù)員進行話務(wù)連接。15年后步進制的交換機問世，這種交換機屬于“直接控制”方式，用戶通過話機撥號脈沖直接控制步進接續(xù)器做升降和旋轉(zhuǎn)動作。從而自動完成用戶間的接續(xù)。這種交換機雖然實現(xiàn)了自動接續(xù)，但存在著速度慢、效率低、雜音大與機械磨損嚴(yán)重等缺點。直到1938年發(fā)明了縱橫制交換機由直接控制過渡到間接控制方式，隨后又出現(xiàn)了電子交換機及程控交換機等。用戶的撥號脈沖不再直接控制接線器動作，而先由記發(fā)器接收、存儲，然后通過標(biāo)志器驅(qū)動接線器，以完成用戶間接續(xù)。由于標(biāo)志器記取的是標(biāo)識信號，所以可以不采用不穩(wěn)定的脈沖計數(shù)方式，逐步發(fā)展到數(shù)字號碼標(biāo)識，所以快速、穩(wěn)定靈活的DTMF音頻

5、撥號方式取代了脈沖計數(shù)方式。1.2 數(shù)字信號簡介信號是反映消息的物理量，是消息的表現(xiàn)形式。由于非電的物理量可以通過各種傳感器較容易地轉(zhuǎn)換成電信號，而電信號又容易傳送和控制，所以使其成為應(yīng)用最廣的信號。電信號是指隨著時間而變化的電壓或電流，因此在數(shù)學(xué)描述上可將它表示為時間的函數(shù)，并可畫出其波形。信息可以通過電信號進行傳送、交換、存儲、提取等。較之傳統(tǒng)的方法，數(shù)字信號處理系統(tǒng)自身有著無可比擬的優(yōu)點：（1）便于加密處理（2）便于存儲、處理和交換（3）設(shè)備便于集成化、微型（4）便于構(gòu)成綜合數(shù)字網(wǎng)和綜合業(yè)務(wù)數(shù)字網(wǎng)（5）占用信道頻帶較寬一路模擬電話的HYPERLINK /view/488112.htm t

6、 _blank頻帶為4kHz帶寬，一路數(shù)字電話約占64kHz，這是模擬通信目前仍有生命力的主要原因。1.3 數(shù)字濾波器數(shù)字濾波器由數(shù)字加法器、單位延時和常數(shù)乘法器組成的一種算法或裝置。數(shù)字濾波器的功能是對輸入離散信號的數(shù)字代碼進行運算處理，即把輸入序列通過一定的運算變換成輸出序列，以達(dá)到改變信號頻譜的目的。數(shù)學(xué)公示表示如下： （1.3-1)通常可以用兩種方法實現(xiàn)數(shù)字濾波器：一種方法是把濾波器所完成的運算編成程序并讓計算機執(zhí)行，也就是采用就算機軟件實現(xiàn)；另一種方法是設(shè)計專用的數(shù)字硬件、專用的數(shù)字信號處理器或采用通用的數(shù)字信號處理器來實現(xiàn)。存在著兩種數(shù)字濾波器的的基本類型：無限長單位沖擊響應(yīng)（II

7、R）濾波器和有限長單位沖擊響應(yīng)（FIR）濾波器。它們的系統(tǒng)函數(shù)分別為： (1.3-2) (1.3-3)式（1.3-2）中的H(Z)稱為N階IIR濾波器函數(shù)式，（1.3-3）中的H(Z)稱為N-1階FIR濾波器函數(shù)。無限長單位沖激響應(yīng)（IIR）濾波器的單位沖激響應(yīng)h（n）是無限長的；系統(tǒng)函數(shù)H(Z)在有限Z平面（0Z）上有極點存在；結(jié)構(gòu)上存在著輸出到輸入的反饋，也就是結(jié)構(gòu)上是遞歸型的。1.4 頻率分析 任何給定信號的頻率分析，都包含著時域信號向其頻率分量的轉(zhuǎn)換。之所以要在頻域內(nèi)對信號進行描述，因為在實現(xiàn)信號處理時，通常都利用以頻率響應(yīng)形式描述的系統(tǒng)。而這些頻域信號處理概念時傅立葉變換是必不可少的

8、。同時Z變換對進行數(shù)字系統(tǒng)分析和實現(xiàn)也極其重要。1.4.1 傅立葉級數(shù)及傅立葉變換任何周期信號X(t)均可以表示為無窮多個諧波調(diào)諧相關(guān)的正弦波和復(fù)指數(shù)之和。具有周期T0(秒)周期信號x(t)，其基本的數(shù)學(xué)表達(dá)式為具有系列定義的傅立葉級數(shù)： (1.4.1-1)式中Ck為傅立葉級系數(shù)，第k個傅立葉系數(shù)Ck可以表示為 (1.4.1-2)用周期信號的傅立葉級數(shù)通過極限的方法導(dǎo)出的非周期信號的頻譜表示式，稱為傅立葉變換：傅立葉逆變換： (1.4.1-3)傅立葉正變換： = (1.4.1-4)1.4.2 Z變換及系統(tǒng)函數(shù)對于離散時間系統(tǒng)來說，與拉普拉斯對應(yīng)的變換是是Z變換，Z變換提供了離散時間信號和系統(tǒng)的

9、頻域描述，并為數(shù)字濾波器的設(shè)計和實現(xiàn)，提供了有力的工具。序列x(n)的Z變換定義為 ： (1.4.2-1)式中z是一個復(fù)變量，定義中對n求和是在和之間求和，稱為雙邊Z變換。還有一種稱為單邊Z變換的定義，如下式： (1.4.2-2)將系統(tǒng)的單位脈沖響應(yīng)h（n）進行Z變換，得到H（z）,一般稱H（z）為系統(tǒng)的系統(tǒng)函數(shù)，它表征了系統(tǒng)的頻域特性。對N階差分方程進行Z變換，得到系統(tǒng)函數(shù)的一般表示式 (1.4.2-3)令Z=ej即得系統(tǒng)的傳輸函數(shù)H(ej)。1.4.3 離散系統(tǒng)傅立葉變換（DFT）離散傅立葉變換是針對有限長序列或周期序列從存在的，相當(dāng)把序列的連續(xù)傅立葉變換加以抽樣，頻率的離散化造成時間函數(shù)

10、也呈周期，故級數(shù)應(yīng)限制在一個周期內(nèi)。令，則，我們即可得到離散傅立葉變換對為： (1.4.31) (1.4.32)其中表示有限長序列的抽樣點數(shù)，或周期序列一個周期的抽樣點數(shù)。2 雙音多頻（DTMF）通信系統(tǒng)2.1 雙音多頻（DTMF）信號過去的電話撥號時靠脈沖計數(shù)確定09這10個數(shù)字的，撥號速度慢，也不能擴展電話上網(wǎng)其他服務(wù)功能?，F(xiàn)在均采用雙音撥號。每一位號碼由兩個不同的單音頻組成，所有的頻率可分成高頻帶和低頻帶兩組。每一位號碼由一個高頻信號和一個低頻信號疊加形成，DTMF信號有16個編碼。一般情況下，聲音信號很難造成對DTMF接收器的錯誤觸發(fā)。電話中的雙音多頻信號有兩個作用：用撥號信號去控制交

11、換機接通被叫的用戶電話機；控制電話機的各種動作，如播放留言、語音信箱等。2.2 雙音多頻（DTMF）信號的產(chǎn)生DTMF編碼器基于兩個二階數(shù)字正弦波振蕩器，一個用于產(chǎn)生行頻，一個用于產(chǎn)生列頻。典型的DTMF信號頻率范圍是7001700Hz，選取8000Hz作為采樣頻率，DTMF雙音頻信號由兩個二階數(shù)字正弦振蕩器產(chǎn)生，一個用來產(chǎn)生行音頻信號，另一個產(chǎn)生列音頻信號。在電話鍵盤上的每次按鍵，會產(chǎn)生式（2.2-1）表示的兩個音調(diào)之和。規(guī)定用8KHz對DTMF信號進行采樣，采樣后得到時域離散信號為 (2.2-1)兩個頻率f1,f2唯一確定了被按壓的鍵。具體數(shù)字對應(yīng)的兩個頻率如表2.2.1所示。表2-1 D

12、TMF撥號頻率（第四列常被省略）頻低頻高1209Hz1336Hz1477Hz1633Hz697Hz123A770Hz456B852Hz789C942Hz*0#D形成上面序列的方法有兩種，即計算法和查表法。用計算法求正弦波的序列值容易，但實際中要占用一些計算時間，影響運行速度。查表法是預(yù)先將正弦波的各序列值計算出來，寄存在存儲器中，運行時只要按順序和一定的速度取出便可。這種方法要占用一定的存儲空間，但是速度快。因為采樣頻率是8000Hz，因此要求每125ms輸出一個樣本，得到的序列再送到D/A變換器和平滑濾波器，輸出便是連續(xù)時間的DTMF信號。DTMF信號通過電話線路送到交換機。2.3雙音多頻（

13、DTMF）信號的檢測在接收端，要對收到的將收到的雙音多頻信號進行檢測，即檢測兩個正弦波的頻率，以判斷其對應(yīng)的十進制數(shù)字或者符號。用數(shù)字方法進行檢測，需要將接收到的模擬音頻信號進行A/D變換，恢復(fù)為數(shù)字信號，然后檢測其中的音頻頻譜來確定所發(fā)送的數(shù)字。檢測方法有兩種，一種是用一組濾波器來提取所需頻率，判斷對應(yīng)的數(shù)字或符號；另一種是用快速傅立葉變換（FFT）算法的DFT對雙音多頻信號進行頻譜分析，有信號的幅度譜，判斷信號的兩個頻率，最后確定對應(yīng)的數(shù)字或符號。以下章節(jié)對兩種方法進行分析并對Goertzel算法進行詳細(xì)的推導(dǎo)。由上面可以知道一個DTMF信號是由兩個不通頻率f1和f2的正弦波組成，它可以用

14、下式表示 (2.3-1)我們可以考慮通過離散傅立葉變換進行信號的頻譜分析來檢測離散的雙音多頻碼。這是因為只有時域及頻域都是離散的情況下，才能適合于在計算機上運算；也就是周期的離散時間信號與離散頻率間的變換對。周期性離散時間信號x（n）的離散傅立葉變換周期性離散頻率函數(shù)X（k）兩相鄰譜線分量之間的角頻率增量與周期之間的關(guān)系可表示為， 代表信號的基頻。取樣頻率與取樣周期T的關(guān)系是 取樣角頻率 T為時域取樣間隔，在一個周期內(nèi)取樣點數(shù)為N。在自變量為t及f的情況下，在一個時域中對函數(shù)進行取樣，兩取樣點間增量的倒數(shù)，必是另一個域中函數(shù)的周期?，F(xiàn)序列的周期為NT，所以對頻譜取樣的間距是。以數(shù)字頻率表示時，

15、則頻譜間距是 k0，1，N1令并稱之為因子。離散傅立葉變換（DFT）可寫成如下形式：， 0 kN1 (2.3-2)DFT的表達(dá)式可以直接寫成：， k0，1，N1 ( 2.3-3) 一般來說，x（n）和都是復(fù)數(shù)，X（k）也是復(fù)數(shù)，每計算一個X（k）值，需要N次復(fù)數(shù)乘法(x（n）與相乘)以及(N1)次復(fù)數(shù)加法。而X（k）一共有N個點（k從0取到N1），所以完成整個DFT運算總共需要次復(fù)數(shù)乘法及N（N1）次復(fù)數(shù)加法。復(fù)數(shù)運算實際上是由實數(shù)運算來完成的，因此上式可寫成： (2.3-4)由此可見，一次復(fù)數(shù)乘法需用四次實數(shù)乘法和兩次實數(shù)加法；一次復(fù)數(shù)加法則需兩次復(fù)數(shù)加法。因而每運算一個需要4N次復(fù)數(shù)乘法和

16、2N2（N-1）=2(2N-1)次實數(shù)加法。 所以，直接計算DFT，乘法次數(shù)和加法次數(shù)都是和成正比的，當(dāng)N很大時，運算量是很可觀的，例如，當(dāng)N8時，DFT需64次復(fù)數(shù)乘法，而當(dāng)N1024時，DFT所需復(fù)數(shù)乘法為1，048，576次，即一百多萬次復(fù)數(shù)乘法運算，這對實時性很強的信號處理來說，對計算速度的要求太高了。因而需要改進對DFT的計算方法，以大大減少運算次數(shù)。仔細(xì)觀察DFT的運算量可以看出，利用系數(shù)的以下固有特性，可以減少DFT的運算量。(1)的共軛對稱性： (2)的周期性： (3) 的可約性：由此得出，這樣，（1）利用這些特性，使DFT運算中有些項可以合并；（2）利用的周期性和對稱性，可以

17、將長序列的DFT分解為短序列的DFT。而前面已經(jīng)說到，DFT的運算量使與成正比的，所以N越小越有利，因而小點數(shù)的DFT比大點數(shù)的DFT的運算量要小?；谶@樣的思路，快速傅立葉變換（FFT）發(fā)展起來了。如果按時間抽取的FFT算法，系統(tǒng)的運算量是，運算量明顯比DFT運算量減少。但它有一個缺點，需要全部抽樣數(shù)據(jù)到達(dá)時，才能進行計算；不適合雙音多頻信號頻譜的計算，滿足不了雙音撥號系統(tǒng)實時性的要求。因此開發(fā)出了適合于DMFT的特殊算法Goertzel算法。3 戈澤爾(Goertzel)算法3.1 戈澤爾(Goertzel)算法原理戈澤爾算法利用DFT中的旋轉(zhuǎn)因子的周期性，將DFT的運算轉(zhuǎn)換成一種線性濾波

18、運算。下面推導(dǎo)戈澤爾算法的計算公式和實現(xiàn)結(jié)構(gòu)。假設(shè)長度為N的序列x（n）的N點DFT用X（k）表示，因為 （3.1-1）因此 k=0,1,12,N-1 (3.1-2)按照上式定義序列 （3.1-3）這個方程可以解釋為有限長序列x(n)，與序列的卷積，因此表示為 （3.1-4）令 （3.1-5）則 （3.1-6）由上式，將看成是序列x（n）通過單位脈沖響應(yīng)為的濾波器的輸出，根據(jù)式（3.1-2）和式（3.1-3），得到 (3.1-7)那么，x（n）的DFT的第k點就是序列x(n)通過濾波器輸出的第n=N點樣值。那么N點DFT就是這N個濾波器分別對序列x(n)的響應(yīng)序列的第N點輸出。對式（3.1-5

19、）進行Z變換，得到濾波器的系統(tǒng)函數(shù) （3.1-8）該濾波器是一個一階純極點濾波器，極點為，極點頻率為。該濾波器的結(jié)構(gòu)圖如圖3-1所示，戈澤爾算法的原理方框圖如圖3-2所示。圖3-1單極點濾波器圖3-2戈澤爾算法的濾波器組在圖3-1中存在一次復(fù)乘算法，為了避免復(fù)乘算法，將一階純極點濾波器變?yōu)槎A濾波器，推導(dǎo)如下 （3.1-9）按照圖3-2所示。再按照該結(jié)構(gòu)圖，可以用兩個差分方程表示該二階濾波器，即 （3.1-10）= （3.1-11）對于音調(diào)檢測，只有的幅值平方是需要的，通過實現(xiàn)這種需求，可以使算法作進一步簡化，由（3.1-11）式，的平方幅值可以計算為： （3.1-12）這樣輸入信號是實序列，

20、用式（3.1-10）計算中間變量，用式（3.1-12）計算輸出信號的幅度，避免了復(fù)數(shù)的運算。由此得到=。因為有8種音頻要檢測，所以需要8個式（3.1-8）表示的濾波器，或者8個式（3.1-9）表示的濾波器。8個濾波器的中心頻率分別對應(yīng)8種音頻。3.2 Goertzel算法程序下面是按（3.1-12）式編寫的Goertzel算法子程序。它根據(jù)輸入序列x（n）和指定的DFT樣本的序號k計算待求的DFT樣本X。其核心語句只有四句，具體如下：function X=gfft(x,k)N=length(x);x1=x,0;d1=2*cos(2*pi*k/N)v=filter(1,1,-d1,1,x1);w

21、=exp(-i*2*pi*k/N);x=v(N+1)-w*v(N);5 雙音多頻（DTMF）系統(tǒng)的仿真DTMF信號的產(chǎn)生與識別仿真實驗在MATLAB環(huán)境下進行，編寫仿真程序，運行程序，送入8位電話號碼，程序自動產(chǎn)生每一位號碼數(shù)字相應(yīng)的DTMF信號，并送出雙頻聲音，再用DFT進行譜分析，顯示每一位號碼數(shù)字的DTMF信號的DFT幅度譜，安照幅度譜的最大值確定對應(yīng)的頻率，再安照頻率確定每一位對應(yīng)的號碼數(shù)字，最后輸出16位電話號碼，程序名為exp8。下面用MATLAB程序?qū)υ撓到y(tǒng)進行模擬，程序分四段：第一段（第27行）設(shè)置參數(shù)，并讀入8位電話號碼；第二段（第920行）根據(jù)輸入的8位電話號碼產(chǎn)生時域離散

22、DTMF信號，并連續(xù)發(fā)出8位號碼對應(yīng)的雙音頻聲音；第三段（第2225行）對時域離散DTMF信號進行頻率檢測，畫出幅度譜；第四段（第2633行）根據(jù)幅度譜的兩個峰值，分別查找并確定輸入的8位電話號碼。程序如下：% DTMF雙頻撥號信號的生成程序：%clear all;clc;tm=1,2,3,65;4,5,6,66;7,8,9,67;42,0,35,68; % DTMF信號代表的16個數(shù)N=205;K=18,20,22,24,31,34,38,42;%16個基頻對應(yīng)的16個k值f1=697,770,852,941; % 行頻率向量f2=1209,1336,1477,1633; % 列頻率向量TN

23、=input(鍵入16位電話號碼= ); % 輸入16位數(shù)字TNr=0; %接收端電話號碼初值為零for l=1:8; %分別對每位號碼數(shù)字處理：產(chǎn)生信號，發(fā)聲，檢測 d=fix(TN/10(8-l); %計算出第l位號碼數(shù)字 TN=TN-d*10(8-l); for p=1:4; for q=1:4; if tm(p,q)=abs(d); break,end % 檢測與l位號碼相符的列號q endif tm(p,q)=abs(d); break,end % 檢測與l位號碼相符的行號pEnd n=0:1023; % 為了發(fā)聲，加長序列 x = sin(2*pi*n*f1(p)/8000) +

24、sin(2*pi*n*f2(q)/8000);% 構(gòu)成雙音頻信號 sound(x,8000); % 發(fā)出聲音 pause(0.1) %相鄰信號響聲之間加0.1秒停頓 % 接收檢測端的程序 X=goertzel(x(1:N),K+1); % 用Goertzel算法計算16點DFT樣本 val = abs(X); % 列出八點DFT的摸 subplot(4,2,l); stem(K,val,.);grid;xlabel(k);ylabel(|X(k)|) % 畫出16點DFT(k)的幅度 axis(10 50 0 120) limit = 80; %基頻檢測門限為80 for s=5:8; if

25、val(s) limit, break, end % 查找列號 end for r=1:4; if val(r) limit, break, end % 查找行號end18 TNr=TNr+tm(r,s-4)*10(16-l); %將16位號碼表示成一個16位數(shù)，以便顯示enddisp(接收端檢測到的號碼為：) % 顯示接收到的16位電話號碼disp(TNr)編譯程序，彈出操作界面，點擊操作界面的數(shù)字按鍵，程序開始運行，數(shù)據(jù)經(jīng)連接函數(shù)與接收端程序連接輸入鍵值，經(jīng)Goertzel算法的處理將頻譜疊加，在頻譜圖顯示窗口顯示輸入頻率（axes2）和疊加頻率(axes1),同時接收端程序檢根據(jù)合成頻率

26、檢測出鍵值，在editor窗口顯示檢測到輸入的數(shù)字。結(jié)束語本文采用GEORTZEL算法，實現(xiàn)了對雙音多頻信號的檢測，在MATLAB環(huán)境下，對雙音多頻信號的解碼檢測過程進行了模擬仿真。本文只是完成了雙音多頻撥號系統(tǒng)的軟件仿真，如何實現(xiàn)高性價比的硬件系統(tǒng)是下一步需要研究和完成的工作。在做MATLAB界面，做好以后同時會產(chǎn)一個M文件，只要把代碼編寫在對應(yīng)的控件中即可以方便的實現(xiàn)仿真。試驗結(jié)果表明基于MATLAB 的DTMF 編解碼具有速度快、精度高、穩(wěn)定性好以及便于靈活應(yīng)用等特點。而利用MATLAB語言，可以較好地、快速地進行算法的仿真。參考文獻1 高西全、丁玉美、闊永紅.數(shù)字信號處理原理、實現(xiàn)及應(yīng)

27、用.電子工業(yè)出版社,20062 程佩青.數(shù)字信號處理教程.清華大學(xué)出版社,2001年, 第二版3 Joyce van de Vegte,侯正信,王國安譯.數(shù)字信號處理基礎(chǔ).電子科工業(yè)出版社,2003年4 趙霞.雙音多頻信號產(chǎn)生及解碼研究.微電子學(xué),2001年,第6期5 陳懷琛.數(shù)字信號處理教程MATLAB釋疑與實現(xiàn).電子工業(yè)出版社,2004附錄資料：不需要的可以自行刪除各類濾波器的MATLAB程序理想低通濾波器IA=imread(lena.bmp);f1,f2=freqspace(size(IA),meshgrid);Hd=ones(size(IA);r=sqrt(f1.2+f2.2);Hd(

28、r0.2)=0;Y=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 二、理想高通濾波器IA=imread(lena.bmp);f1,f2=freqspace(size(IA),meshgrid);Hd=ones(size(IA);r=s

29、qrt(f1.2+f2.2);Hd(r0.2)=0;Y=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); Butterworth低通濾波器IA=imread(lena.bmp);f1,f2=freqspace(size(IA),

30、meshgrid);D=0.3;r=f1.2+f2.2;n=4;for i=1:size(IA,1) for j=1:size(IA,2) t=r(i,j)/(D*D); Hd(i,j)=1/(tn+1); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,non

31、e,Facelighting,phong); Butterworth高通濾波器IA=imread(lena.bmp);f1,f2=freqspace(size(IA),meshgrid);D=0.3;r=f1.2+f2.2;n=4;for i=1:size(IA,1) for j=1:size(IA,2) t=(D*D)/r(i,j); Hd(i,j)=1/(tn+1); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(ui

32、nt8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 高斯低通濾波器IA=imread(lena.bmp);IB=imread(babarra.bmp);f1,f2=freqspace(size(IA),meshgrid);D=100/size(IA,1);r=f1.2+f2.2;Hd=ones(size(IA);for i=1:size(IA,1) for j=1:size(IA,2) t=r(i,j)/(D*D); Hd(i,j)=

33、exp(-t); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 高斯高通濾波器IA=imread(lena.bmp);IB=imread(babarra.bmp);f1,f2=freqspace(size(I

34、A),meshgrid);%D=100/size(IA,1);D=0.3;r=f1.2+f2.2;for i=1:size(IA,1) for j=1:size(IA,2) t=r(i,j)/(D*D); Hd(i,j)=1-exp(-t); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,in

35、terp,Edgecolor,none,Facelighting,phong); 梯形低通濾波器IA=imread(lena.bmp);IB=imread(babarra.bmp);f1,f2=freqspace(size(IA),meshgrid);%D=100/size(IA,1);D0=0.1;D1=0.4;r=sqrt(f1.2+f2.2);Hd=zeros(size(IA);Hd(r=D0 & r(i,j)=D1 Hd(i,j)=(D1-r(i,j)/(D1-D0); end endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifft

36、shift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 梯形高通濾波器IA=imread(lena.bmp);IB=imread(babarra.bmp);f1,f2=freqspace(size(IA),meshgrid);%D=100/size(IA,1);D0=0.1;D1=0.4;r=sqrt(f1.2+f2.2);

37、Hd=ones(size(IA);Hd(r=D0 & r(i,j)=D1 Hd(i,j)=(D0-r(i,j)/(D0-D1); end endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 用其他方法編寫的理想低通、理想高通、Butterworth低通、同態(tài)濾波程序理想低通i1=imread(lena.bmp);i2=imnoise(i1,salt & pepper,0.1);f=double(i2);k=fft2(f);g=fftshift(k);N1,N2=size(g);d0=50;u0=floor(N1/2)+1;v0=floor(N2/2)+1;for i=1:N1 for j=1:N2 d=sqrt(i-u0)2+(j-v0)2); if d=d0 h=1;