圖像頻率域處理程序設計

1、武漢理工大學專業(yè)綜合課程設計說明書目錄摘要21 MATLAB軟件功能簡介32 數字圖像處理簡介42.1 數字圖像處理的目的42.2 數字圖像處理的方法52.3 圖像頻域處理的概述53 二維傅里葉變換73.1 二維連續(xù)傅里葉變換73.2 二維離散傅里葉變換83.3 二維離散傅里葉變換的性質83.4 周期延拓在卷積中的作用94 圖像頻率域處理程序設計步驟134.1 找出兩幅大小不一的256級的灰度圖像134.2 頻率域處理程序設計175 運行結果及結果分析206課程設計心得體會22參考文獻24摘要圖像的頻域處理是指根據一定的圖像模型，對圖像頻譜進行不同程度修改的技術。二維正交變換是圖像處理中常用的

2、變換，其特點是變換結果的能量分布向低頻成份方向集中，圖像的邊緣、線條在高頻成份上得到反映，因此正交變換在圖像處理中得到廣泛運用。傅里葉作為一種典型的正交變換，在數學上有比較成熟和快速的處理方法。卷積特性是傅里葉變換性質之一，由于它在通信系統(tǒng)和信號處理中的重要地位應用最廣。在用頻域方法進行卷積過程中尤其要注意傅里葉變換的周期性，注意周期延拓的重要作用，本次課設將對此作詳細的介紹。關鍵字：頻域處理，二維傅里葉變換，卷積，周期延拓1 MATLAB軟件功能簡介MATLAB的名稱源自Matrix Laboratory,1984年由美國Mathworks公司推向市場。它是一種科學計算軟件，專門以矩陣的形式

3、處理數據。MATLAB將高性能的數值計算和可視化集成在一起，并提供了大量的內置函數，從而被廣泛的應用于科學計算、控制系統(tǒng)和信息處理等領域的分析、仿真和設計工作。MATLAB軟件包括五大通用功能，數值計算功能（Nemeric）、符號運算功能（Symbolic）、數據可視化功能（Graphic）、數字圖形文字統(tǒng)一處理功能（Notebook）和建模仿真可視化功能（Simulink）。其中，符號運算功能的實現(xiàn)是通過請求MAPLE內核計算并將結果返回到MATLAB命令窗口。該軟件有三大特點，一是功能強大；二是界面友善、語言自然；三是開放性強。目前，Mathworks公司已推出30多個應用工具箱。 MAT

4、LAB在線性代數、矩陣分析、數值及優(yōu)化、數值統(tǒng)計和隨機信號分析、電路與系統(tǒng)、系統(tǒng)動力學、圖像處理、控制理論分析和系統(tǒng)設計、過程控制、建模和仿真、通信系統(tǒng)以及財政金融等眾多領域的理論研究和工程設計中得到了廣泛應用。MATLAB 的應用范圍非常廣，包括信號和圖像處理、通訊、控制系統(tǒng)設計、測試和測量、財務建模和分析以及計算生物學等眾多應用領域。附加的工具箱（單獨提供的專用 MATLAB 函數集）擴展了 MATLAB 環(huán)境，以解決這些應用領域內特定類型的問題。 MATLAB應用：MATLAB 產品族可以用來進行以下各種工作：數值分析 數值和符號計算 工程與科學繪圖 控制系統(tǒng)的設計與仿真 數字

5、圖像處理 數字信號處理 通訊系統(tǒng)設計與仿真 財務與金融工程 2 數字圖像處理簡介2.1 數字圖像處理的目的數字圖像處理(Digital Image Processing)是通過計算機對圖像進行去除噪聲、增強、復原、分割、提取特征等處理的方法和技術。數字圖像處理的產生和迅速發(fā)展主要受三個因素的影響：一是計算機的發(fā)展；二是數學的發(fā)展（特別是離散數學理論的創(chuàng)立和完善）;三是廣泛的農牧業(yè)、林業(yè)、環(huán)境、軍事、工業(yè)和醫(yī)學等方面應用需求的增長。一般來講，對圖像進行處理（或加工、分析）的主要目的有三個方面：（1）提高圖像的視感質量，如進行圖像的亮度、彩色變換，增強、抑制某些成分，對圖像進行幾何變換等，以改善圖

