DSP学习笔记（四）| Twiddle Factor

Taeriyakki
120 阅读2分钟

Twiddle Factor

What are twiddle factors?

Twiddle Factor（用字母 W 表示）是一组用于加快 DFT 和 IDFT 计算速度的数值。
对于离散序列 x(n)，我们可以用以下公式计算其DFT和IDFT。

image.png

说得迂腐一点，DFT 需要 NxN 次复数乘法和 N(N-1) 次复数加法。如果没有Twiddle Factor，DFT 的计算复杂度为O(n2)O(n^2)。如果使用捻系数，计算复杂度则为N2log⁡2N\frac { N }{ 2 } \log _{ 2 }{ N }。

Twiddle Factor的数学表达为：

WN=e−j2πNW_N = e^{-j \dfrac{2\pi}{N}}

使用Twiddle Factor重写 DFT 和 IDFT 的计算公式，我们就得到了：

image.png

Why do we use twiddle factors?

我们使用Twiddle Factor来降低计算 DFT 和 IDFT 的计算复杂度。

twiddle factor是一个旋转矢量。这意味着，在给定的 N 点 DFT 或 IDFT 计算中，可以观察到twiddle factor的值每 N 个周期重复一次。在每(N-1)步预期有一组熟悉的值，使得计算变得稍微容易一些。(N-1 是因为第一个序列是 0）

另外，我们也可以说 "twiddle factor"具有周期性/循环性。下面我们将用图形来说明这一点。

Cyclic property of twiddle factors

推导一个4点DFT的twiddle factor值

对于n = 0 和 k = 0,

WNnk=W40=e0=1W_N^{nk} = W_4^0 = e^0 = 1

对其余数值进行类似计算，得出以下值：

W40=1W41=−jW42=−1W43=jW44=1W_4^0 = 1\\ W_4^1 = -j\\ W_4^2 = -1\\ W_4^3 = j\\ W_4^4 = 1

可以看到，数值在第 4 个时刻开始重复。这种周期特性可以如下图所示。

image.png

推导一个8点DFT的twiddle factor值

W80=1W81=0.707−0.707jW82=−jW83=−0.707−0.707jW84=−1W85=−0.707+0.707jW86=jW87=0.707+0.707jW88=1W_8^0 = 1\\ W_8^1 = 0.707-0.707j\\ W_8^2 = -j\\ W_8^3 = -0.707-0.707j\\ W_8^4 = -1\\ W_8^5 = -0.707+0.707j\\ W_8^6 = j\\ W_8^7 = 0.707+0.707j\\ W_8^8 = 1\\

image.png

Note:
WNN=1WNN2=−1W_N^N = 1\\ W_N^{\frac{N}{2}} = -1

Matrix method of calculating DFT and IDFT with twiddle factors

image.png

同样，IDFT 也可以用矩阵形式计算，公式如下

X(n)=WN∗NX(n)=\frac{W_N^*}{N}

WN∗{ W }_{ N }^{ * } 是twiddle factors的复共轭。要得到复共轭的值，只需求反twiddle factors复分量的符号即可。例如 0.707+0.707j 的复共轭将变为 0.707-0.707j。

DFT as linear transform

WN{ W }_{ N }的矩阵称为线性变换矩阵。请看下面 DFT 中线性变换的场景。

使用矩阵计算 DFT 的公式为

X(k)=WNx(n)X(n)=WN−1x(k)X(k)={W_N}x(n)\\ X(n)={W_N^{-1}}x(k)\\

使用矩阵计算 IDFT 的公式为

X(n)=WN∗xNX(n)=\frac{W_N^{*}x}{N}\\
WN−1=WN∗NW_N^{-1} = \frac{W_N^{*}}{N}\\
WN∗WN=NINW_N^{*}W_N = NI_N\\

Mixed-Radix Cooley-Tukey FFT

REFERENCE

  1. 快速傅里叶变换（FFT）之一：Radix-2 DIT FFT

  2. MIXED-RADIX COOLEY-TUKEY FFT

  3. Fast Fourier Transform (FFT) (ntu.edu.tw)

Cooley-Tukey FFT 的两个基本类型是时间抽取（DIT）和频域抽取（DIF）。

Decimation in time

回顾上文，DFT定义为：

X(k)=∑n=0N−1x(n)Wnkn,0≤k≤N−1X(k) = \sum_{n=0}^{N-1}x(n)W_n^{kn}, \quad 0\leq k\leq N-1

twiddle factor定义为：

WN≜e−j2πNW_N \triangleq e^{-j \dfrac{2\pi}{N}}

以及性质

WN2=WN2W_N^2 = W_{\frac{N}{2}}
WNk+N2=−WNkW_N^{k+\frac{N}{2}} = -W_N^k

当 N 为偶数时(即N=2MN=2^M, M为整数)，DFT 求和可分为

  • 输入信号偶数索引 x1(n)=x(2n),n=0,1,⋯ ,N2−1x_1(n)=x(2n), \quad n=0, 1, \cdots, \frac{N}{2}-1
  • 输入信号奇数索引 x2(n)=x(2n+1),n=0,1,⋯ ,N2−1x_2(n)=x(2n+1), \quad n=0, 1, \cdots, \frac{N}{2}-1

的求和：

X(k)=∑n=0N−1x(n)WNnk,k=0,1,⋯ ,N−1=∑n=0N2−1x(2n)WN2nk+∑n=0N2−1x(2n+1)WN(2n+1)k=∑n=0N2−1x1(n)WN2nk+WNk∑n=0N2−1x2(n)WN2nk=X1(k)+WNk⋅X2(k)≜DFTN2,k(Downsample2(x))+WNk⋅DFTN2,k(Downsample2[Shift1(x)])\begin{align} X(k) &= \sum_{n=0}^{N-1}x(n)W_N^{nk}, \qquad k = 0, 1, \cdots, N-1\\ &= \sum_{n=0}^{\frac{N}{2}-1}x(2n)W_N^{2nk} + \sum_{n=0}^{\frac{N}{2}-1}x(2n+1)W_N^{(2n+1)k}\\ &= \sum_{n=0}^{\frac{N}{2}-1}x_1(n)W_{\frac{N}{2}}^{nk} + W_N^k\sum_{n=0}^{\frac{N}{2}-1}x_2(n)W_{\frac{N}{2}}^{nk}\\ &= X_1(k)+ W_N^k\cdot X_2(k)\\ &\triangleq DFT_{\frac{N}{2}, k}(Downsample_2(x))+W_N^k\cdot DFT_{\frac{N}{2}, k}(Downsample_2[Shift_1(x)]) \end{align}

如果 N 不满足这个条件，可以通过补零，使之达到这个要求。

Again, x1(n)≜x(2n)x_1(n) \triangleq x(2n) 和 x2(n)≜x(2n+1)x_2(n) \triangleq x(2n+1)表示来自 xx 的偶数和奇数索引的采样点，长度为 N/2N/2。
X1(k)X_1(k)和X2(k)X_2(k)分别为x1(n)x_1(n)和x2(n)x_2(n)的 N/2N/2 点DFT，且具有周期性，周期为 N/2N/2。 因此我们有 X1(k+N2)=X1(k)X_1(k+\frac{N}{2}) = X_1(k) 和 X2(k+N2)=X2(k)X_2(k+\frac{N}{2}) = X_2(k) ，旋转因子 WNk+N2=−WNkW_N^{k+\frac{N}{2}} = -W_N^k ，从而上式可以表示成

X(k)=X1(k)+WNk⋅X2(k),k=0,1,⋯ ,N2−1X(k+N2)=X1(k)−WNk⋅X2(k),k=0,1,⋯ ,N2−1X(k) = X_1(k) + W_N^k\cdot X_2(k), \qquad k = 0, 1, \cdots, \frac{N}{2}-1\\ X(k+\frac{N}{2}) = X_1(k) - W_N^k\cdot X_2(k), \qquad k = 0, 1, \cdots, \frac{N}{2}-1\\

因此，长度为 NN 的 DFT 可用两个长度为 N/2N/2 的 DFT 来计算。

可以用下面的示意图来表示上述过程：

image.png

图中的蝶形结构称为Radix-2 DIT蝶形结构，可以用下图表示：

image.png

Radix 2 FFT

从上图可以看出每一次蝶形运算需要1次复数乘法，两次复数加法。通过将x(n)x(n)在时间上抽取后，N点的DFT变换，现在变成了两个N/2点的DFT，再加上右边的N/2个蝶形运算，因此计算的复杂度为：O(N2/4+N2/4+N/2)O(N^2/4+N^2/4+N/2)，即约为 O(N2/2)O(N^2/2)。计算量相比之前已经减半了。

我们可以按照这种分而治之(divide and conquer)的思路继续将 x1(n)x_1(n) 和 x2(n)x_2(n) 分别抽取为两个长度为N/4N/4的子序列，然后再进行蝶形运算。我们将 x1(n)x_1(n) 的抽取过程画出来：

image.png

我们可以继续用上述的方法将每个子序列继续划分，直到最后序列长度为2。2点的DFT计算如下：

image.png

image.png

总而言之，当 NN 是 22 的幂时，例如 N=2MN=2^M ，其中 M>1M>1 是整数，那么上述 DIT 分解可以执行 M−1M-1 次，直到每个 DFT 长度为 22。长度为 22 的 DFT 不需要乘法。结果称为radix 2 FFT。如果是对频率抽取，可以得到radix 2 DIF FFT。

通过这种思想，以8点DFT计算为例。8点DFT计算以树形阶段进行，首先计算四个2点 DFT，然后是两个4点 DFT，最后是一个8点 DFT。

image.png

或者

企业微信截图_17068547602284.png

需要关注的是，输入数据序列经过 (M-1) 次抽取后的顺序。以8点DFT为例，我们知道第一次抽取后生成序列x(0) x(4) x(2) x(6) x(1) x(3) x(5) x(7)x(0)\ x(4)\ x(2)\ x(6)\ x(1)\ x(3)\ x(5)\ x(7)，第二次抽取后生成序列x(0) x(4) x(2) x(6) x(1) x(5) x(3) x(7)x(0)\ x(4)\ x(2)\ x(6)\ x(1)\ x(5)\ x(3)\ x(7)。通过观察下图可以确定，输入数据序列的抽取有明确的顺序。

image.png

总结就是，Radix-2 DIT FFT算法，信号时域的输入为倒位序，频域的输出为自然顺序。 倒位序指的是将输入序列的二进制位反转。

从理论上讲，split radix FFT 比纯radix 2 算法更高效，因为它能最大限度地减少实际算术运算。然而，在现代通用处理器上，计算时间往往不能通过最小化算术运算次数来实现。术语 "split radix"是指一种 DIT 分解，它将一个radix 2 和两个radix 4 的 FFT 部分结合起来。

In-Place Computation

在Radix-2 DIT FFT算法中，我们可以发现，每个阶段的蝶形运算都是在前一阶段的序列上进行的，因此蝶形运算的输出能够被存储在之前用于蝶形运算输入数据的同样的硬件存储器内，不需要中间存储器。这种运算方式称为同址运算(in-place computation)。同址运算可以节省内存空间。

image.png

Fixed-point and NFFTS

回想一下，IDFT 需要除以 NN，而DFT 则不需要。在定点运算中，当 NN 是 2 的幂时，除以 NN 可以通过在二进制字中右移 log⁡2(N)\log_2(N) 位来实现。IFFT的定点实现通常是通过在每个Butterfly stage后右移。然而，通过在每隔一个Butterfly stage后进行右移，可以获得更好的整体数值性能，这相当于将FFT和IFFT 都除以 N\sqrt{N}（除以N\sqrt{N}等同于将右移次数除以 NN 后的一半）。因此，在定点中，数值性能可以得到改善，不需要额外的工作，归一化操作（右移）在正向变换和反向变换中平均分配，而不是将所有 NN个右移都集中在反变换中。因此，NDFT 对定点实现非常有吸引力。

由于从一个Butterfly stage到下一个Butterfly stage，信号幅度（amplotude）可能增长 2 倍，因此后续每一对 NDFT butterfly stage都需要一个额外的保护位（guard bit）。

还要注意的是，如果 DFT 长度 N=2KN=2^K 不是 44 的幂次方，则正向变换和反向变换中的右移次数必须相差一个（因为 K=log⁡2(N)K=\log_2(N) 是奇数而不是偶数）。

avatar
文章
阅读
粉丝