FFT算法介紹及OAI中FFT的實作

• OAI中FFT的實作
• • FFT介紹
  • DIT radix-2 FFT（時域抽取基-2 FFT）
  • DIF radix-2 FFT (頻域抽取基-2 FFT)
  • DIF radix-4 FFT (頻域抽取基-4 FFT)
  • OAI中FFT的實作

OAI中FFT的實作

FFT介紹

FFT是DFT（離散傅立葉變換）的一種快速算法。

DFT是一種常用的信号處理方式，可以将時域信号轉換到頻域。它的原理在《信号與系統》課程中有詳細闡述。

但是因為DFT的運算複雜度太高，在實際應用中都是采用FFT的快速算法。

快速傅立葉變換算法(FFT)最開始是由庫利（T.W.Cooley）和圖基（J.W.Tukey）于1965年提出。其基本思想就是分治。從時域或者頻域序列中抽取出不同的子序列進而降低計算複雜度。

根據抽取對象的不同又細分為不同的算法：基于時域的資料做抽取的稱為DIT(Decimation in Time), 基于頻域的抽取就稱為DIF(Decimation in Frequency)。

另外一個區分不同FFT算法的标準稱為基(radix)。嗯，基無處不在。

基決定FFT算法的最小組成機關，這個最小組成機關一般是用被稱為蝶形圖的方式來表示的。

圖1

圖中，(a)表示一個基-2的基本蝶形圖，(b)的左半部分表示一個基-4的基本蝶形圖。

可以看到，基數決定來基本蝶形圖的輸入資料數。除了2和4之外，當然也有可能是其他的基數，或者混合使用幾種不同的基數。

庫利和圖基當時提出的就是時域的基-2算法，是以這個算法又被稱為庫利-圖基算法。

下面我們舉幾個FFT的例子，來形象的了解一下這個偉大的算法。

DIT radix-2 FFT（時域抽取基-2 FFT）

先給出離散傅立葉變換的一般形式：

(1) X [ k ] = ∑ n = 0 N − 1 x [ n ] W N k n , 0 ≤ k ≤ N − 1 X[k]=\sum_{n=0}^{N-1} x[n]W_N^{kn},\qquad \scriptsize 0 \leq k \leq N-1 \normalsize \tag{1} X[k]=n=0∑N−1​x[n]WNkn​,0≤k≤N−1(1)

其中， N N N是時域離散序列 x [ n ] x[n] x[n]的長度，并且是2的指數倍。 W N = e − j 2 π / N W_N=e^{-j2\pi /N} WN​=e−j2π/N稱為旋轉因子，有如下性質：

1)周期性

W N ( k + N ) n = W N k n = W N k ( n + N ) W_N^{(k+N)n}=W_N^{kn}=W_N^{k(n+N)} WN(k+N)n​=WNkn​=WNk(n+N)​

2)對稱性

W N n k + N / 2 = − W N n k W_N^{nk+N/2}=-W_N^{nk} WNnk+N/2​=−WNnk​

3)可約性

W m N n m k = W N n k W_{mN}^{nmk}=W_N^{nk} WmNnmk​=WNnk​

所謂radix-2 DIT就是把時域的序列按照奇偶抽取為2個子序列，這2個子序列的長度都是原序列的一半長度，使得計算複雜度減小到原序列的一半。把這個分治過程遞歸應用到這兩個子序列，直到最後剩下2點的DFT。（其實2點DFT還可以繼續遞歸這個過程到1點）

用公式表示，就是把時域序列按照奇偶劃分為兩個子序列，如下：

