从DFT到FFT
关键词：离散傅里叶变换 快速傅里叶变换
摘要：快速傅里叶变换算法是现代工程技术中广泛使用的一种算法。文章从离散傅里叶变换与复指数函数引入，分析了快速傅里叶变换的基本思路，最后介绍了实际运算用到的蝶形网络。
傅里叶变换广泛应用于信号的频谱分析中，在现代工程中有着重大的意义。但是，作为信息时代核心的计算机系统只能够处理离散的数值，为了充分利用计算机的性能进行工程分析，很有必要把傅里叶变换从连续空间带到离散空间上来。
一、离散傅里叶变换
离散傅里叶变换（Discrete Fourier Transform，缩写为DFT），是傅里叶变换在时域和频域上都呈离散的形式，将信号的时域采样变换为其DTFT的频域采样。从连续傅里叶变换可以推导到DFT的正变换式为：
为了分析方便，转换成下面的矩阵运算：
def dft(x: Vector[Complex]): Vector[Complex] = {
val n = x.size
(0 until n).map { i =>
(0 until n).foldLeft(Complex.zero) { (sum, j) => sum + x(j) * w(j * i, n) }
}.toVector
}
def idft(y: Vector[Complex]): Vector[Complex] = {
val n = y.size
(0 until n).map { i =>
(0 until n).foldLeft(Complex.zero) { (sum, j) => sum + y(j) * w(-j * i, n) } / n
}.toVector
}
def w(k: Double, n: Double): Complex = Complex(Math.E, 0) pow -2 * Math.PI * Complex.i * k / n
用到的 Complex 类
// from Breeze: https://github.com/scalanlp/breeze/blob/master/math/src/main/scala/breeze/math/Complex.scala
case class Complex(real: Double, imag: Double) {
override def toString: String = s"$real + ${imag}i"
def +(that: Complex) =
Complex(this.real + that.real, this.imag + that.imag)
def +(that: Int) =
Complex(this.real + that, this.imag)
def +(that: Long) =
Complex(this.real + that, this.imag)
def +(that: Float) =
Complex(this.real + that, this.imag)
def +(that: Double) =
Complex(this.real + that, this.imag)
def -(that: Complex) =
Complex(this.real - that.real, this.imag - that.imag)
def -(that: Int) =
Complex(this.real - that, this.imag)
def -(that: Long) =
Complex(this.real - that, this.imag)
def -(that: Float) =
Complex(this.real - that, this.imag)
def -(that: Double) =
Complex(this.real - that, this.imag)
def *(that: Complex) =
Complex(this.real * that.real - this.imag * that.imag, this.real * that.imag + this.imag * that.real)
def *(that: Int) =
Complex(this.real * that, this.imag * that)
def *(that: Long) =
Complex(this.real * that, this.imag * that)
def *(that: Float) =
Complex(this.real * that, this.imag * that)
def *(that: Double) =
Complex(this.real * that, this.imag * that)
def /(that: Complex) = {
val denom = that.real * that.real + that.imag * that.imag
Complex(
(this.real * that.real + this.imag * that.imag) / denom,
(this.imag * that.real - this.real * that.imag) / denom)
}
def /(that: Int) =
Complex(this.real / that, this.imag / that)
def /(that: Long) =
Complex(this.real / that, this.imag / that)
def /(that: Float) =
Complex(this.real / that, this.imag / that)
def /(that: Double) =
Complex(this.real / that, this.imag / that)
def %(that: Complex) = {
val div = this./(that)
this - (Complex(floor(div.re()), floor(div.im())) * div)
}
def %(that: Int): Complex = this.%(Complex(that, 0))
def %(that: Long): Complex = %(Complex(that, 0))
def %(that: Float): Complex = %(Complex(that, 0))
def %(that: Double): Complex = %(Complex(that, 0))
def unary_- =
Complex(-real, -imag)
def abs =
math.sqrt(real * real + imag * imag)
def conjugate =
Complex(real, -imag)
def log =
Complex(math.log(abs), math.atan2(imag, real))
def exp = {
val expreal = math.exp(real)
Complex(expreal * math.cos(imag), expreal * math.sin(imag))
}
def pow(b: Double): Complex = pow(Complex(b, 0))
def pow(b: Complex): Complex = {
if (b == Complex.zero) Complex.one
else if (this == Complex.zero) {
if (b.imag != 0.0 || b.real < 0.0) Complex.nan
else Complex.zero
} else {
val c = log * b
val expReal = math.exp(c.real)
Complex(expReal * math.cos(c.imag), expReal * math.sin(c.imag))
}
}
override def equals(that: Any): Boolean = that match {
case that: Complex => this.real == that.real && this.imag == that.imag
case real: Double => this.real == real && this.imag == 0
case real: Int => this.real == real && this.imag == 0
case real: Short => this.real == real && this.imag == 0
case real: Long => this.real == real && this.imag == 0
case real: Float => this.real == real && this.imag == 0
case _ => false
}
// ensure hashcode contract is maintained for comparison to non-Complex numbers
// x ^ 0 is x
override def hashCode(): Int = real.## ^ imag.##
}
object Complex {
/** Constant Complex(0,0). */
val zero = new Complex(0, 0)
/** Constant Complex(1,0). */
val one = new Complex(1, 0)
/** Constant Complex(NaN, NaN). */
val nan = new Complex(Double.NaN, Double.NaN)
/** Constant Complex(0,1). */
val i = new Complex(0, 1)
implicit def realToComplex(re: Double): Complex = Complex(re, 0)
}
显然，为了计算 ，需要进行 N 次乘法（忽略 WNk 本身值的计算）与 N-1
次加法。在计算机中，一般乘法运算要明显慢于加法，所以下面都仅仅考虑乘法的次数。那么完成整个 DFT 需要进行
N
这种朴素方法易于理解与计算，但是 O(N^2)
的时间复杂度意味着计算时间随序列长度增加而增加的速度是相当高的平方级。这严重阻碍了 DFT 的实际应用。
二、WNk 的性质
我们用 WNk 表示 e2πjk/N 以简化公式， 通过对复指数 e2πjk/N 的研究，可以发现 WNk 的一些性质。 通过欧拉公式可以发现其表示复平面中的单位圆，那么显然有：
在单位圆上转整数圈回到原处（周期性）：
在单位圆上转半圈矢量方向相反：
利用这些性质可以消去朴素的DFT算法中的重复步骤。
三、快速傅里叶变换
快速傅里叶变换（Fast Fourier Transform，缩写为 FFT），是快速计算串行的离散傅里叶变换或其逆变换的方法[1]。快速傅里叶变换广泛的应用于工程、科学和数学领域。这里的基本思想在1965年才得到普及，但早在 1805 年就已推导出来[2]。库利-图基算法是最常用的快速傅里叶变换算法。这一算法利用分治法把 DFT 求解优化到了 O(NlogN)。在 1965 年 J.W. Cooley 和 J. W. Tukey 合作发表的 An algorithm for the machine calculation of complex Fourier series 中介绍这种方法以及 FFT 的基本思路。虽然后来发现，实际上高斯在 1805 年就已经提出了这样的算法。
其基本思路是先将偶数与奇数样本点序列分开，
方便起见，令：
由消去引理易得：
这意味着 X(k) 具备着一定的周期性，我们只需要计算一半的 X(k)。但是这还没有结束， 留意 G(k)、H(k) 与 X(k) 的形式非常相似，令：
则：
由此可得，G(k) 可以像 X(k) 一样继续奇偶分解，H(k) 同理。
def fft(x: Vector[Complex]): Vector[Complex] = {
val n = x.size
if (n == 1) return x
val g = fft(x.grouped(2).map(_.head).toVector) // 偶数索引信号
val h = fft(x.grouped(2).map(_.last).toVector) // 奇数索引信号
(0 until n).map { k => g(k % (n / 2)) + w(k, n) * h(k % (n / 2)) }.toVector
}
// Defferences between fft and ifft:
// 1. X(k) = G(k) + W(k,n)H(k) -> X(k) = G(k) + W(-k,n)H(k)
// 2. 递归结果缩放到原来的 1/N
def ifft(y: Vector[Complex]): Vector[Complex] = {
def backtrace(y: Vector[Complex]): Vector[Complex] = {
val n = y.size
if (n == 1) return y
val g = backtrace(y.grouped(2).map(_.head).toVector) // 偶数索引信号
val h = backtrace(y.grouped(2).map(_.last).toVector) // 奇数索引信号
(0 until n).map { k => g(k % (n / 2)) + w(-k, n) * h(k % (n / 2)) }.toVector
}
backtrace(y).map(_ / y.size)
}
// w(k,n) remains the same as dft
def w(k: Double, n: Double): Complex = Complex(Math.E, 0) pow -2 * Math.PI * Complex.i * k / n
实际计算时正向进行，为了方便递推引出蝶形运算：
其中输入点列按位反转（Bit Reversal）排序：首先依据二进制进行编号，然后对各个编号按位倒置并按此重新排序。 例如，对于一个8点变换：
001倒置以后变成 100
000 → 000
001 → 100
010 → 010
011 → 110
100 → 001
101 → 101
110 → 011
111 → 111
倒置后的编号为 {0,4,2,6,1,5,3,7}。依此顺序输入蝶形网络计算出的变换序列就是正序的。
def fft(input: Vector[Complex]): Vector[Complex] = {
val x = butterflyTransform(input)
val len = x.size
for (h <- Iterator.iterate(2)(_ * 2) takeWhile (_ <= len)) {
val wn = Complex(cos(-2 * Math.PI / h), sin(-2 * Math.PI / h))
x.indices by h foreach { j =>
var w = Complex.one
j until j + h / 2 foreach { k =>
val (u, t) = (x(k), w * x(k + h / 2))
x(k) = u + t
x(k + h / 2) = u - t
w = w * wn
}
}
}
x.toVector
}
def butterflyTransform(x: Vector[Complex]): mutable.Buffer[Complex] = {
val res = mutable.ArrayBuffer(x: _*)
val n = x.size
var (j, k) = (n / 2, 0)
1 until n - 1 foreach { i =>
if (i < j) {
res(i) -> res(j) match {
case (a, b) =>
res(j) = a
res(i) = b
}
}
k = n / 2
while (j >= k) {
j = j - k
k = k / 2
}
if (j < k) j += k
}
res
}
相较于DFT，FFT 在 N 很大时，计算量节省相当可观。不同点数时 DFT 与 FFT 的复数乘法次数和比较见下表[3]
| N | 8 | 16 | 32 | 64 | 128 | 256 | 512 | 1024 | 2048 |
|---|---|---|---|---|---|---|---|---|---|
| DFT | 64 | 256 | 1024 | 4096 | 16384 | 65536 | 262144 | 1048576 | 4194304 |
| FFT | 12 | 32 | 80 | 192 | 448 | 1024 | 2034 | 5120 | 11264 |
参考文献：
[1] 杨毅明. 数字信号处理（第2版）. 北京: 机械工业出版社. 2017年: 第95页. ISBN 9787111576235.
[2] Heideman, M. T.; Johnson, D. H.; Burrus, C. S. Gauss and the history of the fast Fourier transform. IEEE ASSP Magazine. 1984, 1 (4): 14–21. doi:10.1109/MASSP.1984.1162257.
[3] 李勇. 数字信号处理原理与应用 西安 西北工业大学出版社2016年: 第96页. ISBN 9787561248287
最后附上 GitHub：https://github.com/gonearewe
