第4章頻率域濾波著名的法國數(shù)學(xué)家傅里葉在其著作《熱分析理論》中指出：任何周期函數(shù)都可以分解為不同頻率的正弦或余弦級數(shù)的形式，即傅里葉級數(shù)。該方法從本質(zhì)上完成了空間信息到頻域信息的變換，通過變換將空域信號處理問題轉(zhuǎn)化為頻域信號處理。20世紀(jì)50年代之后，隨著計算機的發(fā)展以及快速傅里葉變換算法的出現(xiàn)，使得頻率域濾波相關(guān)技術(shù)得到了廣泛的發(fā)展。形象的講，傅里葉變換可以看作是“數(shù)學(xué)中的棱鏡”，可以將任何周期函數(shù)分解為不同頻率的信號成分。頻率域變換為信號處理提供了不同的思路，有時在空域無法處理的問題，通過頻率域變換卻變得非常容易。傅里葉變換是十九世紀(jì)數(shù)學(xué)界和工程界最輝煌的成果之一。一直是信號處理領(lǐng)域應(yīng)用最廣泛、實踐效果也比較好的分析手段。頻域變換也是眾多圖像變換中的一種。為了更加有效的對數(shù)字圖像進行處理，常常需要將原始圖像以某種方式變換到另外一個空間，并利用圖像在變換空間中的特有性質(zhì)對圖像信息進行加工，然后再轉(zhuǎn)換回圖像空間就可以得到所需的效果。類似此類的轉(zhuǎn)換方法稱為圖像變換技術(shù)。圖像變換是雙向的，一般將從圖像空間轉(zhuǎn)換到其它空間的操作稱為正變換，由其它空間轉(zhuǎn)換回圖像空間稱為反變換或者逆變換。圖4-1圖像變換傅里葉變換中將圖像看作二維信號，其水平方向和垂直方向作為二維空間的坐標(biāo)軸，將圖像本身所在的域稱為空間域。圖像灰度值隨空間坐標(biāo)變換的節(jié)奏可以通過頻率來度量，稱為空間頻率或者頻率域。針對數(shù)字圖像的傅里葉變換是將原始圖像通過傅里葉變換轉(zhuǎn)換到頻率域，然后在頻率域當(dāng)中對圖像進行處理的方法?；诟道锶~變換的數(shù)字圖像頻域處理過程如圖4-2所示。首先通過正向傅里葉變換將原始圖像從空域轉(zhuǎn)換到頻率域。然后使用頻域濾波器將某些頻率成分過濾掉，保留某些特定頻率。最后使用傅里葉逆變換將濾波后的頻率域圖像重新轉(zhuǎn)換到空間域，得到處理后的圖像。圖4-2頻域圖像濾波過程相較于圖像空間域處理，頻率域圖像處理有如下優(yōu)點：一、可以通過頻率域成分的特殊性質(zhì)完成一些空間域圖像處理難以完成的任務(wù)。二、頻率域處理更有利于信號處理的解釋，它可以對濾波過程中產(chǎn)生的某些效果做出比較直觀的解釋。三、頻率域濾波器可以作為空間濾波器設(shè)計的指導(dǎo)，通過傅里葉反變換可以將頻率域濾波器轉(zhuǎn)換為空間域變換的操作。通過頻率域濾波做前期設(shè)計，然后實施階段用空間域濾波實現(xiàn)。本章主要基于空間域濾波對數(shù)字圖像進行處理，主要以傅里葉變換做案例，包含以下內(nèi)容：5.1節(jié)講解二維傅里葉變換，5.2節(jié)介紹二維傅里葉變換的性質(zhì)，5.3節(jié)講解快速傅里葉變換，5.4節(jié)對常用的圖像頻域濾波器進行介紹。4.1二維傅里葉變換傅里葉變換是常見正交數(shù)學(xué)變換的一種，可以將一維信號或函數(shù)分解為具有不同頻率不同幅度的正弦信號或余弦信號的組合。著名的法國數(shù)學(xué)家傅里葉在其1822年出版的《熱分析理論》一書中指出：任何周期信號都可以表示為不同頻率的正弦波和余弦波的組合，即所謂的傅里葉級數(shù)。隨著快速傅里葉變換算法的發(fā)明及數(shù)字計算機的出現(xiàn)，使得傅里葉變換得到了廣泛的應(yīng)用。傅里葉分析中最重要的結(jié)論就是幾乎“所有”的信號或函數(shù)都可以分解成簡單的正弦波和余弦波之和。從而提供了一種具有物理意義的函數(shù)表達方式。傅里葉變換的核心貢獻在于一、如何將每種正弦波和余弦波的比例(頻率)求出來。二、給定每種正弦波和余弦波的比例，可以恢復(fù)出原始信號。一種簡單的傅里葉變換如圖4-3所示。本節(jié)首先對一維傅里葉變換進行簡單介紹，并給出傅里葉變換的簡單實現(xiàn)。然后將一維傅里葉變換擴展到二維空間，并給出二維傅里葉變換的簡單實現(xiàn)。圖4-3一種簡單的傅里葉變換示意4.1.1一維傅里葉變換傅里葉變換中一般要求函數(shù)滿足迪利克雷條件(即在周期內(nèi)存在有限個間斷點)、有限極值條件、絕對可積條件(即)三個條件，則函數(shù)的傅里葉變換是存在的。一個函數(shù)的傅里葉變換可以表示為其對應(yīng)傅里葉逆變換表示為為其中傅里葉變換中基函數(shù)的物理意義非常明確，每一個基函數(shù)都是一個單頻率諧波，而對應(yīng)的系數(shù)(又稱頻譜)則表明了原函數(shù)在此基函數(shù)上投影的大小，或者也可以看作原函數(shù)中此種頻率諧波成分的比重。實際應(yīng)用中，需要求解的更多問題是離散信號的處理。定義離散情況的傅里葉變換公式：，其中，則其傅里葉正變換為傅里葉逆變換為觀察傅里葉逆變換，通過歐拉公式得可以看到，空域函數(shù)可表示為個正弦(余弦)函數(shù)的累加，其中F(u)/M為對應(yīng)頻率分量的幅度(系數(shù))。覆蓋的域(值)稱為頻率域。令=0，可得，對應(yīng)的均值，又可稱為直流分量。其余值對應(yīng)的則稱為的交流分量。通過變量可以確定變換后的頻率成分，而的取值范圍稱為頻率域。對每一個值，其對應(yīng)稱為傅里葉變換的頻率分量(或稱振幅)。我們可以注意到傅里葉變換后的函數(shù)是在復(fù)數(shù)域內(nèi),又可以表示為，或者以極坐標(biāo)的形式表示為。這里我們把量稱為傅里葉變換的幅度(Magnitude)或者譜(Spectrum)。通過譜可以表示原函數(shù)(或圖像)對某一頻譜分量的貢獻。而稱為變換的相位角或者相位譜,用來表示原函數(shù)中某一頻譜分量的起始位置。譜的平方稱為功率譜(有時也叫能量譜、譜密度)，表示為。下面通過示例講解傅里葉變換例：是門限函數(shù)，求其傅里葉變換。解：該門限函數(shù)及其傅里葉變換后函數(shù)圖像如圖4-4所示圖4-4門限函數(shù)及其傅里葉變換圖像下面給出一組函數(shù)及其傅里葉變換的示例代碼importmatplotlib.pyplotasplt