(2) X [ k ] = ∑ n = 0 N − 1 x [ n ] W N k n , k = 0 , 1 , . . . , N − 1 = ∑ n   e v e n x [ n ] W N k n + ∑ n   o d d x [ n ] W N k n = ∑ m = 0 N / 2 − 1 x [ 2 m ] W N 2 m k + ∑ m = 0 N / 2 − 1 x [ 2 m + 1 ] W N k ( 2 m + 1 ) X[k] =\sum_{n=0}^{N-1}x[n]W_N^{kn},\qquad \scriptsize k=0,1,...,N-1 \normalsize \\ =\sum_{n \, even}x[n]W_N^{kn} + \sum_{n \, odd}x[n]W_N^{kn} \\ =\sum_{m=0}^{N/2-1}x[2m]W_N^{2mk}+\sum_{m=0}^{N/2-1}x[2m+1]W_N^{k(2m+1)} \tag{2} X[k]=n=0∑N−1​x[n]WNkn​,k=0,1,...,N−1=neven∑​x[n]WNkn​+nodd∑​x[n]WNkn​=m=0∑N/2−1​x[2m]WN2mk​+m=0∑N/2−1​x[2m+1]WNk(2m+1)​(2)

(公式沒有對齊，LaTex的對齊文法老是提示錯誤。)

上式中等号右邊第一部分是偶數部分，第二部分是奇數部分，兩者合起來正好是原時域離散數列。

如果我們令：

f 1 [ m ] = x [ 2 m ] f 2 [ m ] = x [ 2 m + 1 ] f_1[m]=x[2m] \\ f_2[m]=x[2m+1] f1​[m]=x[2m]f2​[m]=x[2m+1]

再令：

F 1 [ k ] = D F T ( f 1 [ m ] ) F 2 [ k ] = D F T ( f 2 [ m ] ) F_1[k]=DFT(f_1[m]) \\ F_2[k]=DFT(f_2[m]) F1​[k]=DFT(f1​[m])F2​[k]=DFT(f2​[m])

那麼利用前面提到的旋轉因子的性質， ( 2 ) (2) (2)式就可以表示為：

X [ k ] = ∑ m = 0 N / 2 − 1 f 1 [ m ] W N / 2 k m + W N k ∑ m = 0 N / 2 − 1 f 2 [ m ] W N / 2 k m = F 1 [ k ] + W N k F 2 [ k ] , k = 0 , 1 , . . . , N − 1 X[k]=\sum_{m=0}^{N/2-1}f_1[m]W_{N/2}^{km}+W_N^k\sum_{m=0}^{N/2-1}f_2[m]W_{N/2}^{km} \\ =F_1[k]+W_N^kF_2[k], \qquad \scriptsize k=0,1,...,N-1 X[k]=m=0∑N/2−1​f1​[m]WN/2km​+WNk​m=0∑N/2−1​f2​[m]WN/2km​=F1​[k]+WNk​F2​[k],k=0,1,...,N−1

注意 F 1 [ k ] 和 F 2 [ k ] F_1[k]和F_2[k] F1​[k]和F2​[k]都是 N / 2 N/2 N/2為周期的周期函數，并且 W N k + N / 2 = − W N k W_N^{k+N/2}=-W_N^k WNk+N/2​=−WNk​, 于是可以得到：

X [ k ] = F 1 [ k ] + W N k F 2 [ k ] , k = 0 , 1 , . . . , N / 2 − 1 X [ k + N / 2 ] = F 1 [ k ] − W N k F 2 [ k ] , k = 0 , 1 , . . . , N / 2 − 1 X[k]=F_1[k]+W_N^kF_2[k], \qquad \scriptsize k=0,1,...,N/2-1 \normalsize \\ X[k+N/2]=F_1[k]-W_N^kF_2[k], \qquad \scriptsize k=0,1,...,N/2-1 X[k]=F1​[k]+WNk​F2​[k],k=0,1,...,N/2−1X[k+N/2]=F1​[k]−WNk​F2​[k],k=0,1,...,N/2−1

如果把這個結論應用到8點DFT，我們可以得到下面這張圖：

圖2

圖2中，第一個4點DFT得出的是 F 1 [ k ] F_1[k] F1​[k]，第二個4點DFT得出的是 F 2 [ k ] F_2[k] F2​[k]。注意，旋轉因子隻在第二部分有，并且與 k k k值相關。

可以繼續對這兩個4點DFT繼續遞歸使用上述分治方法，進而得到下面的結果：

