傅里叶变换 | FFT 与 DFT 原理及算法

注：本文为 “傅里叶变换 | FFT 与 DFT” 相关合辑。
英文引文，机翻未校。
中文引文，略作重排。
图片清晰度受引文原图所限。
如有内容异常，请看原文。

Fast Fourier Transform (FFT)

快速傅里叶变换（FFT）

In this section we present several methods for computing the DFT efficiently. In view of the importance of the DFT in various digital signal processing applications, such as linear filtering, correlation analysis, and spectrum analysis, its efficient computation is a topic that has received considerable attention by many mathematicians, engineers, and applied scientists.
本节介绍几种高效计算离散傅里叶变换（DFT）的方法。鉴于 DFT 在各种数字信号处理应用（如线性滤波、相关分析和频谱分析）中的重要性，其高效计算问题已受到众多数学家、工程师和应用科学家的广泛关注。

From this point, we change the notation that X ( k ) X (k) X(k), instead of y ( k ) y (k) y(k) in previous sections, represents the Fourier coefficients of x ( n ) x (n) x(n).
从此处开始，我们改变符号表示： X ( k ) X (k) X(k)（而非前几节中的 y ( k ) y (k) y(k)）代表 x ( n ) x (n) x(n) 的傅里叶系数。

Basically, the computational problem for the DFT is to compute the sequence { X ( k ) X (k) X(k)} of N N N complex-valued numbers given another sequence of data { x ( n ) x (n) x(n)} of length N N N, according to the formula
基本上，DFT 的计算问题是：根据给定长度为 N N N 的数据序列 { x ( n ) x (n) x(n)}，按照以下公式计算 N N N 个复数值序列 { X ( k ) X (k) X(k)}

X ( k ) = ∑ n = 0 N − 1 x ( n ) W N k n , 0 ≤ k ≤ N − 1 X(k) = \sum_{n=0}^{N-1} x(n) W_N^{kn}, \quad 0 \leq k \leq N-1 X(k)=n=0∑N−1​x(n)WNkn​,0≤k≤N−1

where
其中
W N = e − j 2 π / N W_N = e^{-j2\pi/N} WN​=e−j2π/N.

In general, the data sequence x ( n ) x (n) x(n) is also assumed to be complex valued. Similarly, The IDFT becomes
通常，数据序列 x ( n ) x (n) x(n) 也被假定为复数值。类似地，逆离散傅里叶变换（IDFT）为

x ( n ) = 1 N ∑ k = 0 N − 1 X ( k ) W N − n k , 0 ≤ n ≤ N − 1 x(n) = \frac{1}{N} \sum_{k=0}^{N-1} X(k) W_N^{-nk}, \quad 0 \leq n \leq N-1 x(n)=N1​k=0∑N−1​X(k)WN−nk​,0≤n≤N−1

Since DFT and IDFT involve basically the same type of computations, our discussion of efficient computational algorithms for the DFT applies as well to the efficient computation of the IDFT.
由于 DFT 和 IDFT 涉及基本相同的计算类型，我们对 DFT 高效计算算法的讨论同样适用于 IDFT 的高效计算。

We observe that for each value of k k k, direct computation of X ( k ) X (k) X(k) involves N N N complex multiplications ( 4 N 4N 4N real multiplications) and N − 1 N-1 N−1 complex additions ( 4 N − 2 4N-2 4N−2 real additions). Consequently, to compute all N N N values of the DFT requires N 2 N^2 N2 complex multiplications and N 2 − N N^2-N N2−N complex additions.
我们观察到，对于每个 k k k 值，直接计算 X ( k ) X (k) X(k) 涉及 N N N 次复数乘法（ 4 N 4N 4N 次实数乘法）和 N − 1 N-1 N−1 次复数加法（ 4 N − 2 4N-2 4N−2 次实数加法）。因此，计算 DFT 的全部 N N N 个值需要 N 2 N^2 N2 次复数乘法和 N 2 − N N^2-N N2−N 次复数加法。

Direct computation of the DFT is basically inefficient primarily because it does not exploit the symmetry and periodicity properties of the phase factor W N W_N WN​. In particular, these two properties are :
直接计算 DFT 基本上效率低下，主要是因为它没有利用相位因子 W N W_N WN​ 的对称性和周期性。具体而言，这两个性质为：

Symmetry:  W N k + N / 2 = − W N k Periodicity:  W N k + N = W N k \text {Symmetry: } W_N^{k+N/2} = -W_N^k\\[1em] \text {Periodicity: } W_N^{k+N} = W_N^k Symmetry: WNk+N/2​=−WNk​Periodicity: WNk+N​=WNk​

The computationally efficient algorithms described in this sectio, known collectively as fast Fourier transform (FFT) algorithms, exploit these two basic properties of the phase factor.
本节描述的计算高效算法统称为快速傅里叶变换（FFT）算法，它们利用了相位因子的这两个基本性质。

Radix-2 FFT Algorithms

基 - 2 FFT 算法

Let us consider the computation of the N = 2 v N = 2^v N=2v point DFT by the divide-and conquer approach. We split the N N N-point data sequence into two N / 2 N/2 N/2-point data sequences f 1 ( n ) f_1 (n) f1​(n) and f 2 ( n ) f_2 (n) f2​(n), corresponding to the even-numbered and odd-numbered samples of x ( n ) x (n) x(n), respectively, that is,
让我们考虑用分治法计算 N = 2 v N = 2^v N=2v 点 DFT。我们将 N N N 点数据序列分成两个 N / 2 N/2 N/2 点数据序列 f 1 ( n ) f_1 (n) f1​(n) 和 f 2 ( n ) f_2 (n) f2​(n)，分别对应 x ( n ) x (n) x(n) 的偶数编号和奇数编号样本，即

f 1 ( n ) = x ( 2 n ) f 2 ( n ) = x ( 2 n + 1 ) , n = 0 , 1 , … , N 2 − 1 \begin{aligned} f_1 (n) &= x (2n)\\ f_2 (n) &= x (2n+1), \quad n = 0, 1, \ldots, \frac {N}{2}-1 \end{aligned} f1​(n)f2​(n)​=x(2n)=x(2n+1),n=0,1,…,2N​−1​

Thus f 1 ( n ) f_1 (n) f1​(n) and f 2 ( n ) f_2 (n) f2​(n) are obtained by decimating x ( n ) x (n) x(n) by a factor of 2, and hence the resulting FFT algorithm is called adecimation-in-time algorithm.
因此， f 1 ( n ) f_1 (n) f1​(n) 和 f 2 ( n ) f_2 (n) f2​(n) 是通过以因子 2 对 x ( n ) x (n) x(n) 进行抽取（decimation）得到的，因此所得到的 FFT 算法称为按时间抽取算法（decimation-in-time algorithm）。

Now the N-point DFT can be expressed in terms of the DFT’s of the decimated sequences as follows:
现在，N 点 DFT 可以用抽取序列的 DFT 表示如下：
X ( k ) = ∑ n = 0 N − 1 x ( n ) W N k n , k = 0 , 1 , … , N − 1 = ∑ n n  even x ( n ) W N k n + ∑ n n  odd x ( n ) W N k n = ∑ m = 0 ( N / 2 ) − 1 x ( 2 m ) W N 2 m k + ∑ m = 0 ( N / 2 ) − 1 x ( 2 m + 1 ) W N k ( 2 m + 1 ) \begin{aligned} X(k) &= \sum_{n=0}^{N-1} x(n) W_N^{kn}, \quad k = 0,1,\dots,N-1 \\ &= \sum_{\substack{n \\ n \text{ even}}} x(n) W_N^{kn} + \sum_{\substack{n \\ n \text{ odd}}} x(n) W_N^{kn} \\ &= \sum_{m=0}^{(N/2)-1} x(2m) W_N^{2mk} + \sum_{m=0}^{(N/2)-1} x(2m+1) W_N^{k(2m+1)} \end{aligned} X(k)​=n=0∑N−1​x(n)WNkn​,k=0,1,…,N−1=nn even​∑​x(n)WNkn​+nn odd​∑​x(n)WNkn​=m=0∑(N/2)−1​x(2m)WN2mk​+m=0∑(N/2)−1​x(2m+1)WNk(2m+1)​​

