天天看點
傳回

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

       之前的博文主要介紹了空間域内的濾波器，本文主要從頻域的角度進行分析。主要使用傅裡葉變換，将空間域的圖像轉換到頻域内，在頻域内進行數字圖像處理。這部分的内容及其重要，頻域内的處理可以解決空間域内無法完成的圖像增強。本文首先從數學角度，對圖像的頻域内的性質進行分析，然後在着重介紹濾波器在頻域内的性質。

        1.傅裡葉變換與頻域

        在之前的文中，我們已經進行過一些基本的圖像處理。比如，使用低通濾波可以将圖像模糊，也有些許降噪的作用。這些都是在空間域内進行的濾波處理，這個處理主要是依靠卷積來進行計算的。首先，從連續的一維卷積入手，如下所示。

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

       将上式進行傅裡葉變換，可以得到如下結果。

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

        從這個式子，我們可以得到一個重要的結論。也就是，函數

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

與

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

卷積的傅裡葉變換所得到的結果，是函數

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

與

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

的傅裡葉變換

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

與

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

的乘積。再将其總結得簡單易懂一些，有如下結論。

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

        在将其擴充到二維的形況下，假設尺寸為MxN的圖像，如下關系是成立的。

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

       其實到這，基本的原理就明了的。我們所看到的圖像，均為空間域内的表現形式，我們無法辨識出頻域内的圖像。要進行頻域内的濾波器處理，首先就需要進行傅裡葉變換，然後直接進行濾波處理，最後再用反傅裡葉變換倒回到空間域内。

       到此，已經可以開始空間域内的濾波處理了。但是，還有一點需要注意的地方。使用某個一維信号來舉例子，一維信号的傅裡葉變換是以2π為周期的函數。是以，我們常常使用的範圍[-π,π]來表示這個信号的傅裡葉變換，如下所示。

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

        這樣做的好處是，靠近0的成分就是低頻，靠近-π與π的成分就表示高頻。而對于圖像而言，在Matlab中，我們使用fft2()這個函數來求取圖像的傅裡葉變換。

g = fft2(f);

        上面這個代碼求取的，其實是 範圍[0,π]内的傅裡葉變換。為了友善了解，下圖畫出了本行代碼所求取的圖像的傅裡葉變換的範圍（右）和與其等效的一維傅裡葉變換的範圍（左）。

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

       很顯然，這并不是希望的範圍，下面這個代碼可以求取[0,2π]内的傅裡葉變換。

P = 2*M;
Q = 2*N;
F = fft2(f,P,Q);

       下圖畫出了本行代碼所求取的圖像的傅裡葉變換的範圍（右）和與其等效的一維傅裡葉變換的範圍 （左） 。所得到的圖像F(u,v)的尺寸為PxQ。

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

         我們需要對其移動一下，如下圖所示，我們需要的是粉色範圍的區域。

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

下面，從數學上分析一下，如何獲得這個部分的頻譜。對于傅裡葉變換，有如下性質。

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

這個特性稱為平移特性，粉色部分的頻譜，将

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

帶入上式，我們可以得到如下式子。

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

為次，我們已經得到了粉色範圍的頻譜。越靠近傅裡葉頻譜圖像中間的成分，代表了低頻成分。其Matlab代碼如下所示。

[M,N] = size(f);
P = 2*M;
Q = 2*N;
fc = zeros(M,N);

for x = 1:1:M
    for y = 1:1:N
        fc(x,y) = f(x,y) * (-1)^(x+y);
    end
end

