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逆变换差异较大。