But W N 2 = W N / 2 W_N^2 = W_{N/2} WN2​=WN/2​. With this substitution, the equation can be expressed as
由于 W N 2 = W N / 2 W_N^2 = W_{N/2} WN2​=WN/2​，代入后该式可表示为
X ( k ) = ∑ m = 0 ( N / 2 ) − 1 f 1 ( m ) W N / 2 k m + W N k ∑ m = 0 ( N / 2 ) − 1 f 2 ( m ) W N / 2 k m = F 1 ( k ) + W N k F 2 ( k ) , k = 0 , 1 , … , N − 1 \begin{aligned} X(k) &= \sum_{m=0}^{(N/2)-1} f_1(m) W_{N/2}^{km} + W_N^k \sum_{m=0}^{(N/2)-1} f_2(m) W_{N/2}^{km} \\ &= F_1(k) + W_N^k F_2(k), \quad k = 0,1,\dots,N-1 \end{aligned} X(k)​=m=0∑(N/2)−1​f1​(m)WN/2km​+WNk​m=0∑(N/2)−1​f2​(m)WN/2km​=F1​(k)+WNk​F2​(k),k=0,1,…,N−1​

where F 1 ( k ) F_1 (k) F1​(k) and F 2 ( k ) F_2 (k) F2​(k) are the N / 2 N/2 N/2-point DFTs of the sequences f 1 ( m ) f_1 (m) f1​(m) and f 2 ( m ) f_2 (m) f2​(m), respectively.
其中 F 1 ( k ) F_1 (k) F1​(k) 和 F 2 ( k ) F_2 (k) F2​(k) 分别是序列 f 1 ( m ) f_1 (m) f1​(m) 和 f 2 ( m ) f_2 (m) f2​(m) 的 N / 2 N/2 N/2 点 DFT。

Since F 1 ( k ) F_1 (k) F1​(k) and F 2 ( k ) F_2 (k) F2​(k) are periodic, with period N / 2 N/2 N/2, we have F 1 ( k + N / 2 ) = F 1 ( k ) F_1 (k+N/2) = F_1 (k) F1​(k+N/2)=F1​(k) and F 2 ( k + N / 2 ) = F 2 ( k ) F_2 (k+N/2) = F_2 (k) F2​(k+N/2)=F2​(k). In addition, the factor W N k + N / 2 = − W N k W_N^{k+N/2} = -W_N^k WNk+N/2​=−WNk​. Hence the equation may be expressed as
由于 F 1 ( k ) F_1 (k) F1​(k) 和 F 2 ( k ) F_2 (k) F2​(k) 具有周期性，周期为 N / 2 N/2 N/2，因此有 F 1 ( k + N / 2 ) = F 1 ( k ) F_1 (k+N/2) = F_1 (k) F1​(k+N/2)=F1​(k) 和 F 2 ( k + N / 2 ) = F 2 ( k ) F_2 (k+N/2) = F_2 (k) F2​(k+N/2)=F2​(k)。此外，因子 W N k + N / 2 = − W N k W_N^{k+N/2} = -W_N^k WNk+N/2​=−WNk​。因此该方程可表示为

X ( k ) = F 1 ( k ) + W N k F 2 ( k ) , k = 0 , 1 , … , N 2 − 1 X ( k + N 2 ) = F 1 ( k ) − W N k F 2 ( k ) , k = 0 , 1 , … , N 2 − 1 \begin{aligned} X (k) &= F_1 (k) + W_N^k F_2 (k), \quad k = 0, 1, \ldots, \frac {N}{2}-1\\ X\left (k+\frac {N}{2}\right) &= F_1 (k) - W_N^k F_2 (k), \quad k = 0, 1, \ldots, \frac {N}{2}-1 \end{aligned} X(k)X(k+2N​)​=F1​(k)+WNk​F2​(k),k=0,1,…,2N​−1=F1​(k)−WNk​F2​(k),k=0,1,…,2N​−1​

We observe that the direct computation of F 1 ( k ) F_1 (k) F1​(k) requires ( N / 2 ) 2 (N/2)^2 (N/2)2 complex multiplications. The same applies to the computation of F 2 ( k ) F_2 (k) F2​(k). Furthermore, there are N / 2 N/2 N/2 additional complex multiplications required to compute W N k F 2 ( k ) W_N^k F_2 (k) WNk​F2​(k). Hence the computation of X ( k ) X (k) X(k) requires 2 ( N / 2 ) 2 + N / 2 = N 2 / 2 + N / 2 2 (N/2)^2 + N/2 = N^2/2 + N/2 2(N/2)2+N/2=N2/2+N/2 complex multiplications. This first step results in a reduction of the number of multiplications from N 2 N^2 N2 to N 2 / 2 + N / 2 N^2/2 + N/2 N2/2+N/2, which is about a factor of 2 for N N N large.
我们观察到，直接计算 F 1 ( k ) F_1 (k) F1​(k) 需要 ( N / 2 ) 2 (N/2)^2 (N/2)2 次复数乘法。 F 2 ( k ) F_2 (k) F2​(k) 的计算同样需要这么多。此外，计算 W N k F 2 ( k ) W_N^k F_2 (k) WNk​F2​(k) 还需要额外的 N / 2 N/2 N/2 次复数乘法。因此，计算 X ( k ) X (k) X(k) 需要 2 ( N / 2 ) 2 + N / 2 = N 2 / 2 + N / 2 2 (N/2)^2 + N/2 = N^2/2 + N/2 2(N/2)2+N/2=N2/2+N/2 次复数乘法。这第一步将乘法次数从 N 2 N^2 N2 减少到 N 2 / 2 + N / 2 N^2/2 + N/2 N2/2+N/2，对于较大的 N N N，约为因子 2 的缩减。

Figure TC.3.1 First step in the decimation-in-time algorithm.
图 TC.3.1 按时间抽取算法的第一步。

By computing N / 4 N/4 N/4-point DFTs, we would obtain the N / 2 N/2 N/2-point DFTs F 1 ( k ) F_1 (k) F1​(k) and F 2 ( k ) F_2 (k) F2​(k) from the relations
通过计算 N / 4 N/4 N/4 点 DFT，我们可以从以下关系式得到 N / 2 N/2 N/2 点 DFT F 1 ( k ) F_1 (k) F1​(k) 和 F 2 ( k ) F_2 (k) F2​(k)

F 1 ( k ) = F { f 1 ( 2 n ) } + W N / 2 k F { f 1 ( 2 n + 1 ) } , F 1 ⁣ ( k + N 4 ) = F { f 1 ( 2 n ) } − W N / 2 k F { f 1 ( 2 n + 1 ) } , F 2 ( k ) = F { f 2 ( 2 n ) } + W N / 2 k F { f 2 ( 2 n + 1 ) } , F 2 ⁣ ( k + N 4 ) = F { f 2 ( 2 n ) } − W N / 2 k F { f 2 ( 2 n + 1 ) } ( k = 0 , 1 , … , N 4 − 1 ;   n = 0 , 1 , … , N 4 − 1 ) \begin{align*} F_1(k) &= \mathrm{F}\{f_1(2n)\} + W_{N/2}^k \mathrm{F}\{f_1(2n+1)\}, \\ F_1\!\left(k+\frac{N}{4}\right) &= \mathrm{F}\{f_1(2n)\} - W_{N/2}^k \mathrm{F}\{f_1(2n+1)\}, \\ F_2(k) &= \mathrm{F}\{f_2(2n)\} + W_{N/2}^k \mathrm{F}\{f_2(2n+1)\}, \\ F_2\!\left(k+\frac{N}{4}\right) &= \mathrm{F}\{f_2(2n)\} - W_{N/2}^k \mathrm{F}\{f_2(2n+1)\} \end{align*} \quad (k = 0,1,\dots,\frac{N}{4}-1;\ n = 0,1,\dots,\frac{N}{4}-1) F1​(k)F1​(k+4N​)F2​(k)F2​(k+4N​)​=F{f1​(2n)}+WN/2k​F{f1​(2n+1)},=F{f1​(2n)}−WN/2k​F{f1​(2n+1)},=F{f2​(2n)}+WN/2k​F{f2​(2n+1)},=F{f2​(2n)}−WN/2k​F{f2​(2n+1)}​(k=0,1,…,4N​−1; n=0,1,…,4N​−1)

The decimation of the data sequence can be repeated again and again until the resulting sequences are reduced to one-point sequences. For N = 2 v N = 2^v N=2v, this decimation can be performed v = log ⁡ 2 N v = \log_2 N v=log2​N times. Thus the total number of complex multiplications is reduced to ( N / 2 ) log ⁡ 2 N (N/2)\log_2 N (N/2)log2​N. The number of complex additions is N log ⁡ 2 N N\log_2 N Nlog2​N.
数据序列的抽取可以反复进行，直到结果序列缩减为单点序列。对于 N = 2 v N = 2^v N=2v，这种抽取可以进行 v = log ⁡ 2 N v = \log_2 N v=log2​N 次。因此，复数乘法的总次数减少到 ( N / 2 ) log ⁡ 2 N (N/2)\log_2 N (N/2)log2​N。复数加法的次数为 N log ⁡ 2 N N\log_2 N Nlog2​N。