importnumpyasnp

defshow(ori_func,ft,sampling_period=5):

n=len(ori_func)

interval=sampling_period/n

#繪制原始函數(shù)

plt.subplot(2,1,1)

plt.plot(np.arange(0,sampling_period,interval),ori_func,'black')

plt.xlabel('Time'),plt.ylabel('Amplitude')

#繪制變換后的函數(shù)

plt.subplot(2,1,2)

frequency=np.arange(n/2)/(n*interval)

nfft=abs(ft[range(int(n/2))]/n)

plt.plot(frequency,nfft,'red')

plt.xlabel('Freq(Hz)'),plt.ylabel('Amp.Spectrum')

plt.show()

#生成頻率為1（角速度為2*pi）的正弦波

time=np.arange(0,5,.005)

x=np.sin(2*np.pi*1*time)

y=np.fft.fft(x)

show(x,y)結(jié)果如圖4-5所示圖4-5單一正弦波傅里葉變換結(jié)果繼續(xù)上述代碼，進行三個正弦波形疊加，然后進行傅里葉變換：#將其與頻率為20和60的波疊加起來

x2=np.sin(2*np.pi*20*time)

x3=np.sin(2*np.pi*60*time)

x+=x2+x3

y=np.fft.fft(x)

show(x,y)結(jié)果如圖4-6所示圖4-6三個正弦波形疊加的傅里葉變換下面進行方波傅里葉變換#生成方波，振幅是1，頻率為10Hz

