cooly-Turkey FFT详解

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

---

前言

`这篇文章是写给自己的备忘录，最近需要用FPGA ZU5ev来实现毫米波雷达FFT算法，发现以前对FFT的理解都在遗忘。所以用此篇文章记录下来，以备自己以后查看

一、离散傅里叶变换

​ DFT是时域输入向量x(n)到输出向量X(K)的线性频域变化，其公式如下：
$$X[k]=\sum _{n=0}^{N-1}x[n]e^{-j(2\pi \frac{nk}{N})}$$
向量x(n)和X(K)的长度均为N，k=0,1,…,N-1，也称之为进行了N点DFT，其中$W_N=e^{\frac{-j2\pi }{N}}$为旋转因子(Twiddle factor)，则DFT公式也可以重写为：
$$X[k]=\sum _{n=0}^{N-1}x[n]W_N^{nk}$$
为了分析方便，本例中采用8点FFT为例，采用矩阵形式，公式如下：
$$\left(\begin{array}{c}X[0]\\ X[1]\\ X[2]\\ X[3]\\ X[4]\\ X[5]\\ X[6]\\ X[7]\end{array}\right)=\left[\begin{array}{cccccccc}1& 1& 1& 1& 1& 1& 1& 1\\ 1& W_8& W_8^2& W_8^3& W_8^4& W_8^5& W_8^6& W_8^7\\ 1& W_8^2& W_8^4& W_8^6& W_8^8& W_8^{10}& W_8^{12}& W_8^{14}\\ 1& W_8^3& W_8^6& W_8^9& W_8^{12}& W_8^{15}& W_8^{18}& W_8^{21}\\ 1& W_8^4& W_8^8& W_8^{12}& W_8^{16}& W_8^{20}& W_8^{24}& W_8^{28}\\ 1& W_8^5& W_8^{10}& W_8^{15}& W_8^{20}& W_8^{25}& W_8^{30}& W_8^{35}\\ 1& W_8^6& W_8^{12}& W_8^{18}& W_8^{24}& W_8^{30}& W_8^{36}& W_8^{42}\\ 1& W_8^7& W_8^{14}& W_8^{21}& W_8^{28}& W_8^{35}& W_8^{42}& W_8^{49}\end{array}\right]\left(\begin{array}{c}x[0]\\ x[1]\\ x[2]\\ x[3]\\ x[4]\\ x[5]\\ x[6]\\ x[7]\end{array}\right)\text{\hspace{1em}(1)}$$

二、Cooly-Tukey FFT物理意义

我们以8点FFT为例，输入数据为$x_0\dots x_7$，Cooly算法对输入（时域序列）或者输出（频域序列）进行重新分组。在时域序列上重新分组称为时域抽取（DIT），在频域上重新分组称之为频域抽取（DIF），DIT其中一种经典的抽取图如下：

上图中，$x_0\dots x_7$的FFT变换的结果是X(K)，$x_0,x_2,x_4,x_6$的FFT变换结果是$F_1(K)$，$x_1,x_3,x_5,x_7$的FFT变换结果是$F_2(K)$，$x_0，x_4$的FFT变换结果是$G_1(K)$，$x_2，x_6$的FFT变换结果是$G_2(K)$，$x_1，x_5$的FFT变换结果是$G_3(K)$，$x_3，x_7$的FFT变换结果是$G_4(K)$。
我们将$x_0\dots x_7$，按照下标奇偶分组为$x_0,x_2,x_4,x_6$和$x_1,x_2,x_4,x_7$，

$X[k]={\sum }_{n=0}^{\frac{N}{2}-1}x[2n]W_N^{2nk}+{\sum }_{n=0}^{\frac{N}{2}-1}x[2n+1]W_N^{(2n+1)k}=\underset{x_0{x}_2{x}_4{x}_6}{\underbrace{\sum _{n=0}^{\frac{N}{2}-1}x[2n]W_{\frac{N}{2}}^{nk}}}+W_N^k\underset{x_1,x_3,x_5,x_7}{\underbrace{\sum _{n=0}^{\frac{N}{2}-1}x[2n+1]W_{\frac{N}{2}}^{nk}}}k=0,1,...,N-1$，
由于FFT是线性变换，所以上述公式是有物理意义的，假设8点输入序列$x_0\dots x_7$的采样间隔为T，$x_1,x_3,x_5,x_7$相对于$x_0,x_2,x_4,x_6$在时钟上延迟了一个时钟节拍T，所以对于频率K，那么$X[k]=F_1[k]+W_N^k{F}_2[k]k=0,1,...,N-1$中$W_N^k$的物理意义就是频率K延迟一一个采样周期T，由于FFT时线性变化，所以可以X[K]由$F_1[k]$和$F_2[k]$的延迟T的线性相加。
所以输入序列第一次奇偶拆分后的矩阵表示为
$$\left(\begin{array}{c}X[0]\\ X[1]\\ X[2]\\ X[3]\\ X[4]\\ X[5]\\ X[6]\\ X[7]\end{array}\right)=\left[\begin{array}{cccccccc}1& 0& 0& 0& 1& 0& 0& 0\\ 0& 1& 0& 0& 0& W_8& 0& 0\\ 0& 0& 1& 0& 0& 0& W_8^2& 0\\ 0& 0& 0& 1& 0& 0& 0& W_8^3\\ 1& 0& 0& 0& -1& 0& 0& 0\\ 0& 1& 0& 0& 0& -W_8& 0& 0\\ 0& 0& 1& 0& 0& 0& -W_8^2& 0\\ 0& 0& 0& 1& 0& 0& 0& -W_8^3\end{array}\right]\left(\begin{array}{c}{F}_1[0]\\ F_1[1]\\ F_1[2]\\ F_1[3]\\ F_2[0]\\ F_2[1]\\ F_2[2]\\ F_2[3]\end{array}\right)\text{\hspace{1em}(2)}$$
这里$F_1[k]$和$F_2[k]$的周期都是4，但是X[K]周期时8，大家都是周期循环序列相加，因为输入序列是时域离散的，所以频域里都是周期的。$W_8^7=W_8^4\ast W_8^3=-W_8^3$，可以有效解释上述矩阵。
接着我们对$x_0,x_2,x_4,x_6$继续按照奇偶拆分为$x_0,x_4$和$x_2,x_6$，对于$x_0,x_4$我们在物理意义上需要理解为一个两个点输入序列，采样间隔为（4-0）=4T的输入序列。矩阵表达式如下
$$\left(\begin{array}{c}{F}_1[0]\\ F_1[1]\\ F_1[2]\\ F_1[3]\end{array}\right)=\left[\begin{array}{cccc}1& 0& 1& 0\\ 0& 1& 0& W_8^2\\ 1& 0& -1& 0\\ 0& 1& 0& -W_8^2\end{array}\right]\left(\begin{array}{c}{G}_1[0]\\ G_1[1]\\ G_2[0]\\ G_2[1]\end{array}\right)\text{\hspace{1em}(3)}$$
我们详细解释一下上述矩阵表达式的物理意义，由于$G_1[K]$和$G_2[K]$是周期为2的循环序列，所以$G_1[1]=G_1[3]$，其中$x_2,x_6$相对于$x_0,x_4$在时域上延迟了2T，所以$F_1[1]=G_1[1]+G_1[1]\ast W_{8T}^{2T}$,所以T抵消掉。对于$F_1[3]=G_1[3]+G_1[3]\ast W_{8T}^{2T\ast 3}=G_1[3]+G_1[3]\ast W_{8T}^{2T\ast 1}\ast W_{8T}^{2T\ast 2}=G_1[3]-G_1[3]\ast W_{8T}^{2T}=G_1[1]-G_1[1]\ast W_{8T}^{2T}$
同理 $F_2$可以用相同方法拆分为$G_3$和$G_4$:
$$\left(\begin{array}{c}{F}_2[0]\\ F_2[1]\\ F_2[2]\\ F_2[3]\end{array}\right)=\left[\begin{array}{cccc}1& 0& 1& 0\\ 0& 1& 0& W_8^2\\ 1& 0& -1& 0\\ 0& 1& 0& -W_8^2\end{array}\right]\left(\begin{array}{c}{G}_3[0]\\ G_3[1]\\ G_4[0]\\ G_4[1]\end{array}\right)\text{\hspace{1em}(4)}$$
所以
$$\left(\begin{array}{c}{F}_1[0]\\ F_1[1]\\ F_1[2]\\ F_1[3]\\ F_2[0]\\ F_2[1]\\ F_2[2]\\ F_2[3]\end{array}\right)=\left[\begin{array}{cccccccc}1& 0& 1& 0& 0& 0& 0& 0\\ 0& 1& 0& W_8^2& 0& 0& 0& 0\\ 1& 0& -1& 0& 0& 0& 0& 0\\ 0& 1& 0& -W_8^2& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 1& 0\\ 0& 0& 0& 0& 0& 1& 0& W_8^2\\ 0& 0& 0& 0& 1& 0& -1& 0\\ 0& 0& 0& 0& 0& 1& 0& -W_8^2\end{array}\right]\left(\begin{array}{c}{G}_1[0]\\ G_1[1]\\ G_2[0]\\ G_2[1]\\ G_3[0]\\ G_3[1]\\ G_4[0]\\ G_4[1]\end{array}\right)\text{\hspace{1em}(5)}$$

三、矩阵公式和蝶形计算对应关系

记住$x_1,x_3,x_5,x_7$相对于$x_0,x_2,x_4,x_6$在时钟上延迟了一个时钟节拍T，所以频域需要补偿$W_8^{KT}$，对所有频率补偿一个周期T的延迟，其次$F_1(K)$和$F_2(K)$都是周期性的，本利中X[K]周期为8，$F_1(K)$和$F_2(K)$周期都为4，通过补偿$W_8^{KT}$由周期4变为周期8，所以蝶形变换中Twiddle因子$W_8^{KT}$将周期限制为4，然后通过$W_8^4=-1$延伸扩展到4-7。语言表达有点绕圈，要是蝶形运算理解困难，直接理解矩阵运算就行。对应我们矩阵运算（2）（3）（4）（5）的蝶形图为：

是不是还是矩阵更好理解，线性代数真是个好工具。

总结

写这篇文章真是费劲，最大的收获时学会了如何在markdown里输入LaTex公式，手动狗头一个！
先写到这，让我吐会。。。。。