图像处理之离散傅里叶变换（DFT）

上学期修了数字图像处理这门课程，想着正好趁这个机会写（shui）几篇文章，告诉自己没有白学。傅里叶变换，是图像处理中的一个重要内容，频率域处理的操作都要建立在傅里叶变换的基础上，所以作为这个专栏的开篇，不如就简单介绍和实现一下离散傅里叶变换（DFT）。

一、傅里叶级数和变换

法国数学家傅里叶提出，任何周期函数都可表示为不同频率的正弦函数和/或余弦函数之和，其中每个正弦函数和/或余弦函数都要乘以不同的系数，这个和就称为傅里叶级数。按照这个思想，周期为 $T$ 的连续变量 $t$ 的周期函数 $f(t)$ ，可表示为乘以适当系数的正弦函数和余弦函数之和，即

$f(t)=\sum_{n=-\infty }^{\infty }c_{n}e^{j\tfrac{2\pi n}{T}t}$ ，系数为

$c_{n}=\frac{1}{T}\int_{-T/2}^{T/2}f(t)e^{-j\tfrac{2\pi n}{T}t}dt, n=0,\pm 1,\pm 2,...$

另外根据欧拉公式有 $e^{j\theta }=cos\theta+jsin\theta$ ，所以上述式子便可展开为正弦函数和余弦函数之和（其中涉及到复数的知识请读者自行了解，这里不作过多说明）。

而一些（曲线下方面积有限的）非周期函数也能用正弦函数和/或余弦函数乘以加权函数的积分来表示。这种情况下的公式就是傅里叶变换，其在许多理论和应用科学中起到非常大的作用。连续单变量函数的傅里叶变换和反变换表示为