#我們的間隔是0.05s，每秒有200個點

#所以需要每隔20個點設(shè)為1

x=np.zeros(len(time))

x[::20]=1

y=np.fft.fft(x)

show(x,y)結(jié)果如圖4-7所示圖4-7方波的傅里葉變換結(jié)果下面進行脈沖波的傅里葉變換#生成脈沖波

x=np.zeros(len(time))

x[380:400]=np.arange(0,1,.05)

x[400:420]=np.arange(1,0,-.05)

y=np.fft.fft(x)

show(x,y)圖4-8脈沖波的傅里葉變換結(jié)果4.1.2二維傅里葉變換二維傅里葉變換本質(zhì)上是一維傅里葉變換情形向二維進行簡單擴展對應(yīng)二維傅里葉變換的逆變換可以表示為離散情形完全與之類似。設(shè)是一幅尺寸為M×N的圖像函數(shù)，相應(yīng)的二維離散傅里葉變換可以表示為：均為頻率分量。通過傅里葉變換失去了空間關(guān)系，只保留了頻率關(guān)系。其中空間域是由所張成的坐標(biāo)系，和是變量。而頻率域則是由所張成的坐標(biāo)系，是變量。定義的矩形區(qū)域稱為頻率矩形，其大小與圖像的大小相同。是傅里葉系數(shù)。該圖像函數(shù)對應(yīng)傅里葉逆變換可以表示為針對正方形圖像，如果圖像是的圖像，則其離散傅里葉正變換可以簡化如下對應(yīng)傅里葉逆變換可簡化如下：使用歐拉公式對離散傅里葉公式進行展開可得進一步可以表示為其中又成為頻譜或幅度譜，譜圖像就是將作為亮度進行可視化。而稱為相位譜。下面給出二維傅里葉變換的對應(yīng)實現(xiàn)代碼及結(jié)果fromskimageimportdata

importnumpyasnp

frommatplotlibimportpyplotasplt

img=data.camera()

f=np.fft.fft2(img)#快速傅里葉變換算法得到頻率分布

fshift=np.fft.fftshift(f)#默認(rèn)結(jié)果中心點位置是在左上角，轉(zhuǎn)移到中間位置

fimg=np.log(np.abs(fshift))#fft結(jié)果是復(fù)數(shù)，求絕對值結(jié)果才是振幅

#展示結(jié)果

plt.subplot(121),plt.imshow(img,'gray'),plt.title('Original')

