Matlab實現基於頻域對二維信號的低通濾波

基於頻域的低通濾波（二維信號——圖像）

算法分析

1. 傅里葉變換，將灰度圖由f(x,y)->F(u,v)(空域轉頻域),得到圖像在頻域中的頻譜（在頻譜中低頻信號分布在頻譜的四個角落，其余部分為高頻信號，這樣的分布難以濾出高頻或低頻信號）

2. 中心化，將頻譜F(u,v)中心化，將低頻點移到頻譜中心（這樣就可以通過設置一個截止頻率D0，來過濾信號）

3. 遍歷頻譜圖，使用巴特沃茲低通濾波器和高斯高斯低通濾波器進行低通濾波，計算濾波器函數h(u,v)與F(u,v)的乘積G(u,v)，直到頻譜圖遍歷完
   巴特沃茲低通濾波器和高斯低通濾波器，計算公式
   

4. 反中心化， 對濾波后的頻譜G(u,v)的低頻點移回到頻譜的四周

5. 傅里葉反變換，將第4步的結果傅里葉反變換G(u,v)->g(x,y)，取g(x,y)實數部分作為最后濾波之后的結果，虛數部分是浮點運算存在的誤差造成的（虛數部分的絕對值很小，可以忽略不計），得到濾波后的空域圖像

偽代碼

x = 讀入灰度圖，若是對彩色圖處理，將彩色圖先轉為灰度圖
y = 圖像加噪
f = 將y的數據類型轉換成double，便於運算
F = 二維傅里葉變換（這里可以直接調用fft2函數），得到圖像在頻域里的頻譜
F1 = 對頻譜你中心化，將低頻頻譜由四周轉換到中心，距離中心越遠的頻率越高
[M, N] = 求得頻譜的大小
n = 巴特沃斯低通濾波器的階數
D0 = 設置截止頻率
m = 頻譜中心行坐標
n = 頻譜中心列坐標
遍歷頻譜
        D = 計算頻譜中當前點與頻譜中心點的距離
        h1 = 根據步驟3的公式1計算巴特沃茲低通濾波器
        h2 = 根據步驟3的公式2計算高斯低通濾波器
        G1(i,j) = 計算巴特沃茲低通濾波器當前像素點灰度值的乘積
        G2(i,j) = 計算高斯低通濾波器當前像素點灰度值的乘積
G1 = 將巴特沃茲低通濾波后的頻譜反中心化
G11 = 將頻譜G1二維傅里葉反變換后，取實部
G2 = 將高斯低通濾波后的頻譜反中心化
G22 = 將頻譜G2二維傅里葉反變換后，取實部
%將巴特沃茲低通濾波后的空域圖像G11和高斯低通濾波后的空域圖像G22顯示出來

代碼

A = imread('lena.png');
x = rgb2gray(A);
y=imnoise(x,'gaussian',0,0.02);  %加入均值為0，方差為0.02的高斯噪聲
f = double(y); %將數據類型轉換成double
F = fft2(f); %二維傅里葉變換
F1 = fftshift(F); %中心化，將低頻頻譜由四周轉換到中心，距離中心越遠的頻率越高
[M, N] = size(F);
n = 2;  %巴特沃斯低通濾波器的階數
D0 = 50;  %設置截止頻率為50
m = fix(M/2);
n = fix(N/2);
for i=1:M  
    for j=1:N
        D = sqrt((i-m)^2+(j-n)^2);
        h1 = 1 / (1 + (D / D0)^(2 * n));  %計算巴特沃茲低通濾波器
        h2 = exp(-(D.^2)./(2 * (D0^2)));  %計算高斯低通濾波器
        G1(i,j) = h1 * F1(i,j);  
        G2(i,j) = h2 * F1(i,j);
    end
end
G1 = ifftshift(G1);  %將巴特沃茲低通濾波后的頻譜反中心化
G11 = uint8(real(ifft2(G1)));  %取二維傅里葉反變換后的實部
G2 = ifftshift(G2);  %將高斯低通濾波后的頻譜反中心化
G22 = uint8(real(ifft2(G2)));  %取二維傅里葉反變換后的實部
subplot(2,4,1),imshow(x),title('原圖');
subplot(2,4,2),imshow(y),title('加入高斯噪聲后');
subplot(2,4,3),imshow(G11),title('二階巴特沃斯低通濾波后圖像');
subplot(2,4,4),imshow(G22),title('高斯低通濾波后圖像');
subplot(2,4,5),imshow(log(1+abs(F)),[ ]),title('加噪圖像的頻譜');
subplot(2,4,6),imshow(log(1+abs(F1)),[ ]),title('頻譜中心化');
subplot(2,4,7),imshow(log(1+abs(G1)),[ ]),title('巴特沃茲低通濾波器頻譜反中心化');
subplot(2,4,8),imshow(log(1+abs(G2)),[ ]),title('高斯低通濾波器頻譜反中心化');

實驗結果

實驗分析

• 由圖所示，直接調用fft2函數將空域轉到頻域時，低頻點分布在頻譜的四周，不管是低通還是高通濾波都比較難限制，所以我們可以通過調用fftshift將頻譜中心化，把低頻信號移動到頻譜中心，這樣我們可以通過指定D0（截止頻率），通過設置低通濾波器計算頻譜上的一點到頻譜中心的距離D(U,V) 與D0比較，大於D0舍棄，只保留小於D0的低頻點

• 經過以上在頻域中低通濾波，將頻域反變換回空域時要將低頻信號和高頻信號反中心化回去，再進行反變換

• 原理總結：傅里葉變換后可以得到信號中有哪些頻率的信號，將需要濾除的頻率成分的幅值置為0，然后再經過傅里葉反變化就得到濾波的目的

• 為什么要將時域轉換到頻域？在時域中無法轉換的數學函數，轉換到頻域能簡單的實現