基2和基4FFT

1.1 FFT的必要索引变换

基2算法需要位顺序的反转位逆序，而基4算法需要首先构成一个2位的数字，再反转这些数字，称为数字逆序。

1.1 位逆序和数字逆序

1.2 FFT的复数乘法转实数乘法

\[X(K)=∑_{n=0}^{N-1}x(n)W_N^{kn}, k=0,1,2…..N-1 \]

\[=∑_{n=0}^{N-1}(Re[x(n)]+jIm[x(n)] )∙(Re[W_N^{kn} ]+jIm[W_N^{kn}]) \]

\[=∑_{n=0}^{N-1}(Re[x(n)]Re[W_N^{kn}]-Im[x(n)]Im[W_N^{kn}])+j(Re[x(n)]Im[W_N^{kn} ]+Im[x(n)]Re[W_N^{kn}]) \]

由式知，计算一个X(K)需要N次复数乘法和N-1次复数加法，计算N点X(K)需要\(N^2\)次复数乘法和N(N-1)次复数加法，计算一次复数乘法需要4次实数乘法和2次实数加法，一次复数的加法需要2次实数的加法。所以每计算一个X(K)需要4N次实数乘法和2N+2(N-1)=2(2N-1)次实数加法, 计算N点X(K)需要\({4N}^2\)次实数乘法和2N(2N-1)次复实加法。

3.基2FFT算法

设输入序列长度为\(N=2^M\)（M为整数），将该序列按照时间顺序的奇偶分解为越来越短的子序列，称为基2按时间抽取的FFT算法。也成为Coolkey-Tukey算法。其中即2表示：\(N=2^M\),M为整数。若不满足这个条件，可以认为地加上若干零值（加零补长）使其达到\(N=2^m\)。

3.1 基2FFT算法步骤:

首先分组，变量置换:

\[X(K)=∑_{n=0}^{N-1}x(n)W_N^{kn}, k=0,1,2,... N-1 \]

将x(n)按n地奇偶分为两组，变量置换：当n=偶数时，令n=2r,\(x(2r)=x_1(r),=0,1,2,... N/2-1\);当n=奇数时，令n=2r+1,\(x(2r+1)=x_2(r),=0,1,2,... N/2-1\);

\[X(K)=∑_{r=0}^{r=N/2-1}x(2r)W_N^{k2r}+∑_{r=0}^{r=N/2-1}x(2r+1)W_N^{k(2r+1)} \]

由旋转因子的性质\(W_N^{k2r}=e^{-j\frac{2π}N*2kr}=e^{-j\frac{2π}{N/2}*kr}=W_{N/2}^{kr}\)，所以，

\[X(K)=∑_{r=0}^{r=N/2-1}x(2r)W_{N/2}^{kr} +W_N^k ∑_{r=0}^{r=N/2-1}x(2r+1)W_{N/2}^{kr} \]

令\(X_0 (k)=∑_{r=0}^{r=N/2-1}x(2r)W_{N/2}^{kr}\)，\(X_1 (k)=∑_{r=0}^{r=N/2-1}x(2r+1)W_{N/2}^{kr}\),
\(X(K)= X_0 (k)+W_N^k X_1 (k) ， k=0,1,2,... N/2-1\)，
又根据\(W_N^{k+N/2}=W_N^{N/2}*W_N^{k}=-W_N^{k}\),将K+N/2代替K带入X(K)得到，

\[X(K+ N/2)= X_0 (k)-W_N^k X_1 (k)， k=0,1,2,... N/2-1， \]

其中X(K)和X(K+ N/2)的式子称之为蝶形运算，蝶形运算流图符号如下所示

对蝶形运算进行说明：（1）左边两路为输入；（2）右边两路为输出；（3）中间以一个小圆表示加、减运算（右上为相加输出，右下为相减输出）。一次蝶形运算需要1次复乘，2次复加。
此时一个N点的DFT可转换为2个N/2DFT，再经过第二次分解，又将N/2点DFT分解成2个N/4DFT，以此类推，经过M次分解，最后将N点FFT分解成N个1点DFT和M级蝶形运算，而1点DFT就是时域序列本身。
对于基2的时间抽取，分解到最后一级，就是两点的时域序列，两点序列求DFT，用下面的表达式：
\(X[0]=x[0]+ w_2^0x[1]\),\(X[10]=x[0]- w_2^0x[1]\),即\(X[0]=x[0]+ x[1]\),\(X[10]=x[0]- x[1]\)。

3.1 4点的基2FFT算法

\[X(K)=∑_{n=0}^{N-1}x(n)W_N^{kn}=∑_{n=0}^3x(n)W_4^{kn}=∑_{r=0}^{r=1}x(2r)W_N^{k2r}+∑_{r=0}^{r=1}x(2r+1)W_N^{k(2r+1)}=∑_{r=0}^{r=1}x(2r)W_2^{kr} +W_4^k ∑_{r=0}^{r=1}x(2r+1)W_2^{kr} \]

\[X_0 (k)=∑_{r=0}^{r=1}x(2r)W_2^{kr}， X_1 (k)=∑_{r=0}^{r=1}x(2r+1)W_2^{kr} \]

\[X(K)= X_0 (k)+W_N^k X_1 (k)， k=0,1 \]

\[X(K+ 2)= X_0 (k)-W_N^k X_1 (k)， k=0,1 \]

\[X_0 (0)= x(0)+x(2)，X_0 (1)= x(0)-x(2) \]

