前言

在野外數(shù)據(jù)采集中，雖然單個儀器采集的是一維信號，但是當把多臺儀器數(shù)據(jù)匯總并生成做二維剖面的圖像時，噪聲可不只有一維的，更有x,y兩個方差同時存在的"二維噪聲"！我們已經(jīng)知道一維噪聲可以用一維傅里葉變換到頻域濾波，同理二維噪聲也可以用二維傅里葉變換到"頻率濾波"。

二維傅里葉正變換的原理

筆者很討厭一上來就看到一連串復雜的公式！因此當我看懂一個原理后，我就會用最好理解的方式來重述它，畢竟我更偏重于應用。二維傅里葉變化，只用一個公式：

公式1：二維傅里葉變換公式(其實不用看它)

公式中參數(shù)說明：

• M、N分別是圖像的長和寬；
• u、x范圍從1到M-1；v、y范圍從1到N-1。

公式看上去不難，但其實還是不太明確到底怎么用??！它其實可以矩陣相乘的形式表示：

公式2：實際用這個公式，矩陣好編程(就看它就行)

公式2中f是原始二維數(shù)據(jù)矩陣，和分別是如下：

公式3：G1矩陣的內容

公式4：G2矩陣的內容

這用Matlab太好實現(xiàn)了叭！和這么有規(guī)律，很容易就編程出來了~ 然后3個矩陣做個乘積，就換到頻率域了！二維傅里葉變換就結束了！

二維傅里葉反變換

定義的公式為公式5(不用看)，矩陣形式的公式為公式6(看這個即可)：

公式5：二維傅里葉反變換公式

公式6：矩陣形式(看這個)

其中和矩陣為：

公式7：G3矩陣內容

公式8：G4公式內容

到此，二維傅里葉逆變換也結束了！整個二維傅里葉變換就都結束了！真的很簡單！下面我們就在Matlab中手寫實現(xiàn)正、逆這兩個過程。

Matlab程序實現(xiàn)

首先實現(xiàn)正變換程序，對應Matlab自帶函數(shù)為：fft2

clc; clear;

data = imread('zxc.jpg');  % 數(shù)據(jù)——最好比卷積核的尺寸大
data = im2double(data); 
data = rgb2gray(data);     % rgb轉為灰度圖像
subplot(1,3,1);
imshow(data);
title('原始圖像')

zidai = fft2(data);   % matlab自帶函數(shù)，來用對比
subplot(1,3,2);
imshow(real(zidai));  % 一般只要實部
title('自帶的fft2生成的"頻域"圖像');

size_data = size(data);
M = size_data(1);  % 圖(原始數(shù)據(jù)矩陣)的長
N = size_data(2);  % 圖(原始數(shù)據(jù)矩陣)的寬

% 下面是傅里葉正變換必備的一些矩陣:
Wm = exp(-j*2*pi/M);
Wn = exp(-j*2*pi/N); % 不同G中用不同的W
Em = zeros(M);
En = zeros(N);     % E是輔助計算矩陣
Gm = zeros(M)+Wm;
Gn = zeros(N)+Wn;  % G是計算時要用的矩陣
F = zeros(M,N);    % F是轉換到頻域的結果

% 對Gm的計算: 循環(huán)長度為M
fprintf('二維離散傅里葉變換開始:\n');
for row = 0:M-1
    for col = 0:M-1
        Em(row+1,col+1) = row * col;
        Gm(row+1,col+1) = Gm(row+1,col+1)^Em(row+1,col+1);
    end
end
% 對Gn的計算: 循環(huán)長度為N
for row = 0:N-1
    for col = 0:N-1
        En(row+1,col+1) = row * col;
        Gn(row+1,col+1) = Gn(row+1,col+1)^En(row+1,col+1);
    end
end

F = real(Gm*data*Gn);  % F = Gm*f*Gn是計算公式，一般只要實部
subplot(1,3,3);
imshow(F);
title('手寫的myfft2生成的"頻域"圖像');

error = sum(sum((real(F)-real(zidai)).^2));
if error < 10^(-10)
    fprintf('自帶與手寫結果一致!\n');
else
    fprintf('不一致!\n');
end

接著正變換結果(把頻域結果當輸入)逆變換程序如下，對應Matlab自帶函數(shù)：ifft2

% 鑒于正向fft2手寫與自帶結果一致；
% ifft2的輸入就直接用自帶的fft2的結果。
clc; clear;

data = imread('zxc.jpg');  % 數(shù)據(jù)——最好比卷積核的尺寸大
data = im2double(data); 
data = rgb2gray(data);     % rgb轉為灰度圖像
subplot(1,3,1);
imshow(data);
title('原始圖像')

F = fft2(data);
subplot(1,3,2);
imshow(real(F));  % 一般畫圖只要實部, 作為輸入時實虛都要??！
title('自帶的fft2生成的"頻域"圖像');

% s = ifft2(F);
% subplot(1,3,3);
% imshow(s);
% return;

size_data = size(F);
M = size_data(1);  % 圖(原始數(shù)據(jù)矩陣)的長
N = size_data(2);  % 圖(原始數(shù)據(jù)矩陣)的寬

% 下面是傅里葉逆變換必備的一些矩陣:
Wm = exp(-j*2*pi/M);
Wn = exp(-j*2*pi/N);  % 不同G中用不同的W
Em = zeros(M);
En = zeros(N);        % E是輔助計算矩陣
Gm = zeros(M)+Wm;
Gn = zeros(N)+Wn;  % G是計算時要用的矩陣
f = zeros(M,N);    % F是轉換到頻域的結果

% 對Gm的計算: 循環(huán)長度為M
fprintf('二維離散反傅里葉變換開始:\n');
for row = 0:M-1
    for col = 0:M-1
        Em(row+1,col+1) = -row * col;
        Gm(row+1,col+1) = Gm(row+1,col+1)^Em(row+1,col+1);
    end
end
Gm = Gm/M;
% 對Gn的計算: 循環(huán)長度為N
for row = 0:N-1
    for col = 0:N-1
        En(row+1,col+1) = -row * col;
        Gn(row+1,col+1) = Gn(row+1,col+1)^En(row+1,col+1);
    end
end
Gn = Gn/N;   % 注意：這個/N和上面的/M都是算完G之后才除以的！因為上面計算的時候是冪項變化！

f = real(Gm*F*Gn);  % f = Gm*F*Gn是計算公式，一般只要實部
subplot(1,3,3);
imshow(f);
title('手寫的myidft2生成的"原始"圖像');

error = sum(sum((real(f)-real(data)).^2));
if error < 10^(-10)
    fprintf('反變換后與原圖一致!\n');
else
    fprintf('不一致!\n');
end