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逆變換差別較大。

---