$\left\{\begin{matrix} F(\mu )=\int_{-\infty }^{\infty }f(t)e^{-j2\pi \mu t}dt\\ \\ f(t) =\int_{-\infty }^{\infty }F(\mu)e^{j2\pi \mu t}d\mu \end{matrix}\right.$

上述两个式子共同构成傅里叶变换对，可以表示为 $f(t)\Leftrightarrow F(\mu)$ 。

同样，令 $f(t,z)$ 是两个连续变量 $t$ 和 $z$ 的连续函数，则其二维连续傅里叶变换对为

$\left\{\begin{matrix} F(\mu,\nu )=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(t,z)e^{-j2\pi (\mu t+\nu z)}dtdz \\ \\f(t,z)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}F(\mu,\nu)e^{j2\pi (\mu t+\nu z)}d\mu d\nu \end{matrix}\right.$

二、离散傅里叶变换（DFT）

有了上一节中所讲的数学基础，我们这里直接给出离散傅里叶的变换对

$\left\{\begin{matrix} F(u)=\sum_{x=0}^{M-1}f(x)e^{-j2\pi ux/M}, u=0,1,2,...,M-1\\ \\ f(x)=\frac{1}{M} \sum_{x=0}^{M-1}F(u)e^{j2\pi ux/M},x=0,1,2,...,M-1 \end{matrix}\right.$

将式子与连续形式的傅里叶变换进行对比，应该是容易理解的。有时，我们也会发现有的公式把 $1/M$ 放在了第一个式子中，这并不会影响两个公式形成一个傅里叶变换对。此外在实现的时候，为了方便我们通常将DFT公式写成下面这种形式

$\left\{\begin{matrix} F(u)=\sum_{x=0}^{M-1}f(x)W_{M}^{ux}, u=0,1,2,...,M-1\\ \\f(x)=\frac{1}{M} \sum_{x=0}^{M-1}F(u)W_{M}^{-ux},x=0,1,2,..,M-1 \end{matrix}\right.$

$W_{M}=e^{-j2\pi /M}$

若用矩阵计算的形式表达上面这个过程，就是

$\begin{bmatrix} F(0)\\ F(1)\\ F(2)\\ ...\\ F(M-1) \end{bmatrix}=\begin{bmatrix} W_{M}^{0} &W_{M}^{0} &W_{M}^{0} &... &W_{M}^{0} \\ W_{M}^{0}&W_{M}^{1} &W_{M}^{2} &... &W_{M}^{M-1} \\ W_{M}^{0}& W_{M}^{2} &W_{M}^{4} &... &W_{M}^{2(M-1)} \\ ...& ... &... & ... &... \\ W_{M}^{0}&W_{M}^{M-1} & W_{M}^{2(M-1)} & ...& W_{M}^{(M-1)^{2}} \end{bmatrix}\begin{bmatrix} f(0)\\ f(1)\\ f(2)\\ ...\\ f(M-1) \end{bmatrix}$

综上，给出一维DFT的实现过程如下。

#一维离散傅里叶变换
def dft(f):#得到长度M = f.shape[0]#计算变换矩阵x, y = np.mgrid[0:M, 0:M]w = x * ybase = np.exp(-1j*2*np.pi/M)W = base**w#矩阵相乘计算结果F = np.dot(W, f)return F

三、二维DFT及其可分离性

类似于一维DFT，我们也可得到如下的二维离散傅里叶变换对

$\left\{\begin{matrix} F(u,v)=\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)e^{-j2\pi (ux/M+vy/N)},u=0,1,2,...,M-1,v=0,1,2,...,N-1\\ \\ f(x,y)=\frac{1}{MN}\sum_{u=0}^{M-1}\sum_{v=0}^{N-1}F(u,v)e^{j2\pi (ux/M+vy/N)} ,x=0,1,2,...,M-1,y=0,1,2,...,N-1\end{matrix}\right.$

接下来我们对正变换进行化简，可得

$F(u,v)=\sum_{x=0}^{M-1}e^{-j2\pi ux/M}\sum_{x=0}^{N-1}f(x,y)e^{-j2\pi vy/N}=\sum_{x=0}^{M-1}F(x,v)e^{-j2\pi ux/M}$

$F(x,v)=\sum_{y=0}^{N-1}f(x,y)e^{-j2\pi xy/N}$

上述过程可以看到， $f(x,y)$ 的二维DFT可通过计算 $f(x,y)$ 的每一行的一维变换，然后沿计算结果的每一列计算一维变换来得到。也就是说二维DFT通过多次一维DFT计算即可得到，参考下面实现过程。

#二维离散傅里叶变换
def dft_2d(f):#得到输入的行数和列数M, N = f.shape[0], f.shape[1]#初始化两个复数类型数组，用于保存结果F, F_x = np.zeros(shape=(M,N), dtype=np.complex128), np.zeros(shape=(M,N), dtype=np.complex128)#逐行逐列进行两次一维傅里叶变换for i in range(M):F_x[i,:] = dft(f[i,:])for i in range(N):F[:,i] = dft(F_x[:,i])return F

四、使用DFT算法计算IDFT

我们将二维离散傅里叶反变换过程两边取复共轭，并将得到的结果乘以 $MN$ 得

$MNf^{*}(x,y)=\sum_{u=0}^{M-1}\sum_{v=0}^{N-1}F^{*}(u,v)e^{-j2\pi (ux/M+vy/N)}$

然后，我们可以发现上式的右侧是 $F^{*}(u,v)$ 的DFT。这告诉我们，若把 $F^{*}(u,v)$ 代入计算二维傅里叶正变换的算法中，则结果将是 $MNf^{*}(x,y)$ 。所以将这个结果取复共轭并乘以 $1/MN$ 就可以得到 $f(x,y)$ ，就是我们想要的反变换，实现过程如下。

#二维离散傅里叶反变换
def idft_2d(F):#得到行数和列数M, N = F.shape[0], F.shape[1]#先对F的共轭进行正向离散傅里叶变换，除以常数后再取共轭f = np.conjugate(dft_2d(np.conjugate(F)))/M/Nreturn f

五、傅里叶变换中心化

可以证明，一维正离散变换和反离散变换都是无限周期的，周期为 $M$ ，即

$\left\{\begin{matrix} F(u)=F(u+kM)\\ \\ f(x) = f(x+kM) \end{matrix}\right.$

如一维情况那样，二维傅里叶变换及其反变换在 $u$ 方向和 $v$ 方向是无限周期的，即

$\left\{\begin{matrix} F(u,v)=F(u+k_{1}M,v)=F(u,v+k_{2}N)=F(u+k_{1}M,v+k_{2}N)\\ \\ f(x,y)=f(x+k_{1}M,y)=f(x,y+k_{2}N)=f(x+k_{1}M,y+k_{2}N) \end{matrix}\right.$

周期性在实现基于DFT的算法中很重要。区间 $0$ 到 $M-1$ 上的变换数据由在点 $M/2$ 处相遇的两个半周期组成，但该周期的较低部分出现在较高的频率处。针对显示和滤波的目的，在该区间内有一个数据连续并且正确排序的完整变换周期更为方便。也就说我们需要先将 $F(u,v)$ 变换到 $F(u-M/2,v-N/2)$ ，方便我们进一步操作，下图显示了一维和二维情况下的这个过程。

下面用代码实现了中心化和逆中心化的过程。

#中心化
def fshift(F):#两个轴平移半个周期M, N = F.shape[0], F.shape[1]F = np.roll(F, int(N/2), axis=0)F = np.roll(F, int(M/2), axis=1)return F#去中心化
def ifshift(F):#两个轴平移半个周期M, N = F.shape[0], F.shape[1]F = np.roll(F, -int(N/2), axis=0)F = np.roll(F, -int(M/2), axis=1)return F

六、效果展示

首先，我们使用DFT得到图片的频率域，然后展示其频谱图。方法是先得到变换后结果的模长，然后使用log函数将其映射到图片可表示的范围。

#图像频谱
def fimage(img):return np.log(np.abs(img))

图2 (a)原图；(b)使用自己完成的DFT函数变换得到频谱图；
(c)使用自己完成的DFT函数和中心化函数变换得到的频谱图；
(d)使用numpy包中FFT函数和中心化函数变换得到的频谱图。

接下来是将自己完成的函数和numpy包中函数混用的情况，从结果上看基本没有差别，可以认为实现效果还不错。

图3 (a)原图；
(b)使用自己实现的函数进行正变换后使用numpy包中的函数进行反变换；
(c)使用numpy包中的函数进行正变换后使用自己实现的函数进行反变换。

七、多说两句

这篇文章参考《冈萨雷斯数字图像处理（第四版）》，只是对书中核心内容进行了总结并加以实现，未涉及到的细节可见原著。