F 1 [ k ] = F { f 1 [ 2 n ] } + W N / 2 k F { f 1 [ 2 n + 1 ] } , k = 0 , 1 , . . . , N / 4 − 1 n = 0 , 1 , . . . , N / 4 − 1 F 1 [ k + N / 4 ] = F { f 1 [ 2 n ] } − W N / 2 k F { f 1 [ 2 n + 1 ] } , k = 0 , 1 , . . . , N / 4 − 1 n = 0 , 1 , . . . , N / 4 − 1 F 2 [ k ] = F { f 2 [ 2 n ] } + W N / 2 k F { f 2 [ 2 n + 1 ] } , k = 0 , 1 , . . . , N / 4 − 1 n = 0 , 1 , . . . , N / 4 − 1 F 2 [ k + N / 4 ] = F { f 2 [ 2 n ] } − W N / 2 k F { f 2 [ 2 n + 1 ] } , k = 0 , 1 , . . . , N / 4 − 1 n = 0 , 1 , . . . , N / 4 − 1 F_1[k]=F\{f_1[2n]\}+W_{N/2}^kF\{f_1[2n+1]\}, \quad \scriptsize k=0,1,...,N/4-1 \quad n=0,1,...,N/4-1 \normalsize \\ F_1[k+N/4]=F\{f_1[2n]\}-W_{N/2}^kF\{f_1[2n+1]\}, \quad \scriptsize k=0,1,...,N/4-1 \quad n=0,1,...,N/4-1 \normalsize \\ F_2[k]=F\{f_2[2n]\}+W_{N/2}^kF\{f_2[2n+1]\}, \quad \scriptsize k=0,1,...,N/4-1 \quad n=0,1,...,N/4-1 \normalsize \\ F_2[k+N/4]=F\{f_2[2n]\}-W_{N/2}^kF\{f_2[2n+1]\}, \quad \scriptsize k=0,1,...,N/4-1 \quad n=0,1,...,N/4-1 \normalsize \\ F1​[k]=F{f1​[2n]}+WN/2k​F{f1​[2n+1]},k=0,1,...,N/4−1n=0,1,...,N/4−1F1​[k+N/4]=F{f1​[2n]}−WN/2k​F{f1​[2n+1]},k=0,1,...,N/4−1n=0,1,...,N/4−1F2​[k]=F{f2​[2n]}+WN/2k​F{f2​[2n+1]},k=0,1,...,N/4−1n=0,1,...,N/4−1F2​[k+N/4]=F{f2​[2n]}−WN/2k​F{f2​[2n+1]},k=0,1,...,N/4−1n=0,1,...,N/4−1

其中 F { ∗ } F\{*\} F{∗}表示傅立葉變換，我們把這個結論再應用到上述8點DFT的流程圖得到下面的結果：

圖3

基-2的FFT算法，最終目的就是要把輸入資料劃分成2點的蝶形圖。2點蝶形圖替換上圖中的2點DFT得到最終的8點DIT radix-2 FFT蝶形圖。

圖4

圖4表示了完整的8點FFT流程，但是這其中有一個問題：由于對時域資料不斷應用分治政策，導緻最後要使用的時域資料的順序是被打亂的。

實際上，被打亂後的序号與最開始的序号是bit-reverse的關系。如下圖：

圖5

是以，在實際編碼時，需要預先把輸入資料按照bit-reverse後的順序重新排列再輸入計算。

DIF radix-2 FFT (頻域抽取基-2 FFT)

另一種重要的FFT算法是基于對頻域序列的抽取的，稱為DIF FFT。

DIF radix-2 FFT的算法與DIT radix-2FFT的差別是，它把頻域資料分成兩個子序列。

先将式 ( 1 ) (1) (1)中的 X ( k ) X(k) X(k)劃分成如下兩部分：

X [ k ] = ∑ n = 0 N / 2 − 1 x [ n ] W N k n + W N k N / 2 ∑ n = 0 N / 2 − 1 x [ n + N / 2 ] W N k n X[k]=\sum_{n=0}^{N/2-1}x[n]W_N^{kn}+W_N^{kN/2}\sum_{n=0}^{N/2-1}x[n+N/2]W_N^{kn} X[k]=n=0∑N/2−1​x[n]WNkn​+WNkN/2​n=0∑N/2−1​x[n+N/2]WNkn​

