离散傅里叶变换代码实现

在 重磅推出离散傅里叶变换 这篇文章中给出了周期为 $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)$ 后面有机会的话再介绍吧。