\[X_1 (0)= x(1)+x(3)，X_1 (1)= x(1)-x(3) \]

\[X(0)= X_0 (0)+W_4^0 X_1 (0) = X_0 (0)+X_1 (0)；X(1)= X_0 (1)+W_4^1 X_1 (1)=X_0 (0)-jX_1 (1) \]

\[X(2)= X_0 (0)-W_4^0 X_1 (0) = X_0 (0)-X_1 (0)；X(3)= X_0 (1)-W_4^1 X_1 (1)=X_0 (0)+jX_1 (1) \]

由\(X(K)= X_0 (k)+W_N^k X_1 (k)， k=0,1\)和\(X(K+ 2)= X_0 (k)-W_N^k X_1 (k)， k=0,1\),可得下图，

进一步根据，\(X_0 (0)= x(0)+x(2)，X_0 (1)= x(0)-x(2)\)和\(X_1 (0)= x(1)+x(3)，X_1 (1)= x(1)-x(3)\)，可得下图，只是本公式X_0和X_1用的下标为0和1，而图中用的1和2，对应公式的0和1。

在上图中，只标注了-1的地方，凡是1的地方都没有标注，这样使得图更简洁。

3.2 8点的基2FFT算法

按照基2时间抽取FFT算法，先将N=8点的DFT分解成2个4点DFT，
\(X(K)=∑_{n=0}^{N-1}x(n)W_N^{kn}=∑_{n=0}^7x(n)W_8^{kn}=∑_{r=0}^{r=3}x(2r)W_N^{k2r}+∑_{r=0}^{r=3}x(2r+1)W_N^{k(2r+1)}\)
\(X(K)=∑_{r=0}^{r=3}x(2r)W_2^{kr}+W_8^k ∑_{r=0}^{r=3}x(2r+1)W_2^{kr}=X_0 (k)+W_8^k X_1 (k), k=0,1,2,3,4\) 所以，

\[X(K)=X_0 (k)+W_8^k X_1 (k) \]

其中\(X_0 (k)=∑_{r=0}^{r=3}x(2r)W_2^{kr}，X_1 (k)=∑_{r=0}^{r=3}x(2r+1)W_2^{kr}\)
根据蝶形银子性质\(W_N^{k+N/2}=W_N^{N/2}*W_N^{k}=-W_N^{k}\),计算后半部分的傅里叶变换，

\[X(K+4)= X_0 (k)-W_8^k X_1 (k) \]

接着再将4点的DFT分解成2个2点的DFT，首先将\(X_0 (k)\)分解

\[X_0 (k)=∑_{r=0}^{r=3}x(2r) W_4^{kr}=∑_{l=0}^1x(4l) W_2^{kl}+W_4^k ∑_{l=0}^1x(4l+2) W_2^{kl}=X_2 (k)+W_4^k X_3 (k)， k=0,1 \]

\[X_0 (k+2)=X_2 (k)-W_4^k X_3 (k)，k=0,1 \]

\[X_2 (k)=∑_{l=0}^1x(4l) W_2^{kl}，X_2 (0)=x(0)+x(4) \]

同理计算\(X_3 (k)\),

\[X_3 (k)=∑_{l=0}^1x(4l+2) W_2^{kl},X_3 (0)= x(2)+x(6) \]

继续对\(X_1 (k)\)分解

\[X_1 (k)=∑_{r=0}^{r=3}x(2r+1)W_4^{kr}=∑_{l=0}^1x(4l+1)W_2^{kl}+W_4^k ∑_{l=0}^1x(4l+3) W_2^{kl}= X_4 (k)+W_4^k X_5 (k)， k=0,1 \]

\[X_1 (k+2)=X_4 (k)-W_4^k X_5 (k) k=0,1 \]

\[X_4 (k)=∑_{l=0}^1x(4l+1)W_2^{kl}，X_5 (k)=∑_{l=0}^1x(4l+3)W_2^{kl} \]

\[X_4 (0)= x(1)+x(5)，X_5 (0)= x(3)+x(7) \]

按照基2时间抽取FFT算法，先将N=8点的DFT分解成2个4点DFT，再将4点的DFT分解成2个2点的DFT，再继续分解。
首先分解为两个4点的DFT，有如下分解图：

再分解为2点的DFT，有如下示意图：

最后分解为1点的DFT，有如下示意图：

对于8点的序列来说，一共分为3级；

第一级是把时域的序列转化为它对应的DFT的值；第二级是两个2点序列DFT合成一个4点序列DFT；第三级是两个4点序列DFT合成一个8点序列DFT；
第一级的时候，蝶形图的间隔都是1；第二级的时候，蝶形图的间隔都是2；第二级的时候，蝶形图的间隔就变成4；
旋转因子的规律：旋转因子到最后一级一定是N0，N1，N2，… N/2 - 1。

序列的关系：

将原本的下标0,1,2,3,4,5,6,7二进制表示之后，将这个二进制序列倒一下个，就得到它对应的所要的下标。所以FFT在实现的时候，我们都用这样倒序的方式来实现。

FFT还有一个非常奇妙的地方，它可以实现原位运算。

就是后面一级的数据会自动覆盖前面一级的数据。输入序列是8点序列，在整个运算过程中始终占着8个位置，无需存储中间计算结果。
原位运算：原位运算是指当数据输入到存储器后，每级运算的结果仍然存储在原来的位置，直到最后的输出，无需存储中间的运算结果。在N=8的FFT运算中，输入序列x[0]、x[4]、x[2]、x[6]、...x[7]可以分别存入A[0]、A[4]、A[2]、A[6]、...A[7]的8个存储单元中。