For illustrative purposes, Figure TC.3.2 depicts the computation of N = 8 N = 8 N=8 point DFT. We observe that the computation is performed in tree stages, beginning with the computations of four two-point DFTs, then two four-point DFTs, and finally, one eight-point DFT. The combination for the smaller DFTs to form the larger DFT is illustrated in Figure TC.3.3 for N = 8 N = 8 N=8.
为便于说明，图 TC.3.2 描绘了 N = 8 N = 8 N=8 点 DFT 的计算。我们观察到计算分三个阶段进行：首先计算四个两点 DFT，然后计算两个四点 DFT，最后计算一个八点 DFT。图 TC.3.3 说明了 N = 8 N = 8 N=8 时较小 DFT 组合形成较大 DFT 的过程。

Figure TC.3.2 Three stages in the computation of an N = 8 N = 8 N=8-point DFT.
图 TC.3.2 计算 N = 8 N = 8 N=8 点 DFT 的三个阶段。

Figure TC.3.3Eight-point decimation-in-time FFT algorithm.
图 TC.3.3八点按时间抽取 FFT 算法。

Figure TC.3.4 Basic butterfly computation in the decimation-in-time FFT algorithm.
图 TC.3.4 按时间抽取 FFT 算法中的基本蝶形运算。

An important observation is concerned with the order of the input data sequence after it is decimated ( v − 1 ) (v-1) (v−1) times. For example, if we consider the case where N = 8, we know that the first decimation yeilds the sequence x ( 0 ) , x ( 2 ) , x ( 4 ) , x ( 6 ) , x ( 1 ) , x ( 3 ) , x ( 5 ) , x ( 7 ) x (0), x (2), x (4), x (6), x (1), x (3), x (5), x (7) x(0),x(2),x(4),x(6),x(1),x(3),x(5),x(7), and the second decimation results in the sequence 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) x(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7). Thisshufflingof the input data sequence has a well-defined order as can be ascertained from observing Figure TC.3.5, which illustrates the decimation of the eight-point sequence.
一个重要的观察点是输入数据序列经过 ( v − 1 ) (v-1) (v−1) 次抽取后的顺序。例如，若考虑 N = 8 N = 8 N=8 的情况，我们知道第一次抽取产生序列 x ( 0 ) , x ( 2 ) , x ( 4 ) , x ( 6 ) , x ( 1 ) , x ( 3 ) , x ( 5 ) , x ( 7 ) x (0), x (2), x (4), x (6), x (1), x (3), x (5), x (7) x(0),x(2),x(4),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) x(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7)。这种输入数据序列的混洗（shuffling）具有明确的顺序，可从图 TC.3.5 中看出，该图说明了八点序列的抽取过程。

Figure TC.3.5 Shuffling of the data and bit reversal.
图 TC.3.5 数据的混洗与位反转。

Another important radix-2 FFT algorithm, called the decimation-in-frequency algorithm, is obtained by using the divide-and-conquer approach. To derive the algorithm, we begin by splitting the DFT formula into two summations, one of which involves the sum over the first N / 2 N/2 N/2 data points and the second sum involves the last N/2 data points. Thus we obtain
另一种重要的基 - 2 FFT 算法称为按频率抽取算法（decimation-in-frequency algorithm），通过分治法得到。为推导该算法，我们首先将 DFT 公式分成两个求和式，其中一个涉及前 N / 2 N/2 N/2 个数据点的求和，第二个涉及后 N / 2 N/2 N/2 个数据点的求和。因此我们得到

X ( k ) = ∑ n = 0 N / 2 − 1 x ( n ) W N n k + ∑ n = N / 2 N − 1 x ( n ) W N n k = ∑ n = 0 N / 2 − 1 [ x ( n ) + x ⁣ ( n + N 2 ) W N N k / 2 ] W N n k \begin{align*} X(k) &= \sum_{n=0}^{N/2-1} x(n) W_N^{nk} + \sum_{n=N/2}^{N-1} x(n) W_N^{nk} \\ &= \sum_{n=0}^{N/2-1} \left[ x(n) + x\!\left(n+\frac{N}{2}\right) W_N^{Nk/2} \right] W_N^{nk} \end{align*} X(k)​=n=0∑N/2−1​x(n)WNnk​+n=N/2∑N−1​x(n)WNnk​=n=0∑N/2−1​[x(n)+x(n+2N​)WNNk/2​]WNnk​​

Now, let us split (decimate) X (k) into the even- and odd-numbered samples. Thus we obtain
现在，让我们将 X ( k ) X (k) X(k) 分成偶数编号和奇数编号的样本。因此我们得到
X ( 2 k ) = ∑ n = 0 N / 2 − 1 [ x ( n ) + x ⁣ ( n + N 2 ) ] W N / 2 n k , k = 0 , 1 , … , N 2 − 1 , X ( 2 k + 1 ) = ∑ n = 0 N / 2 − 1 [ x ( n ) − x ⁣ ( n + N 2 ) ] W N n W N / 2 n k , k = 0 , 1 , … , N 2 − 1. \begin{align*} X(2k) &= \sum_{n=0}^{N/2-1} \left[ x(n) + x\!\left(n+\frac{N}{2}\right) \right] W_{N/2}^{\,nk}, && k = 0, 1, \ldots, \frac{N}{2}-1, \\ X(2k+1) &= \sum_{n=0}^{N/2-1} \left[ x(n) - x\!\left(n+\frac{N}{2}\right) \right] W_N^{\,n} W_{N/2}^{\,nk}, && k = 0, 1, \ldots, \frac{N}{2}-1. \end{align*} X(2k)X(2k+1)​=n=0∑N/2−1​[x(n)+x(n+2N​)]WN/2nk​,=n=0∑N/2−1​[x(n)−x(n+2N​)]WNn​WN/2nk​,​​k=0,1,…,2N​−1,k=0,1,…,2N​−1.​

where we have used the fact that W N 2 = W N / 2 W_N^2 = W_{N/2} WN2​=WN/2​
其中我们利用了 W N 2 = W N / 2 W_N^2 = W_{N/2} WN2​=WN/2​ 这一事实

The computational procedure above can be repeated through decimation of the N / 2 N/2 N/2-point DFTs X ( 2 k ) X (2k) X(2k) and X ( 2 k + 1 ) X (2k+1) X(2k+1). The entire process involves v = log ⁡ 2 N v = \log_2 N v=log2​N stages of decimation, where each stage involves N / 2 N/2 N/2 butterflies of the type shown in Figure TC.3.7. Consequently, the computation of the N-point DFT via the decimation-in-frequency FFT requires ( N / 2 ) log ⁡ 2 N (N/2)\log_2 N (N/2)log2​N complex multiplications and N log ⁡ 2 N N\log_2 N Nlog2​N complex additions, just as in the decimation-in-time algorithm. For illustrative purposes, the eight-point decimation-in-frequency algorithm is given in Figure TC.3.8.
上述计算过程可以通过对 N / 2 N/2 N/2 点 DFT X ( 2 k ) X (2k) X(2k) 和 X ( 2 k + 1 ) X (2k+1) X(2k+1) 的抽取来重复。整个过程涉及 v = log ⁡ 2 N v = \log_2 N v=log2​N 个抽取阶段，每个阶段涉及 N / 2 N/2 N/2 个如图 TC.3.7 所示类型的蝶形运算。因此，通过按频率抽取 FFT 计算 N 点 DFT 需要 ( N / 2 ) log ⁡ 2 N (N/2)\log_2 N (N/2)log2​N 次复数乘法和 N log ⁡ 2 N N\log_2 N Nlog2​N 次复数加法，与按时间抽取算法相同。为便于说明，八点按频率抽取算法示于图 TC.3.8。

Figure TC.3.6 First stage of the decimation-in-frequency FFT algorithm.
图 TC.3.6 按频率抽取 FFT 算法的第一阶段。

Figure TC.3.7 Basic butterfly computation in the decimation-in-frequency.
图 TC.3.7 按频率抽取中的基本蝶形运算。

Figure TC.3.8 N = 8 N = 8 N=8-point decimation-in-frequency FFT algorithm.
图 TC.3.8 N = 8 N = 8 N=8 点按频率抽取 FFT 算法。

