因为 FFT 是基于 DFT 的，所以这里不得不提一下它。
抛开一些意义来看，所谓的 DFT 其实就是矩阵运算罢了，通过原向量 X 和矩阵 F 来获得目标矩阵 Y。DFT 的定义是：
\[y_k=\sum\limits_{n=0}^{N-1}x_ne^{-i2\pi k\frac nN}=(1, e^{-i2\pi k\frac1N},e^{-i2\pi k\frac2N} \dots, e^{-i2\pi k\frac{N-1}{N}}) \left(\begin{matrix} x_0 \\x_1\\x_2\\\vdots\\x_{N-1}\end{matrix}\right) k=0,1,\dots, N-1\]
它的反变换是：\(x_n = \frac1N \sum\limits_{k=0}^{N-1}e^{2\pi ik\frac nN}y_k \quad (n=0, 1, 2, \dots, N-1)\)
令 \( X=(x_0, x_1, x_2, \dots, x_{N-1})^H , Y=(y_0, y_1, y_2, \dots, y_{N-1})^H \)
以及 \( \phi_k = (1, e^{-i2\pi k\frac1N}, e^{-i2\pi k\frac2N}, \dots, e^{-i2\pi k\frac{N-1}{N}})^H \)，易知 \( \phi_k \) 两两正交。有 \( y_k = \phi_k \cdot X = \phi_k^HX \)
现在令 \( F = (\phi_0, \phi_1, \dots, \phi_{N-1})^H \)，从而有 \(Y=FX\)。若令 \( W_N^{n} = e^{-i2\pi \frac nN } \)，则整个矩阵乘式可表示为：
\[ \left(\begin{matrix}y_0\\y_1\\y_2\\\vdots\\y_{N-1}\end{matrix}\right) = \left(\begin{matrix}1&1&1&\dots&1\\1&W_N^1&W_N^2&\dots&W_N^{N-1}\\1&W_N^2&W_N^4&\dots&W_N^{2(N-1)}\\\vdots&\vdots&\vdots&\cdots&\vdots\\1&W_N^{N-1}&W_N^{2(N-1)}&\dots&W_N^{(N-1)^2}\end{matrix}\right) \left(\begin{matrix}x_0\\x_1\\x_2\\\vdots\\x_{N-1}\end{matrix}\right) \]
反变换为：
\[ \left(\begin{matrix}x_0\\x_1\\x_2\\\vdots\\x_{N-1}\end{matrix}\right) = \frac1N\left(\begin{matrix}1&1&1&\dots&1\\1&W_N^{-1}&W_N^{-2}&\dots&W_N^{-(N-1)}\\1&W_N^{-2}&W_N^{-4}&\dots&W_N^{-2(N-1)}\\\vdots&\vdots&\vdots&\cdots&\vdots\\1&W_N^{1-N}&W_N^{2(1-N)}&\dots&W_N^{-(1-N)^2}\end{matrix}\right) \left(\begin{matrix}y_0\\y_1\\y_2\\\vdots\\y_{N-1}\end{matrix}\right) \]
若令 \(w = e^{-i2\pi\frac1N}=W_N^1\)，则可将（正）变换表达的更加简单：
\[ \left(\begin{matrix}y_0\\y_1\\y_2\\\vdots\\y_{N-1}\end{matrix}\right) = \left(\begin{matrix}1&1&1&\dots&1\\1&w&w^2&\dots&w^{N-1}\\1&w^2&w^4&\dots&w^{2(N-1)}\\\vdots&\vdots&\vdots&\cdots&\vdots\\1&w^{N-1}&w^{2(N-1)}&\dots&w^{(N-1)(N-1)}\end{matrix}\right) \left(\begin{matrix}x_0\\x_1\\x_2\\\vdots\\x_{N-1}\end{matrix}\right) \]
这个矩阵也被叫做傅里叶矩阵，不过这里我没有写正规化因数 \( \frac{1}{\sqrt{N}} \)。根据情况或惯例，\(w\) 也可取为 \( e^{i2\pi\frac1N} \)。\(w\) 的取值就决定了变换的正逆，如果再加上正规化因数的话，正变换和逆变换只需要调整一下 \(w\) 即可区分。
自然，N 阶阵乘 N 阶列向量的时间复杂度是 \( O(n^2) \)，这一点可以通过朴素乘法实现体现出来：
// 为了方便，这里用一维数组存储 N 阶方阵，x[a][b] == *(x+a*n+b)
void N_m_n(double *res, double *Nm, double*nm, int n)
{
for (int i = 0; i < n; i++)
{
res[i] = 0.0;
for (int j = 0; j < n; j++)
{
res[i] += Nm[i*n+j] * nm[j];
}
}
return;
}这里我们使用信号 \( y = cos(2x) + cos(3x) + cos(5x) \)进行实验。根据奈奎斯特定理，采样角频率至少为两倍最大角频率，即 \( 10 rad/s \)。取采样角频率为 \( 12 rad/s \)，采样区间为 \((0, 2\pi)\)，则采样点个数为 \(N = 12\)。
首先使用 CL 实现：
;; 采样信号
(defun yy-f (x)
(+ (cos (* 2 x))
(cos (* 3 x))
(cos (* 5 x))))
;;傅里叶矩阵
(defun make-fourier (N sgn)
(let ((a (make-array `(,N ,N))))
(do ((i 0 (+ i 1)))
((= i N) a)
(do ((j 0 (+ j 1)))
((= j N))
(setf (aref a i j)
(exp (* sgn #C(0 -2) pi i j (/ 1.0 N))))))))
;;矩阵向量乘法
(defun mat-mul-vec (m v)
(let* ((N (length v))
(vec (make-array N)))
(do ((i 0 (+ i 1)))
((= i N) vec)
(setf (aref vec i) 0.0)
(do ((j 0 (+ j 1)))
((= j N))
(incf (aref vec i)
(* (aref m i j) (aref v j)))))))
(let* ((N 12) ;; 取 12 个采样点
(x (make-array N)))
(loop for i below N
do (setf (aref x i) (yy-f (/ (* i pi 2) N)))) ;; 获取采样信号值
(let ((res1 (mat-mul-vec (make-fourier N 1) x))) ;; 获取变换
(print res1 t) ;; 打印变换
;; 打印比对原值和逆变换得到的值
(print
(map 'vector '-
(map 'vector #'(lambda (x) (* x N)) x) ;乘 N 用来抵消逆变换未乘 1/N
(mat-mul-vec (make-fourier N -1) res1))
t)))
;;为方便观察手动做了一些近似处理
#(#C(0 0);0
#C(0 0);1
#C(6 0);2
#C(6 0);3
#C(0 0);4
#C(6 0);5
#C(0 0);6
#C(6 0);7 = 12 - 5
#C(0 0);8
#C(6 0);9 = 12 - 3
#C(6 0);10 = 12 - 2
#C(0 0));11
;;都是接近 0 的值，说明逆变换与原值相近（毕竟有舍入误差）
#(#C(-3.0695160191385185d-6 1.6852850822087755d-6)
#C(-2.413088492581039d-6 1.6852830708727321d-6)
#C(-1.7509792851200245d-6 1.6852822343317603d-6)
#C(-1.1659329626212411d-6 1.6852817821595204d-6)
#C(-3.4761838731147854d-7 1.6852815497630256d-6)
#C(5.560783016278492d-7 1.6852815206765511d-6)
#C(1.3139656651617315d-7 1.685281628337206d-6)
#C(1.4592174473193609d-6 1.6852819000591794d-6)
#C(1.598376088542303d-6 1.6852824023629203d-6)
#C(2.2046287107002627d-6 1.6852833358663046d-6)
#C(4.087005883945949d-6 1.6852854214532373d-6)
#C(1.0166025507629683d-5 1.6852919591023401d-6))上面的结果符合预期。如果你会 python 的话，下面的代码是给你准备的：
import cmath
def yy_f (x):
return cmath.cos(2*x) + cmath.cos(3*x) + cmath.cos(5*x)
# 嵌套列表推导
def make_fourier(N, sgn):
return [[cmath.exp(-2*cmath.pi*1j/N*i*j) for i in range(0, N)] for j in range(0, N)]
# 矩阵向量乘法运算
def mat_mul_vec (m, v):
a = [0 for i in range(0, len(v))]
for i in range(0, len(v)):
for j in range(0, len(v)):
a[i] += m[i][j] * v[j]
return a
# 信号个数
N = 12
# 获取采样信号
x = [yy_f(cmath.pi * 2 * i / N) for i in range(0, N)]
# 获得 DFT 结果
y = mat_mul_vec(make_fourier(N, 1), x)
# 进行 IDFT，即 DFT 逆变换
x1 = mat_mul_vec(make_fourier(N, -1), y)
x1 = list(map(lambda x: x/N, x1))
# 打印 DFT 结果
for i in range(0, N):
print(y[i])
print ('--------')
# 打印 IDFT 与原信号的差值
for i in range(0, N):
print(x[i] - x1[i])
# 结果数据，与 CL 一致（未手工近似处理）
(-5.329070518200751e-15+0j)
(-4.274358644806853e-15-1.609823385706477e-15j)
(6.000000000000001-6.994405055138486e-15j)
(6.0000000000000036+6.106226635438361e-16j)
(-3.0253577421035516e-15+0j)
(6.000000000000004-7.438494264988549e-15j)
(6.217248937900877e-15-6.88324302725924e-15j)
(6.000000000000008+1.2323475573339238e-14j)
(3.885780586188048e-15-3.3306690738754696e-15j)
(5.999999999999999+8.93729534823251e-15j)
(6+1.715294573045867e-14j)
(-1.4876988529977098e-14-1.0297318553398327e-14j)
--------
-2.0586549626730855e-16j
(6.6058269965196814e-15+5.488559255325794e-16j)
(4.440892098500626e-16+4.019214628772775e-15j)
(7.771561172376096e-16+2.0131296001480145e-15j)
(-7.358223809900208e-15+1.458762778231005e-15j)
(-1.1102230246251565e-15-2.4917462186266056e-15j)
(-5.551115123125783e-16+6.17037270246773e-17j)
(3.1086244689504383e-15-7.728298151441925e-15j)
(9.079069498069204e-15+1.5319004943530082e-15j)
(1.5543122344752192e-15+3.1805645120342776e-15j)
(-2.3314683517128287e-15+4.1663672132031236e-15j)
(-5.551115123125783e-15-3.9483364539879407e-16j)这里有一篇文章可以解释出现上面数据的原理：如何通俗地解释什么是离散傅里叶变换。
FFT，即“快速傅里叶变换”，它可以将时间复杂度为 \(O(n^2)\) 的 DFT 变为 \(O(nlog(n))\)。需要说明的是它并不是一个算法，而是一类算法。我这里介绍的只是一种思路而已，而且只给出了最简实现。
当看到 \(O(nlog(n))\) 时，我最先想到的就是归并排序，想必 FFT 也用到了分治的思路。分治的思想在于如何把问题分解为多个更小的子问题，且原问题可以由子问题的合并解决。
FFT 能进行的根本也许就在于傅里叶矩阵不是个普通的矩阵，在于 \(e^{-2i\pi k\frac nN}\) 的良好性质。当 N 为偶数时，\(y_k = \sum\limits_{i=0}^{N-1}W_N^{ki}x_{n}\) 可化为：
\[\begin{align} y_k &= \sum\limits_{i=0}^{\frac N2-1}W_N^{2ki}x_{2i} + \sum\limits_{i=0}^{\frac N2 -1}W_N^{k(2i+1)}x_{2i+1} \notag \\ &= \sum\limits_{i=0}^{\frac N2-1}W_{\frac N2}^{ki}x_{2i} + W_N^k \sum\limits_{i=0}^{\frac N2-1}W_{\frac N2}^{ki}x_{2i+1} \notag \end{align} \quad k=0, 1, 2, \dots, N - 1\]
据此易知： \(Y_{half} = \left(\begin{matrix}I_{\frac N2}&D_{\frac N2}^N\end{matrix}\right) \left(\begin{matrix}F_{\frac N2}&0\\0&F_{\frac N2}\end{matrix}\right)\left(\begin{matrix}X_{even}\\X_{odd}\end{matrix}\right)\)，其中
\[D_n^N= \left(\begin{matrix}1\\&W_N^1\\&&W_N^2\\&&&\ddots\\&&&&W_N^{n-1}\end{matrix}\right)\]
\(Y_{half}\) 表示 Y 向量的前一半。
接着，由 \(W_N^{k + \frac N2} = e^{-2i\pi k\frac1N -i\pi}\ = -W_N^k\) 和 \(W_N^{k+N} = W_N^k\)，从而有
\[y_{k+\frac N2} = \sum\limits_{i=0}^{\frac N2-1}W_{\frac N2}^{ki}x_{2i} - W_N^k \sum\limits_{i=0}^{\frac N2-1}W_{\frac N2}^{ki}x_{2i+1}\]
由上式体现出的 \(y_{k+\frac N2}\) 和 \(y_k\) 之间的关系，我们可以补完向量 Y 的另一半：
\[Y = \left(\begin{matrix}I_{\frac N2}&D_{\frac N2}^N\\I_{\frac N2}&-D_{\frac N2}^N\end{matrix}\right) \left(\begin{matrix}F_{\frac N2}&0\\0&F_{\frac N2}\end{matrix}\right)\left(\begin{matrix}X_{even}\\X_{odd}\end{matrix}\right)\]
乘掉最左边的形式主义矩阵，可以得到：
\(y_k = y_{k(even)} + W_N^ky_{k(odd)}, \quad k=0, 1, \dots, \frac N2-1 \)
\(y_{k+\frac N2} = y_{k(even)} - W_N^ky_{k(odd)}, \quad k=0, 1, \dots, \frac N2 -1\)
\(Y_{(even)} = F_{\frac N2}X_{even}, \quad Y_{(odd)} = F_{\frac N2}X_{odd} \)
这应该就是最直接的递归公式了。另外，拿掉奇偶分组的 X 向量，我们就得到了 \(F_{2n}\) 与 \(F_n\) 之间的关系：
\[F_{2n} = \left(\begin{matrix}I_n&D_n^{2n}\\I_n&-D_n^{2n}\end{matrix}\right) \left(\begin{matrix}F_n&0\\0&F_n\end{matrix}\right)P_{2n}\]
其中：
\[P_{2n}=\left(\begin{matrix}1\\&&1\\&&&&1\\&&&&&&\ddots\\&&&&&&&&1&\\&1\\&&&1\\&&&&&1\\&&&&&&&\ddots\\&&&&&&&&&1\end{matrix}\right) P_{ij} = \begin{cases}1 ,j = 2i + 1, i \le n \\ 1, j=2(i-n), i\ge n+1\\ 0 \ otherwise\end{cases} \]
傅里叶矩阵能够分解应该是快速傅里叶变换可行的根本原因。
我没有系统地学习过计算理论，所以无法给出一个严谨的分析过程，这里仅对采样点个数为 \(n=2^k \quad k=1, 2, 3, \dots\) 的情况进行分析。
假设 n = 64，根据上面的公式，共有 32 个来自偶数项的 FFT 结果和 32 个来自奇数项的 FFT 结果需要组合得到最终结果，要得到 64 个 y 值就需要 64 次（也许更多）运算。同理，对于 n = 32，从 n = 16 的两个奇偶 FFT 结果合并到 n = 32 需要 32 次运算。
这也就是说，从子节点向父节点合并时，所需要的计算次数仅与父节点的 n 值有关，顶层点 n 值为 64，第二层所有父节点 n 值之和为 \(32 + 32\)，第三层为 \(16 \cdot 4\) …… 每一层所有父节点的 n 值之和都相同，故每一层的合并操作都需要 64 次操作。
由二分特性，把 64 分到不可再分只需要 6 次，其时间复杂度自然就是 \(O(nlog(n))\)。
上面我给出了对信号 \(y=cos(2x)+cos(3x)+cos(5x)\) 的 DFT 实现，并且我们取 \(N=12\) 来作为区间采样点数。根据推导出的递归公式，我们可以很自然地写出 FFT 的递归实现。
整个算法描述大致如下：
（仅针对 2 次幂个数点）
- X 为待变换向量，它的长度是 N
- fft(X, N) 为 FFT 变换函数，它会返回变换后的向量 Y
- W(n, N) 是虚指数函数 e^(-2*j*pi*n/N)
fft(X,N):
if N is 1 then
return X
else
X_even <- X 中的偶数项（序号从 0 开始）
X_odd <- X 中的奇数项
even_f <- fft(X_even, N/2)
odd_f <- fft(X_odd, N/2)
result <- 长度为 N 的空向量
for i <- 0; i < N/2; i++
a <- even_f[i]
b <- W(i, N) * odd_f[i]
result[i] <- a + b
result[i + N/2] <- a - b
return result可以看到，这是个非常简单的递归实现。我在上面给出的递推公式无法处理奇数阶的情况，下面的具体实现仅对 2 的幂次阶有效。因此下面我将采样点个数设为 16 而非 12，16 是 2 的 4 次幂。
(defun yy-f (x)
(+ (cos (* 2 x))
(cos (* 3 x))
(cos (* 5 x))))
(defun yy-fft (v N)
(if (= N 1) v
(let* ((half (/ N 2))
(evens (make-array half))
(odds (make-array half)))
(do ((k 0 (+ k 1))
(i 0 (+ i 2)))
((= i N))
(setf (aref evens k) (aref v i)
(aref odds k) (aref v (+ i 1))))
(let ((evens-r (yy-fft evens half))
(odds-r (yy-fft odds half))
(res (make-array N)))
(loop for i below half
do (let ((a (aref evens-r i))
(b (* (exp (* #C(0 -2) pi i (/ 1.0 N)))
(aref odds-r i))))
(setf (aref res i) (+ a b)
(aref res (+ i half)) (- a b))))
res))))
(let* ((N 16) ;; 取 16 个采样点
(x (make-array N)))
(loop for i below N
do (setf (aref x i) (yy-f (/ (* i pi 2) N)))) ;; 获取采样信号值
(yy-fft x N))
;; result
#(#C(-7.105427357601002d-15 0.0d0)
#C(-4.379167817162229d-15 -1.990128627232129d-15)
#C(7.9999999999999964d0 -5.10702591327572d-15)
#C(7.999999999999999d0 -4.595312355070685d-15)
#C(-1.7763568394002505d-15 -1.3322676295501878d-15)
#C(7.999999999999997d0 -1.1034605897896593d-14)
#C(6.217248937900877d-15 -6.217248937900877d-15)
#C(7.931881495962731d-15 8.964512367932808d-16)
#C(1.7763568394002505d-15 0.0d0)
#C(7.93188149596273d-15 -1.3405404466433406d-15)
#C(6.661338147750939d-15 5.995204332975845d-15)
#C(7.999999999999997d0 1.1478695107746655d-14)
#C(-1.7763568394002505d-15 1.3322676295501878d-15)
#C(7.999999999999999d0 5.0394015649207474d-15)
#C(7.999999999999997d0 5.329070518200751d-15)
#C(-4.37916781716223d-15 1.5460394173820636d-15))以下是 Python 实现：
import cmath
def yy_f (x):
return cmath.cos(2*x) + cmath.cos(3*x) + cmath.cos(5*x)
def yy_fft (v, N):
if N == 1:
return [v[0]]
else:
half = N // 2;
# 获取奇偶向量
evens = [v[i] for i in range(0, N, 2)]
odds = [v[i] for i in range(1, N, 2)]
# 子 fft
evens_r = yy_fft(evens, half)
odds_r = yy_fft(odds, half)
# 结果向量
res = [0 for i in range(0, N)]
for i in range(0, half):
# 递归公式
a = evens_r[i]
b = cmath.exp(-2j * cmath.pi * i / N) * odds_r[i]
res[i] = a + b
res[i + half] = a - b
return res
# 获取信号
M = 16
x = [yy_f(cmath.pi * 2 * i / M) for i in range(0, M)]
r = yy_fft(x, M)
for i in range(0, M):
print(r[i])
# result:
(-6.8833827526759706e-15+0j)
(-4.294195025611368e-15-1.7849861614373374e-15j)
(7.9999999999999964-4.6629367034256575e-15j)
(7.999999999999999-4.6802851466215464e-15j)
(-1.7763568394002505e-15-1.5543122344752192e-15j)
(7.999999999999997-1.1119578689447454e-14j)
(7.105427357601002e-15-6.217248937900877e-15j)
(7.84690870441187e-15+1.101593702588072e-15j)
(1.5543122344752192e-15+0j)
(7.84690870441187e-15-1.5456829124381322e-15j)
(6.661338147750939e-15+5.551115123125783e-15j)
(7.999999999999997+1.1563667899297517e-14j)
(-1.7763568394002505e-15+1.5543122344752192e-15j)
(7.999999999999999+5.124374356471609e-15j)
(7.9999999999999964+5.329070518200751e-15j)
(-4.294195025611368e-15+1.3408969515872725e-15j)如你所见，本文仅仅介绍了 FFT 的最简单的一种算法，这也是我写这篇文章的目的所在，即学习 FFT 的基本原理。如果有时间的话，后续我会介绍非 2 次幂个数点的处理方式，以及各种各样的 FFT 算法。时间精力所限，本文就到这里了。
关于傅里叶变换和 FFT 的文章可谓不可胜数，既然已经有了这么多好文章，那我还有什么写的必要呢？好文章是普适的，但是它的编排和范围不可能适合每一个人，对于比较重要的知识，进行一定的整理是必要的。要说的话，写文章也算是自我实现的一种方式，它可以作为我曾经存在于网上的记录（这说法好怪）。
以下是我在写这篇文章中发现的不错的文章，由于太过麻烦，下面全部都是中文内容（或自带中文翻译）。