互質因子算法

互質因子算法（Prime-factor FFT algorithm, PFA），又稱為Good-Thomas算法^{[1] [2]}，是一種快速傅立葉變換（FFT），把N = N₁N₂大小的離散傅立葉變換重新表示為N₁ * N₂大小的二維離散傅立葉變換，其中N₁與N₂需互質。變成N₁和N₂大小的傅立葉變換後，可以繼續遞迴使用PFA，或用其他快速傅立葉變換算法來計算。

較流行的Cooley-Tukey算法經由mixed-radix一般化後，也是把N = N₁N₂大小的離散傅立葉變換分割為N₁和N₂大小的轉換，但和互質因子算法 (PFA)作法並不相同，不應混淆。Cooley-Tukey算法的N₁與N₂不需互質，可以是任何整數。然而有個缺點是比PFA多出一些乘法，和單位根twiddle factors相乘。相對的，PFA的缺點則是N₁與N₂需互質 (例如N 是2次方就不適用)，而且要藉由中國剩餘定理來進行較複雜的re-indexing。互質因子算法 (PFA)可以和mixed-radix Cooley-Tukey算法相結合，前者將N 分解為互質的因數，後者則用在重複質因數上。

PFA也與nested Winograd FFT算法密切相關，後者使用更為精巧的二維摺積技巧分解成N₁ * N₂的轉換。因而一些較古老的論文把Winograd算法稱為PFA FFT。

儘管PFA和Cooley-Tukey算法並不相同，但有趣的是Cooley和Tukey在他們1965年發表的有名的論文中，沒有發覺到高斯和其他人更早的研究，只引用Good在1958年發表的PFA作為前人的FFT結果。剛開始的時候人們對這兩種作法是否不同有點困惑。

算法

離散傅立葉變換（DFT）的定義如下:

${\displaystyle X_{k}=\sum _{n=0}^{N-1}x_{n}e^{-{\frac {2\pi i}{N}}nk}\qquad k=0,\dots ,N-1}$

PFA將輸入和輸出re-indexing，代入DFT公式後轉換成二維DFT。

Re-indexing

設N = N₁N₂，N₁與N₂兩者互質，然後把輸入n 和輸出k 一一對應到

${\displaystyle n=n_{1}N_{2}+n_{2}N_{1}\mod N}$

