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

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

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

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

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

: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 = N1N2，N1與N2兩者互質，然後把輸入n 和輸出k 一一對應到

:n = n_1 N_2 + n_2 N_1 \mod N

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

:k = k_1 \mod N_1
:k = k_2 \mod N_2

舉例來說:

如果N=15, N_1=5, N_2=3, n=0,1,2,...,12, 13,14, 對於任一 n 都可以對應到

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

0=0\centerdot N_2+0\centerdot N_1 \mod 15

1=2\centerdot N_2+2\centerdot N_1 \mod 15

2=4\centerdot N_2+1\centerdot N_1 \mod 15

3=1\centerdot N_2+0\centerdot N_1 \mod 15

4=3\centerdot N_2+2\centerdot N_1 \mod 15

5=0\centerdot N_2+1\centerdot N_1 \mod 15

6=2\centerdot N_2+0\centerdot N_1 \mod 15

7=4\centerdot N_2+2\centerdot N_1 \mod 15

8=1\centerdot N_2+1\centerdot N_1 \mod 15

9=3\centerdot N_2+0\centerdot N_1 \mod 15

10=0\centerdot N_2+2\centerdot N_1 \mod 15

11=2\centerdot N_2+1\centerdot N_1 \mod 15

12=4\centerdot N_2+0\centerdot N_1 \mod 15

13=1\centerdot N_2+2\centerdot N_1 \mod 15

14=3\centerdot N_2+1\centerdot N_1 \mod 15

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

:k = k_1 N_2^{-1} N_2 + k_2 N_1^{-1} N_1 \mod N

其中N1-1表示N1在模N2之下的反元素，N2-1反之。

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

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

DFT re-expression
表示方法一:

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

: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} k n_1} e^{-\frac{2\pi i}{N_2} k n_2} = e^{-\frac{2\pi i}{N_1} k_1 n_1} e^{-\frac{2\pi i}{N_2} k_2 n_2}

( 因為e2πi = 1，所以兩個指數的k 部份可以分別模N1與N2 )。剩下的部分變成

: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 }.

則內部和外部的總和分別轉換成大小為N2與N1的DFT。

表示方法二:

如果令 k=k_1 N_2 + k_2 N_1 \quad for\quad k=0,1,...,N-1,

令 n = ((n_1 N_2 + n_2 N_1))_N，(\cdot)_N相當於取 N的餘數，n_1 = 0,\dots,N_1-1 , n_2 = 0,\dots,N_2-1

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)}

=\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)}

=\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)}

=\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)}.

對於每一個 n_2 都要做一個 N_1 點的 DFT，而因為 n_2 = 0,\dots,N_2-1 有 N_2個，所以需要 N_2 個 N_1 點 DFT,

對於每一組((k_1N_2))_{N_1}都要做一個 N_2 點的 DFT，而因為 N_2為常數，k_1 = 0,\dots,N_1-1 有 N_1 個，所以需要 N_1 個 N_2 點 DFT，

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

假設N_1 點的 DFT 需要 M_1個乘法器,

假設N_2 點的 DFT 需要 M_2個乘法器,

則總共需要 N_2 M_1 + N_1 M_2個乘法器。

範例
以N = 6為例，有兩種可能，N1 = 2, N2 = 3或N1 = 3, N2 = 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的定義

:X_k = \sum_{n=0}^{N-1} x_n e^{-\frac{2\pi i}{N} nk }
\qquad
k = 0,\dots,N-1

中代入n = n1 + n2N1 , k = k1N2 + k2，則

: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}