plt.subplot(122),plt.imshow(fimg,'gray'),plt.title('Fourier’)

plt.show()結(jié)果如圖4-9所示圖4-9一般圖像的傅里葉變換結(jié)果將代碼修改一下，看一下棋盤圖像對應(yīng)的傅里葉變換fromskimageimportdata

importnumpyasnp

frommatplotlibimportpyplotasplt

img=data.checkerboard()

f=np.fft.fft2(img)#快速傅里葉變換算法得到頻率分布

fshift=np.fft.fftshift(f)#默認(rèn)結(jié)果中心點位置是在左上角，轉(zhuǎn)移到中間位置

fimg=np.log(np.abs(fshift))#fft結(jié)果是復(fù)數(shù)，求絕對值結(jié)果才是振幅

#展示結(jié)果

plt.subplot(121),plt.imshow(img,'gray'),plt.title('Original')

plt.subplot(122),plt.imshow(fimg,'gray'),plt.title('Fourier')

plt.show()圖4-10棋盤圖像的傅里葉變換圖像經(jīng)傅里葉變換后，直流分流正比于圖像的均值，高頻分量則表明了圖像中目標(biāo)邊緣的強度及方向。下面給出二維傅里葉變換的極坐標(biāo)表示如下其中稱為頻率譜，稱為相位譜，而稱為功率譜。4.2二維傅里葉變換的性質(zhì)二維傅里葉變換也需要用到一維傅里葉變換的理論，所以也滿足一維傅里葉變換所滿足的基本性質(zhì)，同時也存在二維傅里葉變換的一些獨特性質(zhì)。下面就從兩個側(cè)面來描述二維傅里葉變換的性質(zhì)：1、傅里葉變換基本性質(zhì)。2、二維傅里葉變換性質(zhì)。4.2.1傅里葉變換基本性質(zhì)線性性質(zhì)：傅里葉變換的線性特性可以表示為若則其中為任意常數(shù)，利用傅氏變換的線性特性，可以將待求信號分解為若干基本信號之和。證：時延性：傅里葉變換的時延(移位)特性可以表示為若則時延（移位）性說明波形在時間軸上時延，并不會改變信號振幅頻譜，僅使信號增加線性相位。證：圖4-11時延移位對傅里葉頻譜的影響如圖4-11展示了時延移位對傅里葉頻譜的影響頻移性：傅里葉變換的頻移(調(diào)制)特性表示為若則頻移(調(diào)制)特性表明信號在時域中與復(fù)因子相乘，則在頻域中將使整個頻譜搬移證：尺度變換：傅里葉變換的尺度變換特性表示為若則尺度特性說明，信號在時域中壓縮，頻域中就擴展；反之，信號在時域中擴展，在頻域中就一定壓縮；即信號的脈寬與頻寬成反比。一般來說時寬有限的信號，其頻寬無限，反之亦然?？梢岳斫鉃樾盘柌ㄐ螇嚎s（擴展）倍，信號隨時間變化加快（慢）倍，所以信號所包含的頻率分量增加（減少）倍，頻譜展寬（壓縮）倍。又因能量守恒原理，各頻率分量分量的大小減?。ㄔ黾樱┍?。針對門限函數(shù)的尺度變換如圖4-12所示圖4-12針對門限函數(shù)的尺度變換及其傅里葉變換時域微分特性：傅里葉變換的時域微分特性表示為若則證：交換微、積分運算次序所以同理，可推廣到高階導(dǎo)數(shù)的傅里葉變換其中是微分因子。頻域微分特性：傅里葉變換的頻域微分特性表示為若則一般頻域微分特性的實用形式為對頻譜函數(shù)的高階導(dǎo)數(shù)亦成立或?qū)ΨQ（偶）性：傅里葉變換的對稱特性表示為若則或時域卷積定理：傅里葉變換的時域卷積定理表示為若則證：交換積分次序利用時延性由這個性質(zhì)，我們可將兩個時間函數(shù)的卷積運算變?yōu)閮蓚€頻譜函數(shù)的相乘(代數(shù))運算。由此我們可以用頻域法求解信號通過系統(tǒng)的響應(yīng)。頻域卷積定理：傅里葉變換的頻域卷積定理表示為若則證：交換積分次序利用移頻性通過該定理可以指導(dǎo)空間域的乘法操作對應(yīng)頻域卷積操作4.2.2二維傅里葉變換性質(zhì)可分離性：二維DFT可視為由沿方向的兩個一維傅里葉變換所構(gòu)成。這一性質(zhì)可以有效降低二維傅里葉變換的計算復(fù)雜性。如傅里葉逆變換也可以進行分離這樣原本在或平面需要時間復(fù)雜度才可以完成的操作，經(jīng)過分離之后可以通過方向的兩次時間復(fù)雜度為的一維傅里葉變換操作代替。逆變換同理。fromskimageimportdata,color

importnumpyasnp

frommatplotlibimportpyplotasplt

img=data.coffee()

img1=img

img=color.rgb2gray(img)

#在x方向?qū)崿F(xiàn)傅里葉變換

m,n=img.shape

fx=img

forxinrange(n):

