Matlab圖像處理系列4———傅立葉變換和反變換的圖像

注意：這一系列實驗的圖像處理程序，使用Matlab實現最重要的圖像處理算法

---

1.Fourier兌換

（1）頻域加強

除了在空間域內可以加工處理圖像之外，咱們還可以將圖像變換到其它空間後進行處理。這些方法稱爲變換域方法，最多見的變換域是頻域。

使用Fourier變換把圖像從空間域變換到頻域。在頻域內作對應加強處理，再從頻域變換到空間域獲得處理後的圖像。

咱們這裏主要學習Fourier變換和FFT變換的算法，沒有學過通訊原理，我對信號、時域分析也不是很是清楚。

---

2.FFT算法

（1）離散Fourier變換，DFT

函數 f(x) 的DFT F(v) 的計算式爲：

F(v)=∑x=1Nf(x)e−i2πvxN,v=1,2,...,N

很是顯然求出所有的長度爲N的信號的DFT需要 N×O(N)=O(N2) 的時間，比較慢，所以咱們需要引入時間複雜度爲 O(NlogN) 的FFT算法。

（2）高速Fourier變換，FFT

高速傅立葉變換FFT是利用單位複數根的特殊性質（消去引理和折半引理，見《算法導論》（第3版中文版）P532具體論述）。在 O(NlogN) 時間內計算出DFT的一種高速算法。

FFT有兩種基本實現方式：

• 遞歸FFT
• 迭代FFT。也叫FFT蝶式運算

遞歸FFT由於遞歸棧開銷大且容量有限等缺陷（但理解easy），通常計算採取迭代FFT實現。

（3）迭代FFT

直接給出算法導論版本號的迭代FFT算法：

當中求逆序數拷貝的函數爲：

FFT採用折半迭代的思想，所以速度能減小到 O(NlogN) 。

（4）迭代FFTMatlab實現

Matlab有fft函數，咱們也可以本身實現：

function [ fft_m ] = IterativeFFT( vec )
    clear i;

    n = length(vec);

    fft_m = BitReverseCopy(vec);
    for s = 1 : log2(n)
        m = power(2, s);
        wm = exp(- 2 * pi * i / m);

        for k = 0 : m : n - 1
            w = 1;
            for j = 0 : m / 2 - 1
                t = w * fft_m(k + j + m / 2 + 1);
                u = fft_m(k + j + 1);
                fft_m(k + j + 1) = u + t;
                fft_m(k + j + m / 2 + 1) = u - t;
                w = w * wm;
            end
        end
    end
end

BitReverseCopy函數例如如下：

function [ copy ] = BitReverseCopy( vec )
    n = length(vec);
    copy = zeros(1, n);

    bitn = log2(n);

    for i = 0 : n - 1
        revi = bin2dec(fliplr(dec2bin(i, bitn)));
        copy(revi + 1) = vec(i + 1);
    end
end

需要特別注意的是:

• 通常給出的FFT算法僞代碼都是基於下標從零開始的數組。寫在Matlab需要考慮映射關係
• clear i是爲了怕以前有變量i和複數記號i混淆，清楚Matlab workspace中的緩存
• 默認vec是double類型的。由於中間採用不少double類型運算

---

3.圖像的二維Fourier變換

（1）二維DFT

二維DFT定義公式爲：

F(u,v)=∑x=1N∑y=1Nf(x,y)e−i2π(ux+vy)N,u,v=1,2,...,N

計算一個頻域點需要 O(N2) 的時間，那麼整個二維頻域計算需要 O(N4) 的時間，效率很是低。

（2）將二維DFT分解爲兩個一維DFT

Fourier變換的變換核（公式中和 f(x,y) 無關的部分）具備對稱的性質（詳見《圖像project（上冊）圖像處理》P81），所以利用對稱性可以將二維DFT分解爲兩個一維DFT：
先對二維矩陣的每一列作一維DFT：

F(x,v)=∑y=1Nf(x,y)e−i2πvyN,v=1,2,...,N

再對變換後的矩陣的每一行作一維DFT：

F(u,v)=∑x=1NF(x,v)e−i2πuxN,u=1,2,...,N

最後以 2×N×O(N2)=O(N3) 的時間完畢二維傅立葉變換。

（3）用一維FFT實現二維FFT

相同的咱們可以用兩個一維FFT實現二維FFT，時間複雜度 O(N2logN) ：

function [ mfft2 ] = JCGuoFFT2( data )
    h = size(data, 1);
    w = size(data, 2);
    mfft2 = data;

    if power(2, log2(h)) ~= h || power(2, log2(w)) ~= w
        disp('JCGuoFFT2 exit: h and w must be the power of 2!')
    else
        for i = 1 : h
            mfft2(i, :) = IterativeFFT(mfft2(i, :));
        end

        for j = 1 : w
            mfft2(:, j) = IterativeFFT(mfft2(:, j));
        end
    end
end

代碼很是簡單，先作FFT行變換再作FFT列變換。以前忘記提到。我這裏實現的都是長度必須是2的次方的Fourier變換。所以有時候會作一些長度規範檢查。

（4）變換結果

通過JCGuoFFT2的二維傅立葉變換，咱們可以獲得複平面內各個點的複數值，那麼顯示在圖像中需要先求出幅值：

pic1_fft = JCGuoFFT2(double(pic1));
pic1_fft_amp = abs(pic1_fft);

在對幅值作一次log變換獲得較好的頻域圖像：

pic1_fft_amp_log = log(1 + pic1_fft_amp);

