我們知道，給定一個空（時）間濾波器，既可以在空（時）域直接完成數字信號的濾波，也可以在頻域完成。空域濾波的數學運算為卷積/相關，對應頻域則為點乘/頻域數據的共軛（G*）與濾波器(H)的乘積。

簡單起見，以一個1D數字信號為例加以說明。

例如：給定信號x = [1 2 2 4 4]，濾波器h = [1 2 3]；

空域濾波（卷積）為：

y = conv(x,h,’same’); % 此處未考慮h系數的歸一化問題

則濾波結果為：

y = [4 9 14 18 20]

下面介紹頻域方式的濾波過程。由于時域濾波屬于有限序列的線性卷積，頻域濾波方式實際上是利用離散傅里葉變換（DFT）求時域線性卷積的過程，而DFT本質上是對應時域濾波中針對周期序列的循環卷積。要滿足循環卷積與線性卷積計算結果一致，時域信號x(m點)與濾波器h(n點)必須等長。這可以采用補零的方法，使x與h的長度均為L≥m+n–1。這樣做的目的，也是為了避免周期函數卷積中因周期靠近引起所謂的頻率纏繞錯誤（混疊）。

一般做法是，對空（時）域數據采用后端補0方式延拓至2倍數據長度。

信號延拓： xp = [1 2 2 4 4 0 0 0 0 0];

由于濾波器h與延拓后的信號xp應保持相同的長度，故也需要補0做延拓。有兩種方式：

(1) 雙邊延拓，濾波器h居中

hp = [0 0 0 0 1 2 3 0 0 0];

(2) 單邊延拓，濾波器h居左

hp = [1 2 3 0 0 0 0 0 0 0];

然后，分別對xp,hp做DFT，完成頻域濾波并做反變換，即

y = real(ifft(fft(xp).*fft(hp)));

最后，剪裁掉補0多出的后半部分數據，即保留主值序列。則兩種方式得到的濾波結果分別為：

y1 = y(1:5) = [12 0 0 0 1]; 和 y2 = y(1:5) = [1 4 9 14 18];

可以看出，第一種延拓，濾波結果與空域結果完全不一致。第二種延拓與空域結果基本一致，但在左右短點處（邊界部分）是不一致的。這里，第二種延拓也等價于：

y = real(ifft(fft(x,10).*fft(h,10)));

從傅立葉變換的時域性質知道，兩種延拓的頻譜是一樣的，但相位會發生變化。因此，第一種延拓（兩端補零方式，h居中）是不可取的。第二種延拓與空域濾波結果基本一致，但邊界上有差異。那么，如何消除這種邊界差異，達到與空域考濾波完全一致的結果呢？

我們只要在濾波器h延拓方式上稍作改動，即把空域濾波器中心元素放到最前段（起始點），左端被擠出的元素順序放在尾部，即所謂的循環移位（circularly shift）法，則有：

hp = [2 3 0 0 0 0 0 0 0 1];

再按以上的相同步驟進行濾波處理，可得：

y = [4 9 14 18 20];

頻域濾波與空域濾波的結果就可完全保持一致。

假如，空域做的是相關運算來完成濾波，那么，只要濾波器旋轉180度后，采取同樣的0填充方式。

顯然，以上方式很容易推廣到二維情況。如果有以下二維濾波器，即一個計算y方向梯度的Sobel算子。

由于圖像的模板運算默認為相關運算，而頻域的乘積對應空域卷積。要達到一致性，以上濾波器需要上下顛倒，即旋轉180度后(若為對稱濾波器，可省此步)，即

以上模板在頻域進行點乘可對應h的相關運算。

若需要填充0方式擴大到10×10（根據濾波圖像的尺寸而定），則填充方式為：

以上兩種方式中，hp的填充方式會出現與空域濾波結果邊界上的不一致，而hp’填充方式則可以解決這個問題。

道理很簡單，因為空域模板運算的當前像素（原點）一般是在濾波模板的中心像素，即

頻域乘積（點乘）對應空域卷積的情況下，是按以下公式：

即原點（見橢圓標記處）是從左上角開始的。

下面的MATLAB代碼可以簡單實現以上方法，保持h中心像素位于填充區域的左上角。

center_h = ceil((size(h) + 1)/2); % 確定空域濾波器h中心像數坐標hp = zeros(P, Q); % 生成P×Q的全零矩陣hp(1:size(h,1), 1:size(h,2)) = h % 左上角填充h，形成延拓濾波器hprow_indices = [center_h(1): P, 1: (center_h(1)-1)]';col_indices = [center_h(2): Q, 1: (center_h(2)-1)]';hp = hp(row_indices, col_indices); % 原點在左上角的延拓濾波器hp其中，P, Q是擴充后的濾波器尺寸，與待處理圖像（M×N）有關（一般取，P = 2M,Q = 2N）。以上代碼也就是實現前面提到的濾波器h補零延拓+循環移位過程。

濾波器h補零延拓做FFT，雖然可以簡單調用Hp = fft2(h,P,Q)來完成（MATLAB內部處理方式），但并未做循環移位，僅僅是右下部補零后計算FFT。因此，濾波結果的邊界處并不能保證與空域濾波結果的一致性。空域濾波器h尺寸越大，這種邊界差異越明顯。例如用25×25，方差為2的空域高斯低通波器進行實驗，結果如下。


