注：本系列来自于图像处理课程实验，用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由于递归栈开销大且容量有限等缺陷（但理解容易），一般计算采取迭代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
endBitReverseCopy函数如下：
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)无关的部分）具有对称的性质（详见《图像工程（上册）图像处理》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文档发现它是采用第二种方法，对图像做以下子矩阵交换：
那么我们可以很容易的写出自己的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逆变换区别较大。