4.DIT-FFT的运算规律及编程思想

4.1 运算规律

为了最终写出DIT-FFT运算程序或设计出硬件实现电路，下面介绍DIT-FFT的运算规律和编程思想。
1.原位计算
由图4.2.4可以看出，DIT-FFT的运算过程很有规律。N＝2^M点的FFT共进行M级运算，每级由N/2个叠形运算组成。同一级中，每个螺形的两个输入数据只对计算本蝶形有用，而且每个蝶形的输人、输出数据结点又同在一条水平线上，这就意味着计算完1个蝶形后，所得输出数据可立即存人原输入数据所占用的存储单元(数组元素)这样，经过M级运算后，原来存放输入序列数据的N个存储单元(数组A)中便依次存放X(k)的N个值。这种利用同一存储单元存储蝶形计算输入、输出数据的方法称为原位(址)计算，原位计算可节省大量内存，从而使设备成本降低。
2.旋转因子的变化规律
如上所述，N点 DIT-FFT运算流图中，每级都有N/2个蝶形。每个蝶形都要乘以因子W_N^{P，称其为旋转因子，p为旋转因子的指数。但各级的旋转因子和循环方式都有所不同，为了编写计算程序，应先找出旋转因子W_N}P与运算级数的关系，用L表示从左到右的运算级数(L＝1，2・…，M)。观察基2的8FFT运算图不难发现，第L级共有2^{(L-1)个不同的旋转因子，N＝2}3＝8时的各级旋转因子表示如下，

L=1时，\(W_N^P=W_{N/4}^J=W_{2^L}^J,J=0\)

L=2时，\(W_N^P=W_{N/2}^J=W_(2^L)^J,J=0,1\)

L=3时，\(W_N^P=W_N^J=W_{2^L}^J,J=0,1,2,3\)
对于\(N=2^M\)的一般情况，第L级的旋转因子为\(W_N^P=W_{2^L}^J，J=0,1,2,.... 2^{L-1}-1\).因为\(2^L=2^M×2^{L-M}=N∙2^{L-M}\),所以:

\(W_N^P=W_{N∙2^{L-M}}^J=W_N^{J∙2^{M-L} }, J=0,1,2,.... 2^{L-1}-1, (1)\)

\(P= J∙2^{M-L}, (2)\)
这样就可以按照（1）和（2）式确定第L级运算的的旋转因子（实际编程时L为最外层循环变量）。
3.蝶形运算规律
设序列x(n)经时域抽选(倒序)后，按图4.2.4所示的次序(倒序)存人数组A中，如果蝶形运算的两个输入数据相距B个点，应用原位计算，则蝶形运算可表示成如下形式:
\(A_L(J)= A_{L-1}(J)+ A_{L-1}(J+B) W_N^P\)

\(A_L(J+B)= A_{L-1}(J)- A_{L-1}(J+B) W_N^P,P= J∙2^(M-L), J=0,1,2,.... 2^(L-1)-1;L=1,2,3,…M.\)
下标L表示第L级运算、\(A _L(J)\)则表示第L级运算后的数组元素A(J)的值(即第L级蝶形的输出数据)。而\(A _{L-1}(J)\)表示第L级运算前A(J)的值(即第L级蝶形的输入数据).

4.2编程思想及程序框图

仔细观察基2的8点DITFFT运算图，还可以归纳出一些对编程有用的运算规律:第L级中，每个蝶形的两个输入数据相距B=2^{(L-1)的力个点，每级有B个不同的旋转因子;同一旋转因子对应着间隔为2}L点的2^(M-L)个蝶形。
总结上述运算规律，便可采用下述运算方法：先从输入端(第1级)开始，逐级进行，共进行M级运算。在进行第L级运算时，依次求出B个不同的旋转因子,每求出一个旋转因子，计算完它对应的所有2^(M-L)个蝶形，这样，我们可用三重循环程序实现 DIT FFT运算，程序框图如上图所示。
另列FFT算法运算流图的输出X(K)为自然顺序但为了适应原位计算、其输入序列不是按x(n)的自然顺序排列，这种经过M次偶奇抽选后的排序称为序列x(n)的倒序(倒位)，因此，在运算M级蝶形之前先对序列x(n)进行倒序，下面介绍倒序算法。
显而易见，只要将顺序数(n2n1n)的二进制位倒置，则得对应的二进制倒序值。按这一规律，用硬件电路和汇编语言程序产生倒序数很容易。但用有些高级语言程序实现时，直接倒置二进制数位是不行的，因此必须找出产生倒序数的十进制运算规律。自然顺序数增加1，是在顺序数的二进制数最低位加1，逢2向高位进位，而倒序数则是在M位二进数最高位加1，逢2向低位进位。

5.基4FFT算法

5.1基4FFT算法推导过程

\(X(K)=∑_{n=0}^{N-1}x(n)W_N^{kn}\)

\(=∑_{r=0}^{r=N/4-1}x(4r)W_N^{k4r} +∑_{r=0}^{r=N/4-1}x(4r+1)W_N^{k(4r+1)}\)

\(+∑_{r=0}^{r=N/4-1}x(4r+2)W_N^{k(4r+2)} +∑_{r=0}^{r=N/4-1}x(4r+3)W_N^{k(4r+3)}\)