上一排圖分別為原圖，空域濾波結果及模板置于左上角的頻域濾波結果。可以看出，頻域濾波結果有明顯的邊界效應。下圖為模版中心像素置左上角的測試結果（右下圖），則與空域濾波結果是完全一致

二

前文（如何保持空域與頻域濾波結果的一致性）中談到，從空域小模板轉化到頻域濾波時，為了保持空頻濾波的一致性，需要合理處理空域濾波器的延拓和布局問題。一般做法是，將小模板的中心通過循環移位后置于補零延拓矩陣的左上角。然后做傅里葉變換，得到對應的頻域濾波器，再與延拓后的圖像在頻率域與濾波器做乘法運算。具體的原理和方法可參看前述博文。

設有一幅M×N的圖像f(x,y)，m×n的空域濾波器為h(x,y)，則頻域濾波的處理步驟如下：

(1) 消除折疊現象的填充（Zero padding）。即分別對f(x,y)，h(x,y)的右下部補零至P×Q得到fp(x,y)和hp (x,y)，其中h(x,y)需要做循環移位，以使小模板h(x,y)的中心像素置于hp (x,y)的左上角。一般取：P=2M，Q=2N 。

(2) fp(x,y)，hp (x,y)分別做傅里葉變換產生Fp(u,v)，Hp(u,v)。

(3) 中心變換（頻譜中心化）。此步也可以不變換，則Hp(u,v)要改變(針對直接在頻域生成對稱濾波器情況)。

(4) 頻域濾波：Hp(u,v)點乘Fp(u,v)。

(5) 傅里葉反變換。

(6) 取實數部分。絕對值很小的虛數部分是浮點運算存在誤差造成的。

(7) 空域中心還原變換（反中心化）。若Fp(u,v)未做中心化，此步可省。

(8) 截取有效數據，即左上角的原始圖像尺寸M×N部分數據。

以上步驟的濾波，僅限于空域濾波的邊界處理為零填充方式。如果空域濾波的邊界處理為其他方式，如對稱邊界（’symmetric’）重復邊界（’replicate’）和周期邊界（’circular’）等，則依然會存在空頻域濾波結果在邊界上的差異性。如圖1所示，是一個方差為4的25×25高斯低通濾波器對cameraman圖像分別在空域和頻域濾波結果的對比。可以看出，頻域濾波與空域濾波在邊界上是不一致的。

圖1 僅在空域考慮了邊界因素的濾波結果

那么，如何有效解決這個問題呢？其實也很簡單，只要在步驟(1)和(8)上稍稍改進，就可以保持空/頻域濾波結果邊界上的一致性。

第(1)步，根據h(x,y)的尺寸對f(x,y)先做重復邊界的擴充，在此基礎上做消除折疊的補零延拓，即得到擴大至（2M+行重復邊界數）×（2N+列重復邊界數）的fp(x,y),如圖2所示。對h(x,y)右下部補零至與fp(x,y)的相同尺寸，并循環移位后得到hp (x,y)。

圖2 濾波前圖像邊界擴充及補零延拓

第(8)步，截取有效數據時，應除去左上角單邊邊界數后的原始圖像尺寸（M×N）部分數據(即圖2中的紅框區域)。其他步驟不變，即可得到與空域濾波完全一致的結果，如圖3所示。

MATLAB參考代碼如下：

%=============================================================================inimg = imread('cameraman.tif');subplot(131)imshow(inimg), title('Original image')[M,N] = size(inimg); % Original image size%====================================================================h = fspecial('gaussian',25,4); % Gaussian filter%====================================================================% 空域濾波gx = imfilter(inimg,h,'same','replicate'); % 空域圖像濾波subplot(132)imshow(gx,[]);title('Spatial domain filtering')%====================================================================% 頻域濾波%====================================================================h_hf = floor(size(h)/2); % 空域濾波器半高/寬imgp = padarray(inimg, [h_hf(1),h_hf(2)],'replicate'); % Padding boundary with copying pixels% PQ = paddedsize(size(imgp)); % Gonzalez DIP教材提供的函數，非MATLAB內部函數PQ = 2*size(imgp);Fp = fft2(double(imgp), PQ(1), PQ(2)); % 延拓圖像FFT% h = rot90(h,2); % Mask旋轉180度，非對稱h需此步驟！因頻域乘積對應空域卷積，而空域濾波為相關。P = PQ(1); Q = PQ(2);center_h = h_hf+1; % 空域小模板h中心位置hp = zeros(P,Q); % 預分配內存，產生P×Q零矩陣hp(1:size(h,1),1:size(h,2)) = h; % h置于hp左上角hp = circshift(hp,[-(center_h(1)-1),-(center_h(2)-1)]); % 循環移位，h中心置于hp左上角%====================================================================Hp = fft2(double(hp)); % hp濾波器做FFT%====================================================================Gp = Hp.*Fp; % 頻域濾波gp = real(ifft2(Gp)); % 反變換，取實部gf = gp(h_hf(1)+1:M+ h_hf(1), h_hf(2)+1:N + h_hf(2)); % 截取有效數據subplot(133)imshow(uint8(gf),[]), title('Frequency domain filtering')% 注：以上處理中，頻域圖像Fp與濾波器Hp均未中心化，因此，返回空域時無需反中心化。% 另外，直接調用Hp = freqz2(h,P,Q)獲得的2D頻域響應，則是中心化的。

總結

以上是生活随笔為你收集整理的如何保持空域与频域滤波结果的一致性的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。