因N₁與N₂ 互質，故根據最大公因數表現定理，對每個n 都存在滿足上式的整數n₁與n₂，且在同餘N 之下n₁可以調整至0～N₁ –1之間，n₂可以調整至0～N₂ –1之間。並根據同餘理論易知滿足上式且在以上範圍內的整數n₁與n₂是唯一的。這稱為Ruritanian 映射 (或Good's 映射)，

${\displaystyle k=k_{1}\mod N_{1}}$

${\displaystyle k=k_{2}\mod N_{2}}$

舉例來說:

如果${\displaystyle N=15,N_{1}=5,N_{2}=3,n=0,1,2,...,12,13,14,}$對於任一 ${\displaystyle n}$ 都可以對應到

${\displaystyle n=n_{1}N_{2}+n_{2}N_{1}\mod N,n_{1}=0,1,...,N_{1}-1,n_{2}=0,1,...,N_{2}-1}$

${\displaystyle 0=0\centerdot N_{2}+0\centerdot N_{1}\mod 15}$

${\displaystyle 1=2\centerdot N_{2}+2\centerdot N_{1}\mod 15}$

${\displaystyle 2=4\centerdot N_{2}+1\centerdot N_{1}\mod 15}$

${\displaystyle 3=1\centerdot N_{2}+0\centerdot N_{1}\mod 15}$

${\displaystyle 4=3\centerdot N_{2}+2\centerdot N_{1}\mod 15}$

${\displaystyle 5=0\centerdot N_{2}+1\centerdot N_{1}\mod 15}$

${\displaystyle 6=2\centerdot N_{2}+0\centerdot N_{1}\mod 15}$

${\displaystyle 7=4\centerdot N_{2}+2\centerdot N_{1}\mod 15}$

${\displaystyle 8=1\centerdot N_{2}+1\centerdot N_{1}\mod 15}$

${\displaystyle 9=3\centerdot N_{2}+0\centerdot N_{1}\mod 15}$

${\displaystyle 10=0\centerdot N_{2}+2\centerdot N_{1}\mod 15}$

${\displaystyle 11=2\centerdot N_{2}+1\centerdot N_{1}\mod 15}$

${\displaystyle 12=4\centerdot N_{2}+0\centerdot N_{1}\mod 15}$

${\displaystyle 13=1\centerdot N_{2}+2\centerdot N_{1}\mod 15}$

${\displaystyle 14=3\centerdot N_{2}+1\centerdot N_{1}\mod 15}$

因N₁與N₂ 互質，故根據中國剩餘定理，對於每組 ( k₁ , k₂ ) (其中k₁在0～N₁ – 1之間, k₂在0～N₂ – 1之間)，都有存在且唯一的k 在0～N - 1之間且滿足上兩式。這稱為 CRT 映射。 CRT 映射的另一種表示法如下

${\displaystyle k=k_{1}N_{2}^{-1}N_{2}+k_{2}N_{1}^{-1}N_{1}\mod N}$

其中N₁^-1表示N₁在模N₂之下的反元素，N₂^-1反之。

( 也可以改成對輸入n 用 CRT 映射以及對輸出k 用Ruritanian 映射)

對於有效re-indexing (理想上是達到原地)的方法有許多研究^[3]，以減少耗費時間的模運算。

DFT re-expression

表示方法一:

將以上的re-indexing代入DFT公式裡指數部分的nk 之中，

${\displaystyle e^{-{\frac {2\pi i}{N}}nk}=e^{-{\frac {2\pi i}{N}}(n_{1}N_{2}+n_{2}N_{1})k}=e^{-{\frac {2\pi i}{N_{1}}}kn_{1}}e^{-{\frac {2\pi i}{N_{2}}}kn_{2}}=e^{-{\frac {2\pi i}{N_{1}}}k_{1}n_{1}}e^{-{\frac {2\pi i}{N_{2}}}k_{2}n_{2}}}$

( 因為e^2πi = 1，所以兩個指數的k 部份可以分別模N₁與N₂ )。剩下的部分變成

${\displaystyle X_{k_{1},k_{2}}=\sum _{n_{1}=0}^{N_{1}-1}\left(\sum _{n_{2}=0}^{N_{2}-1}x_{n_{1}N_{2}+n_{2}N_{1}}e^{-{\frac {2\pi i}{N_{2}}}n_{2}k_{2}}\right)e^{-{\frac {2\pi i}{N_{1}}}n_{1}k_{1}}.}$

則內部和外部的總和分別轉換成大小為N₂與N₁的DFT。

表示方法二:

如果令 ${\displaystyle k=k_{1}N_{2}+k_{2}N_{1}\quad for\quad k=0,1,...,N-1,}$

令 ${\displaystyle n=((n_{1}N_{2}+n_{2}N_{1}))_{N}}$，${\displaystyle (\cdot )_{N}}$相當於取 ${\displaystyle N}$的餘數，${\displaystyle n_{1}=0,\dots ,N_{1}-1}$, ${\displaystyle n_{2}=0,\dots ,N_{2}-1}$

${\displaystyle X[((k_{1}N_{2}+k_{2}N_{1}))_{N}]=\sum _{n=0}^{N-1}x[((n_{1}N_{2}+n_{2}N_{1}))_{N}]e^{-j{\frac {2\pi }{N_{2}N_{1}}}(k_{1}N_{2}+k_{2}N_{1})(n_{1}N_{2}+n_{2}N_{1})}}$

${\displaystyle =\sum _{n=0}^{N-1}x[((n_{1}N_{2}+n_{2}N_{1}))_{N}]e^{-j{\frac {2\pi }{N_{2}N_{1}}}(k_{1}n_{1}N_{2}N_{2}+k_{2}n_{2}N_{1}N_{1}+k_{1}n_{2}N_{2}N_{1}+k_{2}n_{1}N_{1}N_{2})}}$

${\displaystyle =\sum _{n=0}^{N-1}x[((n_{1}N_{2}+n_{2}N_{1}))_{N}]e^{-j{\frac {2\pi }{N_{1}}}(k_{1}n_{1}N_{2})}e^{-j{\frac {2\pi }{N_{2}}}(k_{2}n_{2}N_{1})}}$

${\displaystyle =\sum _{n_{2}=0}^{N_{2}-1}\{\sum _{n_{1}=0}^{N_{1}-1}x[((n_{1}N_{2}+n_{2}N_{1}))_{N}]e^{-j{\frac {2\pi }{N_{1}}}(k_{1}n_{1}N_{2})}\}e^{-j{\frac {2\pi }{N_{2}}}(k_{2}n_{2}N_{1})}.}$

對於每一個 ${\displaystyle n_{2}}$ 都要做一個 ${\displaystyle N_{1}}$ 點的 ${\displaystyle DFT}$，而因為 ${\displaystyle n_{2}=0,\dots ,N_{2}-1}$有 ${\displaystyle N_{2}}$個，所以需要 ${\displaystyle N_{2}}$ 個 ${\displaystyle N_{1}}$ 點 ${\displaystyle DFT}$,

對於每一組${\displaystyle ((k_{1}N_{2}))_{N_{1}}}$都要做一個 ${\displaystyle N_{2}}$ 點的 ${\displaystyle DFT}$，而因為 ${\displaystyle N_{2}}$為常數，${\displaystyle k_{1}=0,\dots ,N_{1}-1}$有 ${\displaystyle N_{1}}$ 個，所以需要 ${\displaystyle N_{1}}$ 個 ${\displaystyle N_{2}}$ 點 ${\displaystyle DFT}$，

因此如果要計算複雜度，可以乘法器的數量當作考量,

假設${\displaystyle N_{1}}$ 點的 ${\displaystyle DFT}$ 需要 ${\displaystyle M_{1}}$個乘法器,

假設${\displaystyle N_{2}}$ 點的 ${\displaystyle DFT}$ 需要 ${\displaystyle M_{2}}$個乘法器,

則總共需要 ${\displaystyle N_{2}M_{1}+N_{1}M_{2}}$個乘法器。

範例

以N = 6為例，有兩種可能，N₁ = 2, N₂ = 3或N₁ = 3, N₂ = 2。

N₁ = 2, N₂ = 3

N₁ = 3, N₂ = 2

第一種情形所產生的流程圖如左圖所示。先做2次3點DFT後再做3次2點DFT。

第二種情形所產生的流程圖如右圖所示。先做3次2點DFT後再做2次3點DFT。

其中2點DFT的部份因構造單純，皆以交錯的蝴蝶圖來顯示。

可以看出即使在這個簡單的例子中，輸入和輸出的index也都經過有點複雜的重新排列。

與Cooley-Tukey算法的比較

如首段所述，Cooley-Tukey算法和互質因子算法 (PFA)曾被誤認為很類似。兩者皆有各自優點可適用於不同狀況，因此分辨它們的不同是很重要的。在1965年著名的論文中發表的Cooley-Tukey算法，是在DFT的定義

${\displaystyle X_{k}=\sum _{n=0}^{N-1}x_{n}e^{-{\frac {2\pi i}{N}}nk}\qquad k=0,\dots ,N-1}$

中代入n = n₁ + n₂N₁ , k = k₁N₂ + k₂，則

${\displaystyle e^{-{\frac {2\pi i}{N}}nk}=e^{-{\frac {2\pi i}{N}}(n_{1}+n_{2}N_{1})(k_{1}N_{2}+k_{2})}=e^{-{\frac {2\pi i}{N_{1}}}n_{1}k_{1}}e^{-{\frac {2\pi i}{N}}n_{1}k_{2}}e^{-{\frac {2\pi i}{N_{2}}}n_{2}k_{2}}}$

${\displaystyle X_{k_{1}N_{2}+k_{2}}=\sum _{n_{1}=0}^{N_{1}-1}\left(\sum _{n_{2}=0}^{N_{2}-1}x_{n_{1}+n_{2}N_{1}}e^{-{\frac {2\pi i}{N_{2}}}n_{2}k_{2}}\right)e^{-{\frac {2\pi i}{N}}n_{1}k_{2}}e^{-{\frac {2\pi i}{N_{1}}}n_{1}k_{1}}}$

比PFA多了一些要乘的因子${\displaystyle e^{-{\frac {2\pi i}{N}}n_{1}k_{2}}}$ (稱為twiddle factors )，但index較為簡單，且適用於任何N₁、N₂。在J. Cooley稍後發表的關於FFT歷史探討的論文^[4]中使用N = 24點FFT為例，顯示兩種作法在index結構上的不同。

N1 與 N2 不互質情況

${\displaystyle N=N_{1}*N_{2}}$，若${\displaystyle N_{1}}$ 與 ${\displaystyle N_{2}}$ 並不互質，仍可透過拆解得到運算所需的乘法數量。 離散傅立葉變換（DFT）的定義如下: ${\displaystyle F[m]=\sum _{n=0}^{N-1}f[n]e^{-{\frac {2\pi i}{N}}mn}\qquad m,n=0,\dots ,N-1}$ 可以拆解成${\displaystyle N_{2}}$ 個 ${\displaystyle N_{1}-point}$ DFTs 和 ${\displaystyle N_{1}}$ 個 ${\displaystyle N_{2}-point}$ ${\displaystyle DFT_{s}}$，以及twiddle factor。

算法

令 ${\displaystyle n=n_{1}N_{1}+n_{2}}$，${\displaystyle m=m_{1}+m_{2}N_{2}}$。 ${\displaystyle n_{1},m_{1}=0,1,2,....,N_{2}-1}$, ${\displaystyle n_{2},m_{2}=0,1,2,....,N_{1}-1}$

可得出下列式子

${\displaystyle F[m_{1}+m_{2}N_{2}]=\sum _{n=0}^{N-1}f[n_{1}N_{1}+n_{2}]e^{-{\frac {2\pi i}{N}}(m_{1}+m_{2}N_{2})(n_{1}N_{1}+n_{2})}\qquad }$

${\displaystyle \qquad \qquad \qquad \ \ =\sum _{n=0}^{N-1}f[n_{1}N_{1}+n_{2}]e^{-{\frac {2\pi i}{N_{1}N_{2}}}(m_{1}n_{1}N_{1}+m_{1}n_{2}+m_{2}n_{1}N_{1}N_{2}+m_{2}n_{2}N_{2})}}$

${\displaystyle \qquad \qquad \qquad \ \ =\sum _{n=0}^{N-1}f[n_{1}N_{1}+n_{2}]e^{-{\frac {2\pi i}{N_{2}}}(m_{1}n_{1})}e^{-{\frac {2\pi i}{N_{1}}}(m_{2}n_{2})}e^{-{\frac {2\pi i}{N_{1}N_{2}}}(m_{1}n_{2})}}$

${\displaystyle \qquad \qquad \qquad \ \ =\sum _{n_{2}=0}^{P_{1}-1}{\Bigl (}\sum _{n_{1}=0}^{N_{2}-1}f[n_{1}N_{1}+n_{2}]e^{-{\frac {2\pi i}{N_{2}}}(m_{1}n_{1})}{\Bigr )}e^{-{\frac {2\pi i}{N_{1}}}(m_{2}n_{2})}e^{-{\frac {2\pi i}{N_{1}N_{2}}}(m_{1}n_{2})}}$

上述的 ${\displaystyle e^{-{\frac {2\pi i}{N_{1}N_{2}}}(m_{1}n_{2})}}$ 被稱為twiddle factor，需要額外的乘法去做運算。

由於${\displaystyle n_{1},m_{1}=0,1,2,....,N_{2}-1}$, ${\displaystyle n_{2},m_{2}=0,1,2,....,N_{1}-1}$

twiddle factors的數量為 ${\displaystyle N=N_{1}*N_{2}}$，移去 ${\displaystyle m_{1}=0}$ 以及 ${\displaystyle n_{2}=0}$ 的情況，至多需要 ${\displaystyle (N_{1}-1)*(N_{2}-1)}$個twiddle factors。

分階段解析

Step1

令 ${\displaystyle g[n_{1},n_{2}]=f[n_{1}N_{1}+n_{2}]}$

Step2

固定 ${\displaystyle n_{2}}$，對 ${\displaystyle n_{1}}$ 做 ${\displaystyle N_{2}-point}$ ${\displaystyle DFT}$。

${\displaystyle G_{1}[m_{1},n_{2}]=\sum _{n_{1}=0}^{N_{2}-1}g[n_{1},n_{2}]e^{-{\frac {2\pi i}{N_{2}}}(m_{1}n_{1})}}$

由於${\displaystyle n_{2}=0,1,2,....,N_{1}-1}$，所以有 ${\displaystyle N_{1}}$ 個 ${\displaystyle N_{2}-point}$ ${\displaystyle DFT_{s}}$

Step3

${\displaystyle G_{2}[m_{1},n_{2}]=G_{1}[m_{1},n_{2}]e^{-{\frac {2\pi i}{N_{1}N_{2}}}(m_{1}n_{2})}}$ (twiddle factors)

Step4

固定 ${\displaystyle n_{2}}$，對 ${\displaystyle n_{2}}$ 做 ${\displaystyle N_{1}-point}$ ${\displaystyle DFT}$。

${\displaystyle G_{3}[m_{1},m_{2}]=\sum _{n_{2}=0}^{N_{1}-1}g[m_{1},n_{2}]e^{-{\frac {2\pi i}{N_{1}}}(m_{2}n_{2})}}$

由於${\displaystyle m_{1}=0,1,2,....,N_{2}-1}$，所以有 ${\displaystyle N_{2}}$ 個 ${\displaystyle N_{1}-point}$ ${\displaystyle DFT_{s}}$

Step5

${\displaystyle F[m_{1}+m_{2}N_{2}]=G_{3}[m_{1},m_{2}]}$

實際乘法量計算

${\displaystyle N=N_{1}*N_{2}}$

假設 ${\displaystyle N_{1}-point}$ DFT 的乘法量為 ${\displaystyle B_{1}}$， ${\displaystyle N_{2}-point}$ ${\displaystyle DFT}$ 的乘法量為 ${\displaystyle B_{2}}$

且 ${\displaystyle m_{1}n_{2}}$ 中，

有 ${\displaystyle D_{1}}$ 個值不為 ${\displaystyle N/12}$ 以及 ${\displaystyle N/8}$ 的倍數

有 ${\displaystyle D_{2}}$ 個值為 ${\displaystyle N/12}$ 或 ${\displaystyle N/8}$ 的倍數，但不為 ${\displaystyle N/4}$ 的倍數

則 ${\displaystyle N-point}$ ${\displaystyle DFT}$ 的乘法量數量為 ${\displaystyle N_{2}B_{1}+N_{1}B_{2}+3D_{1}+2D_{2}}$ 個。

補充說明:

${\displaystyle a}$為一複數，則 ${\displaystyle a*e^{i\theta }}$ 一般需要 3個乘法，但若 ${\displaystyle \theta }$ = ${\displaystyle \pi /4}$ 或 ${\displaystyle \theta }$ = ${\displaystyle \pi /3}$ 則只需要2個乘法。

簡單範例

${\displaystyle 16-point}$ ${\displaystyle DFT}$, ${\displaystyle 16=4*4}$

${\displaystyle 4-point}$ needs ${\displaystyle 0\ MUL_{s}}$

${\displaystyle m_{1}n_{2}}$ = ${\displaystyle 0,0,0,0,0,1,2,3,0,2,4,6,0,3,6,9}$

${\displaystyle D_{1}}$ = 4 (1, 3, 3, 9), ${\displaystyle D_{2}}$ = 4 (2, 2, 6, 6)

總計需要 ${\displaystyle 4*0+4*0+3*4+2*4=20MUL_{s}}$

相關條目

• 快速傅立葉變換
• 中國剩餘定理
• Bézout引理