頻域濾波基礎

1、頻域濾波與空域濾波的關系

傅立葉變換可以將影像從空域變換到頻域，而傅立葉反變換則可以將影像的頻譜逆變換為空域影像，這樣一來，我們可以利用空域影像與頻域之間的對應關系，嘗試將空域卷積濾波變換為頻域濾波，而后再將頻域濾波處理后的影像再反變換回空間域，從而達到影像增強的目的，這樣做的一個最主要的吸引力在于頻域濾波的直觀性特點，

根據著名的卷積定理，兩個二維連續函式在空間域中的卷積可由其相應的兩個傅立葉變換乘積的反變換而得到；反之，在頻域中的卷積可由在空間域中乘積的傅立葉變換而得到，即：
f ( x , y ) ? h ( x , y ) ?? ? ?? F ( u , v ) H ( u , v ) f(x,y)*h(x,y) \iff F(u,v)H(u,v) f(x,y)?h(x,y)?F(u,v)H(u,v)

f ( x , y ) h ( x , y ) ?? ? ?? F ( u , v ) ? H ( u , v ) f(x,y)h(x,y) \iff F(u,v)*H(u,v) f(x,y)h(x,y)?F(u,v)?H(u,v)

其中，F(u,v) 和 H(u,v) 分別表示 f(x,y) 和 h(x,y) 的傅立葉變換，而符號表示傅立葉變換對，

式(1) 構成了整個頻域濾波的基礎，

2、頻域濾波的基本步驟

根據式(1) 進行頻域濾波通常應遵循以下步驟：

計算原始影像 f(x,y) 的 DFT，得到頻譜 F(u,v)；
將頻譜 F(u,v) 的零頻點移動到頻譜圖的中心位置；
計算濾波器函式 H(u,v) 與 F(u,v) 的乘積 G(u,v)；
將頻譜 G(u,v) 的零頻點移回到頻譜圖的左上角位置；
計算第 (4) 步計算結果的傅立葉反變換 g(x,y)；
取 g(x,y) 的實部作為最終濾波后的結果影像，

由上面敘述可知，濾波能否取得理想的關鍵取決于頻域濾波函式 H(u,v)，常常稱之為濾波器，或濾波器傳遞函式，因為它在濾波中抑制或濾除了頻譜中某些頻率的分量，而保留其他的一些頻率不受影響，在這里只關心其值為實數的濾波器，這樣濾波程序中 H 的每一個實數元素分別乘以 F 中對于位置的復數元素，從而使 F 中元素的實部和虛部等比例的變化，不會改變 F 的相位譜，這種濾波器也因此稱為 “零相位” 濾波器，這樣，最終反變換回空間域得到的濾波結果影像 g(x,y) 理論上也應當為實函式，然而由于計算舍入誤差等原因，可額能會有非常小的虛部，通常將虛部直接忽略，

3、頻域濾波的 Matlab 實作

撰寫了一個 imfreqfilt 函式：

function out = imfreqfilt(I,ff)
%imfreqfilt函式     對灰度影像進行頻域濾波
%引數I              輸入的空域影像
%引數ff             應用的與原影像等大的頻域濾鏡

if (ndims(I)==3) && (size(I,3)==3) % RGB影像
    I = rgb2gray(I);
end

if (size(I) ~= size(ff))
    msg1 = sprintf('%s: 濾鏡與原影像不等大，檢查輸入', mfilename);
    msg2 = sprintf('%s: 濾波操作已經取消', mfilename);
    eid = sprintf('Images:%s:ImageSizeNotEqual',mfilename);
    error(eid,'%s %s',msg1,msg2);
end

% 快速傅立葉變換
f = fft2(double(I));

% 移動原點
s = fftshift(f);

% 應用濾鏡及反變換
out = s .* ff; %對應元素相乘實作頻域濾波
out = ifftshift(out);
out = ifft2(out);

% 求模值
out = abs(out);

% 歸一化以便顯示
out = out/max(out(:));

end

4、頻率低通濾波器

在頻譜中，低頻主要對應影像在平滑區域的總體灰度級分布，而高頻對應影像的細節部分，如邊緣和噪聲，因此，影像平滑可以通過衰減影像頻譜中的高頻部分來實作，這就建立了空間域影像平滑和頻域低通濾波之間對應關系，

