通过矩阵从整体角度搞懂快速傅里叶变换原理

离散傅里叶变换公式

公式

$$f[k]=\sum _{n=0}^{N-1}g[n]e^{-i(2\pi /N)kn},其中(0<=n<N)$$

逆变换公式

$$g[n]=\frac{1}{N}\sum _{k=0}^{N-1}f[k]e^{i(2\pi /N)kn},其中(0<=k<N)$$

快速傅里叶变换

从以上公式看，如果直接按照公式来求离散傅里叶变换，其时间复杂度是O(N^2)

快速傅里叶变换就是一种能在O(n*log(n))时间复杂度内进行傅里叶变换及其逆变换的算法

离散傅里叶变换公式矩阵表示

令

$$\begin{gathered} G=\left[\begin{array}{c}g[0]\\ g[1]\\ ⋮\\ g[n-1]\end{array}\right] \\ F=\left[\begin{array}{c}f[0]\\ f[1]\\ ⋮\\ f[n-1]\end{array}\right] \end{gathered}$$

$$E[i]=\left[\begin{array}{cccc}{e}^{-i(2\pi /N)\ast 0\ast 0}& e^{-i(2\pi /N)\ast 0\ast 1}& \cdots & e^{-i(2\pi /N)\ast 0\ast (n-1)}\\ e^{-i(2\pi /N)\ast 1\ast 0}& 95& \cdots & e^{-i(2\pi /N)\ast 1\ast (n-1)}\\ ⋮& ⋮& \ddots & \ddots \\ ⋮& ⋮& \ddots & e^{-i(2\pi /N)\ast (n-1)\ast (n-1)}\end{array}\right]$$

则离散傅里叶公式可改写为

$$F=E[i]\ast G$$

逆变换可改写为

$$G=\frac{1}{N}E[-i]\ast F$$

注意: E[i]里的 i 是一个变量，i 即正虚数单位，代表正变换，-i 代表逆变换，下同。

递归求解

令

$$\begin{gathered} E[i][n]=\left[\begin{array}{cccc}{e}^{-i(2\pi /N)\ast k\ast 0}& e^{-i(2\pi /N)\ast k\ast 1}& \dots & e^{-i(2\pi /N)\ast k\ast (n-1)}\end{array}\right] \\ E[i][n{]}^{0n}=\left[\begin{array}{cccc}{e}^{-i(2\pi /N)\ast k\ast 0}& e^{-i(2\pi /N)\ast k\ast 2}& \dots & e^{-i(2\pi /N)\ast k\ast (n-2)}\end{array}\right] \\ E[i][n{]}^{1n}=\left[\begin{array}{cccc}{e}^{-i(2\pi /N)\ast k\ast 1}& e^{-i(2\pi /N)\ast k\ast 3}& \dots & e^{-i(2\pi /N)\ast k\ast (n-1)}\end{array}\right] \end{gathered}$$

则

$$E[i]=\left[\begin{array}{cc}E[0][n{]}^{0n}& E[0][n{]}^{1n}\\ ⋮& ⋮\\ E[n-1][n{]}^{0n}& E[n-1][n{]}^{1n}\end{array}\right]$$

将E[i]矩阵竖着再切一刀，平分两半，用数学语言描述就是，

令

$$\begin{gathered} {E}_{00(n/2)}=\left[\begin{array}{c}E[0][n{]}^{0n}\\ ⋮\\ E[n/2-1][n{]}^{0n}\end{array}\right] \\ {E}_{01(n/2)}=\left[\begin{array}{c}E[0][n{]}^{1n}\\ ⋮\\ E[n/2-1][n{]}^{1n}\end{array}\right] \\ {E}_{10(n/2)}=\left[\begin{array}{c}E[n/2][n{]}^{0n}\\ ⋮\\ E[n-1][n{]}^{0n}\end{array}\right] \\ {E}_{11(n/2)}=\left[\begin{array}{c}E[n/2][n{]}^{1n}\\ ⋮\\ E[n-1][n{]}^{1n}\end{array}\right] \end{gathered}$$

于是

$$E[i]=\left[\begin{array}{cc}{E}_{00(n/2)}& E_{01(n/2)}\\ E_{10(n/2)}& E_{11(n/2)}\end{array}\right]$$

再令

$$w[k]=e^{-i(2\pi /N)\ast k}$$

