FFT具体实现原理（FFT硬件结构）

前面文章中讲解的FFT基础知识，从傅里叶级数——DFT——FFT的整个推导过程，详细内容参考：

本文将正式讲解FFT钻石吊坠实现流程（倾向于底层硬件的实现），主要包括：算法的简介、蝶形运算单元的结构、旋转因子的生成、存储结构（输入数据的存储、中间操作数的存储、旋转因子的存储等）、相关的寻址规律和硬件架构。

一、算法的简介输入数据量=输出的数据量。如八点FFT每次输入8个数据，那么对应输出8个数据。输入和输出是倒序关系。这里的倒序指的是二进制的倒序，如下图，输入的顺序为（用前四个数举例）：000、100、010、110（0、4、2、6）那么对应的输出就是000、001、010、011（1、2、3、4）。图1

3. 本次主要讲解DIT的FFT实现原理。

二、蝶形运算单元的结构

通过前面文章中的推导，可以得出实际的每个蝶形运算单元的结构： \begin{aligned} &Y A=X A+W^{*} X B \\ &Y B=X A-W^{*} X B \end{aligned}

图2

上面的公式中，XA和XB是输入数据；YA和YB是输出数据。上面的两个公式就构成了一对基2的DIT蝶形结构。

具体的硬件单元为：

首先XB和W相乘，输出的结构进行加减运算，得出YA和YB。

下面针对蝶形运算做几点相关的说明。

N点基2FFT要经历的级数M。 M=\log _{2} N 如8点FFT需要3级，如图1所示男士卫衣品牌。N点基2FFT要G个蝶形运算结构。 G=\frac{N}{2} \log _{2} N 如8点基2FFT需要3级，每一级需要8/2=4个蝶形结构，共需要8/2*3=12个。需要经历C个复数乘法运算由基2的蝶形结构来看，每个蝶形运算需要一个复数乘法，所以和蝶形运算单元结构数目相同： C=\frac{N}{2} \log _{2} N 需要经历A个复数乘法运算由基2的蝶形结构来看，每个蝶形运算需要一个加法、一个减法运算，所以是蝶形运算单元结构数目的两倍： A=N*log _{2} N 三、旋转因子的生成

旋转因子用 W_{N}^{k} 来表示，k的取值是0~N-1，用公式来表示就是 :

W_{N}^{k}=e^{-j 2 \pi k/ N}=cos(2\pi k/N)-jsin(2\pi k/N)

我们已知N（FFT的点数），k的取值（0~N-1），我们就可以计算出 W_{N}^{k} 的数值了。

下面说明一下蝶形运算单元中用到旋转因子和输入XB的乘法运算：

设 XB=a+bj ， XB*W_{N}^{k}=(a+b邓典j)*[cos(2\pi k/N)-jsin(2\pi k/N)]

把实部和虚部拆开来计算，化简得：

实部： a*cos(2\pi k/N)+b*sin(2\pi k/N)