6、像的質量。（2）提取圖像中所包含的某些特征或特殊信息，這些被提取的特征或信息往往為計算機分析圖像提供便利。提取特征或信息的過程是模式識別或計算機視覺的預處理。提取的特征可以包括很多方面，如頻域特征、灰度或顏色特征、邊界特征、區(qū)域特征、紋理特征、形狀特征、拓撲特征和關系結構等。（3）圖像數據的變換、編碼和壓縮，以便于圖像的存儲和傳輸。不管是何種目的的圖像處理，都需要由計算機和圖像專用設備組成的圖像處理系統(tǒng)對圖像數據進行輸入、加工和輸出。圖像增強的目的是采用某種技術手段，改善圖像的視覺效果，或將圖像轉換到更適合人眼觀察和機器分析、識別的形式，以便從圖像中獲取更有用的信息。圖像增強的基本方法可分為兩

7、大類：空間域和頻域方法。空間域是指圖像平面自身，這類方法是以對圖像的像素直接處理為基礎的；而頻率域處理技術是以修改圖像的傅里葉變換為基礎的。兩者的具體方法包括以下內容：（1）空間域處理：點處理，模板處理即鄰域處理（2）頻率域處理：高、低通濾波，同態(tài)濾波等。2.2 數字圖像處理的方法傅里葉變換能將滿足一定條件的某個函數表示成三角函數（正弦和/或余弦函數）或者它們的積分的線性組合。在不同的研究領域，傅里葉變換具有多種不同的變體形式，如連續(xù)傅里葉變換和離散傅里葉變換。最初傅里葉分析是作為熱過程的解析分析的工具被提出的1。傅里葉變換屬于諧波分析。傅里葉變換的逆變換容易求出，而且形式與正變換非常類似。正

8、弦基函數是微分運算的本征函數，從而使得線性微分方程的求解可以轉化為常系數的代數方程的求解。在線性時不變的物理系統(tǒng)內，頻率是個不變的性質，從而系統(tǒng)對于復雜激勵的響應可以通過組合其對不同頻率正弦信號的響應來獲取。卷積定理指出：傅里葉變換可以化復雜的卷積運算為簡單的乘積運算，從而提供了計算卷積的一種簡單手段。離散形式的傅里葉變換可以利用數字計算機快速的實現(xiàn)（其算法稱為快速傅里葉變換算法（FFT）。MATLAB中提供的變換函數（1）fft2:用于計算二維快速傅立葉變換，語句格式：B=fft2(I,m,n) 按指定的點數計算m，返回矩陣B的大小為m×n，不寫默認為原圖像大小（2）fftn:用于

9、計算n維快速傅立葉變換（3）fftshift:用于將變換后的圖像頻譜中心從矩陣的原點移到矩陣地中心，語法格式B=fftshift(I) （4）ifft2:用于計算圖像的二維傅立葉反變換，語法格式：B=ifft2(i)（5）ifftn:用于計算n維傅立葉變換，快速卷積實驗：傅立葉變換一個重要特性是可以實現(xiàn)快速卷積設A為M×N矩陣，B為P×Q的矩陣，快速卷積方法如下： （1）對A和B補0，使其大小都為(M+P-1)×(N+Q-1) （2）利用fft2對矩陣A和B進行二維變換 （3）將兩個FFT結果相乘，利用ifft2對得到的乘積進行傅立葉反變換2.3 圖像頻域處理的概