fx[:,x]=np.fft.fft(img[:,x])

foryinrange(m):

fx[y,:]=np.fft.fft(img[y,:])

fshift=np.fft.fftshift(fx)#默認(rèn)結(jié)果中心點位置是在左上角，轉(zhuǎn)移到中間位置

fimg=np.log(np.abs(fshift))#fft結(jié)果是復(fù)數(shù)，求絕對值結(jié)果才是振幅

#展示結(jié)果

plt.subplot(121),plt.imshow(img1,'gray'),plt.title('Original')

plt.subplot(122),plt.imshow(fimg,'gray'),plt.title('Fourier')

plt.show()圖4-13使用兩次一維傅里葉變換代替二維傅里葉變換平移特性：二維傅里葉變換的平移特性表示如下在空間平移了，相當(dāng)于把傅里葉變換與一個指數(shù)相乘。在空間與一個指數(shù)項相乘相當(dāng)于平移其傅里葉變換。當(dāng)u0=M/2,v0=N/2時通常在變換前用乘以輸入圖像函數(shù)，實現(xiàn)頻域中心化變換：旋轉(zhuǎn)特性：對旋轉(zhuǎn)一定角度，相當(dāng)于將其傅里葉變換旋轉(zhuǎn)一定角度。如圖4-14所示，圖像旋轉(zhuǎn)一定角度之后，頻譜圖像也做出相應(yīng)旋轉(zhuǎn)。圖4-14傅里葉變換旋轉(zhuǎn)特性公式表述如下如通過以上性質(zhì)可以發(fā)現(xiàn)，在空間域所做的全局變化，也會反映到頻率域中。通過傅里葉變換可以在空間域和頻率域之間建立良好的關(guān)聯(lián)關(guān)系。4.3快速傅里葉變換離散傅立葉變換已成為數(shù)字信號處理的重要工具，然而，它的計算量大，運算時間長，在某種程度上卻限制了它的使用范圍?？焖偎惴ù蟠筇岣吡诉\算速度，在某些應(yīng)用場合已能作實時處理，并且應(yīng)用在控制系統(tǒng)中?？焖俑盗⑷~變換不是一種新的變換，它是離散傅立葉變換的一種算法，它是在分析離散傅立葉變換中的多余運算的基礎(chǔ)上，進而消除這些重復(fù)工作的思想指導(dǎo)下得到的。4.3.1快速傅里葉變換原理對于一個有限長序列，它的傅立葉變換由下式表示：令傅里葉變換可以表示為將進行展開可得從上式可以看出，要得到每一個頻率分量，需進行次乘法和次加法運算。要完成整個變換需要次乘法和次加法運算。當(dāng)序列較長時，必然要花費大量的時間。矩陣表示為其系數(shù)矩陣為觀察上面的系數(shù)矩陣，發(fā)現(xiàn)是以為周期的，即，當(dāng)=8時，其周期性如圖4-15所示圖4-15周期由于所以當(dāng)=8時，可得通過推導(dǎo)可以發(fā)現(xiàn)，離散傅立葉變換中的乘法運算有許多冗余操作。庫利-圖基于1965年提出了原始N點序列可以依次分解成一系列短序列，然后，求出這些短序列的離散傅立葉變換，以此來減少乘法運算，進一步降低計算復(fù)雜度。例設(shè)：由此，離散傅立葉變換可改寫成下面的形式：因為所以其中，和分別是和的N/2點的傅立葉變換。由于由于和均是以N/2為周期，所以這說明當(dāng)時，上式也是成立的，因此下式成立：由上面的分析可見，一個點的離散傅立葉變換可由兩個點的傅立葉變換得到。離散傅立葉變換的計算時間主要由乘法決定，分解后所需的乘法次數(shù)大為減少。第一項為次，第二項為＋次，總共為2×＋次運算即可完成，而原來需要次運算，可見分解后的乘法計算次數(shù)減少了近一半。當(dāng)為2的整數(shù)冪時，則上式中的和還可以再分成兩個更短的序列，因此計算時間會更短。由此可見，利用的周期性和分解運算，從而減少乘法運算次數(shù)是實現(xiàn)快速運算的關(guān)鍵。4.3.2快速傅里葉變換實現(xiàn)快速傅里葉變換的基本思想：快速傅里葉變換(FFT)基于逐次倍乘法(Successivedoublingmethod)。這個方法的主要思想是利用傅里葉變換(基底)的性質(zhì)，將個數(shù)據(jù)的傅里葉變換轉(zhuǎn)化為2組個數(shù)據(jù)的傅里葉變換。這樣，原來4的運算量就降低為2的運算量了。依次類推，便可得到快速算法。以一維為例作簡單介紹，二維情形可以通過兩次一維計算實現(xiàn)。將傅里葉變換改寫為其中僅考慮具有2的冪次方的形式，為正整數(shù)。故還可以表示為，也是正整數(shù)。從而注意：，,且，故有對，定義可知由于對任意的正整數(shù)，有以及，最后得到這樣就可以將原來計算比較復(fù)雜的傅里葉運算分解為兩個計算較為簡單傅里葉運算，且還可以繼續(xù)分解，如此循環(huán)推導(dǎo)下去，直到最后剩下若干組兩個點對。4.4圖像頻率域濾波圖像變換是對圖象信息進行變換，使能量保持但重新分配。以利于加工、處理（濾除不必要信息(如噪聲)，加強/提取感興趣的部分或特征）。傅里葉變換在圖像分析、濾波、增強、壓縮等處理中都有著非常重要的應(yīng)用。本節(jié)主要介紹基于傅里葉變換的圖像頻率域濾波。假定原圖像，經(jīng)傅立葉變換為，頻域增強就是選擇合適的濾波器函數(shù)對的頻譜成分進行調(diào)整，然后經(jīng)逆傅立葉變換得到增強的圖像。該過程可以通過下面流程描述：其中，=·，稱為傳遞函數(shù)或濾波器函數(shù)?？梢酝ㄟ^選擇合適的頻率傳遞函數(shù)來突出的某一方面的特征，從而得到需要的圖像。如:利用傳遞函數(shù)突出高頻分量，以增強圖像的邊緣信息，即高通濾波；如果突出低頻分量，就可以使圖像顯得比較平滑，即低通濾波。頻率域濾波基本步驟如下：一、對原始圖像進行傅里葉變換得到。二、將與傳遞函數(shù)進行卷積運算得到。三、將進行傅里葉逆變換得到增強圖。頻域濾波的核心在于如何確定傳遞函數(shù)，即。二維圖像頻域濾波的基本原理如圖4-16所示圖4-16二維圖像的頻率域濾波目前基于傅里葉變換的頻域濾波主要包括低通濾波、高通濾波、帶通濾波及同態(tài)濾波四類，后續(xù)將對低通濾波、高通濾波進行重點講述。4.4.1低通濾波圖像從空間域變換到頻率域后，其低頻分量對應(yīng)圖像中灰度值變化比較緩慢的區(qū)域，高頻分量則表征圖像中物體的邊緣和隨機噪聲等信息。低通濾波是指保留低頻分量，而通過濾波器函數(shù)減弱或抑制高頻分量的在頻域進行的濾波。低通濾波與空域中的平滑濾波器一樣可以消除圖像中的隨機噪聲，減弱邊緣效應(yīng)，起到平滑圖像的作用。下面介紹幾種常用的頻率域低通濾波器：一、理想低通濾波器二維理想低通濾波器的傳遞函數(shù)如下函數(shù)圖像如圖4-17所示圖4-17理想低通濾波器及其圖像截斷頻率是一個非負(fù)整數(shù)，是從點到頻率平面原點的距離即。理想低通濾波器的含義是指小于的頻率，即以為半徑的圓內(nèi)的所有頻率分量可以完全無損地通過，而圓外的頻率，即大于的頻率分量則完全被除掉。理想低通濾波器的平滑作用非常明顯，但由于變換有一個陡峭的波形，它的反變換有強烈的振鈴特性，使濾波后圖像產(chǎn)生模糊效果。因此這種理想低通濾波實用中并不采用。典型及其逆變換圖像如圖4-18所示圖4-18理想低通濾波器的傅里葉逆變換理想低通濾波器相關(guān)代碼及結(jié)果如下#導(dǎo)入相關(guān)庫