虚部： b*cos(2\pi k/N)-a*sin(2\pi k/N

单纯从乘法器的数目上来看，计算过程需要四个乘法器，两个加法器，我们可以通过下面的计算方式将乘法器的数目减少到三个，加法器的数目保持两个不变：

实部： (a+b)*cos(2\pi k/N)+b*[sin(2\pi k/N)-cos(2\pi k/N)]

虚部： (a+b)*cos(2\pi k/N)金属雕塑-a*[sin(光纤2\pi k/N)+cos(2\pi k/N)]

四、存储结构1. 操作数和中间计算结果的存储

这中乒乓操作的存储结构是针对蝶形运算产生的中间数据和输出结果进行数据缓冲使用，举个列子：以厦门亲子1024点基2 FFT为例：首先把待操作的1024个数存储到M0中，下个阶段开始向蝶形单元中输入数据，假设我们每次取16个数输入到8个蝶形运算单元中，计算的16个结果我们存储到M1中，同样下次再取16个数输入到8个蝶形运算单元中，计算的16个结果我们存储到M1中......等到一级中的所有运算都结束后，从M1中读取数据送入蝶形运算单元中，M0用来存储计算结果......直到9级蝶形运算全部结束。这样，就可以保证蝶形运算计算结束后，可同时的读取数据和写入数据，使得蝶形运算单元的计算效率充分发挥。

2. 旋转因子的存储

其实，旋转因子的存储根据FFT点数会有所不同，如果做小点数的FFT，如8、16、32、64等，我们可以选择将 W_{N}^{0}-W_{N}^{N-1} 全部存储，这样只需要用到哪个直接读取对应位置的地址；对于中等规模的FFT，会选择存储1/4的旋转因子，即存储 W_{N}^{0}-W_{N}^{(N/4-1)} ，其它的3/4只是sin、cos前面的符号有所不同；对于大规模的FFT，达到数十K的规模时，采用旋转因子的压缩算法，具体实现方法：

（1） 存储1/4的旋转因子

// wn_re = cos(-2pi*n/64) wn_im = sin(-2pi*n/64) n:0--63// 0 1.000 -0.000// 1 0.995 -0.098// 2 0.981 -0.195// 3 0.957 -0.290// 4 0.924 -0.383// 5 0.882 -0.471// 6 0.831 -0.556// 7 0.773 -0.634// 8 0.707 -0.707// 9 0.634 -0.773// 10 0.556 -0.831// 11 0.471 -0.882// 12 0.383 -0.924// 13 0.290 -0.957// 14 0.195 -0.981// 15 0.098 -0.995// 16 0.000 -1.000// 17 -0.098 -0.995// 18 -0.195 注册投资分析师-0.981// 19 -0.290 -0.957// 20 -0.383 -0.924// 21 -0.471 -0.882// 22 -0.556 -0.831// 23 -0.634 -0.7体检报告单73// 24 -0.707 -0.707// 25 -0.773 -0.634// 26 -0.831 -0.556// 27 -0.882 -0.471// 28 -0.924 -0.383// 29 -0.957 -0.290// 30 -0.981 -0.195// 31 -0.995 -0.098// 32 -1.000 -0.000// 33 -0.995 0.098// 34 -0.981 0.195// 35 -0.957 0.290// 36 -0.924 0.383// 37 -0.882 0.471// 38 -0.831 0.556// 39 -0.773 0.634// 40 -0.707 0.707// 41 -0.634 0.773// 42 -0.556 0.831// 43 -0.471 0.882// 44 -0.383 0.924// 45 -0.290 0.957// 46 -0.195 0.981// 47 -0.098 0.995// 48 -0.000 1.000// 49 0.098 0.995// 50 0.195 0.981// 51 0.290 0.95滑雪场7// 52 0.383 0.924// 53 0.471 0.882// 54 0.556 0.831// 55 0.634 0.773// 56 0.707 0.707// 57 0.773 0.634// 58 0.831 0.556// 59 0.882 0.471// 60 0.924 0.383// 61 0.957 0.290// 62 0.981 0.195// 63 0.995 0.098

上面部分64点FFT旋转因子的数值（cos、sin），但是只需要存储16个即可孟军，通过读取的旋转因子的不同地址来区分：

2'd0 : {mx_re, mx_im} <= { tw_im, -tw_re}; //16-312'd1 : {mx_re, mx_im} <= { -tw_im, tw_re}; //48-632'd2 : {mx_re, mx_im} <= {-tw_re, -tw_im}; //32-47

这样可以节省3/4的存储空间，但是通过一个数据选择器来实现，速度不如直接存储全部的旋转因子。

（2）旋转因子的压缩算法

W_{N}^{k}=W_{N}^{(\alpha+\beta)}=W_{N}^{\alpha} * W_{N}^{\beta} （指数函数的运算规则牛仔裤美臀）

通过将指数k分解为 \alpha和\beta 两部分之和继而转换为同底数幂相乘的形式 ，如进行32K点的FFT变换，N=256×128， k=256 \alpha+\beta\ \alpha=0,1, \ldots \ldots, 1恒河猴实验27 ; \beta=0,1, \ldots \l新生儿奶粉排行榜dots, 255 ，这样可以把32K个所有的旋转因子全部遍历，把k带入，可得：

W_{256 \times 128}^{k}=W_{256 \times 128}^{(256 \alpha+\beta)}=W_{256 \times 128}^{256 \alpha} * W_{256 \times 128}^{\beta}=W_{128}^{\alpha} * W_{256 \times 128}^{\beta} （256和256约分是根据旋转因子的性质），这样操作，只需要存储128+256=384的存储空间，但是需要额外开销复数乘法器，如果我们以16并行度向蝶形运算单元输入数据，则需要额外增加16个复数乘法器（复数乘法器的结构和 XB*W_{N}^{k} 的结构相同）。

五、相关的寻址规律1. 旋转因子的寻址规律

还是用 W_{N}^{k} 来代表FFT的旋转因子，N代表FFT的点数，各级之间的旋转因子并不是完全不同，我们用 P来代表旋转因子的级数（M=\log _{2} N ，所以P的取值为0~M-1），则旋转因子的个数为 2^M 个，如8点的FFT旋转因子共有 2^3=8 个，则每一级不同的旋转因子个数为： 2^P 个。如基2 8点FFT第一级数只用到 W_8^0 ，第二级用到 W_8^0、W_8^2 （两个不同），第三级 W_8^0-W_8^3 （四个不同），如果用J来表示 2^P 个旋转因子的排列次序（ J=0,1 \ldots \ldots, 2^{P}-1 ），则：

k=J \cdot 2^{M-P-1}

下面介绍旋转因子的具体寻址方法：

实际上，如果FFT的点数是N点，我们使用到的旋转因子数目为N/2，也就是 W_N^0-W_N^{N/2-1} 。接下来为了区分不同蝶形单元对应的旋转因子，对蝶形单元进行分组， M=\log _{2} N 级运算，第n（ \mathrm{n}=0,1, \cdots \cdots, \quad \mathrm{M}-1 ）级FFT运算共有 2^{M-n-1}=p 组。例如，基2 8点FFT第一级将蝶形运算excel图表分为： 2^{3-0-1}=4 组，第二级： 2^{3-1-1}=2 ，第三级..惠普硒鼓加粉....（由图一也可以看出）。

设基2 8点FFT的旋转因子的地址分别为：00 01 10 11（分别对应 W_8^0-W_8^3 ），称之为选装银子的物理地址。

第0级分为4组，每组一个蝶形运算单元，共需要一个不同的旋转因子，即： W_8^0

第1级分为两组，每一组两个蝶形运算单元，共需要两个不同的旋转因子，两个旋转因子的寻址地址分别为00 10（00 01——物理地址的倒序），即为（ W_8^0、W_8^2 ）

第2级分为一组，每一组四个个蝶形运算单元，共需要四个不同的旋转因子，四个个旋转因子的寻址地址分别为00 10 01 11（00 01 10 11——物理地址的倒序），即为（ W_8^0、W_8^2、W_8^1、W_8^3 ）

2. 操作数、中间结果、最终结果的存储与寻址方式

假设FFT输入倒序，输出正序，下面来介绍原始数据的存储、中间结果的存储、输出结果读取的过程：

首先输入按照正常的物理地址来存储，每次蝶形运算计算的结果仍然存储到输入到蝶形运算中数据的地址，但是在寻址时需要按照蝶形运算的组进行选择。还是以基2 8点FFT为例，存储地址和对应的数据如下：

// 000 x(0)// 001 x(4)// 010 x(2)// 011 x(6)// 100 x(1)// 101 x(5)// 110 x(3)// 111 x(7)

每一级蝶形运算的取址有固定的间隔， M=\log _{2} N 级，用S来代表蝶形运算所在（S=0, 1, ... M-1）的级数，取值间隔D为： \mathrm{D}=2^{S}

如，第0级取值间隔为1，就是按序取值，并依此组合四个蝶形运算单元

// 000 x(0) 001 x(4)// 010 六一儿童节礼物 x(2) 011 x(6)// 100 x(1) 101 x(5)// 110 x(3) 111黄玉价格 x(7)

第1级的取值间隔为2，并依此组合四个微信小程序二维码蝶形运算单元

// 000 y(0) // 001 y(4)// 010 y(2)// 011 y(6)// 100 y(1)// 101 y(5)// 110 y(3)// 111 y(7)// first group// 000 y(0) 010 y(2)// 001 y(4) 01援交天使1 y(6)// second group// 100 x(1) 110 x(3)// 101 y(5) 111 y(7)

第三级.....

至此，完全可以看武侠故事懂蝶形运算单元的数据流向的，真是清清爽爽！

六、硬件架构

具体实现方式请参考这篇文章。

后续的更新请关注：

具体设计方案可加wx：Jx_G_zhihu获取