\(=∑_{r=0}^{r=N/4-1}x(4r)W_N^{k4r} +W_N^k ∑_{r=0}^{r=N/4-1}x(4r+1)W_N^{k4r}\)

\(+W_N^{2k}∑_{r=0}^{r=N/4-1}x(4r+2)W_N^{k4r} +W_N^{3k} ∑_{r=0}^{r=N/4-1}x(4r+3)W_N^{k4r}\)

\(=∑_{r=0}^{r=N/4-1}x(4r)W_{N/4}^{kr} +W_N^k ∑_{r=0}^{r=N/4-1}x(4r+1)W_{N/4}^{kr}\)

\(+W_N^{2k}∑_{r=0}^{r=N/4-1}x(4r+2)W_{N/4}^{kr} +W_N^3k ∑_{r=0}^{r=N/4-1}x(4r+3)W_{N/4}^{kr}\)

令\(X_0 (k)=∑_{r=0}^{r=N/4-1}x(4r)W_{N/4}^{kr}\) , \(X_1 (k)=∑_{r=0}^{r=N/4-1}x(4r+1)W_{N/4}^{kr}\)

\(X_2 (k)=∑_{r=0}^{r=N/4-1}x(4r+2)W_{N/4}^{kr}\) , \(X_3 (k)=∑_{r=0}^{r=N/4-1}x(4r+3)W_{N/4}^{kr}\)

\(X(k)= X_0 (k)+W_N^k X_1 (k)+W_N^{2k} X_2 (k)+W_N^{3k} X_3 (k) k=0,1,2,... N/4-1\)

\(X(k+N/4)= X_0 (k)+W_N^{k+N/4} X_1 (k) +W_N^{2(k+N/4)} X_2 (k)+W_N^{3(k+N/4)} X_3 (k)\)

\(=X_0 (k)-jW_N^k X_1 (k)-W_N^{2k} X_2 (k)+jW_N^{3k} X_3 (k)\)

\(X(k+N/2)= X_0 (k)+W_N^{k+N/2} X_1 (k) +W_N^{2(k+N/2)} X_2 (k)+W_N^{3(k+N/2)} X_3 (k)\)

\(=X_0 (k)-W_N^k X_1 (k) +W_N^{2k} X_2 (k)-W_N^{3k} X_3 (k)\)

\(X(k+3N/4)= X_0 (k)+W_N^{k+3N/4} X_1 (k) +W_N^{2(k+3N/4)} X_2 (k)+W_N^{3(k+3N/4)} X_3 (k)\)

\(=X_0 (k)+jW_N^k X_1 (k)-W_N^{2k} X_2 (k)-jW_N^{3k} X_3 (k)\)

所以总结如下：

\(X(k)= X_0 (k)+W_N^k X_1 (k)+W_N^{2k} X_2 (k)+W_N^{3k} X_3 (k) k=0,1,2,... N/4-1\)

\(X(k+N/4) =X_0 (k)-jW_N^k X_1 (k)-W_N^{2k} X_2 (k)+jW_N^{3k} X_3 (k)\)

\(X(k+N/2)= X_0 (k)-W_N^k X_1 (k) +W_N^{2k} X_2 (k)-W_N^{3k} X_3 (k)\)

\(X(k+3N/4) =X_0 (k)+jW_N^k X_1 (k)-W_N^{2k} X_2 (k)-jW_N^{3k} X_3 (k)\)

\(U_0(K)= X_0 (k)+W_N^{2k} X_2 (k), U_1(K)= W_N^k X_1 (k)+W_N^{3k} X_3 (k)\)

\(U_2(K)= X_0 (k)-W_N^{2k} X_2 (k), U_3(K)= W_N^k X_1 (k)-W_N^{3k} X_3 (k)\)

\(X(k)= U_0(K)+ U_1(K), X(k+N/4)= U_2(K)-jU_3(K)\)，

\(X(k+N/2)= U_0(K)- U_1(K)，X(k+3N/4) =U_2(K)+jU_3(K)\)

此时一个N点的DFT可转换为N/4DFT。对于一个4点FFT的基4蝶形运算单元为：

5.1 16点FFT的基4运算

对于一个16点FFT的基4运算：

\(X_0 (k)=∑_{r=0}^{r=3}x(4r)W_{N/4}^{kr}, X_1 (k)=∑_{r=0}^{r=3}x(4r+1)W_{N/4}^{kr}\)

\(X_2 (k)=∑_{r=0}^{r=3}x(4r+2)W_{N/4}^{kr}, X_3 (k)=∑_{r=0}^{r=3}x(4r+3)W_{N/4}^{kr}\)

\(X_0 (0)=x(0)+x(4)+x(8)+x(12); X_0 (1)=x(0)-jx(4)-x(8)+jx(12)\)

\(X_0 (2)=x(0)-x(4)+x(8)-x(12);X_0 (3)=x(0)+jx(4)-x(8)-jx(12)\)

\(X(k)= X_0 (k)+W_16^k X_1 (k)+W_16^{2k} X_2 (k)+W_16^{3k} X_3 (k) k=0,1,2,3\)

\(X(k+4)= X_0 (k)-jW_16^k X_1 (k)-W_{16}^{2k} X_2 (k)+jW_{16}^{3k} X_3 (k)\)

\(X(k+8)=X_0 (k)-W_16^k X_1 (k) +W_{16}^{2k} X_2 (k)-W_{16}^{3k} X_3 (k)\)