fromskimageimportdata,color

importnumpyasnp

importmatplotlib.pyplotasplt

D=10

#讀入圖片

new_img=data.coffee()

new_img=color.rgb2gray(new_img)

#numpy中的傅立葉變化

f1=np.fft.fft2(new_img)

f1_shift=np.fft.fftshift(f1)

#np.fft.fftshift()函數(shù)來實現(xiàn)平移,讓直流分量在輸出圖像的重心

#實現(xiàn)理想低通濾波器

rows,cols=new_img.shape

crow,ccol=int(rows/2),int(cols/2)#計算頻譜中心

mask=np.zeros((rows,cols),np.uint8)#生成rows行cols的矩陣，數(shù)據(jù)格式為uint8

foriinrange(rows):

forjinrange(cols):

ifnp.sqrt(i*i+j*j)<=D:

#將距離頻譜中心小于D的部分低通信息設(shè)置為1，屬于低通濾波

mask[crow-D:crow+D,ccol-D:ccol+D]=1

f1_shift=f1_shift*mask

#傅立葉逆變換

f_ishift=np.fft.ifftshift(f1_shift)

img_back=np.fft.ifft2(f_ishift)

img_back=np.abs(img_back)

img_back=(img_back-np.amin(img_back))/(np.amax(img_back)-np.amin(img_back))