:X_{k_1N_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多了一些要乘的因子e^{-\frac{2\pi i}{N} n_1 k_2} (稱為twiddle factors )，但index較為簡單，且適用於任何N1、N2。在J. Cooley稍後發表的關於FFT歷史探討的論文中使用N = 24點FFT為例，顯示兩種作法在index結構上的不同。

N1 與 N2 不互質情況
N = N_1 * N_2，若N_1 與 N_2 並不互質，仍可透過拆解得到運算所需的乘法數量。
離散傅立葉變換（DFT）的定義如下:
F[m] = \sum_{n=0}^{N-1} f[n] e^{-\frac{2\pi i}{N} mn }
\qquad
m,n = 0,\dots,N-1

可以拆解成N_2 個 N_1-point DFTs 和 N_1 個 N_2-point DFT_s，以及twiddle factor。

算法
令 n = n_1N_1 + n_2，m = m_1 + m_2N_2。 n_1, m_1 = 0, 1, 2, ...., N_2-1, n_2, m_2 = 0, 1, 2, ...., N_1-1

可得出下列式子

F[m_1+m_2N_2] = \sum_{n=0}^{N-1} f[n_1N_1+n_2] e^{-\frac{2\pi i}{N} (m_1+m_2N_2)(n_1N_1+n_2) }
\qquad

\qquad\qquad\qquad \ \ = \sum_{n=0}^{N-1} f[n_1N_1+n_2] e^{-\frac{2\pi i}{N_1N_2} (m_1n_1N_1+m_1n_2+m_2n_1N_1N_2+m_2n_2N_2) }

\qquad\qquad\qquad \ \ = \sum_{n=0}^{N-1} f[n_1N_1+n_2] e^{-\frac{2\pi i}{N_2}(m_1n_1) } e^{-\frac{2\pi i}{N_1}(m_2n_2) } e^{-\frac{2\pi i}{N_1N_2}(m_1n_2) }

\qquad\qquad\qquad \ \ = \sum_{n_2=0}^{P_1-1} \Bigl( \sum_{n_1=0}^{N_2-1} f[n_1N_1+n_2] e^{-\frac{2\pi i}{N_2}(m_1n_1) } \Bigr) e^{-\frac{2\pi i}{N_1}(m_2n_2) } e^{-\frac{2\pi i}{N_1N_2}(m_1n_2) }

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

由於n_1, m_1 = 0, 1, 2, ...., N_2-1, n_2, m_2 = 0, 1, 2, ...., N_1-1

twiddle factors的數量為 N = N_1 N_2，移去 m_1 = 0 以及 n_2 = 0 的情況，至多需要 (N_1-1) (N_2-1) 個twiddle factors。

分階段解析
Step1
令 g[n_1, n_2] = f[n_1N_1+n_2]

Step2
固定 n_2，對 n_1
做 N_2-point DFT。

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_1n_1) }

由於n_2 = 0, 1, 2, ...., N_1-1，所以有 N_1 個 N_2-point DFT_s

Step3
G_2[m_1, n_2] = G_1[m_1, n_2]e^{-\frac{2\pi i}{N_1N_2}(m_1n_2) } (twiddle factors)

Step4
固定 n_2，對 n_2
做 N_1-point DFT。

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_2n_2) }

由於m_1 = 0, 1, 2, ...., N_2-1，所以有 N_2 個 N_1-point DFT_s

Step5
F[m_1+m_2N_2] = G_3[m_1, m_2]

實際乘法量計算
N = N_1 * N_2

假設 N_1-point DFT 的乘法量為 B_1
， N_2-point DFT 的乘法量為 B_2

且 m_1n_2 中，

有 D_1
個值不為 N/12
以及 N/8
的倍數

有 D_2
個值為 N/12
或 N/8
的倍數，但不為 N/4
的倍數

則 N-point DFT 的乘法量數量為 N_2B_1 + N_1B_2 + 3D_1 + 2D_2 個。

補充說明
a

為一複數，則 a*e^{i\theta}

一般需要 3個乘法，但若 \theta

= \pi / 4

或 \theta

= \pi / 3

則只需要2個乘法。

簡單範例
16-point DFT, 16=4*4

4-point needs 0 \ MUL_s

m_1n_2 = 0, 0, 0, 0, 0, 1 ,2, 3, 0 ,2, 4, 6, 0, 3, 6, 9

D_1
= 4 (1, 3, 3, 9), D_2
= 4 (2, 2, 6, 6)

總計需要 40 + 40 + 34 + 24 = 20 MUL_s

相關條目
*快速傅立葉變換
*中國剩餘定理
*Bézout引理

注釋
參考文獻
#

Jian-Jiun Ding, Advanced Digital Signal Processing class note,the Department of Electrical Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2025

外部連結
*[http://www.jjj.de/fft/fftnote.txt fft note by Burrus]
*[http://cnx.org/content/m12033/latest/ cnx]