因為 W N k N / 2 = ( − 1 ) k W_N^{kN/2}=(-1)^k WNkN/2​=(−1)k, 是以上式可以寫成：

X [ k ] = ∑ n = 0 N / 2 − 1 [ x [ n ] + ( − 1 ) k x [ n + N / 2 ] ] W N k n X[k]=\sum_{n=0}^{N/2-1}\left[x[n]+(-1)^kx[n+N/2]\right]W_N^{kn} X[k]=n=0∑N/2−1​[x[n]+(−1)kx[n+N/2]]WNkn​

然後，把 X ( k ) X(k) X(k)分成偶數和奇數兩個子序列，因為 W N 2 = W N / 2 W_N^2=W_{N/2} WN2​=WN/2​，是以我們得到：

X [ 2 k ] = ∑ n = 0 N / 2 − 1 [ x [ n ] + x [ n + N / 2 ] ] W N / 2 k n , k = 0 , 1 , . . . , N / 2 X [ 2 k + 1 ] = ∑ n = 0 N / 2 − 1 [ x [ n ] − x [ n + N / 2 ] ] W N / 2 k n , k = 0 , 1 , . . . , N / 2 X[2k]=\sum_{n=0}^{N/2-1}\left[x[n]+x[n+N/2]\right]W_{N/2}^{kn}, \quad \scriptsize k=0,1,...,N/2 \normalsize \\ X[2k+1]=\sum_{n=0}^{N/2-1}\left[x[n]-x[n+N/2]\right]W_{N/2}^{kn}, \quad \scriptsize k=0,1,...,N/2 X[2k]=n=0∑N/2−1​[x[n]+x[n+N/2]]WN/2kn​,k=0,1,...,N/2X[2k+1]=n=0∑N/2−1​[x[n]−x[n+N/2]]WN/2kn​,k=0,1,...,N/2

令 f 1 [ n ] = x [ n ] + x [ n + N / 2 ] f_1[n]=x[n]+x[n+N/2] f1​[n]=x[n]+x[n+N/2], f 2 [ n ] = x [ n ] − x [ n + N / 2 ] f_2[n]=x[n]-x[n+N/2] f2​[n]=x[n]−x[n+N/2], 那麼 X [ 2 k ] X[2k] X[2k]和 X [ 2 k + 1 ] X[2k+1] X[2k+1]可以看成是 f 1 [ n ] f_1[n] f1​[n]和 f 2 [ n ] f_2[n] f2​[n]在 0 ～ N / 2 0～N/2 0～N/2上的傅立葉變換。

将上式應用于8點DFT，得到如下流程圖：

圖6

注意旋轉因子在每組蝶形圖中隻作用于下半部分的蝶形圖，即 n + N / 2 n+N/2 n+N/2部分的 x [ n ] x[n] x[n]需要多乘一個旋轉因子。

将上述分治過程遞歸應用于 f 1 [ n ] f_1[n] f1​[n]和 f 2 [ n ] f_2[n] f2​[n]，即上圖中的兩個4點DFT，并畫出最後一層的2點DFT，可以得到最終的蝶形圖如下：

圖7

DIF的好處在于它的輸入資料不需要bit-reverse，這樣在實際應用中就不需要對輸入資料進行存儲，有利于硬體實作。

DIF radix-4 FFT (頻域抽取基-4 FFT)

最後再來看一個radix-4的DIF算法，這也是為了友善分析OAI中dft()函數的實作。

利用旋轉因子的性質，經過推導可以把頻域序列 X [ m ] X[m] X[m]劃分成如下4個子序列：

