Digital Image processing 數字影象處理最佳陷波濾波器設計（頻域）
最佳陷波濾波器（傅立葉變換，matlab實現）
1基本概念
1.1頻域濾波步驟小結
在頻域中的濾波是簡單明瞭的。它包含如下步驟：
（1）給定一幅大小為M*N的輸入影象f(x,y),選擇填充引數P，Q，典型地，我們選擇P=2M和Q=2N。
（2）對f(x,y)新增必要數量的0，形成大小為P*Q的填充後的影象fp(x,y)。
（3）用(-1)x+y乘以fp(x,y)移到其變換的中心。
（4）計算來自步驟3的影象DFT，得到F(u,v)。
（5）生成一個實的、對稱的濾波函式H(u,v),其大小為P*Q，中心在(P/2,Q/2)處。用陣列相乘形成乘積G(u,v)=H(u,v)F(u,v);
（6）得到處理後的影象：
（7）通過從gp(x,y)的左上象限提取M*N區域，得到最終處理結果g(x,y)。
參照課本p264進行實驗，如圖1所示，詳細可見程式filtPro.m。
1.2帶阻濾波器
理想帶阻濾波器的表示式為：
其中，D(u,v)是到中心化頻率矩形原點的距離，W是頻帶的寬度，D0是頻帶的中心半徑。
同樣，n階的巴特沃斯帶阻濾波器的表示式為：
高斯帶阻濾波器的表示式為：
詳細可見程式lowpfilter.m,參考課本p336利用帶阻濾波器消除週期性噪聲，如圖4所示，具體程式碼可見filtCircle.m.
1.3帶通濾波器
帶通濾波器執行與帶阻濾波器相反的操作。帶通濾波器的傳遞函式Hbp(u,v)是根據相應的帶阻濾波器的傳遞函式Hbr(u,v)並應用下式：
1.4陷波濾波器
陷波濾波器阻止（或通過）事先定義的中心頻率領域內的頻率。假定頻率矩形的中心已經移動到點（M/2，N/2），（u0，v0）的值對應移動中心，則有
程式notchfilter.m為帶阻陷波濾波器，分別為理想型、巴特沃斯、高斯帶阻陷波濾波器，notchfilter_1.m為帶阻陷波濾波器對。
2.演算法思想及實驗結果分析
(1)第一步，提取干擾模式的主頻率成分。
通過觀察傅立葉變換後的頻譜圖，也叫功率圖，我們首先就可以看出，影象的能量分佈，如果頻譜圖中暗的點數更多，那麼實際影象是比較柔和的（因為各點與鄰域差異都不大，梯度相對較小），反之，如果頻譜圖中亮的點數多，那麼實際影象一定是尖銳的，邊界分明且邊界兩邊畫素差異較大的。對頻譜移頻到原點以後，可以看出影象的頻率分佈是以原點為圓心，對稱分佈的。將頻譜移頻到圓心除了可以清晰地看出影象頻率分佈以外，還有一個好處，它可以分離出有周期性規律的干擾訊號，比如正弦干擾，一副帶有正弦干擾，移頻到原點的頻譜圖上可以看出除了中心以外還存在以某一點為中心，對稱分佈的亮點集合，這個集合就是干擾噪音產生的，這時可以很直觀的通過在該位置放置帶阻濾波器消除干擾。
如圖5為“水星6號”飛船拍攝的活性地形影象，圖6為週期干擾的傅立葉頻譜（由美國國家航空航天局提供），觀察頻譜，在每個尖峰處設一陷波帶通濾波器H(u,v),即設定H(u,v)只可通過與干擾模式相關的成分，如圖7所示，中間使用理想型高通濾波器和“+”型帶阻濾波器，四角使用理想型帶阻陷波濾波器。通過下式計算干擾噪聲模式的傅立葉變換：
其中，G(u,v)為被汙染影象的傅立葉變換。
然後通過傅立葉逆變換獲得通過的噪聲干擾，如圖8所示。
圖5“水星6號”飛船拍攝的活性地形影象
圖6 週期干擾的傅立葉頻譜（由美國國家航空航天局提供）
圖7 疊加在圖6的陷波帶通濾波器
圖8 濾波後圖像的傅立葉逆變換，在空間域顯示噪聲模式
（2）從受汙染的影象g(x,y)減去模式得到f(x,y)是很簡單的，但是這個濾波過程通常只會得到真實幹擾模式的近似值。在噪聲
其中，
此過程的目的就是選取加權函式，使得估計值在每一點(x,y)的指定領域的方差最小。考慮點(x,y)的尺寸為a*b的領域，假設w(x,y)在整個領域內保持不變（詳細參考書p200），則w(x,y)解為：
利用上述兩個式子來獲得復原影象，詳細可見程式optimum_notch.m。如果w(x,y)在某一領域內被假定為常量，則不必再像法1一樣在影象中計算對應每對x,y值的函式值，而是在每個非重疊領域的一點（一般為中心點）計算w(x,y)值，然後處理該領域內包含的所有影象點，如法2，3，這樣演算法效率會搞一點。對於邊界處理，可採用對稱映象填充法、零填充法等，這裡採用了對稱映象填充法。
（實驗分別使用理想濾波器、高斯濾波器和巴特沃斯濾波器進行濾波比較，以及a和b設定的相關實驗，具體參看github）
實際天文影象應用成果
原圖（雲南天文臺拍攝）
傅立葉譜
濾波
噪聲
本Markdown編輯器使用[StackEdit][6]修改而來，用它寫部落格，將會帶來全新的體驗哦：