We observe from Figure TC.3.8 that the input data x ( n ) x (n) x(n) occurs in natural order, but the output DFT occurs in bit-reversed order. We also note that the computations are performed in place. However, it is possible to reconfigure the decimation-in-frequency algorithm so that the input sequence occurs in bit-reversed order while the output DFT occurs in normal order. Furthermore, if we abandon the requirement that the computations be done in place, it is also possible to have both the input data and the output DFT in normal order.
从图 TC.3.8 我们观察到，输入数据 x ( n ) x (n) x(n) 按自然顺序出现，但输出 DFT 按位反转顺序出现。我们还注意到计算是原位进行的（in place）。然而，可以重新配置按频率抽取算法，使输入序列按位反转顺序出现，而输出 DFT 按正常顺序出现。此外，如果我们放弃计算必须原位进行的要求，也可以使输入数据和输出 DFT 都按正常顺序出现。

Radix-4 FFT Algorithm

基 - 4 FFT 算法

When the number of data points N N N in the DFT is a power of 4 (i.e., N = 4 v N = 4^v N=4v), we can, of course, always use a radix-2 algorithm for the computation. However, for this case, it is more efficient computationally to employ a radix-4 FFT algorithm.
当 DFT 中的数据点数 N N N 是 4 的幂（即 N = 4 v N = 4^v N=4v）时，我们当然总是可以使用基 - 2 算法进行计算。然而，在这种情况下，采用基 - 4 FFT 算法在计算上更为高效。

Let us begin by describing a radix-4 decimation-in-time FFT algorithm briefly. We split or decimate the N N N-point input sequence into four subsequences, x ( 4 n ) , x ( 4 n + 1 ) , x ( 4 n + 2 ) , x ( 4 n + 3 ) , n = 0 , 1 , … , N / 4 − 1 x (4n), x (4n+1), x (4n+2), x (4n+3), n = 0, 1, \ldots , N/4-1 x(4n),x(4n+1),x(4n+2),x(4n+3),n=0,1,…,N/4−1.
让我们首先简要描述基 - 4 按时间抽取 FFT 算法。我们将 N N N 点输入序列分成四个子序列： x ( 4 n ) , x ( 4 n + 1 ) , x ( 4 n + 2 ) , x ( 4 n + 3 ) x (4n), x (4n+1), x (4n+2), x (4n+3) x(4n),x(4n+1),x(4n+2),x(4n+3)，其中 n = 0 , 1 , … , N / 4 − 1 n = 0, 1, \ldots , N/4-1 n=0,1,…,N/4−1。
X ( k ) = ∑ n = 0 N / 4 − 1 x ( 4 n ) W N 4 n k + ∑ n = 0 N / 4 − 1 x ( 4 n + 1 ) W N ( 4 n + 1 ) k + ∑ n = 0 N / 4 − 1 x ( 4 n + 2 ) W N ( 4 n + 2 ) k + ∑ n = 0 N / 4 − 1 x ( 4 n + 3 ) W N ( 4 n + 3 ) k = ∑ l = 0 3 W N l k ∑ n = 0 N / 4 − 1 x ( 4 n + l ) W N / 4 n k . \begin{align*} X(k) &= \sum_{n=0}^{N/4-1} x(4n) W_N^{\,4nk} + \sum_{n=0}^{N/4-1} x(4n+1) W_N^{\,(4n+1)k} \\ &\quad + \sum_{n=0}^{N/4-1} x(4n+2) W_N^{\,(4n+2)k} + \sum_{n=0}^{N/4-1} x(4n+3) W_N^{\,(4n+3)k} \\ &= \sum_{l=0}^{3} W_N^{\,lk} \sum_{n=0}^{N/4-1} x(4n+l) W_{N/4}^{\,nk}. \end{align*} X(k)​=n=0∑N/4−1​x(4n)WN4nk​+n=0∑N/4−1​x(4n+1)WN(4n+1)k​+n=0∑N/4−1​x(4n+2)WN(4n+2)k​+n=0∑N/4−1​x(4n+3)WN(4n+3)k​=l=0∑3​WNlk​n=0∑N/4−1​x(4n+l)WN/4nk​.​

Thus the four N / 4 N/4 N/4-point DFTs F ( l , q ) F (l, q) F(l,q) obtained from the above equation are combined to yield the N-point DFT. The expression for combining the N / 4 N/4 N/4-point DFTs defines a radix-4 decimation-in-time butterfly, which can be expressed in matrix form as
因此，从上述方程得到的四个 N / 4 N/4 N/4 点 DFT F ( l , q ) F (l, q) F(l,q) 被组合起来产生 N 点 DFT。组合 N / 4 N/4 N/4 点 DFT 的表达式定义了一个基 - 4 按时间抽取蝶形运算，可以用矩阵形式表示为

[ X ( k ) X ( k + N / 4 ) X ( k + N / 2 ) X ( k + 3 N / 4 ) ] = [ 1 1 1 1 1 − j − 1 j 1 − 1 1 − 1 1 j − 1 − j ] [ W N 0 F ( 0 , k ) W N k F ( 1 , k ) W N 2 k F ( 2 , k ) W N 3 k F ( 3 , k ) ] \begin {bmatrix} X (k) \\ X (k+N/4) \\ X (k+N/2) \\ X (k+3N/4) \end {bmatrix} = \begin {bmatrix} 1 & 1 & 1 & 1 \\ 1 & -j & -1 & j \\ 1 & -1 & 1 & -1 \\ 1 & j & -1 & -j \end {bmatrix} \begin {bmatrix} W_N^0 F (0,k) \\ W_N^k F (1,k) \\ W_N^{2k} F (2,k) \\ W_N^{3k} F (3,k) \end {bmatrix} ​X(k)X(k+N/4)X(k+N/2)X(k+3N/4)​​=​1111​1−j−1j​1−11−1​1j−1−j​​​WN0​F(0,k)WNk​F(1,k)WN2k​F(2,k)WN3k​F(3,k)​​

The radix-4 butterfly is depicted in Figure TC.3.9a and in a more compact form in Figure TC.3.9b. Note that each butterfly involves three complex multiplications, since W N 0 = 1 W_N^0 = 1 WN0​=1, and 12 complex additions.
基 - 4 蝶形运算示于图 TC.3.9a，更紧凑的形式示于图 TC.3.9b。注意，每个蝶形运算涉及三次复数乘法（因为 W N 0 = 1 W_N^0 = 1 WN0​=1）和 12 次复数加法。

Figure TC.3.9 Basic butterfly computation in a radix-4 FFT algorithm.
图 TC.3.9 基 - 4 FFT 算法中的基本蝶形运算。

By performing the additions in two steps, it is possible to reduce the number of additions per butterfly from 12 to 8. This can be accomplished ty expressing the matrix of the linear transformation mentioned previously as a product of two matrices as follows:
通过分两步执行加法，可以将每个蝶形运算的加法次数从 12 减少到 8。这可以通过将前述线性变换的矩阵表示为两个矩阵的乘积来实现，如下所示：

[ X ( 0 , q ) X ( 1 , q ) X ( 2 , q ) X ( 3 , q ) ] = [ 1 0 1 0 0 1 0 − j 1 0 − 1 0 0 1 0 j ] [ 1 0 1 0 1 0 − 1 0 0 1 0 1 0 1 0 − 1 ] [ W N 0 F ( 0 , q ) W N 1 F ( 1 , q ) W N 2 F ( 2 , q ) W N 3 F ( 3 , q ) ] \begin{bmatrix} X(0,q) \\ X(1,q) \\ X(2,q) \\ X(3,q) \end{bmatrix} = \begin{bmatrix} 1 & 0 & 1 & 0 \\ 0 & 1 & 0 & -j \\ 1 & 0 & -1 & 0 \\ 0 & 1 & 0 & j \end{bmatrix} \begin{bmatrix} 1 & 0 & 1 & 0 \\ 1 & 0 & -1 & 0 \\ 0 & 1 & 0 & 1 \\ 0 & 1 & 0 & -1 \end{bmatrix} \begin{bmatrix} W_N^0 F(0,q) \\ W_N^1 F(1,q) \\ W_N^2 F(2,q) \\ W_N^3 F(3,q) \end{bmatrix} ​X(0,q)X(1,q)X(2,q)X(3,q)​​=​1010​0101​10−10​0−j0j​​​1100​0011​1−100​001−1​​​WN0​F(0,q)WN1​F(1,q)WN2​F(2,q)WN3​F(3,q)​​