繪製結果例如如下圖：

（5）低頻信號移到圖像中心點

由於複數運算的週期特性，圖像的Fourier變換在複平面上是全然對稱的。可以想象平面是由無限多個上圖（右）頻域圖像拼接而成的二維陣列。通常研究頻域圖像是把低頻部分，也就是變換後的邊緣部分移到圖像中心點。Matlab提供fftshift函數完畢平移。

平移的思路有兩個

• 經過Fourier變換平移定理先把原始圖像作變換再作FFT
• 先作FFT後再依據頻域圖像的對稱性作對稱變換

查閱Matlab文檔發現它是採用另一種方法，對圖像作下面子矩陣交換：

那麼咱們可以很是easy的寫出本身的fftshift：

function [ after ] = FFTShift( before )
    h = size(before, 1);
    w = size(before, 2);

    after = before;

    if power(2, log2(h)) ~= h || power(2, log2(w)) ~= w
        disp('FFTShift exit: h and w must be the power of 2!')
    else
        hh = h / 2;
        hw = w / 2;
        after(1 : hh, 1 : hw) = before(hh + 1 : h, hw + 1 : w);
        after(1 : hh, hw + 1 : w) = before(hh + 1 : h, 1 : hw);
        after(hh + 1 : h, 1 : hw) = before(1 : hh, hw + 1 : w);
        after(hh + 1 : h, hw + 1 : w) = before(1 : hh, 1 : hw);
    end
end

將低頻部分平移到中心點後結果爲：

---

4.圖像的二維反Fourier變換

（1）一維逆DFT和一維逆FFT

一維離散傅立葉變換的逆變換是將e的指數部分變號，而後整體除以長度N（Fourier變換與逆變換關於符號、係數有很是多種組合的定義。但他們都是等價且對稱的。這裏的定義配合Matlab作fft實際效果。）：

F(v)=1N∑x=1Nf(x)ei2πvxN,v=1,2,...,N

相同咱們可以依據iDFT的定義推演iFFT的算法，其迭代版本號的Matlab實現例如如下：

function [ ifft_m ] = IterativeIFFT( vec )
    clear i;
    n = length(vec);
    ifft_m = BitReverseCopy(vec);
    for s = 1 : log2(n)
        m = power(2, s);
        wm = exp(2 * pi * i / m);

        for k = 0 : m : n - 1
            w = 1;
            for j = 0 : m / 2 - 1
                t = w * ifft_m(k + j + m / 2 + 1);
                u = ifft_m(k + j + 1);
                ifft_m(k + j + 1) = u + t;
                ifft_m(k + j + m / 2 + 1) = u - t;
                w = w * wm;
            end
        end
    end
    ifft_m = ifft_m ./ n;
end

（2）二維逆FFT

二維逆FFT和二維FFT的思路一致。對所有行和列分別做一次iFFT就能夠：

function [ mifft2 ] = JCGuoIFFT2( data )
    h = size(data, 1);
    w = size(data, 2);
    mifft2 = data;

    if power(2, log2(h)) ~= h || power(2, log2(w)) ~= w
        disp('JCGuoIFFT2 exit: h and w must be the power of 2!')
    else
        for i = 1 : h
            mifft2(i, :) = IterativeIFFT(mifft2(i, :));
        end

        for j = 1 : w
            mifft2(:, j) = IterativeIFFT(mifft2(:, j));
        end
    end
end

（3）逆FFT結果

對以前Rect1.bmp用JCGuoFFT2變換的到的Fourier變換（非shift和log以後、非僅幅度部分）做FFT2逆變換可以直接獲得原圖像：

這幅圖有多個結果。可以看title知道每個結果的含義。圖2-1是用JCGuoIFFT2作傅立葉反變換的結果，獲得的圖像和原圖像、和Matlab ifft2反變換後的圖像都是一致的。

---

5.幅頻特性與相頻特性

（1）對振幅和相位單獨作逆FFT

假設咱們僅僅把圖像Fourier變換的振幅部分和相位部分單獨作二維逆FFT，可以直觀的感覺人眼對圖像幅頻特性和相頻特性的敏感度。

複數 z=a+ib 的幅度/振幅定義爲：

|z|=a2+b2‾‾‾‾‾‾‾√

相位角和相位定義爲：

ϕ(z)=arctanbaeiϕ(z)=eiarctanba

對相位反變換需要在乘以一個係數（我是調出來的。針對圖像，確定可以本身主動的作均衡化）：

pic2_fft_angle = angle(pic2_fft);
clear i;
tmp = 10000 * exp(i * pic2_fft_angle);
pic2_ifft_angle = uint8(JCGuoIFFT2(tmp));

對振幅和相位單獨作逆FFT結果例如如下：

（2）人眼敏感度

觀察上圖，幅頻特性主要涵蓋了圖像顏色的分佈，相頻特性主要刻畫了圖像的邊界信息。人眼對圖像的相頻特性更加敏感。看相頻特性可以大概知道圖像內容。

---

6.Fourier變換的旋轉定理

（1）Fourier變換旋轉定理

將f(x,y)旋轉角度θ0。其Fourier變換F(u,v)也旋轉角度θ0

（2）結果

Rect2.bmp是Rect1.bmp旋轉45度的示意圖（注：原Rect2.bmp是二值的。作了預處理，但其自己邊界不平滑。致使效果不太好。但不影響觀察旋轉定理）：

咱們可以看到幅度FFT正變換和相位FFT你變換都是旋轉了45度，但是幅度FFT逆變換差異較大。