10、述圖像的頻率是表征圖像中灰度變化劇烈程度的指標，是灰度在平面空間上的梯度。如大面積的沙漠在圖像中是一片灰度變化緩慢的區(qū)域，對應的頻率值很低；而對于地表屬性變化劇烈的邊緣區(qū)域在圖像中是一片灰度變化劇烈的區(qū)域，對應的頻率值較高。頻域處理是指根據一定的圖像模型，對圖像頻譜進行不同程度修改的技術，通常作如下假設：1)引起圖像質量下降的噪聲占頻譜的高頻段；2)圖像邊緣占高頻段；3)圖像主體或灰度緩變區(qū)域占低頻段?；谶@些假設，可以在頻譜的各個頻段進行有選擇性的修改。為什么要在頻率域研究圖像增強 （1）可以利用頻率成分和圖像外表之間的對應關系。一些在空間域表述困難的增強任務，在頻率域中變得非常普通。（2）

11、濾波在頻率域更為直觀，它可以解釋空間域濾波的某些性質。 （3）可以在頻率域指定濾波器，做反變換，然后在空間域使用結果濾波器作為空間域濾波器的指導。 （4）一旦通過頻率域試驗選擇了空間濾波，通常實施都在空間域進行。3 二維傅里葉變換由于圖像的頻率是表征圖像中灰度變化劇烈程度的指標，是灰度在平面空間上的梯度。傅立葉變換在實際中的物理意義，設f是一個能量有限的模擬信號，則其傅立葉變換就表示f的譜。從純粹的數學意義上看，傅立葉變換是將一個函數轉換為一系列周期函數來處理的。從物理效果看，傅立葉變換是將圖像從空間域轉換到頻率域，其逆變換是將圖像從頻率域轉換到空間域。換句話說，傅立葉變換的物理意義是將圖像的

12、灰度分布函數變換為圖像的頻率分布函數，傅立葉逆變換是將圖像的頻率分布函數變換為灰度分布函數。 3.1 二維連續(xù)傅里葉變換 如果二維連續(xù)函數f(x,y)滿足狄里赫萊條件，則將有下面的傅立葉變換對存在：與一維傅立葉變換類似，二維傅立葉變換的傅立葉譜和相位譜為：3.2 二維離散傅里葉變換一個M×N大小的二維函數f(x,y)，其離散傅立葉變換對為 ：在數字圖像處理中，圖像一般取樣為方形矩陣，即N×N，則其傅立葉變換及其逆變換為 ：3.3 二維離散傅里葉變換的性質 離散傅里葉變換主要有以下性質：1.平移性質2.分配律3.尺度變換（縮放）4.旋轉性5.周期性和共軛對稱性6.平均值7.可

13、分性8.卷積9.相關性。這里主要簡述周期性，卷積相關內容會在下一節(jié)中介紹。離散傅里葉變換有如下周期性性質：反變換也是周期性的：頻譜也是關于原點對稱的：這些等式的有效性是建立在二維離散傅里葉變換公式基礎上的。圖像的周期性在圖像處理中有非常重要的作用，下面會在卷積部分繼續(xù)闡述周期性的相關內容。3.4 周期延拓在卷積中的作用基于卷積理論，頻率域的乘法相當于空間域的卷積，反之亦然。當處理離散變量和傅里葉變換時，要記住不同函數所包含的周期性。雖然可能不太直觀，但周期性是定義離散傅里葉變換對時產生的數學副產品。周期性是處理操作的一部分,不應忽視。圖1列舉了周期性的重要性。圖1 傅里葉變換周期性圖1左邊（a

14、e）：兩個離散函數的卷積，右邊（fj）：相同函數的卷積，考慮DFT周期性的應用。圖的左邊一列是用下式的一維形式計算的卷積：在此詳細地解釋卷積運算的過程。為簡化表示，簡單的數字將代替那些表示函數長度和高度的通用符號。圖1(a)和(b)是兩個要進行卷積的函數。每個函數包含400個點。卷積的第一步是將一個函數關于原點進行鏡像映射(倒轉)，在本例情況下，對第二個函數進行，在圖1(c)中以h(-m)示出。下一步是將h(-m)滑過f(m)。這要增加一個常數x到h(-m)，即變成h(x-m)，如圖1(d)所示。注意只有一個置換值。在第一次遇到時這個簡單步驟通常是引起混亂的根源。而這恰好是卷積計算的全部關鍵。