X [ 4 m ] = ∑ n = 0 N / 4 − 1 [ x [ n ] + x [ n + N / 4 ] + x [ n + N / 2 ] + x [ n + 3 N / 4 ] ] W N / 4 m n X [ 4 m + 1 ] = ∑ n = 0 N / 4 − 1 [ x [ n ] − j x [ n + N / 4 ] − x [ n + N / 2 ] + j x [ n + 3 N / 4 ] ] W N / 4 m n W N n X [ 4 m + 2 ] = ∑ n = 0 N / 4 − 1 [ x [ n ] − x [ n + N / 4 ] + x [ n + N / 2 ] − x [ n + 3 N / 4 ] ] W N / 4 m n W N 2 n X [ 4 m + 3 ] = ∑ n = 0 N / 4 − 1 [ x [ n ] + j x [ n + N / 4 ] − x [ n + N / 2 ] − j x [ n + 3 N / 4 ] ] W N / 4 m n W N 3 n m = 0 , 1 , . . . , N / 4 − 1 X[4m]=\sum_{n=0}^{N/4-1}\left[x[n]+x[n+N/4]+x[n+N/2]+x[n+3N/4]\right]W_{N/4}^{mn} \\ X[4m+1]=\sum_{n=0}^{N/4-1}\left[x[n]-jx[n+N/4]-x[n+N/2]+jx[n+3N/4]\right]W_{N/4}^{mn}W_N^{n} \\ X[4m+2]=\sum_{n=0}^{N/4-1}\left[x[n]-x[n+N/4]+x[n+N/2]-x[n+3N/4]\right]W_{N/4}^{mn}W_N^{2n} \\ X[4m+3]=\sum_{n=0}^{N/4-1}\left[x[n]+jx[n+N/4]-x[n+N/2]-jx[n+3N/4]\right]W_{N/4}^{mn}W_N^{3n} \\ m=0,1,...,N/4-1 X[4m]=n=0∑N/4−1​[x[n]+x[n+N/4]+x[n+N/2]+x[n+3N/4]]WN/4mn​X[4m+1]=n=0∑N/4−1​[x[n]−jx[n+N/4]−x[n+N/2]+jx[n+3N/4]]WN/4mn​WNn​X[4m+2]=n=0∑N/4−1​[x[n]−x[n+N/4]+x[n+N/2]−x[n+3N/4]]WN/4mn​WN2n​X[4m+3]=n=0∑N/4−1​[x[n]+jx[n+N/4]−x[n+N/2]−jx[n+3N/4]]WN/4mn​WN3n​m=0,1,...,N/4−1

如果令：

f 1 [ n ] = [ x [ n ] + x [ n + N / 4 ] + x [ n + N / 2 ] + x [ n + 3 N / 4 ] ] f 2 [ n ] [ x [ n ] − j x [ n + N / 4 ] − x [ n + N / 2 ] + j x [ n + 3 N / 4 ] ] W N n f 3 [ n ] = [ x [ n ] − x [ n + N / 4 ] + x [ n + N / 2 ] − x [ n + 3 N / 4 ] ] W N 2 n f 4 [ n ] = [ x [ n ] + j x [ n + N / 4 ] − x [ n + N / 2 ] − j x [ n + 3 N / 4 ] ] W N 3 n f_1[n]=\left[x[n]+x[n+N/4]+x[n+N/2]+x[n+3N/4]\right] \\ f_2[n]\left[x[n]-jx[n+N/4]-x[n+N/2]+jx[n+3N/4]\right]W_N^{n} \\ f_3[n]=\left[x[n]-x[n+N/4]+x[n+N/2]-x[n+3N/4]\right]W_N^{2n} \\ f_4[n]=\left[x[n]+jx[n+N/4]-x[n+N/2]-jx[n+3N/4]\right]W_N^{3n} f1​[n]=[x[n]+x[n+N/4]+x[n+N/2]+x[n+3N/4]]f2​[n][x[n]−jx[n+N/4]−x[n+N/2]+jx[n+3N/4]]WNn​f3​[n]=[x[n]−x[n+N/4]+x[n+N/2]−x[n+3N/4]]WN2n​f4​[n]=[x[n]+jx[n+N/4]−x[n+N/2]−jx[n+3N/4]]WN3n​