$$\begin{gathered} W^{0(n/2)}=\left[\begin{array}{ccccc}w[0]& 0& 0& \dots & 0\\ 0& w[1]& 0& \dots & 0\\ 0& 0& w[2]& \dots & 0\\ ⋮& ⋮& ⋮& ⋮& 0\\ 0& 0& 0& \dots & w[n/2-1]\end{array}\right] \\ {W}^{1(n/2)}=\left[\begin{array}{ccccc}w[n/2]& 0& 0& \dots & 0\\ 0& w[n/2+1]& 0& \dots & 0\\ 0& 0& w[n/2+2]& \dots & 0\\ ⋮& ⋮& ⋮& ⋮& 0\\ 0& 0& 0& \dots & w[n-1]\end{array}\right] \end{gathered}$$

于是

$$\begin{gathered} E[i]=\left[\begin{array}{cc}{E}_{00(n/2)}& E_{01(n/2)}\\ E_{10(n/2)}& E_{11(n/2)}\end{array}\right] \\ =\left[\begin{array}{cc}{E}_{00(n/2)}& W^{0(n/2)}\ast E_{00(n/2)}\\ E_{10(n/2)}& W^{1(n/2)}\ast E_{10(n/2)}\end{array}\right] \\ =\left[\begin{array}{cc}{E}_{00(n/2)}& W^{0(n/2)}\ast E_{00(n/2)}\\ -E_{00(n/2)}& -W^{1(n/2)}\ast E_{00(n/2)}\end{array}\right] \\ =\left[\begin{array}{cc}{E}_{00(n/2)}& W^{0(n/2)}\ast E_{00(n/2)}\\ -E_{00(n/2)}& W^{0(n/2)}\ast E_{00(n/2)}\end{array}\right] \end{gathered}$$
上面这一步就是整个转化的核心了。

再往后，令

$$\begin{gathered} G^{0(n/2)}=\left[\begin{array}{c}g[0]\\ g[2]\\ ⋮\\ g[n-2]\end{array}\right] \\ {G}^{1(n/2)}=\left[\begin{array}{c}g[1]\\ g[3]\\ ⋮\\ g[n-1]\end{array}\right] \end{gathered}$$

则

$$\begin{gathered} F=GE[i] \\ F=\left[\begin{array}{cc}{E}_{00(n/2)}& W^{0(n/2)}\ast E_{00(n/2)}\\ -E_{00(n/2)}& W^{0(n/2)}\ast E_{00(n/2)}\end{array}\right]\ast \left[\begin{array}{c}{G}^{0(n/2)}\\ G^{1(n/2)}\end{array}\right] \end{gathered}$$
到这一步，可以看到，我们已经成功将一个 n 阶的离散傅里叶变换降为了 n/2 阶，到这里就实现了O(n*logn)时间复杂度的快速傅里叶算法。
逆变换的和正变换的大同小异，就不在赘述了。

如果想要更简洁的形式，更深刻的理解离散傅里叶变换，可以进一步化简。
令

$$\begin{gathered} M=E_{00(n/2)}\ast G^{0(n/2)} \\ N=E_{00(n/2)}\ast G^{1(n/2)} \\ Z=W^{0(n/2)} \end{gathered}$$

则

$$F=\left[\begin{array}{c}M+Z\ast N\\ -M+Z\ast N\end{array}\right]$$

最后，上代码

public class FFT {/*** 傅里叶变换* @param x* @return*/public static Complex[] fft(Complex[] x){return fftTrans(x, true);}/*** 逆傅里叶变换* @param x* @return*/public static Complex[] ifft(Complex[] x) {int N = x.length;Complex[] y = fftTrans(x, false);for (int n = 0; n < N; n++) {y[n] = y[n].divides(N);}return y;}public static Complex[] fftTrans(Complex[] x, boolean dire) {int N = x.length;Complex[] y = new Complex[N];if (N == 1) {y[0] = x[0];return y;}Complex[] x0 = new Complex[N/2];for (int n = 0; n < N; n += 2) {x0[n/2] = x[n];}Complex[] X0 = fftTrans(x0, dire);Complex[] x1 = new Complex[N/2];for (int n = 1; n < N; n += 2) {x1[n/2] = x[n];}Complex[] X1 = fftTrans(x1, dire);for (int k = 0; k < N / 2; k++) {int ci = -1;if (!dire) {ci=-ci;}Complex wnk = getWnk(N, k, ci);Complex wnkX1 = wnk.multiply(X1[k]);y[k] = X0[k].plus(wnkX1);y[k+N/2] = X0[k].minus(wnkX1);}return y;}private static Complex getWnk(int N, int k, int n) {double T = 2 * Math.PI / N;double kth = T * k * n;Complex wk = new Complex(Math.cos(kth), Math.sin(kth));return wk;}}