DFT和FFT
为什么需要DFT
数字语音信号是离散时间信号，对其进行频域分析可以通过离散时间傅里叶变换（Discrete-time Fourier transform, DTFT）或者离散傅里叶变换（Discrete Fourier transform, DFT）。二者的区别在于，DTFT作用于时域离散的非周期信号，变换到频域后得到的是连续的周期信号；DFT作用于时域离散的周期信号，变换到频域后得到的是离散的周期信号\(^{[1]}\)。
当使用计算机处理信号时，由于其不能无限密集地存储连续时间信号，因此必须将信号+结果进行离散化处理\(^{[2]}\)。由于DFT的输入和输出都是离散信号，因此方便了使用计算机进行频域分析。
DFT公式与代码实现
对于一个周期为\(N\)的序列\(x[n]\)，其可以被表示为：
\(x[n]= \sum_{k=0}^{N-1 } a_{k} e^{jk(2\pi /N)n}\)
我们称上式中的\(a_{k}\)为傅里叶级数系数。这个公式的意义是，该信号可以表示为一系列离散时间复指数信号加权后的和。
经过一系列推导后\(^{[2]}\)，为了得到离散时间信号的傅里叶变换，我们将离散时间信号的傅里叶级数系数乘以其周期\(N\)后可得：
\(X[k] = a_{k} \cdot N = \sum_{n=0}^{N-1} x[n] \cdot e^{-j2\pi kn/N} = \sum_{n=0}^{N-1} x[n] [cos(2\pi kn/N) - i\cdot sin(2\pi kn/N)]\)
这个公式的python代码为：
def dft_1(x):
# Original code, from https://pythonnumericalmethods.berkeley.edu/notebooks/chapter24.02-Discrete-Fourier-Transform.html
N = len(x) # 信号长度
n = np.arange(N)
k = n.reshape((N,1))
e = np.exp(-2j * np.pi * k * n /N)
X = np.dot(e, x)
return X
def dft_2(x):
# codes from https://www.zhihu.com/question/41427782
N = len(x)
WN = np.power(np.e, complex(0,-1) * (2*np.pi) / N) # 即dft_1()里的e
Xk = []
for k in range(N):
value = 0
for n in range(N):
value += x[n] * np.power(WN, n*k)
Xk.append(value)
return np.array(Xk)
第二个函数dft_2()更容易理解，但计算时间很长。这个在接下来讲FFT时会展示。
跑下代码
现在，先演示下如何生成一段信号并将其转换到频域的：
import matplotlib.pyplot as plt
import numpy as np
plt.style.use('seaborn-poster')
# 首先生成一段信号
# sampling rate
sr = 100
ts = 1.0/sr
t = np.arange(0,1,ts)
# add signal
freq = 1.0
x = 3*np.sin(2*np.pi*freq*t)
freq = 4.0
x += np.sin(2*np.pi*freq*t)
freq = 7.0
x += 0.5*np.sin(2*np.pi*freq*t)
freq = 8.0
x += 0.5*np.cos(2*np.pi*freq*t)
freq = 8.0
x += 0.5*np.sin(2*np.pi*freq*t)
# do DFT
X = dft_1(x)
N = len(X)
n = np.arange(N)
T = N / sr
freq = n/T
n_oneside = N//2
f_oneside = freq[:n_oneside]
X_oneside = X[:n_oneside]/n_oneside
plt.figure(figsize = (16, 6))
plt.subplot(131)
plt.plot(t, x, 'r')
plt.ylabel('Amplitude')
plt.subplot(132)
plt.stem(f_oneside, abs(X_oneside), 'b', markerfmt=" ", basefmt="-b")
plt.xlabel('Freq (Hz)')
plt.ylabel('DFT Amplitude |X(freq)|')
plt.subplot(133)
plt.stem(f_oneside, abs(X_oneside), 'b', markerfmt=" ", basefmt="-b")
plt.xlabel('Freq (Hz)')
plt.xlim(0, 10)
plt.tight_layout()
plt.show()
最左边红色的是原始信号，中间的是转换后整个频域内的幅度，右边是在[0,10]上的放大图。可以看到经过多次叠加后，原始的时域信号已经分辨不出是如何得来的了，但经过DFT变换后，能够清楚地看到信号在1， 4， 7， 8Hz处有幅度，这跟我们定义信号时设置的freq相同。
我们下来定义一个只有10个点的短信号，看一下在DFT前后的数值变化：
# sampling rate，只有10个点
sr = 10
ts = 1.0/sr
t = np.arange(0,1,ts)
# add signal
freq = 1.0
x = 3*np.sin(2*np.pi*freq*t)
freq = 2.0
x += np.cos(2*np.pi*freq*t)
freq = 3.0
x += 0.5*np.sin(2*np.pi*freq*t)
freq = 4.0
x += 0.5*np.cos(2*np.pi*freq*t)
freq = 4.0
x += 0.5*np.sin(2*np.pi*freq*t)
此时的x为：
array([ 1.5 , 2.43728514, 1.42924017, 2.38029668, 1.84949989, 1.5 , -2.04048289, -3.68931368, -2.73825716, -2.62826814])
得到的DTF结果X_oneside为：
array([-6.21724894e-16+0.00000000e+00j, 4.44089210e-16-3.00000000e+00j, 1.00000000e+00+7.10542736e-16j, 1.04360964e-15-5.00000000e-01j, 5.00000000e-01-5.00000000e-01j])
可以看到X_oneside在实部上的值为定义信号时cos的系数，在虚部上的值为sin的系数（欧拉公式）。
计算耗时
前文里的两个DFT函数在计算时间上也不相同。在我的电脑上对一段2000点的信号进行DFT，dft_1()的执行时间为226 ms，而尽管dft_2()更方便看懂，它的执行时间达到了18.1 s。如果使用快速傅里叶变换（Fast Fourier transform, FFT）的话，这个数字是19.9 ms。为什么呢？
快速计算 -- FFT
由于DFT得到的信号是离散周期的，因此有：
$ X_{k+N} = \sum_{n=0}^{N-1} x_{n} \cdot e^{-i2\pi (k+N)n/N} = \sum_{n=0}^{N-1} x_{n} \cdot e^{-i2\pi kn/N} \cdot e^{-i2\pi n} $
同时由于 \(e^{-i2\pi n}=1\)，因此有:
$ X_{k+N} = \sum_{n=0}^{N-1} x_{n} \cdot e^{-i2\pi kn/N} = X_{k} $
所以，对任意的整数\(i\)，有:
\(X_{k+i\cdot N} = X_{k}\)
利用这个性质，其实我们在计算DFT时，并不需要将中间步骤中每一次的结果都算出来。例如，可以把结果\(k\)按奇数和偶数分开：
$ X_{k} = \sum_{n=0}^{N-1} x_{n}\cdot e^{-i2\pi kn/N} \ =\sum_{m=0}^{N/2-1} x_{2m}\cdot e^{-i2\pi k(2m)/N} + \sum_{m=0}^{N/2-1}x_{2m+1}\cdot e^{-i2\pi k(2m+1)/N} \ =\sum_{m-0}^{N/2-1}x_{2m}\cdot e^{-i2\pi km/(N/2)} + e^{-i2\pi k/N}\sum_{m=0}^{N/2-1}x_{2m+1}\cdot e^{-i2\pi km/(N/2)}$
这样，对一个\(N\)长的序列进行DFT就变成了对两个\(\frac{N}{2}\)长的序列进行DFT。现在这只是一个简单的示范，实际上我们还可以接着对这两个序列进行同样的分解操作，直到最后剩下两个数为止：
def FFT(x):
N = len(x)
if N == 1:
return x
else:
X_even = FFT(x[::2])
X_odd = FFT(x[1::2])
factor = np.exp(-2j*np.pi*np.arange(N)/N)
X = np.concatenate([X_even + factor[:int(N/2)]*X_odd,
X_even + factor[int(N/2):]*X_odd])
return X
跑下代码
看下这个FFT函数是怎么工作的。首先定义一个\(N=4\)的输入：[ 0.0000000e+00 3.0000000e+00 3.6739404e-16 -3.0000000e+00]
Length = 4, X_even = [0.0000000e+00 3.6739404e-16], X_odd = [ 3. -3.]
1.0 X_even calling FFT([0.0000000e+00 3.6739404e-16])
1.1 Now, X_even_even = [0.], X_even_odd = [3.6739404e-16]
1.2 X_even_even calling FFT([0.]) gets [0.]
1.3 X_even_odd calling FFT([3.6739404e-16]) gets [3.6739404e-16]
1.4 Now X_even_even = [0.], X_even_odd = [3.6739404e-16]
X_1 = np.concatenate([X_even_even + factor[:int(N/2)]*X_even_odd ,
X_even_even + factor[int(N/2):]*X_even_odd ])
= [ 3.6739404e-16+0.00000000e+00j -3.6739404e-16-4.49927935e-32j]
2.0 X_odd calling FFT([ 3. -3.])
2.1 Now, X_odd_even = [3.], X_odd_odd = [-3.]
2.2 X_odd_even calling FFT([3.]) gets [3.]
2.3 X_odd_odd calling FFT([-3.]) gets [-3.]
2.4 Now X_odd_even = [3.], X_odd_odd = [-3.]
X_2 = np.concatenate([X_odd_even + factor[:int(N/2)]*X_odd_odd ,
X_odd_even + factor[int(N/2):]*X_odd_odd ])
= [0.+0.0000000e+00j 6.+3.6739404e-16j]
3.0 X = np.concatenate([X_1 + factor[:int(N/2)]*X_2,
X_1 + factor[int(N/2):]*X_2])
= [ 3.6739404e-16+0.j 3.6739404e-16-6.j 3.6739404e-16+0.j -1.8369702e-15+6.j]
[1] 定性聊聊DFT、DFS和DTFT之间的关系
[2] 重磅推出离散傅里叶变换