\(X(k+12)=X_0 (k)+jW_{16}^k X_1 (k)-W_{16}^{2k} X_2 (k)-jW_{16}^{3k} X_3 (k)\)

5.2 N为16的基4FFT分解过程

一个4点序列的DFT运算流图：

对于基4运算，每个蝶形运算有4点输入和4点输出，因此总共有N/4个蝶形运算单元。由于每一个蝶形运算仅需3次复数乘法和8次复数加法，因此，对每一次分解，蝶式运算所需的复数乘法次数和复数加法次数分别为3N/4次和2N次。
在\(N＝4^M\)(M为正整数)的一般情形下，N点序列分解成一点序列需要作\(M=log_4⁡N\)次分解，因此，基4时分FFT算法所需要的复数乘法次数为\(\frac{3N}4log_4⁡N\)次,复数加法次数为\(2Nlog_4⁡N\).

5.3 DIT-FFT的运算规律及编程思想

为了最终写出DIT-FFT运算程序或设计出硬件实现电路，下面介绍DIT-FFT的运算规律和编程思想。

5.3.1 原位计算

由基4的16点FFT蝶形图可知，DIT-FFT的运算过程很有规律。\(N＝4^M\)点的FFT共进行M级运算，每级由N/4个叠形运算组成。同一级中，每个螺形的四个输入数据只对计算本蝶形有用，而且每个蝶形的输人、输出数据结点又同在一条水平线上，这就意味着计算完1个蝶形后，所得输出数据可立即存人原输入数据所占用的存储单元(数组元素)这样，经过M级运算后，原来存放输入序列数据的N个存储单元(数组A)中便依次存放X(k)的N个值。这种利用同一存储单元存储蝶形计算输入、输出数据的方法称为原位(址)计算，原位计算可节省大量内存，从而使设备成本降低。

5.3.2 旋转因子的变化规律

如上所述，N点 DIT-FFT运算流图中，每级都有N/4个蝶形。每个蝶形都要乘以因子\(W_N^P\)，称其为旋转因子，p为旋转因子的指数。但各级的旋转因子和循环方式都有所不同，为了编写计算程序，应先找出旋转因子W_N^{P与运算级数的关系，用L表示从左到右的运算级数(L＝1，2・…，M)。观察蝶形运算图不难发现，第L级共有$4^{L-1}$个不同的旋转因子，N＝4}3＝64时的各级旋转因子表示如下，

L=1时，\(W_N^P=W_{N/16}^J=W_{4^L}^J\), \(W_N^{2P}=W_{N/16}^{2J}=W_{4^L}^{2J}\), \(W_N^{3P}=W_{N/16}^3J=W_{4^L}^{3J}，J=0\)

L=2时，\(W_N^P=W_{N/2}^J=W_{4^L}^J\), \(W_N^{2P}=W_{N/2}^{2J}=W_{4^L}^{2J}\)，\(W_N^{3P}=W_{N/2}^3J=W_{4^L}^{3J}，J=0,1,2,3\)

L=3时，\(W_N^P=W_N^J=W_{4^L}^J\), \(W_N^{2P}=W_N^{2J}=W_{4^L}^{2J}\)，\(W_N^3P=W_N^{3J}=W_{4^L}^{3J}，J=0,1,2,...15\)

对于\(N=4^M\)的一般情况，第L级的旋转因子为\(W_N^P=W_{4^L}^J\)，\(W_N^{2P}=W_{4^L}^{2J}\)，\(W_N^{3P}=W_{4^L}^{3J}\)，\(J=0,1,2,.... 4^{L-1}-1\).因为\(4^L=4^M×4^{L-M}=N∙4^{L-M}\),所以:

\(W_N^P=W_{N∙4^{L-M}}^J=W_N^{J∙4^{M-L} }, J=0,1,2,.... 4^{L-1}-1, (1)\)

\(P= J∙4^{M-L}, (2)\)
这样就可以按照（1）和（2）式确定第L级运算的的旋转因子（实际编程时L为最外层循环变量）。

5.3.3 蝶形运算规律

设序列x(n)经时域抽选(倒序)后，按图4.2.4所示的次序(倒序)存人数组A中，如果蝶形运算的4个输入数据分别相距B个点，应用原位计算，则蝶形运算可表示成如下形式:
\(A_L(J)= A_{L-1}(J)+ A_{L-1}(J+B) W_N^P\)

\(A_L(J+B)= A _{L-1}(J)- A _{L-1}(J+B) W_N^P,P= J∙2^{M-L}, J=0,1,2,.... 4^(L-1)-1;L=1,2,3,…M.\)

\(X(k)= X_0 (k)+W_{16}^k X_1 (k)+W_{16}^{2k} X_2 (k)+W_{16}^{3k} X_3 (k) k=0,1,2,3\)

\(X(k+4)= X_0 (k)-jW_{16}^k X_1 (k)-W_{16}^{2k} X_2 (k)+jW_{16}^{3k} X_3 (k)\)

\(X(k+8)=X_0 (k)-W_{16}^k X_1 (k) +W_{16}^{2k} X_2 (k)-W_{16}^{3k} X_3 (k)\)

\(X(k+12)=X_0 (k)+jW_{16}^k X_1 (k)-W_{16}^{2k} X_2 (k)-jW_{16}^{3k} X_3 (k)\)

下标L表示第L级运算、\(A_L(J)\)则表示第L级运算后的数组元素A(J)的值(即第L级蝶形的输出数据)。而\(A _{L-1}(J)\)表示第L级运算前A(J)的值(即第L级蝶形的输入数据).