Figure TC.3.10 Sixteen-point radix-4 decimation-in-time algorithm with input in normal order and output in digit-reversed order
图 TC.3.10 输入为正常顺序、输出为数字反转顺序的 16 点基 - 4 按时间抽取算法

A 16-point, radix-4 decimation-in-frequency FFT algorithm is shown in Figure TC.3.11. Its input is in normal order and its output is in digit-reversed order. It has exactly the same computational complexity as the decimation-in-time radex-4 FFT algorithm.
图 TC.3.11 示出了一种 16 点基 - 4 按频率抽取 FFT 算法。其输入为正常顺序，输出为数字反转顺序。其计算复杂度与基 - 4 按时间抽取 FFT 算法完全相同。

Figure TC.3.11 Sixteen-point, radix-4 decimation-in-frequency algorithm with input in normal order and output in digit-reversed order.
图 TC.3.11 输入为正常顺序、输出为数字反转顺序的 16 点基 - 4 按频率抽取算法。

For illustrative purposes, let us re-derive the radix-4 decimation-in-frequency algorithm
为便于说明，让我们重新推导基 4 频率抽取算法

by breaking the N N N-point DFT formula into four smaller DFTs. We have
通过将 N N N 点 DFT 公式分解为四个较小的 DFT。我们有

X ( k ) = ∑ n = 0 N − 1 x ( n ) W N k n = ∑ n = 0 N / 4 − 1 x ( n ) W N k n + ∑ n = N / 4 N / 2 − 1 x ( n ) W N k n + ∑ n = N / 2 3 N / 4 − 1 x ( n ) W N k n + ∑ n = 3 N / 4 N − 1 x ( n ) W N k n = ∑ n = 0 N / 4 − 1 x ( n ) W N k n + W N N k / 4 ∑ n = 0 N / 4 − 1 x ( n + N 4 ) W N k n + W N N k / 2 ∑ n = 0 N / 4 − 1 x ( n + N 2 ) W N k n + W N 3 N k / 4 ∑ n = 0 N / 4 − 1 x ( n + 3 N 4 ) W N k n \begin {aligned} X (k) &= \sum_{n=0}^{N-1} x (n) W_N^{kn} \\ &= \sum_{n=0}^{N/4-1} x (n) W_N^{kn} + \sum_{n=N/4}^{N/2-1} x (n) W_N^{kn} + \sum_{n=N/2}^{3N/4-1} x (n) W_N^{kn} + \sum_{n=3N/4}^{N-1} x (n) W_N^{kn} \\ &= \sum_{n=0}^{N/4-1} x (n) W_N^{kn} + W_N^{Nk/4} \sum_{n=0}^{N/4-1} x\left (n+\frac {N}{4}\right) W_N^{kn} + \\ &\quad W_N^{Nk/2} \sum_{n=0}^{N/4-1} x\left (n+\frac {N}{2}\right) W_N^{kn} + W_N^{3Nk/4} \sum_{n=0}^{N/4-1} x\left (n+\frac {3N}{4}\right) W_N^{kn} \end {aligned} X(k)​=n=0∑N−1​x(n)WNkn​=n=0∑N/4−1​x(n)WNkn​+n=N/4∑N/2−1​x(n)WNkn​+n=N/2∑3N/4−1​x(n)WNkn​+n=3N/4∑N−1​x(n)WNkn​=n=0∑N/4−1​x(n)WNkn​+WNNk/4​n=0∑N/4−1​x(n+4N​)WNkn​+WNNk/2​n=0∑N/4−1​x(n+2N​)WNkn​+WN3Nk/4​n=0∑N/4−1​x(n+43N​)WNkn​​

From the definition of the twiddle factors, we have
根据旋转因子的定义，我们有
W N k N / 4 = ( − j ) k , W N k N / 2 = ( − 1 ) k , W N 3 k N / 4 = ( j ) k W_N^{kN/4} = (-j)^k, \quad W_N^{kN/2} = (-1)^k, \quad W_N^{3kN/4} = (j)^k WNkN/4​=(−j)k,WNkN/2​=(−1)k,WN3kN/4​=(j)k

Thus
因此

X ( k ) = ∑ n = 0 N / 4 − 1 [ x ( n ) + ( − j ) k x ( n + N 4 ) + ( − 1 ) k x ( n + N 2 ) + ( j ) k x ( n + 3 N 4 ) ] W N n k X (k) = \sum_{n=0}^{N/4-1} \left [ x (n) + (-j)^k x\left (n+\frac {N}{4}\right) + (-1)^k x\left (n+\frac {N}{2}\right) + (j)^k x\left (n+\frac {3N}{4}\right) \right] W_N^{nk} X(k)=n=0∑N/4−1​[x(n)+(−j)kx(n+4N​)+(−1)kx(n+2N​)+(j)kx(n+43N​)]WNnk​

The relation is not an N / 4 N/4 N/4-point DFT because the twiddle factor depends on N and not on N / 4 N/4 N/4. To convert it into an N / 4 N/4 N/4-point DFT we subdivede the DFT sequence into four N / 4 N/4 N/4-point subsequences, X ( 4 k ) , X ( 4 k + 1 ) , X ( 4 k + 2 ) X (4k), X (4k+1), X (4k+2) X(4k),X(4k+1),X(4k+2), and X ( 4 k + 3 ) , k = 0 , 1 , … , N / 4 X (4k+3), k = 0, 1, \ldots, N/4 X(4k+3),k=0,1,…,N/4. Thus we obtain the radix-4 decimation-in frequency DFT as
该关系式不是 N / 4 N/4 N/4 点 DFT，因为旋转因子依赖于 N 而非 N / 4 N/4 N/4。为将其转换为 N / 4 N/4 N/4 点 DFT，我们将 DFT 序列分成四个 N / 4 N/4 N/4 点子序列： X ( 4 k ) , X ( 4 k + 1 ) , X ( 4 k + 2 ) X (4k), X (4k+1), X (4k+2) X(4k),X(4k+1),X(4k+2) 和 X ( 4 k + 3 ) X (4k+3) X(4k+3)，其中 k = 0 , 1 , … , N / 4 k = 0, 1, \ldots, N/4 k=0,1,…,N/4。因此我们得到基 - 4 按频率抽取 DFT 为

X ( 4 k ) = ∑ n = 0 N / 4 − 1 [ x ( n ) + x ( n + N 4 ) + x ( n + N 2 ) + x ( n + 3 N 4 ) ] W N / 4 n k X ( 4 k + 1 ) = ∑ n = 0 N / 4 − 1 [ x ( n ) − j x ( n + N 4 ) − x ( n + N 2 ) + j x ( n + 3 N 4 ) ] W N n W N / 4 n k X ( 4 k + 2 ) = ∑ n = 0 N / 4 − 1 [ x ( n ) − x ( n + N 4 ) + x ( n + N 2 ) − x ( n + 3 N 4 ) ] W N 2 n W N / 4 n k X ( 4 k + 3 ) = ∑ n = 0 N / 4 − 1 [ x ( n ) + j x ( n + N 4 ) − x ( n + N 2 ) − j x ( n + 3 N 4 ) ] W N 3 n W N / 4 n k \begin{aligned} X(4k) &= \sum_{n=0}^{N/4-1} \left[ x(n) + x\left(n+\frac{N}{4}\right) + x\left(n+\frac{N}{2}\right) + x\left(n+\frac{3N}{4}\right) \right] W_{N/4}^{nk} \\ X(4k+1) &= \sum_{n=0}^{N/4-1} \left[ x(n) - jx\left(n+\frac{N}{4}\right) - x\left(n+\frac{N}{2}\right) + jx\left(n+\frac{3N}{4}\right) \right] W_N^n W_{N/4}^{nk} \\ X(4k+2) &= \sum_{n=0}^{N/4-1} \left[ x(n) - x\left(n+\frac{N}{4}\right) + x\left(n+\frac{N}{2}\right) - x\left(n+\frac{3N}{4}\right) \right] W_N^{2n} W_{N/4}^{nk} \\ X(4k+3) &= \sum_{n=0}^{N/4-1} \left[ x(n) + jx\left(n+\frac{N}{4}\right) - x\left(n+\frac{N}{2}\right) - jx\left(n+\frac{3N}{4}\right) \right] W_N^{3n} W_{N/4}^{nk} \end{aligned} X(4k)X(4k+1)X(4k+2)X(4k+3)​=n=0∑N/4−1​[x(n)+x(n+4N​)+x(n+2N​)+x(n+43N​)]WN/4nk​=n=0∑N/4−1​[x(n)−jx(n+4N​)−x(n+2N​)+jx(n+43N​)]WNn​WN/4nk​=n=0∑N/4−1​[x(n)−x(n+4N​)+x(n+2N​)−x(n+43N​)]WN2n​WN/4nk​=n=0∑N/4−1​[x(n)+jx(n+4N​)−x(n+2N​)−jx(n+43N​)]WN3n​WN/4nk​​