那麼 X [ 4 m ] , X [ 4 m + 1 ] , X [ 4 m + 2 ] , X [ 4 m + 3 ] X[4m],X[4m+1],X[4m+2],X[4m+3] X[4m],X[4m+1],X[4m+2],X[4m+3]就分别是 f 1 [ n ] , f 2 [ n ] , f 3 [ n ] , f 4 [ n ] f_1[n],f_2[n],f_3[n],f_4[n] f1​[n],f2​[n],f3​[n],f4​[n]在 0   N / 4 − 1 0~N/4-1 0 N/4−1上的傅立葉變換。把上述簡化應用于16點DFT，可以得到如下蝶形圖：

圖8

圖中沒有标出時域資料對應的系數。

Radix-4的FFT基本運算單元是4點DFT，把上圖中的4點DFT用蝶形圖替代，得到最終的16點蝶形圖，如下：

圖9

同樣沒有标出輸入資料的系數。

這就是DIF radix-4的16點FFT完整蝶形圖。

OAI中FFT的實作

下面我們來看OAI中對FFT的實作方式。

5G中因為使用的調制方式仍然是OFDM，是以同樣需要FFT。（OFDM和FFT之間的關系教材上應該有，如果有時間以後可以補一篇直覺的解釋）

隻是5G的帶寬更大，所需要的FFT size也更大，如下是LTE與5G的對比。

LTE：

圖10

5G:

圖11

可以看出來，5G在LTE的基礎上又一次提高了FFT size到4096。但是FFT就是FFT，在算法上還是一樣的。

OAI中實作了很多size的FFT，最基本的實作是16點DFT。4096點的DFT，最終也是通過調用dft16()實作的。是以，我們通過dft16()這個函數來分析OAI實作FFT的算法和代碼。

先看會用到的全局變量：

//計算過程中符号反轉用
const static int16_t conjugatedft[32] __attribute__((aligned(32))) = {-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1};
//旋轉因子*32767
const static int16_t tw16a[24] __attribute__((aligned(32))) = {32767,0,30272,12540,23169 ,23170,12539 ,30273,
                                                  32767,0,23169,23170,0     ,32767,-23170,23170,
                                                  32767,0,12539,30273,-23170,23170,-30273,-12539
                                                 };

const static int16_t tw16b[24] __attribute__((aligned(32))) = { 0,32767,-12540,30272,-23170,23169 ,-30273,12539,
                                                   0,32767,-23170,23169,-32767,0     ,-23170,-23170,
                                                   0,32767,-30273,12539,-23170,-23170,12539 ,-30273
                                                 };
           

dft16()函數的實作，第一階段與上述radix-4的分析保持一緻，可以參考圖9蝶形圖。不同點在于dft16()将輸入資料分成4組，每次計算一組，這樣無疑進一步加快了運算速度。

