基-2DIT-FFT算法的原理及MATLAB简单实现

DFT的矩阵表示
$$\left[\begin{array}{c}X(0)\\ X(1)\\ X(2)\\ ⋮\\ X(k)\\ ⋮\\ X(N-1)\end{array}\right]=[W_N^{ij}]\left[\begin{array}{c}x(0)\\ x(1)\\ x(2)\\ ⋮\\ x(k)\\ ⋮\\ x(N-1)\end{array}\right]$$
忽略复数加法，DFT方法需要N²次复数乘法。

基-2DIT-FFT算法（该算法只适用于$N=2^l$的情况）

将N点序列（$N=2^l$）分解为奇序列和偶序列两个$N/2$序列，分别计算其DFT。
$$\begin{array}{rl}{X}_N(k)& =\sum _{n=0}^{N-1}x(n)W_N^{nk}k=0,1,\dots ,N-1\\ & =\sum _{\text{n=偶整数}}x(n)W_N^{nk}+\sum _{\text{n=奇整数}}x(n)W_N^{nk}\\ & =\sum _{r=0}^{(N/2)-1}x(2r)W_N^{2rk}+\sum _{r=0}^{(N/2)-1}x(2r+1)W_N^{(2r+1)k}\\ & =\sum _{r=0}^{(N/2)-1}x(2r)(W_N^2{)}^{rk}+W_N^k\sum _{r=0}^{(N/2)-1}x(2r+1)(W_N^2{)}^{rk}\\ & =\sum _{r=0}^{(N/2)-1}x(2r)W_{N/2}^{rk}+W_N^k\sum _{r=0}^{(N/2)-1}x(2r+1)W_{N/2}^{rk}\\ & =G_{N/2}(k)+W_N^k{H}_{N/2}(k)\end{array}$$

$$\begin{array}{rl}{X}_N(k+\frac{N}{2})& =G_{N/2}(k)+W_N^{k+\frac{N}{2}}{H}_{N/2}(k)\\ & =G_{N/2}(k)-W_N^k{H}_{N/2}(k)\end{array}$$

方便讨论，记N点DFT计算中的DFT矩阵$[W_N^{ij}]$为$F_N$。

将$F_N$分解，则得到
$$F_N=\left[\begin{array}{cc}{I}_{N/2}& D_{N/2}\\ I_{N/2}& -D_{N/2}\end{array}\right]\left[\begin{array}{cc}{F}_{N/2}& \mathrm{O}\\ \mathrm{O}& F_{N/2}\end{array}\right]P_N$$
其中，
$$D_{N/2}=diag[W_N^k]k=0,1,\dots ,\frac{N}{2}-1$$
$P_N$为置换矩阵，用于将$x(n)$的奇偶序列分开，举例说明，当$N=4$和$N=8$时，$P_4$和$P_8$为
$$P_4=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 1\end{array}\right]$$

$$P_8=\left[\begin{array}{cccccccc}1& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right]$$

方便理解，我写出第一次迭代完整过程。
$$\begin{array}{rl}\left[\begin{array}{c}X(0)\\ X(1)\\ X(2)\\ ⋮\\ X(k)\\ ⋮\\ X(N-1)\end{array}\right]=[W_N^{ij}]\left[\begin{array}{c}x(0)\\ x(1)\\ x(2)\\ ⋮\\ x(k)\\ ⋮\\ x(N-1)\end{array}\right]& =F_{N}x\\ & =\left[\begin{array}{cc}{I}_{N/2}& D_{N/2}\\ I_{N/2}& -D_{N/2}\end{array}\right]\left[\begin{array}{cc}{F}_{N/2}& \mathrm{O}\\ \mathrm{O}& F_{N/2}\end{array}\right]P_{N}x\\ & =\left[\begin{array}{cc}{I}_{N/2}& D_{N/2}\\ I_{N/2}& -D_{N/2}\end{array}\right]\left[\begin{array}{cc}{F}_{N/2}& \mathrm{O}\\ \mathrm{O}& F_{N/2}\end{array}\right]\left[\begin{array}{c}x(0)\\ x(2)\\ ⋮\\ x(N-2)\\ \text{- - - - - - -}\\ x(1)\\ x(3)\\ ⋮\\ x(N-1)\end{array}\right]\\ & =\left[\begin{array}{cc}{I}_{N/2}& D_{N/2}\\ I_{N/2}& -D_{N/2}\end{array}\right]\left[\begin{array}{c}{G}_{N/2}(k)\\ H_{N/2}(k)\end{array}\right]\\ & =\left[\begin{array}{c}{G}_{N/2}(k)+W_N^k{H}_{N/2}(k)\\ G_{N/2}(k)-W_N^k{H}_{N/2}(k)\end{array}\right]\\ & =\left[\begin{array}{c}{X}_N(k)\\ X_N(k+\frac{N}{2})\end{array}\right]k=0,1,\dots ,\frac{N}{2}-1\end{array}$$
$F_{N/2}$还能继续分解，直到分解到$F_2$
$$F_2=\left[\begin{array}{cc}1& 1\\ 1& -1\end{array}\right]$$
最终只需$\frac{N}{2}logN$次复数乘法。

function [Xk] = myfft(xn,N)

if(N == 2)
    Xk = [xn(1)+xn(2);xn(1)-xn(2)];
else
    xk_e = myfft(xn(1:2:N-1),N/2);
    xk_o = myfft(xn(2:2:N),N/2);

    WN = exp(-1i*2*pi/N);
    n = 0:1:N/2-1;
    D = (WN*(ones(1,N/2))).^n;
    tmp = (D.') .* xk_o;
    Xk = [xk_e+tmp;xk_e-tmp];
end

IFFT同理，将$D_N$中的元素改为$W_N^{-k}$，并将最后结果乘以$\frac{1}{N}$。

function [xn] = myifft(Xk,N)

if(N == 2)
    xn = [Xk(1)+Xk(2);Xk(1)-Xk(2)];
else
    xn_e = myifft(Xk(1:2:N-1),N/2);
    xn_o = myifft(Xk(2:2:N),N/2);

    WN = exp(1i*2*pi/N);
    n = 0:1:N/2-1;
    D = (WN*(ones(1,N/2))).^n;
    tmp = (D.') .* xn_o;
    xn = [xn_e+tmp;xn_e-tmp] / N;
end