plt.figure(figsize=(15,8))

plt.subplot(121),plt.imshow(new_img,cmap='gray'),plt.title('inputimage')

plt.subplot(122),plt.imshow(img_back,cmap='gray'),plt.title('outputimage')

plt.show()圖4-19二維圖像的理想低通濾波可以發(fā)現(xiàn)當(dāng)距離設(shè)置為10時，出現(xiàn)了較為明顯的振鈴現(xiàn)象，讀者可以嘗試對賦予其它不同值，看一下效果。低通濾波的能量和的關(guān)系，能量在變換域中集中在低頻區(qū)域。以理想低通濾波作用于的數(shù)字圖像為例，其總能量當(dāng)理想低通濾波的變化時，通過的能量和總能量比值必然與有關(guān)，而可表示的通過能量百分?jǐn)?shù)。一個以頻域中心為原點，為半徑的圓就包含了百分之的能量根據(jù)對保留能量的要求來確定濾波器的截止頻率，半徑與包含能量的關(guān)系如表4-1所示半徑包含總能量%590.01196.02298.03699.05399.59899.9表4-1濾波器半徑與包含能量之間的關(guān)系理想低通濾波器在數(shù)學(xué)上定義得很清楚，在計算機模擬中也可實現(xiàn)，但在截斷頻率處直上直下的理想低通濾波器不能用實際的電子器件實現(xiàn)，物理上可實現(xiàn)的是巴特沃斯低通濾波器。巴特沃斯（Butterworth）低通濾波器：巴特沃斯低通濾波器的傳遞函數(shù)為為截止頻率，為函數(shù)的階。一般取使最大值下降到最大值的0.5時的為截止頻率。Butterworth低通濾波器的截面如圖4-20所示圖4-20Butterworth低通濾波器的截面與理想低通濾波器相比，高低頻之間過渡較為平滑，用此濾波后的輸出圖像振鈴現(xiàn)象不明顯。=1時，過渡最為平滑，即尾部包含有大量的高頻成分，所以1階巴特沃斯低通濾波器沒有振鈴現(xiàn)象。但隨著的增加，振鈴現(xiàn)象會越來越明顯。下面給出巴特沃斯低通濾波器代碼及結(jié)果importnumpyasnp

importmatplotlib.pyplotasplt

importnumpyasnp

fromskimageimportdata,color

img=data.coffee()

img=color.rgb2gray(img)#直接讀為灰度圖像

f=np.fft.fft2(img)

fshift=np.fft.fftshift(f)

#取絕對值：將復(fù)數(shù)變化成實數(shù)

#取對數(shù)的目的為了將數(shù)據(jù)變化到0-255