where we have used the property W N 4 k n = W N / 4 k n W_N^{4kn} = W_{N/4}^{kn} WN4kn​=WN/4kn​. Note that the input to each N / 4 N/4 N/4-point DFT is a linear combination of four signal samples scaled by a twiddle factor. This procedure is repeated v times, where v = log ⁡ 4 N v = \log_4 N v=log4​N.
其中我们利用了性质 W N 4 k n = W N / 4 k n W_N^{4kn} = W_{N/4}^{kn} WN4kn​=WN/4kn​。注意，每个 N / 4 N/4 N/4 点 DFT 的输入是四个信号样本经旋转因子加权后的线性组合。此过程重复 v 次，其中 v = log ⁡ 4 N v = \log_4 N v=log4​N。

Split-Radix FFT Algorithms

分裂基 FFT 算法

An inspection of the radix-2 decimation-in-frequency flowgraph shown in Figure TC.3.8 indicates that the even-numbered pints of the DFT can be computed independently of the odd-numbered points. This suggests teh possibility of using different computational methods for independent parts of the algorithm, with the objective of reducing the number of computations. The split-radix FFT (SRFFT) algorithms exploit this idea by using both a radix-2 and a radix-4 decomposition in the same FFT algorithm.
观察图 TC.3.8 所示的基 - 2 按频率抽取流程图可知，DFT 的偶数编号点可以独立于奇数编号点进行计算。这提示我们可以在算法的独立部分使用不同的计算方法，以减少计算次数。分裂基 FFT（SRFFT）算法利用这一思想，在同一个 FFT 算法中同时使用基 - 2 和基 - 4 分解。

First, we recall that in the radix-2 decimation-in-frequency FFT algorithm, the even-numbered samples of the N N N-point DFT are given as
首先，我们回顾在基 - 2 按频率抽取 FFT 算法中，N 点 DFT 的偶数编号样本由下式给出

X ( 2 k ) = ∑ n = 0 N / 2 − 1 [ x ( n ) + x ( n + N 2 ) ] W N / 2 n k , k = 0 , 1 , … , N 2 − 1 X (2k) = \sum_{n=0}^{N/2-1} \left [x (n) + x\left (n+\frac {N}{2}\right)\right] W_{N/2}^{nk}, \quad k = 0, 1, \ldots, \frac {N}{2}-1 X(2k)=n=0∑N/2−1​[x(n)+x(n+2N​)]WN/2nk​,k=0,1,…,2N​−1

A radix-2 suffices for this computation.
此计算使用基 - 2 已足够。

The odd-numbered samples { X ( 2 k + 1 ) X (2k+1) X(2k+1)} of the DFT require the pre-multiplication of the input sequence with the twiddle factors W N n W_N^n WNn​. For these samples a radix-4 decomposition produces some computational efficiency because the four-point DFT has the largest multiplication-free butterfly. Indeed, it can be shown that using a radix greater than 4 does not result in a significant reduction in computational complexity.

DFT 的奇数编号样本 { X ( 2 k + 1 ) X (2k+1) X(2k+1)} 需要先将输入序列与旋转因子 W N n W_N^n WNn​ 相乘。对于这些样本，基 - 4 分解能产生一定的计算效率，因为四点 DFT 具有最大的无乘法蝶形运算。事实上，可以证明使用大于 4 的基不会显著降低计算复杂度。

If we use a radix-4 decimation-in-frequency FFT algorithm for the odd-numbered samples of the N N N-point DFT, we obtain the following N / 4 N/4 N/4-point DFTs:

如果我们对 N 点 DFT 的奇数编号样本使用基 - 4 按频率抽取 FFT 算法，我们得到以下 N / 4 N/4 N/4 点 DFT：

X ( 4 k + 1 ) = ∑ n = 0 N / 4 − 1 { [ x ( n ) − x ( n + N 2 ) ] − j [ x ( n + N 4 ) − x ( n + 3 N 4 ) ] } W N n W N / 4 n k X ( 4 k + 3 ) = ∑ n = 0 N / 4 − 1 { [ x ( n ) − x ( n + N 2 ) ] + j [ x ( n + N 4 ) − x ( n + 3 N 4 ) ] } W N 3 n W N / 4 n k \begin{aligned} X (4k+1) &= \sum_{n=0}^{N/4-1} \left\{\left [x (n) - x\left (n+\frac {N}{2}\right)\right] - j\left [x\left (n+\frac {N}{4}\right) - x\left (n+\frac {3N}{4}\right)\right]\right\} W_N^n W_{N/4}^{nk} \\ X (4k+3) &= \sum_{n=0}^{N/4-1} \left\{\left [x (n) - x\left (n+\frac {N}{2}\right)\right] + j\left [x\left (n+\frac {N}{4}\right) - x\left (n+\frac {3N}{4}\right)\right]\right\} W_N^{3n} W_{N/4}^{nk} \end{aligned} X(4k+1)X(4k+3)​=n=0∑N/4−1​{[x(n)−x(n+2N​)]−j[x(n+4N​)−x(n+43N​)]}WNn​WN/4nk​=n=0∑N/4−1​{[x(n)−x(n+2N​)]+j[x(n+4N​)−x(n+43N​)]}WN3n​WN/4nk​​

Figure TC.3.12 shows the flow graph for an in-place 32-point decimation-in-frequency SFFT algorithm.
图 TC.3.12 示出了原位 32 点按频率抽取 SFFT 算法的流程图。

Figure TC.3.12 Length 32 split-radix FFT algorithms from paper by Duhamel (1986); reprinted with permission from the IEEE
图 TC.3.12 Duhamel（1986）论文中的长度 32 分裂基 FFT 算法；经 IEEE 许可转载

Figure TC.3.13 Butterfly for SRFFT algorithm.

图 TC.3.13 SRFFT 算法的蝶形运算。

Table TC.3.1 Number of Nontrivial Real Multiplcations and Additions to Compute an N N N-point Complex DFT
表 TC.3.1 计算 N N N 点复数 DFT 所需的非平凡实数乘法和加法次数

一篇文章搞懂 DFT

Sean 编辑于 2024-06-18 13:53

1. DFT 的定义

连续时间傅里叶变换的定义为：
X ( f ) = ∫ − ∞ + ∞ x ( t ) e − j 2 π f t d t X(f) = \int_{-\infty}^{+\infty}x(t)\, e^{-j2\pi f t}\,dt X(f)=∫−∞+∞​x(t)e−j2πftdt

离散傅里叶变换（DFT）由连续傅里叶变换离散化得到：
X ( m ) = ∑ n = 0 N − 1 x ( n ) e − j 2 π m n N X(m) = \sum_{n=0}^{N-1} x(n)\, e^{-j\frac{2\pi m n}{N}} X(m)=n=0∑N−1​x(n)e−jN2πmn​

式中：

• N N N：时域离散信号的点数
• n n n：时域离散信号序号，满足 n ∈ { 0 , 1 , … , N − 1 } n\in\{0,1,\dots,N-1\} n∈{0,1,…,N−1}
• m m m：频域信号序号，满足 m ∈ { 0 , 1 , … , N − 1 } m\in\{0,1,\dots,N-1\} m∈{0,1,…,N−1}

DFT 输入为 N N N 个时域离散样本，输出为 N N N 个复数值频域样本。

由欧拉公式：
e − j ω = cos ⁡ ω − j sin ⁡ ω e^{-j\omega} = \cos\omega - j\sin\omega e−jω=cosω−jsinω

DFT 可展开为实部虚部形式：
X ( m ) = ∑ n = 0 N − 1 x ( n ) [ cos ⁡ ( 2 π m n N ) − j sin ⁡ ( 2 π m n N ) ] X(m) = \sum_{n=0}^{N-1} x(n)\left[ \cos\left(\frac{2\pi m n}{N}\right)- j\sin\left(\frac{2\pi m n}{N}\right) \right] X(m)=n=0∑N−1​x(n)[cos(N2πmn​)−jsin(N2πmn​)]

2. DFT 的实例分析

