信号处理(5) 快速Fourier变换

回忆DFT定义：

$$X(k)=\sum\limits_{n=0}^{N-1}x(n)W_N^{nk}$$

则计算每个$X(k)$，需要$N$次复数乘法、$N-1$次复数加法，总复杂度为$\mathcal{O}(N^2)$

整体复杂度较高，可以采用分治的做法进行优化

快速Fourier变换 (FFT)

基础性质

$W_N$本质是第$-1$个$N$次单位根：

$$W_N=e^{-j\frac{2\pi}{N}}$$

因此，$W_N$有一些显然的性质：

$$W_N^{nk+N}=W_N^{nk}$$

$$W_N^{nk+N/2}=-W_N^{nk}$$

原理

设$N$是偶数，采用分治策略，试图将$N$阶DFT转换为$N/2$阶DFT

于是，将奇数项的计算和偶数项的计算分开：

$$X(k)=\sum\limits_{r=0}^{N/2-1}x(2r)W_N^{(2r)k}+\sum\limits_{r=0}^{N/2-1}x(2r+1)W_N^{(2r+1)k}$$

$$=\sum\limits_{r=0}^{N/2-1}x(2r)W_{N/2}^{rk}+W_N^k\sum\limits_{r=0}^{N/2-1}x(2r+1)W_{N/2}^{rk}$$

令偶数项序列为$g(r)=x(2r)$，奇数项子序列为$h(r)=x(2r+1)$，则

$$X(k)=G(k)+W_N^kH(k)$$

注意，其中$G(k)$和$H(k)$的周期为$N/2$，对$k\geq N/2$的项需要进行周期延拓

于是，我们把$N$阶DFT分为了两个$N/2$阶DFT计算，这一过程是可以递归使用的，之后对其加权求和，加权求和步骤的复杂度为$\mathcal{O}(N)$

则总复杂度的递推式为

$$T(N)=2T\left(\frac{N}{2}\right)+\mathcal{O}(N)$$

根据主定理，有

$$T(N)=\mathcal{O}(N\log N)$$

矩阵解读

考虑一个$4$阶的Fourier矩阵：

$$\mathbf{F}_4=\begin{bmatrix} 1 & 1 & 1 & 1\\ 1 & -j & -1 & j\\ 1 & -1 & 1 & -1\\ 1 & j & -1 & -j \end{bmatrix}$$

交换中间两列，使偶数列和奇数列分离：

$$\mathbf{F}_4=\begin{bmatrix} 1&1&1&1\\ 1&-1&-j&j\\ 1&1&-1&-1\\ 1&-1&j&-j \end{bmatrix}\begin{bmatrix} 1 & & &\\ &&1&\\ &1&&\\ &&&1 \end{bmatrix}$$

$$=\begin{bmatrix} \mathbf{F}_2 & \begin{bmatrix}1&\\&-j\end{bmatrix}\mathbf{F}_2\\ \mathbf{F}_2 & \begin{bmatrix}-1&\\&j\end{bmatrix}\mathbf{F}_2 \end{bmatrix}\begin{bmatrix} 1 & & &\\ &&1&\\ &1&&\\ &&&1 \end{bmatrix}$$

$$=\begin{bmatrix}1&&1&\\&1&&-j\\1&&-1&\\&1&&j\end{bmatrix} \begin{bmatrix}\mathbf{F}_2&\\&\mathbf{F}_2\end{bmatrix}\begin{bmatrix} 1 & & &\\ &&1&\\ &1&&\\ &&&1 \end{bmatrix}$$

$$=\begin{bmatrix}\mathbf{I}_2&\mathbf{D}_2\\\mathbf{I}_2&-\mathbf{D}_2\end{bmatrix} \begin{bmatrix}\mathbf{F}_2&\\&\mathbf{F}_2\end{bmatrix} \mathbf{P}_4$$

于是，我们把$\mathbf{F}_4$拆成了两个$\mathbf{F}_2$和一些简单矩阵

我们推广到$\mathbf{F}_N$，$N$为偶数的情形，则有：

$$\mathbf{F}_N=\begin{bmatrix} \mathbf{I}_{N/2}&\mathbf{D}_{N/2}\\ \mathbf{I}_{N/2}&-\mathbf{D}_{N/2} \end{bmatrix} \begin{bmatrix} \mathbf{F}_{N/2}&\\&\mathbf{F}_{N/2} \end{bmatrix} \mathbf{P}_N$$

---

下面我们将上面的矩阵式解读为算法

$\mathbf{P}_N$是一个置换矩阵，它将所有偶数元素排列在前方，所有奇数元素排列在后方，即

$$\mathbf{P}_N\mathbf{x}_N= \begin{bmatrix} \mathbf{g}_{N/2}\\\mathbf{h}_{N/2} \end{bmatrix}$$

这一步实际上不产生任何运算量（我们并不需要真的在内存上交换奇数项和偶数项的位置）

之后乘以右数第二个矩阵：

$$\begin{bmatrix}\mathbf{F}_{N/2}&\\&\mathbf{F}_{N/2}\end{bmatrix}\mathbf{P}_N\mathbf{x}_N= \begin{bmatrix}\mathbf{G}_{N/2}\\\mathbf{H}_{N/2}\end{bmatrix}$$

这一步为分别对$g(n)$和$h(n)$进行DFT

再乘以最左侧的矩阵，则

$$\mathbf{X}_N=\begin{bmatrix} \mathbf{I}_{N/2}&\mathbf{D}_{N/2}\\ \mathbf{I}_{N/2}&-\mathbf{D}_{N/2} \end{bmatrix}\begin{bmatrix}\mathbf{G}_{N/2}\\\mathbf{H}_{N/2}\end{bmatrix}=\begin{bmatrix} \mathbf{G}_{N/2}+\mathbf{D}_{N/2}\mathbf{H}_{N/2}\\ \mathbf{G}_{N/2}-\mathbf{D}_{N/2}\mathbf{H}_{N/2} \end{bmatrix}$$

其中

$$\mathbf{D}_{N/2}=\begin{bmatrix} W_N^{0} &&&&\\ & W_N^{1} &&&\\ && W_N^2 &&\\ &&& \ddots &\\ &&&& W_N^{N/2-1} \end{bmatrix}$$

由于

$$W_N^{k+N/2}=W_N^{k}$$

上面的矩阵式实际上和之前的非矩阵式等价