// First stage : 4 Radix-4 butterflies without input twiddles
  //第一階段的蝶形圖計算，沒有使用旋轉因子,每一個x128[]的元素順序存儲4個輸入資料的實部和虛部，即8個16位整數
  x02t    = _mm_adds_epi16(x128[0],x128[2]);   //x[0]+x[8], x[1]+x[9], x[2]+x[10], x[3]+x[11]
  x13t    = _mm_adds_epi16(x128[1],x128[3]);   //x[4]+x[12], x[5]+x[13], x[6]+x[14], x[7]+x[15]
  xtmp0   = _mm_adds_epi16(x02t,x13t);         //x[0]+x[4]+x[8]+x[12], x[1]+x[5]+x[9]+x[13], x[2]+x[6]+x[10]+x[14], x[3]+x[7]+x[11]+x[15],相當于上圖中第一個4點DFT的輸入
  xtmp2   = _mm_subs_epi16(x02t,x13t);         //x[0]+x[8]-x[4]-x[12], x[1]+x[9]-x[5]-x[13], x[2]+x[10]-x[6]-x[14], x[3]+x[11]-x[7]-x[15]
  x1_flip = _mm_sign_epi16(x128[1],*(__m128i*)conjugatedft);  //将x[4], x[5], x[6], x[7]的實部和虛的符号按照conjugatedft設定
  x1_flip = _mm_shuffle_epi8(x1_flip,complex_shuffle);        //将x[4], x[5], x[6], x[7]的實部和虛的符号的位置互換,與上一步合起來相當于*(-j)
  x3_flip = _mm_sign_epi16(x128[3],*(__m128i*)conjugatedft);
  x3_flip = _mm_shuffle_epi8(x3_flip,complex_shuffle);   //與上一步合起來相當于x[12], x[13], x[14], x[15]分别*(-j)
  x02t    = _mm_subs_epi16(x128[0],x128[2]);     //x[0]-x[8], x[1]-x[9], x[2]-x[10], x[3]-x[11]
  x13t    = _mm_subs_epi16(x1_flip,x3_flip);     //-jx[4]+jx[12], -jx[5]+jx[13], -jx[6]+jx[14], -jx[7]+jx[15]
  xtmp1   = _mm_adds_epi16(x02t,x13t);  // x0 + x1f - x2 - x3f=(x[0],x[1],x[2],x[3])-j(x[4],x[5],x[6],x[7])-(x[8],x[9],x[10],x[11])+j(x[12],x[13],x[14],x[15])
  xtmp3   = _mm_subs_epi16(x02t,x13t);  // x0 - x1f - x2 + x3f=(x[0],x[1],x[2],x[3])+j(x[4],x[5],x[6],x[7])-(x[8],x[9],x[10],x[11])-j(x[12],x[13],x[14],x[15])
           

dft16()函數，在完成第一階段的計算後，為了使輸出結果保持順序，同時也是為了重複使用第一階段的計算過程，把第一階段的輸出結果進行了重新排序。這裡的實作有些特殊，如果按照正常時域抽取DIT，這個重排序應該放到第一階段之前; 而如果是頻域抽取DIF，這個重排序應該放到最後結果出來之後。特殊實作的結果，就是代碼重用性更佳。

參照radix-4的蝶形圖，重新排序第一階段輸出結果的流程如下：

圖12

可見，重新排序中間結果後，第二階段的計算過程與第一階段完全相同，是以可以複用第一階段的代碼。

對中間結果重新排序的代碼：

//正常輸出應該是bit-reverse的，但是這裡調整了stage1的輸出序列的順序，使得stage2的輸出序列是順序的
  ytmp0   = _mm_unpacklo_epi32(xtmp0,xtmp1);
  ytmp1   = _mm_unpackhi_epi32(xtmp0,xtmp1);
  ytmp2   = _mm_unpacklo_epi32(xtmp2,xtmp3);
  ytmp3   = _mm_unpackhi_epi32(xtmp2,xtmp3);
  xtmp0   = _mm_unpacklo_epi64(ytmp0,ytmp2);
  xtmp1   = _mm_unpackhi_epi64(ytmp0,ytmp2);
  xtmp2   = _mm_unpacklo_epi64(ytmp1,ytmp3);
  xtmp3   = _mm_unpackhi_epi64(ytmp1,ytmp3);
           

然後是第二階段蝶形圖計算過程：

// Second stage : 4 Radix-4 butterflies with input twiddles
  xtmp1 = packed_cmult2(xtmp1,tw16a_128[0],tw16b_128[0]);
  xtmp2 = packed_cmult2(xtmp2,tw16a_128[1],tw16b_128[1]);
  xtmp3 = packed_cmult2(xtmp3,tw16a_128[2],tw16b_128[2]);

  x02t    = _mm_adds_epi16(xtmp0,xtmp2);
  x13t    = _mm_adds_epi16(xtmp1,xtmp3);
  y128[0] = _mm_adds_epi16(x02t,x13t);
  y128[2] = _mm_subs_epi16(x02t,x13t);
  x1_flip = _mm_sign_epi16(xtmp1,*(__m128i*)conjugatedft);
  x1_flip = _mm_shuffle_epi8(x1_flip,_mm_set_epi8(13,12,15,14,9,8,11,10,5,4,7,6,1,0,3,2));
  x3_flip = _mm_sign_epi16(xtmp3,*(__m128i*)conjugatedft);
  x3_flip = _mm_shuffle_epi8(x3_flip,_mm_set_epi8(13,12,15,14,9,8,11,10,5,4,7,6,1,0,3,2));
  x02t    = _mm_subs_epi16(xtmp0,xtmp2);
  x13t    = _mm_subs_epi16(x1_flip,x3_flip);
  y128[1] = _mm_adds_epi16(x02t,x13t);  // x0 + x1f - x2 - x3f
  y128[3] = _mm_subs_epi16(x02t,x13t);  // x0 - x1f - x2 + x3f
           

