离散傅里叶变换DFT

离散傅里叶变换DFT

---

DFT和IDFT

时域中长度为$N$的序列$x[n]$的离散傅里叶变换（DFT）和逆变换（IDFT）

$$X[k] = \sum_{n=0}^{N-1}x[n]\cdot \exp \left ( -j\frac{2\pi kn}{N} \right )$$

$$x[n]=\frac{1}{N}\sum_{n=0}^{N-1}X[k]\cdot \exp \left ( j\frac{2\pi kn}{N} \right )$$

为了方便记忆，常用一个中间变量

$$W_N=\exp \left( -j\frac{2\pi}{N}\right)$$
这样DFT和IDFT变成下面容易记忆的形式：
$$X[k] = \sum_{n=0}^{N-1}x[n] \cdot (W_N)^{kn}$$

$$x[n]=\frac{1}{N}\sum_{n=0}^{N-1}X[k]\cdot (W_N)^{-kn}$$
进一步，正反变换都可以看成是两个矩阵的乘法
$$X=W\cdot x=\left [\begin{matrix} 1 & 1 & \cdots & 1\\ 1 & (W_N)^1 & \cdots & (W_N)^{N-1} \\ 1 & (W_N)^2 & \cdots & (W_N)^{2(N-1)}\\ \vdots & \vdots & \vdots & \vdots\\ 1 & (W_N)^{N-1} & \cdots & (W_N)^{(N-1)^2}\end{matrix} \right ] \left[\begin{matrix}x[0]\\x[1]\\x[2]\\\vdots \\ x[N-1]\end{matrix}\right ]$$

$$x=W^{-1}\cdot X=\frac{1}{N}\left [\begin{matrix} 1 & 1 & \cdots & 1\\ 1 & (W_N)^{-1} & \cdots & (W_N)^{-(N-1)} \\ 1 & (W_N)^{-2} & \cdots & (W_N)^{-2(N-1)}\\ \vdots & \vdots & \vdots & \vdots\\ 1 & (W_N)^{-(N-1)} & \cdots & (W_N)^{-(N-1)^2}\end{matrix} \right ] \left[\begin{matrix}X[0]\\X[1]\\X[2]\\\vdots \\ X[N-1]\end{matrix}\right ]$$
留意到 $W$矩阵求逆变换到$W^{-1}$只需要将每个元素求倒数，再除以 N N <script type="math/tex" id="MathJax-Element-51">N</script>即可，不需要直接求逆，所以理论上是很快的。

---

numpy实现

代码参考自这里。

波形显示函数

def show(ori_func, ft, sampling_period = 5): 
    n = len(ori_func) 
    interval = sampling_period / float(n) 
    print interval
    # 绘制原始函数
    plt.subplot(2, 1, 1) 
    plt.plot(np.arange(0, sampling_period, interval), ori_func, 'black') 
    plt.xlabel('Time'), plt.ylabel('Amplitude') 
    # 绘制变换后的函数
    plt.subplot(2,1,2) 
    frequency = np.arange(n / 2) / (n * interval) 
    nfft = abs(ft[range(int(n / 2))] / n ) 
    plt.plot(frequency, nfft, 'red') 
    plt.xlabel('Freq (Hz)'), plt.ylabel('Amp. Spectrum') 
    plt.show() 

多个谐波叠加波形

time = np.arange(0, 5, 0.005) 
x1 = np.sin(2 * np.pi * 1 * time) # 频率为1
x2 = np.sin(2 * np.pi * 20 * time) # 频率为20
x3 = np.sin(2 * np.pi * 60 * time) # 频率为60
x = x1 + x2 + x3 # 叠加
y = np.fft.fft(x) 
show(x, y) 

验证numpy.fft.fft的计算原理

# 傅里叶变换DFT
x = np.random.random(500) 
N = len(x) 
n = np.arange(N) 
k = n.reshape((N, 1)) 
W = np.exp(-2j * np.pi * k * n / N) # 500*500 矩阵
y = np.dot(W, x) 
print np.allclose(y, np.fft.fft(x)) # 应该是 True

验证numpy.fft.ifft的计算原理

# 傅里叶逆变换IDFT
W_inv = np.exp(2j * np.pi * k * n / N) / N
x = np.dot(y, W_inv) 
print np.allclose(x, np.fft.ifft(y)) # 应该是 True

---

参考资料

《信号与系统（第二版）》Alan V. Oppenheim
《数字信号处理——基于计算机的方法（第四版）》Sanjit K. Mitra
【博客】 NumPy 中的傅里叶分析
【API】numpy.fft.fft