s1=np.log(np.abs(fshift))

"""

巴特沃斯低通濾波器

"""

defbutterworthPassFilter(image,d,n):

f=np.fft.fft2(image)

fshift=np.fft.fftshift(f)

defmake_transform_matrix(d):

transfor_matrix=np.zeros(image.shape)

center_point=tuple(map(lambdax:(x-1)/2,s1.shape))

foriinrange(transfor_matrix.shape[0]):

forjinrange(transfor_matrix.shape[1]):

defcal_distance(pa,pb):

frommathimportsqrt

dis=sqrt((pa[0]-pb[0])**2+(pa[1]-pb[1])**2)

returndis

dis=cal_distance(center_point,(i,j))

transfor_matrix[i,j]==1/(1+(dis/d)**(2*n))

returntransfor_matrix

d_matrix=make_transform_matrix(d)

new_img=np.abs(np.fft.ifft2(np.fft.ifftshift(fshift*d_matrix)))

returnnew_img

plt.subplot(221)

plt.axis("off")

plt.title('Original')

plt.imshow(img,cmap='gray')

plt.subplot(222)

plt.axis('off')

plt.title('Butter1001')

butter_100_1=butterworthPassFilter(img,100,1)

plt.imshow(butter_100_1,cmap='gray')

plt.subplot(223)

plt.axis('off')

plt.title('Butter301')

butter_30_1=butterworthPassFilter(img,30,1)

plt.imshow(butter_30_1,cmap='gray')

plt.subplot(224)

plt.axis('off')

plt.title('Butter305')

butter_30_5=butterworthPassFilter(img,30,5)

plt.imshow(butter_30_5,cmap='gray')

plt.show()圖4-21二維圖像的Butterworth低通濾波結(jié)果4.4.2高通濾波圖像的邊緣、細(xì)節(jié)主要在高頻，圖像模糊是由于高頻成分較弱產(chǎn)生的。為了消除模糊，突出邊緣，可以采用高通濾波的方法，使低頻分量得到抑制，從而達到增強高頻分量，使圖像的邊沿或線條變得清晰，實現(xiàn)圖像的銳化。理想高通濾波器：形狀與低通濾波器的形狀正好相反，其傳遞函數(shù)為圖4-22理想高通濾波器函數(shù)圖像如圖4-22所示，代碼與理想低通濾波器相似，將低通濾波器部分代碼修改如下#實現(xiàn)理想高通濾波器

rows,cols=new_img.shape

crow,ccol=int(rows/2),int(cols/2)#計算頻譜中心

mask=np.zeros((rows,cols),np.uint8)#生成rows行cols的矩陣，數(shù)據(jù)格式為uint8

foriinrange(rows):

forjinrange(cols):

ifnp.sqrt(i*i+j*j)<=D:

#將距離頻譜中心小于D的部分低通信息設(shè)置為1，屬于低通濾波

mask[crow-D:crow+D,ccol-D:ccol+D]=1

mask=1-mask

f1_shift=f1_shift*mask圖4-22二維圖像的理想高通濾波巴特沃斯高通濾波器：形狀與巴特沃斯低通濾波器的形狀正好相反，同樣的，因為高低頻率間平滑過渡，振鈴現(xiàn)象不明顯。其傳遞函數(shù)如下代碼實現(xiàn)只需要將巴特沃斯低通濾波器代碼做少量修改，需修改部分及結(jié)果如下"""

巴特沃斯高通濾波器

"""

defbutterworthPassFilter(image,d,n):

f=np.fft.fft2(image)

fshift=np.fft.fftshift(f)

defmake_transform_matrix(d):

transfor_matrix=np.zeros(image.shape)

center_point=tuple(map(lambdax:(x-1)/2,s1.shape))

foriinrange(transfor_matrix.shape[0]):

forjinrange(transfor_matrix.shape[1]):

defcal_distance(pa,pb):

frommathimportsqrt

dis=sqrt((pa[0]-pb[0])**2+(pa[1]-pb[1])**2)

returndis

dis=cal_distance(center_point,(i,j))

transfor_matrix[i,j]=1/(1+(dis/d)**(2*n))

returnt