15、換言之，為了執(zhí)行卷積，倒轉了一個函數，并將它滑過另一個函數。在每一個置換點(的每一個值)都要計算式的全部總和。這個總和不比在給定位移處f和h乘積的和更太。位移x的范圍為h完全滑過f需要的所有值。圖1(e)顯示了h完全滑過f后的結果，并在x的每個點計算式。在此例中，為使h(x-m)完全滑過f,x值的范圍是從0到799。這幅圖是兩個函數的卷積，要清楚地記住卷積中的變量是x.從上面介紹的卷積理論可知，由F(u)H(u)的傅里葉反變換能得到同樣的準確結果。但是，從前面對周期性的討論又知離散傅里葉變換自動地將輸入函數周期化。換言之，采用DFT允許在頻率域進行卷積計算，但函數必須看做周期性的，且周期等于函

16、數的長度?？梢酝ㄟ^圖1右邊一列考察這種隱含的周期性。圖1(f)同圖1(a)一樣，但同樣的函數在兩個方向上周期性地無限擴展(擴展部分用虛線表示)。從圖1(g)到圖1(i)同樣應用該擴展?，F(xiàn)在，可以通過將h(x-m)滑過f(m)進行卷積。如前面一樣，變化x完成滑動。然而，h(x-m)的周期性擴展產生了圖1左邊的計算中所沒有的值。例如，在圖1(i)中，當x=0時，看到h(x-m)右側第一個擴展周期的一部分進 入圖1(f)中所示的f(m)(從原點開始)的一部分。當h(x-m)向右滑動時，在f(m)中的那部分開始向右側移出，但被h(x-m)左側相同部分所取代。這引起卷積產生一個常量值，如圖1(j)所示的

17、0，100的一段.從100到4OO的一段是正確的，但周期性是周而復始的，這樣就引起卷積函數尾部的一部分丟失，由圖1(j)和圖1(e)實線部分的比較可以看出這一點。在頻率域，該過程需要計算圖1(a)和(b)中函數的傅里葉變換。根據卷積理論，兩個變換要相乘，再計算傅里葉反變換。結果包含40O個點的卷積，如圖1(j)的實線部分所示。簡單的解釋表明當使用傅里葉變換得出卷積函數時，錯誤地處理周期性將得到錯誤的結論。結果，在開頭有錯誤數據，結尾將丟失數據。問題的解決辦法很簡單。假設f和h分別由A和B個點組成。對兩個函數同時添加零，以使它們具有相同的周期，表示為P。這個過程產生擴展的或延拓的函數，如下所示：

18、和可以看出,除非選擇PA+B-1，否則卷積的獨立周期將會混疊。已經在圖1中看到了這種現(xiàn)象的結果，這通常歸于纏繞誤差。若P=A+B-1，周期便會鄰接起來。若PA+B-1，周期將會是分隔開的，分隔的程度等于P與A+B-1的差。 擴展后的卷積結果如圖2所示。在這里，選擇P=A+B-1(799)，即可知卷積周期是相鄰的。遵循與前面的解釋相同的過程，得到如圖2所示的卷積函數。該結果的一個周期與圖1(e)相同，是正確的。這樣，如果要在頻率域計算卷積，應該：(1)得到兩個擴展序列的傅里葉變換(每個序列有8OO個點)；(2)將兩個變換相乘；(3)計算傅里葉反變換。結果便得到正確的8OO個點的卷積函數。見圖2中

19、周期加重的部分。圖2 卷積函數這些概念擴展到二維函數時遵循了相同的前提。假設有f(x,y)和h(x,y)兩幅圖像，大小分別為A×B和C×D。如同一維情況，這些行列必須假定在x方向上有相同的周期P，在y方向上有相同的周期Q。二維卷積的混疊可由選擇如下周期避免：擴展f(x,y)和h(x,y)形成如下周期性序列：為了簡化圖例，假設f和h是方形的，且大小相同。圖3 二維函數周期延拓圖3對二維函數周期延拓的說明。圖3（a）沒有延拓執(zhí)行二維卷積的結果；圖3（b）合格的函數延拓；圖3（c）正確的卷積結果。圖3(a)顯示了圖像沒有延拓時得到的濾波結果。這通常是由于沒有對一幅輸入圖像進行延拓