5.3.4 基4算法的模块实现

FFT模块的构成：
1.旋转因子乘法器；2.基4蝶形运算器；3.地址产生器；4.控制信号产生器；5中间变量存放组。

1．旋转因子乘法器的实现
FFT设计一个旋转因子存储单元可以提高FFT的计算速度，因此需要合理的设计旋转因子的存储单元，在FFT便件设计中占有重要地位。可以利用高效的算法来实现存储旋转因子。
3.地址产生单元生成 RAM读写地址, 写使能信号以及相关模块的启动、 控制信号, 是系统的控制核心;4 点蝶形运算单元的最后一级输出不是顺序的 ;旋转因子产生单元生成复乘运算中的旋转因子的角度数据;旋转因子ROM中预置了每一级运算中所需的旋转因子
1.运算模块为基4运算模块，控制模块产生所有的控制信号，存储器1和2的读写地址，写使能，运算模块的启动信号因子表的读写地址等信号。存储器1作为当前输入标志对应输入N点数据缓冲器，存储器2作为中间变量结果存储器，用于存储运算模块计算出各个通道的结果。
举例：外部输入为N点数据段流和启动信号，一方面，外部数据存入储器1中，同时通过控制模块的控制，读出存储器1中的前段N点数据和ROM表中的因子及相关控制信号送入运算核心模块进行各个Pass运算，每个Pass运算的输出都存入存储器2中，最后一个Pass的计算结果存入存储器2中，启动并在下一个启动头到来后，输出计算结果。
4.旋转因子与数据相乘实现
//通过旋转因子的地址查找表查找相位因子对应的旋转因子与输入
//数据相乘。