4.1、理想低通濾波器及其實作

4.4.1、理論基礎

最容易想到的衰減高頻成份的方法就是咋一個稱為 “截止頻率” 的位置 “截斷” 所有的高頻成份，將影像頻譜中所有高于這一截斷頻率的頻譜成份設為 0，低于截止頻率的成份設為保持不變，能夠達到這種效果的濾波器如圖：

請添加圖片描述

如果影像的寬度為 M，高度為 N，那么理想低通濾波器可形式化地描述為：
H ( u , v ) = { 1 , [ ( u ? M 2 ) 2 + ( v ? N 2 ) 2 ] ≤ D 0 0 , [ ( u ? M 2 ) 2 + ( v ? N 2 ) 2 ] ≥ D 0 H(u,v)=\begin{cases}1,[(u-\frac{M}{2})^2+(v-\frac{N}{2})^2]\leq D_0 \\ 0,[(u-\frac{M}{2})^2+(v-\frac{N}{2})^2]\ge D_0 \end{cases} H(u,v)={1,[(u?2M?)2+(v?2N?)2]≤D0?0,[(u?2M?)2+(v?2N?)2]≥D0??
其中，D₀ 表示理想低通濾波器的截止頻率，濾波器的頻率域原點在頻譜影像的中心處，在以截止頻率為半徑的圓形區域之內的濾鏡元素值全部為 1，而該圓之外的濾鏡元素值全部為 0，理想低通濾波器的頻率特性在截止頻率處十分陡峭，無法用硬體實作，這也是我們稱之為理想的原因，

理想低通濾波器可在一定程度上去除影像噪聲，但由此帶來的影像邊緣和細節的模糊效應也是較為明顯，其效果類似于均值濾波，

4.4.2、Matlab 實作

撰寫了一個 imidealflpf 函式可以得到截止頻率為 freq 的理想低通濾波器：

function ff = imidealflpf(I, freq)
 % imidealflpf函式     構造理想的頻域低通濾波器
 % I引數               輸入的灰度影像
 % freq引數            低通濾波器的截止頻率
 % 回傳值 ff           與I等大的頻域濾鏡   
 
 [M,N] = size(I);
 ff = ones(M, N);
 for i = 1:M
     for j = 1:N
         if (sqrt (((i-M/2)^2 + (j-N/2)^2 ))>freq)  
             ff(i,j)=0; %高于截止頻率 設為0
         end
     end
 end

【例】理想低通濾波器的實作

img = imread('lena_Rician10.bmp');

%生成濾鏡  截止頻率為20
ff = imidealflpf(img,20);
out = imfreqfilt(img,ff);

%生成濾鏡  截止頻率為40
ff = imidealflpf(img,40);
out1 = imfreqfilt(img,ff);

%生成濾鏡  截止頻率為80
ff = imidealflpf(img,80);
out2 = imfreqfilt(img,ff);

figure;
subplot(2,4,1);
imshow(img);
title('原圖');

subplot(2,4,2);
imshow(out);
title('截止頻域為20');

subplot(2,4,3);
imshow(out1);
title('截止頻率為40');

subplot(2,4,4);
imshow(out2);
title('截止頻率為80');

%展示頻域的頻譜
f = fft2(img);
f = abs(fftshift(f)); % 頻譜
f = log(1 + f);  %壓縮頻率的范圍

f1 = fft2(out);
f1 = abs(fftshift(f1)); % 頻譜
f1 = log(1 + f1);  %壓縮頻率的范圍

f2 = fft2(out1);
f2 = abs(fftshift(f2)); % 頻譜
f2 = log(1 + f2);  %壓縮頻率的范圍

f3 = fft2(out2);
f3 = abs(fftshift(f3)); % 頻譜
f3 = log(1 + f3);  %壓縮頻率的范圍


subplot(2,4,5);
imshow(f,[]);
title('原圖對應的頻譜');

subplot(2,4,6);
imshow(f1,[]);
title('freq = 20 的影像對應的頻譜');

subplot(2,4,7);
imshow(f2,[]);
title('freq = 40 的影像對應的頻譜');

subplot(2,4,8);
imshow(f3,[]);
title('freq = 80 的影像對應的頻譜');

在這里插入圖片描述

當截止頻率非常低時，只有非常靠近原點的低頻成份能夠通過，影像模糊嚴重；截止頻率非常高，通過的頻率成份就越多，影像模糊的程度就越小，所獲得的影像也就越接近原影像，但可以看出，理想低通濾波器不能很好地兼顧噪聲濾除與細節保留這兩方面，

4.2、高斯低通濾波器及其實作

4.2.1、理論基礎

高斯低通濾波器的頻域二維形式由下式給出：
H ( u , v ) = e ? [ ( u ? M 2 ) 2 + ( v ? N 2 ) 2 ] / 2 σ 2 H(u,v)=e^{-[(u-\frac{M}{2})^2+(v-\frac{N}{2})^2]/2\sigma^2} H(u,v)=e?[(u?2M?)2+(v?2N?)2]/2σ2
高斯函式具有相對簡單的形式，而且它的傅立葉變換和傅立葉反變換都是實高斯函式，為了簡單起見，下面給出一個一維高斯函式的傅立葉變換和傅立葉反變換作為例子，式(5) 告訴我們一維高斯函式的傅立葉正變換和反變換仍為高斯函式，
H ( u ) = A e ? u 2 2 σ 2 ?? ? ?? h ( x ) = 2 π σ A ? 2 π 2 σ 2 x 2 H(u)=Ae^{-\frac{u^2}{2\sigma^2}} \iff h(x)=\sqrt{2\pi\sigma}A^{-2\pi^2\sigma^2x^2} H(u)=Ae?2σ2u2??h(x)=2πσ ?A?2π2σ2x2
式中：σ 是高斯曲線的標準差，

頻域和與之對應的空間域一維高斯函式的圖形如下：

在這里插入圖片描述

當 σ 增大時，H(u) 的影像傾向于變窄，而 h(x) 的影像傾向于變窄和變高，這也體現了頻率域和空間域的對應關系，頻率域濾波器越窄，濾波的高頻成份越多，影像就越平滑（模糊）；而在空間域，對應的濾波器就越窄，相應的卷積模板越平坦，平滑（模糊）效果就越明顯，

下圖給出了 σ 取值為 3 和 1 時的頻域二維高斯濾波器的三維曲面表示，可以看出和二維高斯濾波器擁有同樣的特點，

在這里插入圖片描述

4.2.2、Matlab 實作

撰寫 imgaussflpf 函式：

function ff = imgaussflpf(I,freq,n)
 % imbutterworthflpf函式     構造布特沃斯的頻域低通濾波器
 % I引數                     輸入的灰度影像
 % freq引數                  低通濾波器的截止頻率
 % n引數                     階數
 % 回傳值 ff                 與I等大的頻域濾鏡   
 
 [M,N] = size(I);
 ff = ones(M,N);
 for i = 1:M
     for j = 1:N
         ff(i,j) = 1/(1 + (sqrt((i - M/2)^2 + (j - N/2)^2)/freq)^(2*n));
     end
 end

【例】實作高斯低通濾波

img = imread('lena_Rician10.bmp');

%生成濾鏡  截止頻率為20,階數為2
ff = imgaussflpf(img,20,2);
out = imfreqfilt(img,ff);

%生成濾鏡  截止頻率為20，階數為5
ff = imgaussflpf(img,20,5);
out1 = imfreqfilt(img,ff);

%生成濾鏡  截止頻率為80,階數為2
ff = imgaussflpf(img,80,2);
out2 = imfreqfilt(img,ff);

%生成濾鏡  截止頻率為80，階數為5
ff = imgaussflpf(img,80,5);
out3 = imfreqfilt(img,ff);

figure;
subplot(2,5,1);
imshow(img);
title('原圖');

subplot(2,5,2);
imshow(out);
title('freq = 20,n = 2');

subplot(2,5,3);
imshow(out1);
title('freq = 20,n = 5');

subplot(2,5,4);
imshow(out2);
title('freq = 80,n = 2');

subplot(2,5,5);
imshow(out3);
title('freq = 80,n = 5');

%展示頻域的頻譜
f = fft2(img);
f = abs(fftshift(f)); % 頻譜
f = log(1 + f);  %壓縮頻率的范圍

f1 = fft2(out);
f1 = abs(fftshift(f1)); % 頻譜
f1 = log(1 + f1);  %壓縮頻率的范圍

f2 = fft2(out1);
f2 = abs(fftshift(f2)); % 頻譜
f2 = log(1 + f2);  %壓縮頻率的范圍

f3 = fft2(out2);
f3 = abs(fftshift(f3)); % 頻譜
f3 = log(1 + f3);  %壓縮頻率的范圍

f4 = fft2(out3);
f4 = abs(fftshift(f4)); % 頻譜
f4 = log(1 + f4);  %壓縮頻率的范圍


subplot(2,5,6);
imshow(f,[]);
title('原圖對應的頻譜');

subplot(2,5,7);
imshow(f1,[]);
title('freq = 20,n = 2的頻譜');

subplot(2,5,8);
imshow(f2,[]);
title('freq = 20,n = 5的頻譜');

subplot(2,5,9);
imshow(f3,[]);
title('freq = 80,n = 2的頻譜');

subplot(2,5,10);
imshow(f4,[]);
title('freq = 80,n = 5的頻譜');

在這里插入圖片描述

5、頻率域高通濾波器

影像銳化可以通過衰減影像頻譜中的低頻成份來實作，這就建立了空間域影像銳化與頻域高通濾波之間對應關系，

5.1、高斯高通濾波器及其實作

對應于高斯低通濾波器：
H ( u , v ) = 1 ? e ? [ ( u ? M 2 ) 2 + ( v ? N 2 ) 2 ] / 2 σ 2 H(u,v)=1-e^{-[(u-\frac{M}{2})^2+(v-\frac{N}{2})^2]/2\sigma^2} H(u,v)=1?e?[(u?2M?)2+(v?2N?)2]/2σ2
在這里插入圖片描述

5.1.1、Matlab 實作

撰寫 imgaussfhpf 函式：

function ff = imgaussfhpf(I,sigma)
 % imgaussfhpf函式     構造高斯的頻域低通濾波器
 % I引數               輸入的灰度影像
 % sigma引數           高斯函式的sigma
 % 回傳值 ff           與I等大的頻域濾鏡   
 
 [M,N] = size(I);
 ff = ones(M,N);
 for i = 1:M
     for j = 1:N
         ff(i,j) = 1 - exp(-((i - M/2).^2 + (j - N/2).^2) / 2 /(sigma.^2));
     end
 end

【例】實作高斯高通濾波器

img = imread('lena_Rician10.bmp');

%生成濾鏡  截止頻率為20
ff = imgaussfhpf(img,20);
out = imfreqfilt(img,ff);

%生成濾鏡  截止頻率為40
ff = imgaussfhpf(img,40);
out1 = imfreqfilt(img,ff);

%生成濾鏡  截止頻率為80
ff = imgaussfhpf(img,80);
out2 = imfreqfilt(img,ff);


subplot(2,4,1);
imshow(img);
title('原圖');

subplot(2,4,2);
imshow(out);
title('sigma = 20');

subplot(2,4,3);
imshow(out1);
title('sigma = 40');

subplot(2,4,4);
imshow(out2);
title('sigma = 80');

%展示頻域的頻譜
f = fft2(img);
f = abs(fftshift(f)); % 頻譜
f = log(1 + f);  %壓縮頻率的范圍

f1 = fft2(out);
f1 = abs(fftshift(f1)); % 頻譜
f1 = log(1 + f1);  %壓縮頻率的范圍

f2 = fft2(out1);
f2 = abs(fftshift(f2)); % 頻譜
f2 = log(1 + f2);  %壓縮頻率的范圍

f3 = fft2(out2);
f3 = abs(fftshift(f3)); % 頻譜
f3 = log(1 + f3);  %壓縮頻率的范圍


subplot(2,4,5);
imshow(f,[]);
title('原圖對應的頻譜');

subplot(2,4,6);
imshow(f1,[]);
title('sigma = 20的頻譜');

subplot(2,4,7);
imshow(f2,[]);
title('sigma = 40的頻譜');

subplot(2,4,8);
imshow(f3,[]);
title('sigma = 80的頻譜');

在這里插入圖片描述

高斯高通濾波器可以較好地提取影像中的邊緣資訊，σ 引數取值越小，邊緣資訊提取越不精確，會包含越多的非邊緣資訊；σ 引數取值越大，邊緣提取越精確，但可能包含不完整的邊緣資訊，

5.2、頻域拉普拉斯濾波器及其實作

5.2.1、理論基礎

頻域拉普拉斯算子的推導可以從一維開始，由傅立葉變換的性質可知：
F F T [ d n f ( x ) d x n ] = ( j u ) n F ( u ) FFT[\frac{d^nf(x)}{dx^n}]=(ju)^nF(u) FFT[dxndnf(x)?]=(ju)nF(u)
因此拉普拉斯算子的傅立葉變換計算如下：
F F T [ ? 2 f ( x , y ) ? x 2 + ? 2 f ( x , y ) ? y 2 ] = ( j u ) 2 F ( u , v ) + ( j v ) 2 F ( u , v ) = ? ( u 2 + v 2 ) F ( u , v ) FFT[\frac{\partial^2f(x,y)}{\partial x^2}+\frac{\partial^2f(x,y)}{\partial y^2}]=(ju)^2F(u,v)+(jv)^2F(u,v)=-(u^2+v^2)F(u,v) FFT[?x2?2f(x,y)?+?y2?2f(x,y)?]=(ju)2F(u,v)+(jv)2F(u,v)=?(u2+v2)F(u,v)
因此有下式成立：
F F T [ ? 2 f ( x , y ) ] = ? ( u 2 + v 2 ) F ( u , v ) FFT[\nabla^2f(x,y)]=-(u^2+v^2)F(u,v) FFT[?2f(x,y)]=?(u2+v2)F(u,v)
也即頻域的拉普拉斯濾波器為：
H ( u , v ) = ? ( u 2 + v 2 ) H(u,v)=-(u^2+v^2) H(u,v)=?(u2+v2)
根據頻域影像頻率原點的平移規律，將上式改寫為：
H ( u , v ) = ? [ ( u ? M 2 ) 2 + ( v ? N 2 ) 2 ] H(u,v)=-[(u-\frac{M}{2})^2+(v-\frac{N}{2})^2] H(u,v)=?[(u?2M?)2+(v?2N?)2]

5.2.2、Matlab 實作

撰寫 imlapf 函式：

function out = imlapf(I)
%imlapf 函式    構造頻域拉普拉斯濾波器
%I 引數         輸入的灰度影像

[M, N] = size(I);
out = ones(M, N);
for i = 1:M
    for j = 1:N
        out(i, j) = -((i-M/2)^2+(j-N/2)^2);
    end
end

end

【例】拉普拉斯濾波實作

img = imread('lena.BMP');

%生成濾鏡  
ff = imlapf(img);
out = imfreqfilt(img,ff);

figure;
subplot(2,2,1);
imshow(img);
title('原圖');

subplot(2,2,2);
imshow(out);
title('laplace濾波影像');

%展示頻域的頻譜
f = fft2(img);
f = abs(fftshift(f)); % 頻譜
f = log(1 + f);  %壓縮頻率的范圍

f1 = fft2(out);
f1 = abs(fftshift(f1)); % 頻譜
f1 = log(1 + f1);  %壓縮頻率的范圍


subplot(2,2,3);
imshow(f,[]);
title('原圖對應的頻譜');

subplot(2,2,4);
imshow(f1,[]);
title('laplace濾波影像對應的頻譜');

在這里插入圖片描述

6、Matlab 綜合案例——利用頻域濾波消除周期噪聲

6.1、頻域帶阻濾波器

所謂 “帶阻” 就是阻止頻譜中某一頻帶范圍的分量通過，其他頻率成份則不受影響，常見的帶阻濾波器有理想帶阻濾波器和高斯帶阻濾波器，

6.1.1、理想帶阻濾波器

理想帶阻濾波器的運算式為：
H ( u , v ) = { 0 , ∣ D ? D 0 ∣ ≤ W 1 , ∣ D ? D 0 ∣ > W H(u,v)=\begin{cases}0,|D-D_0|\leq W \\ 1,|D-D_0|>W \end{cases} H(u,v)={0,∣D?D0?∣≤W1,∣D?D0?∣>W?
式中：D₀ 是阻塞頻帶中心頻率到頻率原點的距離；W 是阻塞頻帶寬度；D 是 (u,v) 點到頻率原點的距離，

在這里插入圖片描述

6.1.2、高斯帶阻濾波器

H ( u , v ) = 1 ? e ? 1 2 [ D 2 ( u , v ) ? D 0 2 D ( u , v ) W ] 2 H(u,v)=1-e^{-\frac{1}{2}[\frac{D^2(u,v)-D_0^2}{D(u,v)W}]^2} H(u,v)=1?e?21?[D(u,v)WD2(u,v)?D02??]2

式中：D₀ 是阻塞頻帶中心頻率到頻率原點的距離；W 是阻塞頻帶寬度；D 是 (u,v) 點到頻率原點的距離，

在這里插入圖片描述

6.1.3、高斯帶阻濾波器的 Matlab 實作

撰寫了 imgaussfbrf 函式：

function ff = imgaussfbrf(I,freq,width)
 % imgaussbrf函式            構造高斯帶阻濾波器
 % I引數                     輸入的灰度影像
 % freq引數                  截止頻率
 % width引數                 阻帶寬度
 % 回傳值 ff                 與I等大的頻域濾鏡   
 
 [M,N] = size(I);
 ff = ones(M,N);
 for i = 1:M
     for j = 1:N
         ff(i,j) = 1 - exp(-0.5*((((i - M/2)^2 + (j - N / 2).^2) - freq^2) / (sqrt((i - M/2)^2 + (j - N / 2)^2) * width))^2); %gauss帶阻濾波器公式
     end
 end

6.2、帶阻濾波消除周期噪聲

6.2.1、得到周期噪聲影像

通常使用正弦平面波來描繪周期性噪聲，

img = imread('lena.BMP');

[M,N] = size(img);
I = img;
for i = 1:M
    for j = 1:N;
        I(i,j) = I(i,j) + 20 * sin(20 * i) + 20 * sin(20 * j); %給原影像添加周期性噪聲
    end
end

figure;
subplot(3,2,1);
imshow(img);
title('原圖');

subplot(3,2,2);
imshow(I);
title('added noise');

%得到頻譜
f = fft2(img);
f = abs(fftshift(f));  %頻譜
f = log(1 + f);   %壓縮頻率范圍

f1 = fft2(I); 
f1 = abs(fftshift(f1));  %頻譜
f1 = log(1 + f1);   %壓縮頻率范圍

subplot(3,2,3);
imshow(f,[]);
title('原圖頻譜');

subplot(3,2,4);
imshow(f1,[]);
title('added noise 頻譜');

%帶阻濾波
ff = imgaussfbrf(I,50,5); % 構造帶阻濾波器，截止頻率為50，阻帶寬度為5

subplot(3,2,5);
imshow(ff,[]);
title('高斯帶阻濾波器');

out = imfreqfilt(I,ff);  % 帶阻濾波

subplot(3,2,6);
imshow(out,[]);
title('帶阻濾波結果');

在這里插入圖片描述

轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/290976.html

標籤：其他

上一篇：【游戲開發實戰】Unity從零開發多人視頻聊天功能，無聊了就和自己視頻聊天（附原始碼 | Mirror | 多人視頻 | 詳細教程）

下一篇：機器視覺邊緣檢測算子

數字影像處理——頻域濾波基礎

頻域濾波基礎

1、頻域濾波與空域濾波的關系

2、頻域濾波的基本步驟

3、頻域濾波的 Matlab 實作

4、頻率低通濾波器

4.1、理想低通濾波器及其實作

4.4.1、理論基礎

4.4.2、Matlab 實作

4.2、高斯低通濾波器及其實作

4.2.1、理論基礎

4.2.2、Matlab 實作

5、頻率域高通濾波器

5.1、高斯高通濾波器及其實作

5.1.1、Matlab 實作

5.2、頻域拉普拉斯濾波器及其實作

5.2.1、理論基礎

5.2.2、Matlab 實作

6、Matlab 綜合案例——利用頻域濾波消除周期噪聲

6.1、頻域帶阻濾波器

6.1.1、理想帶阻濾波器

6.1.2、高斯帶阻濾波器

6.1.3、高斯帶阻濾波器的 Matlab 實作

6.2、帶阻濾波消除周期噪聲

6.2.1、得到周期噪聲影像