20、就進行傅里葉變換，然后又乘上同樣大小的函數(也沒有延拓)，計算傅里葉反變換。結果就是與輸入圖像相同的大小為A×B的圖像，如圖3(a)左上象限所示。如同一維情況，圖像前面邊沿(阻影部分)由于周期性而引入了錯誤數據，而在尾部邊沿將丟失數據。如圖3(b)所示，通過對輸入圖像和函數進行合適的延拓，將得到正確的、大小為P×Q的過濾圖像。這幅圖像在兩個坐標方向上是原始圖像的兩倍大小，有原始圖像4倍數量的像素點。4 圖像頻率域處理程序設計步驟4.1 找出兩幅大小不一的256級的灰度圖像選擇兩幅圖片如圖4圖5 圖4 源圖片1 圖5源圖片2（1）檢查源圖片1和源圖片2格式>>in

21、fo=imfinfo('1.jpg')info=imfinfo('2.jpg')顯示如下info = Filename: '1.jpg' FileModDate: '15-Aug-2014 13:44:02' FileSize: 51822 Format: 'jpg' FormatVersion: '' Width: 440 Height: 737 BitDepth: 24 ColorType: 'truecolor' FormatSignature: '' Numb

22、erOfSamples: 3 CodingMethod: 'Huffman' CodingProcess: 'Sequential' Comment: info = Filename: '2.jpg' FileModDate: '08-Jan-2015 08:38:58' FileSize: 16439 Format: 'jpg' FormatVersion: '' Width: 300 Height: 400 BitDepth: 8 ColorType: 'grayscale' F