可以看出來，除了需要乘上旋轉因子之外，剩下的計算代碼與第一階段的計算過程完全一緻。

OAI中更高階的DFT計算是中dft16()的基礎上采用時域抽取的方法實作的，以dft4096為例：

void dft4096(int16_t *x,int16_t *y,int scale)
{

  simd_q15_t xtmp[1024],ytmp[1024],*tw4096_128p=(simd_q15_t *)tw4096,*x128=(simd_q15_t *)x,*y128=(simd_q15_t *)y,*y128p=(simd_q15_t *)y;
  simd_q15_t *ytmpp = &ytmp[0];
  int i,j;

  for (i=0,j=0; i<1024; i+=4,j++) {
    transpose16_ooff(x128+i,xtmp+j,256);   //輸入資料重排序
  }

  //計算4個dft2014
  dft1024((int16_t*)(xtmp),(int16_t*)(ytmp),1);
  dft1024((int16_t*)(xtmp+256),(int16_t*)(ytmp+256),1);
  dft1024((int16_t*)(xtmp+512),(int16_t*)(ytmp+512),1);
  dft1024((int16_t*)(xtmp+768),(int16_t*)(ytmp+768),1);

  //計算蝶形圖
  for (i=0; i<256; i++) {
    bfly4(ytmpp,ytmpp+256,ytmpp+512,ytmpp+768,
          y128p,y128p+256,y128p+512,y128p+768,
          tw4096_128p,tw4096_128p+256,tw4096_128p+512);
    tw4096_128p++;
    y128p++;
    ytmpp++;
  }

           

注意，這裡的旋轉因子tw4096其實就是在之前我們提到的初始化函數中初始化的：

phy_init_nr_top-》init_dfts-》 init_rad4(4096,tw4096);

這個函數的内容如下：

void init_rad4(int N,int16_t *tw) {

  int16_t *twa = tw;
  int16_t *twb = twa+(N/2);
  int16_t *twc = twb+(N/2);
  int i;

  for (i=0;i<(N/4);i++) { 
    *twa = (int16_t)round(32767.0*cos(2*M_PI*i/N)); twa++;    //旋轉因子都是×32767後的結果
    *twa = -(int16_t)round(32767.0*sin(2*M_PI*i/N)); twa++;
    *twb = (int16_t)round(32767.0*cos(2*M_PI*2*i/N)); twb++;
    *twb = -(int16_t)round(32767.0*sin(2*M_PI*2*i/N)); twb++;
    *twc = (int16_t)round(32767.0*cos(2*M_PI*3*i/N)); twc++;
    *twc = -(int16_t)round(32767.0*sin(2*M_PI*3*i/N)); twc++;
  }
}
           

放一個基-4的DIT FFT蝶形圖，對照代碼直覺上感受一下吧：

參考來源：

http://www.sharetechnote.com/

http://www.cmlab.csie.ntu.edu.tw/cml/dsp/training/coding/transform/fft.html

http://face2ai.com/DIP-2-3-FFT算法了解與c語言的實作/

https://www.semanticscholar.org/paper/The-Hybrid-Architecture-Parallel-Fast-Fourier-Palmer-Nelson/e6867a85e61e1228063b4a6e40d0af8c8f20de08

https://www.keysight.com/upload/cmc_upload/All/5GRadioValidation.pdf

http://www.iraj.in/journal/journal_file/journal_pdf/6-33-139038874964-67.pdf