F = fft2(fc,P,Q);

        代碼所得到的結果，如下圖所示。

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

        接下來，我們總結一下頻域濾波的步驟：

        ①：先将圖像做頻域内的水準移動，然後求原圖像f(x,y)的DFT，得到其圖像的傅裡葉譜F(u,v)。

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

        ②：與頻域濾波器做乘積，

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

        ③：求取G(u,v)的IDFT，然後再将圖像做頻域内的水準移動（移動回去），其結果可能存在寄生的虛數，此時忽略即可。

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

        ④：這裡使用ifft2函數進行IDFT變換，得到的圖像的尺寸為PxQ。切取左上角的MxN的圖像，就能得到結果了。

        2.低通濾波器

        2.1理想的低通濾波器

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

       其中，D0表示通帶的半徑。D(u,v)的計算方式也就是兩點間的距離，很簡單就能得到。

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

       使用低通濾波器所得到的結果如下所示。低通濾波器濾除了高頻成分，是以使得圖像模糊。由于理想低通濾波器的過度特性過于急峻，是以會産生了振鈴現象。

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器
[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

        2.2巴特沃斯低通濾波器

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

       同樣的，D0表示通帶的半徑，n表示的是巴特沃斯濾波器的次數。随着次數的增加，振鈴現象會越來越明顯。

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

       2.3高斯低通濾波器

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

       D0表示通帶的半徑。高斯濾波器的過度特性非常平坦，是以是不會産生振鈴現象的。

[數字圖像處理]頻域濾波(1)--基礎與低通濾波器
[數字圖像處理]頻域濾波(1)--基礎與低通濾波器

       3.實作代碼 

close all;
clear all;

%% ---------Ideal Lowpass Filters (Fre. Domain)------------
f = imread('characters_test_pattern.tif');
f = mat2gray(f,[0 255]);

[M,N] = size(f);
P = 2*M;
Q = 2*N;
fc = zeros(M,N);

for x = 1:1:M
    for y = 1:1:N
        fc(x,y) = f(x,y) * (-1)^(x+y);
    end
end

F = fft2(fc,P,Q);

H_1 = zeros(P,Q);
H_2 = zeros(P,Q);

for x = (-P/2):1:(P/2)-1
     for y = (-Q/2):1:(Q/2)-1
        D = (x^2 + y^2)^(0.5);
        if(D <= 60)  H_1(x+(P/2)+1,y+(Q/2)+1) = 1; end    
        if(D <= 160)  H_2(x+(P/2)+1,y+(Q/2)+1) = 1; end
     end
end

G_1 = H_1 .* F;
G_2 = H_2 .* F;

g_1 = real(ifft2(G_1));
g_1 = g_1(1:1:M,1:1:N);

g_2 = real(ifft2(G_2));
g_2 = g_2(1:1:M,1:1:N);         

for x = 1:1:M
    for y = 1:1:N
        g_1(x,y) = g_1(x,y) * (-1)^(x+y);
        g_2(x,y) = g_2(x,y) * (-1)^(x+y);
    end
end

%% -----show-------
figure();
subplot(1,2,1);
imshow(f,[0 1]);
xlabel('a).Original Image');

subplot(1,2,2);
imshow(log(1 + abs(F)),[ ]);
xlabel('b).Fourier spectrum of a');

figure();
subplot(1,2,1);
imshow(H_1,[0 1]);
xlabel('c).Ideal Lowpass filter(D=60)');

subplot(1,2,2);
h = mesh(1:20:P,1:20:Q,H_1(1:20:P,1:20:Q));
set(h,'EdgeColor','k');
axis([0 P 0 Q 0 1]);
xlabel('u');ylabel('v');
zlabel('|H(u,v)|');

figure();
subplot(1,2,1);
imshow(log(1 + abs(G_1)),[ ]);
xlabel('d).Result of filtering using c');

subplot(1,2,2);
imshow(g_1,[0 1]);
xlabel('e).Result image');

figure();
subplot(1,2,1);
imshow(H_2,[0 1]);
xlabel('f).Ideal Lowpass filter(D=160)');

subplot(1,2,2);
h = mesh(1:20:P,1:20:Q,H_2(1:20:P,1:20:Q));
set(h,'EdgeColor','k');
axis([0 P 0 Q 0 1]);
xlabel('u');ylabel('v');
zlabel('|H(u,v)|');

figure();
subplot(1,2,1);
imshow(log(1 + abs(G_2)),[ ]);
xlabel('g).Result of filtering using e');

subplot(1,2,2);
imshow(g_2,[0 1]);
xlabel('h).Result image');
close all;
clear all;

%% ---------Butterworth Lowpass Filters (Fre. Domain)------------
f = imread('characters_test_pattern.tif');
f = mat2gray(f,[0 255]);

[M,N] = size(f);
P = 2*M;
Q = 2*N;
fc = zeros(M,N);

for x = 1:1:M
    for y = 1:1:N
        fc(x,y) = f(x,y) * (-1)^(x+y);
    end
end

F = fft2(fc,P,Q);

H_1 = zeros(P,Q);
H_2 = zeros(P,Q);

for x = (-P/2):1:(P/2)-1
     for y = (-Q/2):1:(Q/2)-1
        D = (x^2 + y^2)^(0.5);
        D_0 = 100;
        H_1(x+(P/2)+1,y+(Q/2)+1) = 1/(1+(D/D_0)^2);   
        H_2(x+(P/2)+1,y+(Q/2)+1) = 1/(1+(D/D_0)^6);
     end
end


G_1 = H_1 .* F;
G_2 = H_2 .* F;

g_1 = real(ifft2(G_1));
g_1 = g_1(1:1:M,1:1:N);

g_2 = real(ifft2(G_2));
g_2 = g_2(1:1:M,1:1:N);         

for x = 1:1:M
    for y = 1:1:N
        g_1(x,y) = g_1(x,y) * (-1)^(x+y);
        g_2(x,y) = g_2(x,y) * (-1)^(x+y);
    end
end

%% -----show-------
figure();
subplot(1,2,1);
imshow(f,[0 1]);
xlabel('a).Original Image');

subplot(1,2,2);
imshow(log(1 + abs(F)),[ ]);
xlabel('b).Fourier spectrum of a');

figure();
subplot(1,2,1);
imshow(H_1,[0 1]);
xlabel('c)Butterworth Lowpass (D_{0}=100,n=1)');

subplot(1,2,2);
h = mesh(1:20:P,1:20:Q,H_1(1:20:P,1:20:Q));
set(h,'EdgeColor','k');
axis([0 P 0 Q 0 1]);
xlabel('u');ylabel('v');
zlabel('|H(u,v)|');

figure();
subplot(1,2,1);
imshow(log(1 + abs(G_1)),[ ]);
xlabel('d).Result of filtering using c');

subplot(1,2,2);
imshow(g_1,[0 1]);
xlabel('e).Result image');

figure();
subplot(1,2,1);
imshow(H_2,[0 1]);
xlabel('f).Butterworth Lowpass (D_{0}=100,n=3)');

subplot(1,2,2);
h = mesh(1:20:P,1:20:Q,H_2(1:20:P,1:20:Q));
set(h,'EdgeColor','k');
axis([0 P 0 Q 0 1]);
xlabel('u');ylabel('v');
zlabel('|H(u,v)|');

figure();
subplot(1,2,1);
imshow(log(1 + abs(G_2)),[ ]);
xlabel('g).Result of filtering using e');

subplot(1,2,2);
imshow(g_2,[0 1]);
xlabel('h).Result image');
close all;
clear all;
clc;
%% ---------Gaussian Lowpass Filters (Fre. Domain)------------
f = imread('characters_test_pattern.tif');
f = mat2gray(f,[0 255]);

[M,N] = size(f);
P = 2*M;
Q = 2*N;
fc = zeros(M,N);

for x = 1:1:M
    for y = 1:1:N
        fc(x,y) = f(x,y) * (-1)^(x+y);
    end
end

F = fft2(fc,P,Q);

H_1 = zeros(P,Q);
H_2 = zeros(P,Q);

for x = (-P/2):1:(P/2)-1
     for y = (-Q/2):1:(Q/2)-1
        D = (x^2 + y^2)^(0.5);
        D_0 = 60;
        H_1(x+(P/2)+1,y+(Q/2)+1) = exp(-(D*D)/(2*D_0*D_0));   
        D_0 = 160;
        H_2(x+(P/2)+1,y+(Q/2)+1) = exp(-(D*D)/(2*D_0*D_0));
     end
end

G_1 = H_1 .* F;
G_2 = H_2 .* F;

g_1 = real(ifft2(G_1));
g_1 = g_1(1:1:M,1:1:N);

g_2 = real(ifft2(G_2));
g_2 = g_2(1:1:M,1:1:N);         

for x = 1:1:M
    for y = 1:1:N
        g_1(x,y) = g_1(x,y) * (-1)^(x+y);
        g_2(x,y) = g_2(x,y) * (-1)^(x+y);
    end
end


%% -----show-------
close all;

figure();
subplot(1,2,1);
imshow(f,[0 1]);
xlabel('a).Original Image');

subplot(1,2,2);
imshow(log(1 + abs(F)),[ ]);
xlabel('b).Fourier spectrum of a');

figure();
subplot(1,2,1);
imshow(H_1,[0 1]);
xlabel('c)Gaussian Lowpass (D_{0}=60)');

subplot(1,2,2);
h = mesh(1:20:P,1:20:Q,H_1(1:20:P,1:20:Q));
set(h,'EdgeColor','k');
axis([0 P 0 Q 0 1]);
xlabel('u');ylabel('v');
zlabel('|H(u,v)|');

figure();
subplot(1,2,1);
imshow(log(1 + abs(G_1)),[ ]);
xlabel('d).Result of filtering using c');

subplot(1,2,2);
imshow(g_1,[0 1]);
xlabel('e).Result image');

figure();
subplot(1,2,1);
imshow(H_2,[0 1]);
xlabel('f).Gaussian Lowpass (D_{0}=160)');

subplot(1,2,2);
h = mesh(1:20:P,1:20:Q,H_2(1:20:P,1:20:Q));
set(h,'EdgeColor','k');
axis([0 P 0 Q 0 1]);
xlabel('u');ylabel('v');
zlabel('|H(u,v)|');

figure();
subplot(1,2,1);
imshow(log(1 + abs(G_2)),[ ]);
xlabel('g).Result of filtering using e');

subplot(1,2,2);
imshow(g_2,[0 1]);
xlabel('h).Result image');

      原文發于部落格：http://blog.csdn.net/thnh169/ 

數字圖像處理 數字信号處理基礎學習 圖像處理 matlab 頻域濾波 低通濾波器
上一篇: [數字圖像處理]圖像複原--逆濾波
下一篇: 數字圖像處理：12.小波變換

繼續閱讀