23、ormatSignature: '' NumberOfSamples: 1 CodingMethod: 'Huffman' CodingProcess: 'Sequential' Comment: 用MATLAB檢查發(fā)現(xiàn)是圖片1是truecolor 格式，ColorType: 'truecolor'。（2）用MATLAB將其轉換為灰度圖像a=imread('1.jpg');b=rgb2gray(a);imwrite(b,'3.jpg')>> info=imfinfo('3.jpg

24、')顯示如下info = Filename: '3.jpg' FileModDate: '08-Jan-2015 08:54:41' FileSize: 43637 Format: 'jpg' FormatVersion: '' Width: 440 Height: 737 BitDepth: 8 ColorType: 'grayscale' FormatSignature: '' NumberOfSamples: 1 CodingMethod: 'Huffman' Codi

25、ngProcess: 'Sequential' Comment: 在 ColorType: 'grayscale' 屬性行如此顯示，則說明以3.jpg命名的文件為灰度圖像如圖6所示。圖6 圖片3（3）檢查圖像灰度級在命令窗口輸入f=imread('2.jpg')whos fg=imread('3.jpg')whos g顯示如下，表明圖像為256級灰度圖像Name Size Bytes Class f 400x300 120000 uint8 arrayName Size Bytes Class g 737x440 324280 u

26、int8 arrayGrand total is 120000 elements using 120000 bytes4.2 頻率域處理程序設計MATLAB中提供的變換函數（1）fft2:用于計算二維快速傅立葉變換，語句格式：B=fft2(I,m,n)按指定的點數計算m,返回矩陣B的大小為m×n,不寫默認為原圖像大小。（2）ifft2:用于計算圖像的二維傅立葉反變換，語法格式：B=ifft2(i)%直接卷積程序>> I1=imread('2.jpg');I2=imread('3.jpg');I5=conv2(I1,I2);figure(3)

27、;imshow(I5,);title('直接函數卷積得到的圖像 (黃深)')%正確的頻域處理程序I1=imread('2.jpg');I2=imread('3.jpg');m1,n1=size(I1);m2,n2=size(I2);I1(m1+m2-1,n1+n2-1)=0;I2(m1+m2-1,n1+n2-1)=0;I3=ifft2(fft2(I1).*fft2(I2);I3=I3(1:m1+m2-1,1:n1+n2-1);I3=real(I3);figure(1);imshow(I3,);title('正確延拓頻域法得到的卷積圖像(黃

28、深)')%比較頻域方法與直接卷積的結果，顯示差矩陣并且顯示錯誤數據數F=minus(I3,I5);figure(4)imshow(F);title('正確延拓差矩陣的二值圖像(黃深)')s=0;for i=1:m1+m2-1 for j=1:n1+n2-1 if (minus(abs(F(i,j),0.000001)>0) s=s+1; end; end;end;disp(sprintf('差錯 1: %d',s);%補0不夠的頻域處理程序I1=imread('2.jpg');I2=imread('3.jpg');m

29、1,n1=size(I1);m2,n2=size(I2);I1(m1+m2-100,n1+n2-100)=0;I2(m1+m2-100,n1+n2-100)=0;I3=ifft2(fft2(I1).*fft2(I2);I3=I3(1:m1+m2-100,1:n1+n2-100);I3=real(I3);I3(m1+m2-1,n1+n2-1)=0;figure(2);imshow(I3,);title('補0不夠頻域法得到的卷積圖像(黃深)')%比較頻域方法與直接卷積的結果，顯示差矩陣并且顯示錯誤數據數F=minus(I3,I5);figure(5)imshow(F);title

30、('補0不夠的差矩陣的二值圖像(黃深)')s=0;for i=1:m1+m2-100 for j=1:n1+n2-100 if (minus(abs(F(i,j),0.000001)>0) s=s+1; end; end;end;disp(sprintf('差錯 2: %d',s);5 運行結果及結果分析在MATLAB中輸入程序后，顯示的卷積結果如下: 圖7 正確延拓頻域法得到的卷積圖像 圖8補0不夠頻域法得到的卷積圖像 圖9直接函數卷積得到的圖像 圖10正確延拓差矩陣的二值圖像圖11 補0不夠的差矩陣的二值圖像比較圖7和圖9，看不出兩個圖像有任何區(qū)別。通

31、過作差，認為舍入誤差小于0.000001的均可作為0來處理，這里S=785，差值矩陣的二值圖像全為黑，可以認為兩圖幾乎沒有任何區(qū)別，即頻域方法的卷積結果是完全正確的。比較圖8和圖9，表面上也看不出兩個圖像有什么區(qū)別，圖8的靠左和靠上部分有亮度增加，這部分是疊加錯誤，而靠下和靠右部分是兩條黑杠，這是補零的數據，也就是原來丟失的數據。通過檢測差值矩陣，S=99099,錯誤的有很多，即沒有補0的頻域方法計算的結果不正確。值得注意的是這里差值矩陣應該四周都是白色，因為左邊和上邊是混疊錯誤的地方應該為，行數:100，列數100；同理右邊和下邊是數據丟失人為補0的地方也有與混疊相同的行數和列數。但因為這里

32、gyy.jpg周圍為0，因此正確卷積的結果也為0，因此差矩陣得到的相應區(qū)域也為0，顯示的2值圖像就看不到白色地方了。通過以上分析說明，二維圖像或矩陣的線性卷積可以通過補零周期延拓后，經二維傅里葉變換相乘，再做反變換來實現(xiàn)。而不補零或補零不足，用此方法求得卷積圖像靠左靠上會有疊加錯誤和靠下靠右會有數據丟失。6課程設計心得體會通過本次課程設計，我主要了解了用MATLAB處理數字圖像中的問題，特別是數字圖像的頻域分析法，快速傅里葉變換，把數字圖像技術和MATLAB等通信類科目的內容應用到本課程設計中來，進一步鞏固復習數字圖像處理技術，MATLAB等課程，以達到融會貫通的目的。本次課設的主要任務是用MATLAB編程來實現(xiàn)數字圖像的傅里葉計算，卷積原理。開始我對數字圖像技術了解不是很多，通過查閱相關資料，我熟悉了用MATLAB處理數字圖