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] 重磅推出离散傅里叶变换