一、頻率域基礎
頻率域濾波實際上是將圖像進行傅里葉變換,然后在變換域進行處理,然后進行傅里葉反變換轉換回空間域,原理是用傅里葉變換表示的函數特征完全可以通過傅里葉反變換來重建,而且不會丟失任何信息(因為任何周期或非周期函數都可以表示為不同頻率的正弦函數和余弦函數之和的形式)。實際上,空間域濾波和頻率域濾波經常是對應的:

空間濾波實際上是圖像與各種空間濾波器(模板)的卷積,而空間卷積的傅里葉變換是頻率域中相應變換的乘積,因此頻率域濾波可以用圖像的傅里葉變換乘以相應的頻率域濾波器。
連續變量函數的傅里葉變換暫且不提,我們直接來看圖像的二維離散傅里葉變換,令f(x,y)表示一幅大小為[M,N]像素的數字圖像,由F(u,v)表示圖像的二維離散傅里葉變換(DFT):

指數項由歐拉公式獲得,可以將其展開為正弦函數和余弦函數,頻率域是使用u和v作為(頻率)變量,由F(u,v)構成坐標系。由u和v構成的大小為[M,N]的矩形區域稱為頻率矩形,大小和輸入圖像的大小相同。
離散傅里葉反變換(IDFT)的形式為:

頻率域原點處變換的值(F(0,0))稱為傅里葉變換的直流分量,不難看出,F(0,0)等於f(x,y)平均值的M*N倍。還需要注意的是即使f(x,y)是實函數,它的變換通常也是復數。因此,傅里葉變換常用的是F(u,v)的頻譜和相角,並將其顯示為一幅圖像。matlab有專有的快速傅里葉變換(FFT)算法實現,可以直接調用:F=fft2(f)。一般來說,使用傅里葉變換濾波時,需要對輸入數據進行0填充:F=fft2(f,P,Q),結果函數的大小為P*Q。令R(u,v)和I(u,v)表示其實部和虛部,則傅里葉頻譜(幅度)為:

傅里葉譜可以實驗函數abs直接獲得:S=abs(F)。變換的相角定義為:

相角可以用atan2調用:phi=atan2(imag(F),real(F));或者直接使用函數angle:phi=angle(F)。有些地方還會提到功率譜,其實就是頻譜的平方。實函數的傅里葉變換是關於原點共軛對稱的,因此頻譜是關於原點偶對稱的,而相角是關於原點奇對稱的。由此得出的結論是傅里葉變換在u,v方向上是無窮周期的,而通過反變換得到的圖像也是無窮周期的,我們感興趣的是在區間[0,M-1]上的一個完整周期,變換之前,通過將f(x,y)乘以(-1)^(x+y)可以得到期望的周期,調用fftshift可以實現。最重要的是零頻率項與f(x,y)的平均值成正比:

傅里葉變換后,頻率矩陣的原點不在中心,可以使用函數fftshift將變換的原點移動到頻率矩形的中心:Fc=fftshift(F);與中心處亮度值占支配地位的8位顯示相比,譜中的范圍很大,顯示不明顯,通過對數變換可以解決這個問題:S2=log(1+abs(Fc))。具體的例子如下:

matlab代碼如下:
%傅里葉變換 f=imread('G:\數字圖像處理(岡薩雷斯)\DIP3E_CH04_Original_Images\DIP3E_Original_Images_CH04\Fig0424(a)(rectangle).tif'); subplot(231);imshow(f);title('原圖'); F = fft2(f); %傅里葉變換 subplot(232);imshow(F,[]);title('傅里葉變換'); S=abs(F); %傅里葉譜 subplot(233);imshow(S,[]);title('傅里葉頻譜') Fc = fftshift(F); %將變換的原點移動到頻率矩形的中心 Sc=abs(Fc); %居中的譜 subplot(234);imshow(Sc,[]);title('居中的譜') S2 = log(1+abs(Sc)); subplot(235);imshow(S2,[]);title('對數變換后的譜') f2=real(ifft2(F)); %傅里葉反變換 subplot(236);imshow(f2,[]);title('傅里葉反變換') phi = atan2(imag(F),real(F)); %subplot(236);imshow(phi,[]);
頻率域變換中的低頻與圖像中緩慢變化的灰度成分有關,高頻由灰度的尖銳過渡造成,如邊緣和噪聲等,衰減高頻而通過低頻的濾波器(低通濾波器)將模糊一幅圖像,相反的高通濾波器將增強尖銳的細節,但會降低圖像的對比度。
二、有填充和無填充的濾波效果
空間濾波由濾波模板卷積一幅圖像組成,函數相對位移,直到一個函數全部滑過另一個函數為止,在頻率域中,F和H是周期的,表明在離散頻率域中執行的卷積也是周期的,保證空間和頻率域給出相同結果的唯一辦法就是適當的0填充。假設函數f(x,y)和h(x,y)的大小均為M*N,通過對f和g補零,構造大小均為P*Q的填充函數:P>=2M-1,Q>=2N-1。填充函數為:
function PQ = paddedsize(AB,CD,PARAM) if nargin ==1 PQ = 2*AB; elseif nargin ==2&&~ischar(CD) PQ = AB+CD-1; PQ = 2*ceil(PQ/2); elseif nargin ==2 m = max(AB); P=2^nextpow2(2*m); PQ=[P,P]; elseif (nargin==3)&&strcmpi(PARAM,'pwr2') m=max([AB CD]); P=2^nextpow2(2*m); PQ=[P,P]; else error('wrong number of imputs'); end
具體是否填充的例子如下:
當濾波器位於虛線圖像的一側時,他將會遇到鄰近這一側的周期分量中的一個相同區域,因為一個恆定區域的均值是同一常數,因此結果的這一部分沒有被模糊。代碼如下:
f=imread('G:\數字圖像處理(岡薩雷斯)\DIP3E_CH04_Original_Images\DIP3E_Original_Images_CH04\Fig0432(a)(square_original).tif'); subplot(221);imshow(f);title('原圖'); [M,N]=size(f); %f = im2double(f); %轉換為浮點類型再計算 [f,revertclass] = tofloat(f); F = fft2(f); %傅里葉變換 H = lpfilter('gaussian',M, N,10); %頻域高斯低通濾波器 G=H.*F; %濾波(乘積) g=ifft2(G); %傅里葉反變換 %g=im2uint8(g); g=revertclass(g); subplot(222);imshow(g);title('無填充時頻率域低通濾波'); %經過填充后的傅里葉變換及其濾波 PQ = paddedsize(size(f)); %以0填充, Fp = fft2(f,PQ(1),PQ(2)); %用填充的方式計算傅里葉變換 Hp = lpfilter('gaussian',PQ(1),PQ(2),2*10); %頻域高斯低通濾波器 Gp = Hp.*Fp; gp = ifft2(Gp); gpc = gp(1:size(f,1),1:size(f,2)); %只取右上角一個周期的圖像 gpc = im2uint8(gpc); subplot(223);imshow(gpc);title('有填充時頻域低通濾波'); subplot(224);imshow(gp);title('未裁剪的全填充圖像');