互质因子算法

互质因子算法（Prime-factor FFT algorithm, PFA），又称为Good-Thomas算法，是一种快速傅立叶变换（FFT），把N = N1N2大小的离散傅立叶变换重新表示为N1 * N2大小的二维离散傅立叶变换，其中N1与N2需互质。变成N1和N2大小的傅立叶变换后，可以继续递回使用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 * N2的转换。因而一些较古老的论文把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 } = <pre> \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 }. </pre>$$

则内部和外部的总和分别转换成大小为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。

<image>

<image>

第一种情形所产生的流程图如左图所示。先做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} = <pre> \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 } </pre>$$

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

总计需要 $4*0 + 4*0 + 3*4 + 2*4 = 20 MUL_s$

## 参考文献

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