设时域连续信号为：
x ( t ) = sin ⁡ ( 2 π ⋅ 1000 ⋅ t ) + 0.5 sin ⁡ ( 2 π ⋅ 2000 ⋅ t + 3 π 4 ) x(t) = \sin(2\pi\cdot 1000\cdot t)+ 0.5\sin\left(2\pi\cdot 2000\cdot t + \frac{3\pi}{4}\right) x(t)=sin(2π⋅1000⋅t)+0.5sin(2π⋅2000⋅t+43π​)

以采样频率 f s = 8000  Hz f_s=8000\ \text{Hz} fs​=8000 Hz 采样，采样周期 T s = 1 / 8000  s T_s=1/8000\ \text{s} Ts​=1/8000 s，得到离散信号：
x ( n ) = sin ⁡ ( 2 π ⋅ 1000 ⋅ n T s ) + 0.5 sin ⁡ ( 2 π ⋅ 2000 ⋅ n T s + 3 π 4 ) = sin ⁡ ( 2 π n 8 ) + 0.5 sin ⁡ ( 4 π n 8 + 3 π 4 ) \begin{aligned} x(n) &= \sin(2\pi\cdot 1000\cdot n T_s)+ 0.5\sin\left(2\pi\cdot 2000\cdot n T_s + \frac{3\pi}{4}\right) \\ &= \sin\left(\frac{2\pi n}{8}\right)+ 0.5\sin\left(\frac{4\pi n}{8}+\frac{3\pi}{4}\right) \end{aligned} x(n)​=sin(2π⋅1000⋅nTs​)+0.5sin(2π⋅2000⋅nTs​+43π​)=sin(82πn​)+0.5sin(84πn​+43π​)​

取 N = 8 N=8 N=8，即 n = 0 , 1 , … , 7 n=0,1,\dots,7 n=0,1,…,7，得到离散样本：
x ( 0 ) = 0.3535 , x ( 1 ) = 0.3535 , x ( 2 ) = 0.6464 , x ( 3 ) = 1.0607 , x ( 4 ) = 0.3535 , x ( 5 ) = − 1.0607 , x ( 6 ) = − 1.3535 , x ( 7 ) = − 0.3535 \begin{aligned} &x(0)= 0.3535,\quad x(1)= 0.3535,\quad x(2)= 0.6464,\quad x(3)= 1.0607,\\ &x(4)= 0.3535,\quad x(5)=-1.0607,\quad x(6)=-1.3535,\quad x(7)=-0.3535 \end{aligned} ​x(0)=0.3535,x(1)=0.3535,x(2)=0.6464,x(3)=1.0607,x(4)=0.3535,x(5)=−1.0607,x(6)=−1.3535,x(7)=−0.3535​

Fig.1 输入信号

视角一：以 n n n 为自变量计算 DFT

对固定序号 m = i m=i m=i：

• X ( i ) X(i) X(i) 的实部： x ( n ) x(n) x(n) 与 cos ⁡ ( 2 π i n 8 ) \cos\left(\frac{2\pi i n}{8}\right) cos(82πin​) 做点积
• X ( i ) X(i) X(i) 的虚部： x ( n ) x(n) x(n) 与 − sin ⁡ ( 2 π i n 8 ) -\sin\left(\frac{2\pi i n}{8}\right) −sin(82πin​) 做点积

本例 DFT 计算结果：
X ( 0 ) = 0.0 − j 0.0 , X ( 1 ) = 0.0 − j 4.0 , X ( 2 ) = 1.414 + j 1.414 , X ( 3 ) = 0.0 + j 0.0 , X ( 4 ) = 0.0 − j 0.0 , X ( 5 ) = 0.0 − j 0.0 , X ( 6 ) = 1.414 − j 1.414 , X ( 7 ) = 0.0 + j 4.0 \begin{aligned} &X(0)= 0.0-j0.0,\quad X(1)= 0.0-j4.0,\quad X(2)= 1.414+j1.414,\quad X(3)= 0.0+j0.0,\\ &X(4)= 0.0-j0.0,\quad X(5)= 0.0-j0.0,\quad X(6)= 1.414-j1.414,\quad X(7)= 0.0+j4.0 \end{aligned} ​X(0)=0.0−j0.0,X(1)=0.0−j4.0,X(2)=1.414+j1.414,X(3)=0.0+j0.0,X(4)=0.0−j0.0,X(5)=0.0−j0.0,X(6)=1.414−j1.414,X(7)=0.0+j4.0​

Fig. 2 X ( 0 ) X(0) X(0) 的计算方法

Fig. 3 X ( 1 ) X(1) X(1) 的计算方法

Fig. 4 X ( 2 ) X(2) X(2) 的计算方法

Fig. 5 X ( 3 ) X(3) X(3) 的计算方法

Fig. 6 X ( 4 ) X(4) X(4) 的计算方法

Fig. 7 X ( 5 ) X(5) X(5) 的计算方法

Fig. 8 X ( 6 ) X(6) X(6) 的计算方法

Fig. 9 X ( 7 ) X(7) X(7) 的计算方法

视角二：以频率索引 k k k 为自变量观察 DFT

在频域视角下，离散傅里叶变换（DFT）从频率维度拆解时域信号的构成，频率索引 k k k 作为自变量可揭示频域信号的构成规律：对每个 k k k 计算 DFT，本质是将时域信号 x ( n ) x(n) x(n) 与复指数基函数 e − j 2 π k n / N e^{-j2\pi kn/N} e−j2πkn/N 做内积，即对 N N N 个等距频率分量的加权求和。

1. 数学推导基础

依据DFT的标准定义：
X ( k ) = ∑ n = 0 N − 1 x ( n ) e − j 2 π k n N k = 0 , 1 , … , N − 1 X(k)=\sum_{n=0}^{N-1} x(n)\,e^{-j\frac{2\pi kn}{N}}\qquad k=0,1,\dots ,N-1 X(k)=n=0∑N−1​x(n)e−jN2πkn​k=0,1,…,N−1

利用欧拉公式 e − j θ = cos ⁡ θ − j sin ⁡ θ e^{-j\theta}= \cos\theta - j\sin\theta e−jθ=cosθ−jsinθ 将复指数展开为正、负正弦形式，可得：
X ( k ) = ∑ n = 0 N − 1 x ( n ) [ cos ⁡ ⁣ ( 2 π k n N ) − j sin ⁡ ⁣ ( 2 π k n N ) ] = ∑ n = 0 N − 1 x ( n ) cos ⁡ ⁣ ( 2 π k n N ) ⏟ Re ⁡ { X ( k ) } − j ∑ n = 0 N − 1 x ( n ) sin ⁡ ⁣ ( 2 π k n N ) ⏟ Im ⁡ { X ( k ) } . \begin{aligned} X(k) &= \sum_{n=0}^{N-1} x(n)\bigl[\cos\!\bigl(\tfrac{2\pi kn}{N}\bigr)-j\sin\!\bigl(\tfrac{2\pi kn}{N}\bigr)\bigr] \\ &= \underbrace{\sum_{n=0}^{N-1} x(n)\cos\!\bigl(\tfrac{2\pi kn}{N}\bigr)}_{\displaystyle \operatorname{Re}\{X(k)\}} \;-\;j\underbrace{\sum_{n=0}^{N-1} x(n)\sin\!\bigl(\tfrac{2\pi kn}{N}\bigr)}_{\displaystyle \operatorname{Im}\{X(k)\}} . \end{aligned} X(k)​=n=0∑N−1​x(n)[cos(N2πkn​)−jsin(N2πkn​)]=Re{X(k)}n=0∑N−1​x(n)cos(N2πkn​)​​−jIm{X(k)}n=0∑N−1​x(n)sin(N2πkn​)​​.​

2. 实部与虚部的规范表述

3. 要素说明

每个频域点 X ( k ) X(k) X(k) 由 N N N 组不同频率的正余弦信号叠加构造，关键要素如下：

注意：k k k 与 n n n 取值范围相同，但物理意义截然不同；k k k 固定、 n n n 遍历时，求得单一频域点 X ( k ) X(k) X(k) 的值；k k k 取遍 0 , 1 , … , N − 1 0,1,\dots ,N-1 0,1,…,N−1 时，可得到完整的频谱 { X ( 0 ) , X ( 1 ) , … , X ( N − 1 ) } \{X(0),X(1),\dots ,X(N-1)\} {X(0),X(1),…,X(N−1)}。

4. 完整 DFT 表达式

