頻域濾波步驟:
1.給定大小為 M × N M\times N M×N的輸入圖像,将其填充為大小 P × Q P\times Q P×Q的圖像其中 P = 2 M , Q = 2 N P=2M,Q=2N P=2M,Q=2N。填充方法為:在M,N的後面添加0(尾部加0)。設經過此步驟的圖像為 f p ( x , y ) f_p(x,y) fp(x,y)
進行第一步的原因(下面卷積公式應為f(m)):
![](https://img.laitimes.com/img/9ZDMuAjOiMmIsIjOiQnIsIyZuBnLmVGZmRTZ0MDM5IGOlVGO5gTMmRzMjJTN0YGOhBjMwU2Lc52YucWbp5GZzNmLn9Gbi1yZtl2Lc9CX6MHc0RHaiojIsJye.png)
如圖,相當于對h(m)關于y軸對稱之後向右不斷平移x,并求其和f(m)的乘積。如果對稱之後,h(m)的圖像與h(x) or f(x)有重疊,就會産生混疊,影響結果。
2.對圖像× ( − 1 ) x + y (-1)^{x+y} (−1)x+y,将其移到變換中心
進行第二步的原因:
如圖,進行運算之後在 [ 0 , M − 1 ] [0,M-1] [0,M−1]内是完整的一個周期,便于DFT
3.計算步驟二之後獲得的圖像(矩陣)的DFT
就直接套公式就完事了
F ( u , v ) = ∑ x = 0 M − 1 ∑ y = 0 N − 1 f ( x , y ) e − 2 π j ( u x M + v y N ) F(u,v)=\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)e^{-2\pi j(\frac{ux}{M}+\frac{vy}{N})} F(u,v)=∑x=0M−1∑y=0N−1f(x,y)e−2πj(Mux+Nvy)
4.生成一個實的,對稱的濾波函數H(u,v),其大小為P×Q,中心在(P/2,Q/2)處,G(u,v) = H(u,v).*F(u,v)
要注意濾波函數的生成大小,低通濾波器可以用
%D(u,v)<=D0時有值
u = 0:Q-1
v = 0:P-1
%calculate D(u,v)
H = double(D<=D0)
控制濾波器的大小
5.對G進行傅裡葉反變換,得到g(x,y),并對其取實部,乘(-1)^(x+y)
其實取abs實際感覺差别不是很大
6.從g(x,y)上提取M*N的部分得到最終結果。
下面是代碼實作
%matlab調庫
function [DFT_img] = FFT_(I,H)
%調用FFT庫函數實作頻域濾波
%輸入 I為圖像f(x,y),H為濾波器H(u,v)
%輸出濾波後的圖像g(x,y)
%在不嚴謹的情況下我們可以不做擴充0,結果仍然與正确結果差不多
F = fft2(I);
Fs = fftshift(F);
out = H.*Fs;
out=ifftshift(out);
out=ifft2(out);
out=real(out);%用abs基本沒差別
out=out/max(out(:));%除以out中最大值
DFT_img = out;
end
%以低通濾波為例
path = 'your own path';
I = imread(path);
M = size(I,2);
N = size(I,1);
u=0:M-1;
v=0:N-1;
[U,V]=meshgrid(u,v);%基于向量 x 和 y 中包含的坐标傳回二維網格坐标
D=sqrt((U-M/2).^2+(V-N/2).^2);
D0=15;%截止頻率
H=double(D<=D0);%理想低通濾波器
n = FFT_(I,H)
imshow(n,[])%不加[]效果相同
原圖:
低通濾波以後(D0=15):
若D0=5,可以看到更為模糊的圖像:
matlab庫函數解析:
-
Y = fft2(X)——輸出X的二維傅裡葉變換
Y = fft2(X,m,n) 将截斷 X 或用尾随零填充 X,以便在計算變換之前形成 m×n 矩陣
-
fftshift(X)将零頻分量移到螢幕中心,相當于我們上述講的 × ( − 1 ) x + y \times(-1)^{x+y} ×(−1)x+y
實際上fftshift操作的時候就是把矩陣X第一象限與第三象限交換,第二象限與第四象限交換。(相當于拆開再拼起來)
- ifftshift(X) 撤銷fftshift(X)的結果,實際還是第一象限與第三象限交換,第二象限與第四象限交換。
- ifft2(X) 二維傅裡葉逆變換