互質因子算法

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

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

PFA 也與 nested Winograd FFT 算法密切相關，後者使用更為精巧的二維摺積技巧分解成 _N_1 \* _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_ = _N_1_N_2，_N_1 與 _N_2 兩者互質，然後把輸入 _n_ 和輸出 _k_ 一一對應到

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

因 _N_1 與 _N_2 互質，故根據最大公因數表現定理，對每個 _n_ 都存在滿足上式的整數 _n_1 與 _n_2，且在同餘 _N_ 之下 _n_1 可以調整至 0～_N_1 –1 之間，_n_2 可以調整至 0～_N_2 –1 之間。並根據同餘理論易知滿足上式且在以上範圍內的整數 _n_1 與 _n_2 是唯一的。這稱為 _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
$$

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

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

其中 _N_1-1 表示 _N_1 在模 _N_2 之下的反元素，_N_2-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}}}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_1 與 _N_2 )。剩下的部分變成

　
: $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_2 與 _N_1 的 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_{1}N_{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 為例，有兩種可能，_N_1 = 2, _N_2 = 3 或 _N_1 = 3, _N_2 = 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_ = _n_1 + _n_2_N_1 , _k_ = _k_1_N_2 + _k_2，則

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

## 相關條目

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

## 注釋

## 參考文獻

1. P. Duhamel and M. Vetterli, Fast Fourier transforms: a tutorial review and a state of the art, Signal Processing, 1990, **19**: 259–299

## 外部連結

* fft note by Burrus（頁面存檔備份，存於網際網路檔案館）
* cnx（頁面存檔備份，存於網際網路檔案館）