`timescale 1ns/1ps
//name: phase_factor
//function: multiply with input data by phase factor

module phase_factor(line1,line2,line3,line4,
    in_pathway1,in_pathway2,in_pathway3,in_pathway4,
                   lut_phase,CLK20M);               
 input CLK20M;
 input [5:0] lut_phase; // address for lut of phase factor,
//最多2^6=64 phase factor
 input [31:0] line1,line2,line3,line4; // input data 
 output [31:0] in_pathway1,in_pathway2,in_pathway3,in_pathway4;
// result/output data
              
 reg [31:0] in_pathway1,in_pathway2,in_pathway3,in_pathway4;
 reg [31:0] delay_line1;
 
 wire [31:0] dout1,dout2,dout3,dout4,
             dout5,dout6,dout7,dout8,
             dout9,dout10,dout11,dout12;
             
always @(posedge CLK20M) // delay for wait other 3 line delay multiply
  delay_line1 <= line1;
//first16 real,last16image
wire [15:0] temp1_r1 = line1 [31:16];
wire [15:0] temp1_i1 = line1 [15: 0];

wire [15:0] temp1_r2 = line2 [31:16];
wire [15:0] temp1_i2 = line2 [15: 0];

wire [15:0] temp1_r3 = line3 [31:16];
wire [15:0] temp1_i3 = line3 [15: 0];

wire [15:0] temp1_r4 = line4 [31:16];
wire [15:0] temp1_i4 = line4 [15: 0];

wire [15:0] temp2_r2; //= fac2 [31:16];
wire [15:0] temp2_i2; //= fac2 [15: 0];

wire [15:0] temp2_r3; //= fac3 [31:16];
wire [15:0] temp2_i3; //= fac3 [15: 0];

wire [15:0] temp2_r4; //= fac4 [31:16];
wire [15:0] temp2_i4; //= fac4 [15: 0];

//d_rom_fft_phase2 m2(.A(lut_phase),.SPO(fac2));
//d_rom_fft_phase3 m3(.A(lut_phase),.SPO(fac3));
//d_rom_fft_phase4 m4(.A(lut_phase),.SPO(fac4));
d_rom_fft_phase2r m2r(.A(lut_phase),.SPO(temp2_r2));
d_rom_fft_phase2i m2i(.A(lut_phase),.SPO(temp2_i2));

d_rom_fft_phase3r m3r(.A(lut_phase),.SPO(temp2_r3));
d_rom_fft_phase3i m3i(.A(lut_phase),.SPO(temp2_i3));

d_rom_fft_phase4r m4r(.A(lut_phase),.SPO(temp2_r4));
d_rom_fft_phase4i m4i(.A(lut_phase),.SPO(temp2_i4));

 
// m16 is 16 multiply 16 ,has one delay for pipline
// for multiply (a+bj) * (c+dj) = ac-bd+(ad+bc)j
m16 u1(.clk(CLK20M),.a(temp1_r2),.b(temp2_r2),.o(dout1)); //ac
m16 u2(.clk(CLK20M),.a(temp1_i2),.b(temp2_i2),.o(dout2)); //bd
m16 u3(.clk(CLK20M),.a(temp1_r2),.b(temp2_i2),.o(dout3)); //ad
m16 u4(.clk(CLK20M),.a(temp1_i2),.b(temp2_r2),.o(dout4)); //bc
wire [15:0] dot1 =dout1[29:14]; //= dout1[30] ? dout1[30:15] : dout1[30:15] +1'b1;
wire [15:0] dot2 =dout2[29:14];//= dout2[30] ? dout2[30:15] : dout2[30:15] +1'b1;
wire [15:0] dot3 =dout3[29:14];//= dout3[30] ? dout3[30:15] : dout3[30:15] +1'b1;
wire [15:0] dot4 =dout4[29:14];//= dout4[30] ? dout4[30:15] : dout4[30:15] +1'b1;


wire [15:0] temp3_r2; //= dout1 [30:15] + (~dout2[30:15] +1'b1);
wire [15:0] temp3_i2; //= dout3 [30:15] +  dout4 [30:15];

subtract16 sub1(.A(dot1),.B(dot2),.S(temp3_r2));
adder16    add1(.A(dot3),.B(dot4),.S(temp3_i2));

m16 u5(.clk(CLK20M),.a(temp1_r3),.b(temp2_r3),.o(dout5));
m16 u6(.clk(CLK20M),.a(temp1_i3),.b(temp2_i3),.o(dout6));
m16 u7(.clk(CLK20M),.a(temp1_r3),.b(temp2_i3),.o(dout7));
m16 u8(.clk(CLK20M),.a(temp1_i3),.b(temp2_r3),.o(dout8));

wire [15:0] dot5 =dout5[29:14];//= dout5[30] ? dout5[30:15] : dout5[30:15] +1'b1;
wire [15:0] dot6 =dout6[29:14];//= dout6[30] ? dout6[30:15] : dout6[30:15] +1'b1;
wire [15:0] dot7 =dout7[29:14];//= dout7[30] ? dout7[30:15] : dout7[30:15] +1'b1;
wire [15:0] dot8 =dout8[29:14];//= dout8[30] ? dout8[30:15] : dout8[30:15] +1'b1;

wire [15:0] temp3_r3; //= dout5 [30:15] + (~dout6[30:15] +1'b1);
wire [15:0] temp3_i3; //= dout7 [30:15] +  dout8 [30:15];

subtract16 sub2(.A(dot5),.B(dot6),.S(temp3_r3));
adder16    add2(.A(dot7),.B(dot8),.S(temp3_i3));

m16 u9(.clk(CLK20M),.a(temp1_r4),.b(temp2_r4),.o(dout9));
m16 u10(.clk(CLK20M),.a(temp1_i4),.b(temp2_i4),.o(dout10));
m16 u11(.clk(CLK20M),.a(temp1_r4),.b(temp2_i4),.o(dout11));
m16 u12(.clk(CLK20M),.a(temp1_i4),.b(temp2_r4),.o(dout12));
wire [15:0] dot9  = dout9[29:14];//= dout9[30] ? dout9[30:15] : dout9[30:15] +1'b1;
wire [15:0] dot10 =dout10[29:14];//= dout10[30] ? dout10[30:15] : dout10[30:15] +1'b1;
wire [15:0] dot11 =dout11[29:14];//= dout11[30] ? dout11[30:15] : dout11[30:15] +1'b1;
wire [15:0] dot12 =dout12[29:14];//= dout12[30] ? dout12[30:15] : dout12[30:15] +1'b1;

wire [15:0] temp3_r4; //= dout9 [30:15] + (~dout10[30:15] +1'b1);
wire [15:0] temp3_i4; //= dout11 [30:15] + dout12 [30:15];

subtract16 sub3(.A(dot9),.B(dot10),.S(temp3_r4));
adder16    add3(.A(dot11),.B(dot12),.S(temp3_i4));

always @(posedge CLK20M) begin
   in_pathway1 <= delay_line1;
   in_pathway2 <= {temp3_r2,temp3_i2};
   in_pathway3 <= {temp3_r3,temp3_i3};
   in_pathway4 <= {temp3_r4,temp3_i4};
end

endmodule
//function finish

5.蝶形运算单元
`timescale 1ns/1ps
//name: butterfly
//function: the butterfly unit for radix 4 structure,ever step is same
//add reg_rx reg_ix for high speed pipline; 添加了rx和ix寄存器用于高速流水线实现。

module ifft_butterfly (in_pathway1,in_pathway2,in_pathway3,in_pathway4,
                  out_pathway1,out_pathway2,out_pathway3,out_pathway4,
                  CLK20M);

input CLK20M;   
input  [31:0] in_pathway1,in_pathway2,in_pathway3,in_pathway4;//4 pathway input butterfly
output[31:0]out_pathway1,out_pathway2,out_pathway3,out_pathway4;//4pathwayout butterfly

wire [15:0] in_real1,in_real2,in_real3,in_real4,
         in_imag1,in_imag2,in_imag3,in_imag4;
            
reg [15:0] reg_r1,reg_r2,reg_r3,reg_r4;    //中间实部寄存器
reg	[15:0] reg_i1,reg_i2,reg_i3,reg_i4;   //中间虚部寄存器

//高16位存储的是实部，低16位存储的是虚部
assign in_real1 = in_pathway1 [31:16];// first 16 for real and last for image
assign in_imag1 = in_pathway1 [15: 0];

assign in_real2 = in_pathway2 [31:16];
assign in_imag2 = in_pathway2 [15: 0];

assign in_real3 = in_pathway3 [31:16];
assign in_imag3 = in_pathway3 [15: 0];

assign in_real4 = in_pathway4 [31:16];
assign in_imag4 = in_pathway4 [15: 0];

                  
// begin caculate
// first step for butterfly
wire [15:0] temp1_r1; //= in_real1 +  in_real3;
wire [15:0] temp1_i1; //= in_imag1 +  in_imag3;

adder16 add4(.A(in_real1),.B(in_real3),.S(temp1_r1));
adder16 add5(.A(in_imag1),.B(in_imag3),.S(temp1_i1));