X ( k ) = ∑ n = 0 N − 1 x ( n ) e − j 2 π k n N = ∑ n = 0 N − 1 x ( n ) [ cos ⁡ ⁣ ( 2 π k n N ) − j sin ⁡ ⁣ ( 2 π k n N ) ] , k = 0 , 1 , … , N − 1 X(k)=\sum_{n=0}^{N-1} x(n)\,e^{-j\frac{2\pi kn}{N}} =\sum_{n=0}^{N-1} x(n)\Bigl[\cos\!\bigl(\tfrac{2\pi kn}{N}\bigr)-j\sin\!\bigl(\tfrac{2\pi kn}{N}\bigr)\Bigr], \qquad k=0,1,\dots ,N-1 X(k)=n=0∑N−1​x(n)e−jN2πkn​=n=0∑N−1​x(n)[cos(N2πkn​)−jsin(N2πkn​)],k=0,1,…,N−1

该形式清晰展示了频率索引 k k k 作为自变量、时域采样 x ( n ) x(n) x(n) 与基函数相位 2 π k n N \frac{2\pi kn}{N} N2πkn​ 交互产生的加权求和过程，便于在频域分析中系统刻画信号的频率分量。

DFT 结果与视角一一致：
X ( 0 ) = 0.0 − j 0.0 , X ( 1 ) = 0.0 − j 4.0 , X ( 2 ) = 1.414 + j 1.414 , X ( 3 ) = 0.0 + j 0.0 , X ( 4 ) = 0.0 − j 0.0 , X ( 5 ) = 0.0 − j 0.0 , X ( 6 ) = 1.414 − j 1.414 , X ( 7 ) = 0.0 + j 4.0 \begin{aligned} &X(0)= 0.0-j0.0,\quad X(1)= 0.0-j4.0,\quad X(2)= 1.414+j1.414,\quad X(3)= 0.0+j0.0,\\ &X(4)= 0.0-j0.0,\quad X(5)= 0.0-j0.0,\quad X(6)= 1.414-j1.414,\quad X(7)= 0.0+j4.0 \end{aligned} ​X(0)=0.0−j0.0,X(1)=0.0−j4.0,X(2)=1.414+j1.414,X(3)=0.0+j0.0,X(4)=0.0−j0.0,X(5)=0.0−j0.0,X(6)=1.414−j1.414,X(7)=0.0+j4.0​

Fig. 10 DFT 计算的第二视角

3. DFT 结果

本例 DFT 输出：
X ( 0 ) = 0.0 − j 0.0 , X ( 1 ) = 0.0 − j 4.0 , X ( 2 ) = 1.414 + j 1.414 , X ( 3 ) = 0.0 + j 0.0 , X ( 4 ) = 0.0 − j 0.0 , X ( 5 ) = 0.0 − j 0.0 , X ( 6 ) = 1.414 − j 1.414 , X ( 7 ) = 0.0 + j 4.0 \begin{aligned} &X(0)= 0.0-j0.0,\quad X(1)= 0.0-j4.0,\quad X(2)= 1.414+j1.414,\quad X(3)= 0.0+j0.0,\\ &X(4)= 0.0-j0.0,\quad X(5)= 0.0-j0.0,\quad X(6)= 1.414-j1.414,\quad X(7)= 0.0+j4.0 \end{aligned} ​X(0)=0.0−j0.0,X(1)=0.0−j4.0,X(2)=1.414+j1.414,X(3)=0.0+j0.0,X(4)=0.0−j0.0,X(5)=0.0−j0.0,X(6)=1.414−j1.414,X(7)=0.0+j4.0​

分别绘制实部、虚部、幅值、相位：

Fig. 11 DFT 计算结果

4. DFT 的共轭对称性

由计算结果可观察到：

∣ X ( 5 ) ∣ = ∣ X ( 3 ) ∣ , ∣ X ( 6 ) ∣ = ∣ X ( 2 ) ∣ , ∣ X ( 7 ) ∣ = ∣ X ( 1 ) ∣ Re ⁡ { X ( 5 ) } = Re ⁡ { X ( 3 ) } , Re ⁡ { X ( 6 ) } = Re ⁡ { X ( 2 ) } , Re ⁡ { X ( 7 ) } = Re ⁡ { X ( 1 ) } Im ⁡ { X ( 5 ) } = − Im ⁡ { X ( 3 ) } , Im ⁡ { X ( 6 ) } = − Im ⁡ { X ( 2 ) } , Im ⁡ { X ( 7 ) } = − Im ⁡ { X ( 1 ) } \begin{aligned} |X(5)| &= |X(3)|, \quad & |X(6)| &= |X(2)|, \quad & |X(7)| &= |X(1)| \\ \operatorname{Re}\{X(5)\} &= \operatorname{Re}\{X(3)\}, \quad & \operatorname{Re}\{X(6)\} &= \operatorname{Re}\{X(2)\}, \quad & \operatorname{Re}\{X(7)\} &= \operatorname{Re}\{X(1)\} \\ \operatorname{Im}\{X(5)\} &= -\operatorname{Im}\{X(3)\}, \quad & \operatorname{Im}\{X(6)\} &= -\operatorname{Im}\{X(2)\}, \quad & \operatorname{Im}\{X(7)\} &= -\operatorname{Im}\{X(1)\} \end{aligned} ∣X(5)∣Re{X(5)}Im{X(5)}​=∣X(3)∣,=Re{X(3)},=−Im{X(3)},​∣X(6)∣Re{X(6)}Im{X(6)}​=∣X(2)∣,=Re{X(2)},=−Im{X(2)},​∣X(7)∣Re{X(7)}Im{X(7)}​=∣X(1)∣=Re{X(1)}=−Im{X(1)}​

上述关系对应 DFT 的共轭对称性：

X ( m ) = X ∗ ( N − m ) X(m) = X^*(N-m) X(m)=X∗(N−m)

式中 ( ⋅ ) ∗ (\cdot)^* (⋅)∗ 表示复共轭。

受共轭对称性约束，DFT 输出中后 N / 2 − 1 N/2-1 N/2−1 个样本存在信息冗余，有效输出可视为 N / 2 + 1 N/2+1 N/2+1 个独立复数。

5. 序号 m m m 的物理含义

在 N N N 点 DFT 中，序号 m m m 对应信号在 N N N 个点内包含 m m m 个周期。

对采样频率 f s f_s fs​，第 m m m 根谱线对应的频率为：
f m = m ⋅ f s N f_m = m\cdot\frac{f_s}{N} fm​=m⋅Nfs​​

本例中：

• m = 1 ⇒ f 1 = 1000  Hz m=1 \Rightarrow f_1 = 1000\ \text{Hz} m=1⇒f1​=1000 Hz
• m = 2 ⇒ f 2 = 2000  Hz m=2 \Rightarrow f_2 = 2000\ \text{Hz} m=2⇒f2​=2000 Hz

6. DFT 幅值的含义

原始信号中：

• 1000 Hz 分量幅值为 1
• 2000 Hz 分量幅值为 0.5

而 DFT 幅值分别为 4 和 2。
原因在于：傅里叶变换输出为频谱密度，与实际频谱之间相差系数 N / 2 N/2 N/2。本例 N = 8 N=8 N=8，对应比例为 4。

7. IDFT 的定义

离散傅里叶反变换（IDFT）定义为：
x ( n ) = 1 N ∑ m = 0 N − 1 X ( m ) e j 2 π m n N x(n) = \frac{1}{N}\sum_{m=0}^{N-1} X(m)\, e^{j\frac{2\pi m n}{N}} x(n)=N1​m=0∑N−1​X(m)ejN2πmn​

与 DFT 相比有两处区别：

• 多系数 1 / N 1/N 1/N
• 指数项为 + j 2 π m n N +j\frac{2\pi m n}{N} +jN2πmn​

利用欧拉公式展开：
x ( n ) = 1 N ∑ m = 0 N − 1 X ( m ) [ cos ⁡ ( 2 π m n N ) + j sin ⁡ ( 2 π m n N ) ] x(n) = \frac{1}{N}\sum_{m=0}^{N-1} X(m)\left[ \cos\left(\frac{2\pi m n}{N}\right)+ j\sin\left(\frac{2\pi m n}{N}\right) \right] x(n)=N1​m=0∑N−1​X(m)[cos(N2πmn​)+jsin(N2πmn​)]

IDFT 以 N N N 个复频域样本为输入，恢复出 N N N 个时域离散样本。

参考文献

• Fast Fourier Transform (FFT)
  https://www.cmlab.csie.ntu.edu.tw/cml/dsp/training/coding/transform/fft.html
• 一篇文章搞懂 DFT - 知乎
  https://zhuanlan.zhihu.com/p/365872805