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

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

---

一、傅里叶级数和变换

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

，系数为

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

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

上述两个式子共同构成傅里叶变换对，可以表示为。

同样，令是两个连续变量和的连续函数，则其二维连续傅里叶变换对为

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

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

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

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

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

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

三、二维DFT及其可分离性

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

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

上述过程可以看到，的二维DFT可通过计算的每一行的一维变换，然后沿计算结果的每一列计算一维变换来得到。也就是说二维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

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

然后，我们可以发现上式的右侧是的DFT。这告诉我们，若把代入计算二维傅里叶正变换的算法中，则结果将是。所以将这个结果取复共轭并乘以就可以得到，就是我们想要的反变换，实现过程如下。

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

五、傅里叶变换中心化

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

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

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

图1 冈萨雷斯原书中关于一维和二维傅里叶变换中心化的说明

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

#中心化
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包中的函数进行正变换后使用自己实现的函数进行反变换。

 七、多说两句

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