离散傅里叶变换代码实现

在 重磅推出离散傅里叶变换 这篇文章中，给出了周期为 $N$ 的周期信号 $x[n]$ 在一个周期内的截断信号（该截断信号长度为 $N$, 取 $n=0,1,\cdots N-1$）的离散傅里叶变换的定义式：

$$X[k]=\sum _{n=0}^{N-1}x[n]W_N^{nk},k=0,1\cdots ,N-1$$
以及对应的傅里叶反变换式：
$$x[n]=\frac{1}{N}\sum _{k=0}^{N-1}X[k]W_N^{-nk},n=0,1\cdots ,N-1$$

离散傅里叶变换存在的意义是让计算机可以进行时域、频域间的相互转换。实际上，离散周期信号是无限长的，但计算机能处理的离散信号肯定是有限长的。假设要处理的某一离散信号长度为 $N$ , 我们就把这个要处理的信号看做是周期为 $N$ 的信号在一个周期内的截断信号。然后，再利用离散傅里叶变换的定义式进行计算。

根据上述公式，以下给出离散傅里叶变换、离散傅里叶反变换的代码实现，并和numpy科学计算库里的fft计算结果比较后，得到了一致的结果，确保了代码的准确性（注：代码是为了和公式有较好的匹配度而写成的，旨在有好的可读性，程序运行速度还可以优化）。

import numpy as np


def dft(x):
    N = len(x)
    WN = np.power(np.e, complex(0, -1) * (2 * np.pi) / N)
    Xk = []
    for k in range(N):
        value = 0
        for n in range(N):
            value += x[n] * np.power(WN, n * k)
        Xk.append(value)
    return np.array(Xk)


def inverse_dft(X):
    N = len(X)
    WN = np.power(np.e, complex(0, -1) * (2 * np.pi) / N)
    xn = []
    for n in range(N):
        value = 0
        for k in range(N):
            value += 1 / N * X[k] * np.power(WN, -k * n)
        xn.append(value)
    return np.array(xn)


if __name__ == '__main__':
    sig = np.array([1, 2, 3, 2, 1, 4, 5, 6])  # 随便定义一个离散数据序列
    fft_res = np.fft.fft(sig)  # 利用numpy封装的快速傅里叶变换公式求解变换结果
    res = dft(sig)  # 利用离散傅里叶变换公式求解变换结果
    rec_sig = inverse_dft(res)  # 用离散傅里叶变换的结果重构原信号
    
    # 如果下行不报错，说明上述dft函数运行结果和numpy的fft是一致的
    np.testing.assert_array_almost_equal(res, fft_res)  
    # 如果下行不报错，说明上述离散傅里叶反变换结果正确
    np.testing.assert_array_almost_equal(sig, rec_sig)  


 

  
 

从代码中能够清晰的看到，用了两层嵌套的 $for$ 循环，所以离散傅里叶变换算法复杂度是 $O(n^2)$ ，为了提高计算速度，在不改变算法原理的情况下，人们对该算法进行了优化，发展出了快速傅里叶变换（FFT），算法复杂度是 $O(n\log n)$ ，后面有机会的话再介绍吧。