FFT
Discrete Fourier Transform(DFT)
对于一个二维周期函数$f(x,y)$，比如一张大小为$M\times N$并延拓到整个平面的灰度图，其离散傅里叶变换公式如下
$$正变换\quad g(u,v)=\frac{1}{\sqrt{MN}}\sum\limits_{x=0}^{M-1}\sum\limits_{y=0}^{N-1}f(x,y)e^{-2\pi i(\frac{x u}{M}+\frac{y v}{N})}\\
\quad\\
逆变换\quad f(x,y)=\frac{1}{\sqrt{MN}}\sum\limits_{u=0}^{M-1}\sum\limits_{v=0}^{N-1}g(u,v)e^{+2\pi i(\frac{x u}{M}+\frac{y v}{N})}
$$
下面是对一张$400\times400$的图片按照定义进行离散傅里叶变换再进行逆变换得到的图案
这样虽然能够还原最初的形状，但是计算的时间复杂度是$O(N^4)$，随着$N$增大计算耗时迅速增加，也会带来更大的舍入误差。
Fast Fourier Transform(FFT)
如上所见，按照定义来进行离散傅里叶变换很不划算，在实际应用起不到太大的作用，这于是催生了许多高效的离散傅里叶变换算法，统称为快速傅里叶变换，最早的一种基于单位根与多项式性质实现的FFT，称为Cooley-Tukey算法，如下只介绍2-FFT亦即只适用于$M$为2的幂情况下的FFT
首先考虑一维的离散傅里叶变换，为简便起见，令其系数为1
$$g(u)=\sum\limits_{x=0}^{M-1}f(x)e^{-2\pi i \frac{xu}{M}}$$
用$W_M$来表示$M$次单位根，即$e^{-\frac{2\pi i}{M}}$，容易得到如下性质
$$1.周期性\quad W_M^{k+M}=W_M^k\\
\quad\\
2.对称性\quad W_M^{k+\frac{M}{2}}=-W_M^k\\
\quad\\
3.可约性\quad W_M^{mk}=W_{\frac{M}{m}}^k
$$
将$g$用单位根来表示
$$g(u)=\sum\limits_{x=0}^{M-1}f(x)W_{M}^{xu}$$
接下来的变换就是FFT的核心思想：将长序列化为短序列的和。对$g$进行一次如下的奇偶拆分
$$\begin{aligned}g(u)&=\sum\limits_{x=2t}^{M-1}f(x)W_{M}^{xu}+\sum\limits_{x=2t+1}^{M-1}f(x)W_{M}^{xu}\\
\quad\\
&=\sum\limits_{t=0}^{\frac{M}{2}-1}f(2t)W_{\frac{M}{2}}^{tu}+W_{M}^{u}\sum\limits_{t=0}^{\frac{M}{2}-1}f(2t+1)W_{\frac{M}{2}}^{tu}\\
\quad\\
&=g_{even}(u)+W_{M}^{u}g_{odd}(u)
\end{aligned}
$$
因为应用了可约性，周期变为原来的一半，这样我们只能算出$u$序列中前一半的值，这时我们便可利用单位根的周期性与对称性，进行如下操作，让所有的$u$都有相对应的变换
$$
g(u)=g_{even}(u)+W_{M}^{u}g_{odd}(u)\\
\quad\\
g(u+\frac{M}{2})=g_{even}(u)-W_{M}^{u}g_{odd}(u)
$$
于是我们将序列长度为$M$的变换拆成了两个长度为$\frac{M}{2}$的变换，由主定理可知其时间复杂度为$O(Mlog(M))$，如此迭代下去，最终可以让每个DFT的序列长度为1，DFT即其本身，我们只需确定最后一层$x$序列的顺序即可。观察上面的奇偶拆分，如果初始$x$序列中某个位置的二进制末位为零，则将此位置对应的元素移动到前一半序列中，亦即二进制首位变为零，如果倒数第二位又为零，则将其移动到前四分之一序列中，亦即二进制第二位变为零，不难发现，最终的$x$序列中各元素编号即为初始的$x$序列中对应元素编号的二进制反转，这样我们就能由最后一层$x$反推出第一层$x$的傅里叶变换。按照这个思路，我们就可以写出如下一维FFT与IFFT的python代码，为保证对称性，在FFT之前加上一个系数$\frac{1}{\sqrt{M}}$
def ib(i, N): #二进制反转
return int(
('0' * (len(str(bin(N - 1))) - len(str(bin(i))))
+str(bin(i))[2::])
[::-1], 2)
def W(N, k): #单位根
return complex(np.cos(2 * PI * k / N), -np.sin(2 * PI * k / N))
def FFT(f):
N = len(f)
g = [0] * N
for i in range(N):
g[ib(i, N)] = f[i]
for i in range(int(np.log2(N))): #遍历每层
L = [0] * N
n = 1 << i
m = N >> (i + 1)
for j in range(m): #遍历每块
for k in range(n): #遍历各块中的元素
L[j * n * 2 +k] = g[j * n * 2 + k] + W(2 * n, k) * g[j * n * 2 + k + n]
L[j * n * 2 + k +n] = g[j * n * 2 + k] - W(2 * n, k) * g[j * n * 2 + k + n]
g = L
return np.array(g)/np.sqrt(N)
def IFFT(f):
G=np.conjugate(f) #取共轭
return np.conjugate(FFT(G)) #再次取共轭
实现了一维的FFT，二维以及更高维的FFT也就可以分解为许多一维FFT的复合
$$g(u)=FFT(f)=\frac{1}{\sqrt{M}}\sum\limits_{x=0}^{M-1}f(x)e^{-2\pi i \frac{xu}{M}}\\
\quad\\
f(x)=IFFT(g)=\overline{FFT(\bar{g})}\\
\quad\\
\begin{aligned}\quad g(u,v)&=FFT2(f)=\frac{1}{\sqrt{MN}}\sum\limits_{x=0}^{M-1}\sum\limits_{y=0}^{N-1}f(x,y)e^{-2\pi i(\frac{x u}{M}+\frac{y v}{N})}\\
\quad\\
&=\frac{1}{\sqrt{M}}\sum\limits_{x=0}^{M-1}(\frac{1}{\sqrt{N}}\sum\limits_{y=0}^{N-1}f(x,y)e^{-2\pi i \frac{y v}{N}})e^{-2\pi i \frac{x u}{M}}\\
\quad\\
&=FFT_y(FFT_x(f))\\
\end{aligned}\\
\quad\\
f(x,y)=IFFT2(g)=\overline{FFT2(\bar{g})}
$$
写出FFT2的python代码如下
def FFT2(f):
M,N=len(f),len(f[0])
G=np.empty(dtype=complex,shape=(M,N))
for j in range(M):
G[j]=FFT(f[j])
G=G.T
for i in range(N):
G[i]=FFT(G[i])
return G.T
def IFFT2(f):
G=np.conjugate(f)
return np.conjugate(FFT2(G))
利用FFT2可以去除图像中特定空间频率的成分，以下是对图像分别进行高通滤波与低通滤波所得到的结果
还可以达到将铁栅栏变成玻璃栅栏的效果