wire [15:0] temp1_r2; //= in_real2 +  in_real4;
wire [15:0] temp1_i2; //= in_imag2 +  in_imag4;

adder16 add6(.A(in_real2),.B(in_real4),.S(temp1_r2));
adder16 add7(.A(in_imag2),.B(in_imag4),.S(temp1_i2));

wire [15:0] temp1_r3; //= in_real1 + (~in_real3 +1'b1);// 补码计算，取反加1，in_real1-in_real3
wire [15:0] temp1_i3; //= in_imag1 + (~in_imag3 +1'b1);// 

subtract16 sub4(.A(in_real1),.B(in_real3),.S(temp1_r3));
subtract16 sub5(.A(in_imag1),.B(in_imag3),.S(temp1_i3));

wire [15:0] temp1_r4; //= in_real2 + (~in_real4 +1'b1);// in_real2-in_real4 
wire [15:0] temp1_i4; //= in_imag2 + (~in_imag4 +1'b1);// 

subtract16 sub6(.A(in_real2),.B(in_real4),.S(temp1_r4));
subtract16 sub7(.A(in_imag2),.B(in_imag4),.S(temp1_i4));

always @(posedge CLK20M) begin // reg for pipline
   reg_r1 <= temp1_r1;
   reg_r2 <= temp1_r2;
   reg_r3 <= temp1_r3;
   reg_r4 <= temp1_r4;

   reg_i1 <= temp1_i1;
   reg_i2 <= temp1_i2;
   reg_i3 <= temp1_i3;
   reg_i4 <= temp1_i4;
end
// second step for butterfly
// (a+bj)*j = -b+aj 
wire [15:0] temp2_r1; //= temp1_r1 +  temp1_r2; // wire1 = 1 +2;
wire [15:0] temp2_i1; //= temp1_i1 +  temp1_i2;

adder16 add8(.A(reg_r1),.B(reg_r2),.S(temp2_r1));
adder16 add9(.A(reg_i1),.B(reg_i2),.S(temp2_i1));

wire [15:0] temp2_r2; //= temp1_r3 +  temp1_i4; // wire2 = 3 -4*j;
wire [15:0] temp2_i2; //= temp1_i3 + (~temp1_r4 +1'b1); // 

adder16    add10(.A(reg_r3),.B(reg_i4),.S(temp2_r2));
subtract16 sub8 (.A(reg_i3),.B(reg_r4),.S(temp2_i2));

wire [15:0] temp2_r3; //= temp1_r1 + (~temp1_r2 +1'b1); // wire3 = 1 -2;
wire [15:0] temp2_i3; //= temp1_i1 + (~temp1_i2 +1'b1); // 

subtract16 sub9 (.A(reg_r1),.B(reg_r2),.S(temp2_r3));
subtract16 sub10(.A(reg_i1),.B(reg_i2),.S(temp2_i3));

wire [15:0] temp2_r4; //= temp1_r3 + (~temp1_i4 +1'b1); // wire4 = 3 +4*j;
wire [15:0] temp2_i4; //= temp1_i3 +  temp1_r4;

subtract16 sub11(.A(reg_r3),.B(reg_i4),.S(temp2_r4));
adder16    add11(.A(reg_i3),.B(reg_r4),.S(temp2_i4));
// end caculate

assign out_pathway1 = {temp2_r1,temp2_i1};
assign out_pathway2 = {temp2_r2,temp2_i2};
assign out_pathway3 = {temp2_r3,temp2_i3};
assign out_pathway4 = {temp2_r4,temp2_i4};
// function finish
endmodule

//analysis: first 16 for real and last for image
in_real1 = in_pathway1 [31:16]; in_imag1 = in_pathway1 [15: 0];
in_real2 = in_pathway2 [31:16]; in_imag2 = in_pathway2 [15: 0];
in_real3 = in_pathway3 [31:16]; in_imag3 = in_pathway3 [15: 0];
in_real4 = in_pathway4 [31:16]; in_imag4 = in_pathway4 [15: 0];
// first step for butterfly
temp1_r1= in_real1 + in_real3;      temp1_i1= in_imag1 + in_imag3;
temp1_r2= in_real2 + in_real4;      temp1_i2= in_imag2 + in_imag4;
temp1_r3= in_real1 - in_real3 ;      temp1_i3= in_imag1 - in_imag3;
temp1_r4= in_real2 - in_real4;      temp1_i4=in_imag2 - in_imag4;
// second step for butterfly, // (a+bj)*j = -b+aj 
temp2_r1= temp1_r1+temp1_r2= in_real1 + in_real3+ in_real2 + in_real4;
temp2_i1=temp1_i1 + temp1_i2= in_imag1 + in_imag3+in_imag2 + in_imag4;

temp2_r2= temp1_r3 + temp1_i4= in_real1 - in_real3+ in_imag2 - in_imag4;
temp2_i2= temp1_i3-temp1_r4= in_imag1 - in_imag3-(in_real2 - in_real4);
 
temp2_r3= temp1_r1 -temp1_r2= in_real1 + in_real3-( in_real2 + in_real4);
temp2_i3= temp1_i1 -temp1_i2= in_imag1 + in_imag3-( in_imag2 + in_imag4);
 
temp2_r4= temp1_r3-temp1_i4= in_real1 - in_real3-( in_imag2 - in_imag4); 
temp2_i4= temp1_i3 + temp1_r4= in_imag1 - in_imag3+ in_real2 - in_real4;