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

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

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

空域濾波（卷積）為：

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

則濾波結(jié)果為：

y = [4 9 14 18 20]

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

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

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

由于濾波器h與延拓后的信號xp應(yīng)保持相同的長度，故也需要補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多出的后半部分?jǐn)?shù)據(jù)，即保留主值序列。則兩種方式得到的濾波結(jié)果分別為：

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

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

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

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

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

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

再按以上的相同步驟進(jìn)行濾波處理，可得：

y = [4 9 14 18 20];

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

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

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

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

以上模板在頻域進(jìn)行點乘可對應(yīng)h的相關(guān)運算。

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

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

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

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

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

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

center_h = ceil((size(h) + 1)/2); % 確定空域濾波器h中心像數(shù)坐標(biāo)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）有關(guān)（一般取，P = 2M,Q = 2N）。以上代碼也就是實現(xiàn)前面提到的濾波器h補零延拓+循環(huán)移位過程。

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


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

二

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

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

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

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

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

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

(5) 傅里葉反變換。

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

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

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

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

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

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

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

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

第(8)步，截取有效數(shù)據(jù)時，應(yīng)除去左上角單邊邊界數(shù)后的原始圖像尺寸（M×N）部分?jǐn)?shù)據(jù)(即圖2中的紅框區(qū)域)。其他步驟不變，即可得到與空域濾波完全一致的結(jié)果，如圖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教材提供的函數(shù)，非MATLAB內(nèi)部函數(shù)PQ = 2*size(imgp);Fp = fft2(double(imgp), PQ(1), PQ(2)); % 延拓圖像FFT% h = rot90(h,2); % Mask旋轉(zhuǎn)180度，非對稱h需此步驟！因頻域乘積對應(yīng)空域卷積，而空域濾波為相關(guān)。P = PQ(1); Q = PQ(2);center_h = h_hf+1; % 空域小模板h中心位置hp = zeros(P,Q); % 預(yù)分配內(nèi)存，產(chǎn)生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)]); % 循環(huán)移位，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)); % 截取有效數(shù)據(jù)subplot(133)imshow(uint8(gf),[]), title('Frequency domain filtering')% 注：以上處理中，頻域圖像Fp與濾波器Hp均未中心化，因此，返回空域時無需反中心化。% 另外，直接調(diào)用Hp = freqz2(h,P,Q)獲得的2D頻域響應(yīng)，則是中心化的。

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。