离散傅里叶变换DFT

DFT和IDFT

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

X[k]=∑n=0N−1x[n]⋅exp(−j2πknN) X [ k ] = ∑ n = 0 N − 1 x [ n ] ⋅ exp ⁡ ( − j 2 π k n N )

x[n]=1N∑n=0N−1X[k]⋅exp(j2πknN) x [ n ] = 1 N ∑ n = 0 N − 1 X [ k ] ⋅ exp ⁡ ( j 2 π k n N )

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

WN=exp(−j2πN) W N = exp ⁡ ( − j 2 π N )

这样DFT和IDFT变成下面容易记忆的形式：

X[k]=∑n=0N−1x[n]⋅(WN)kn X [ k ] = ∑ n = 0 N − 1 x [ n ] ⋅ ( W N ) k n

x[n]=1N∑n=0N−1X[k]⋅(WN)−kn x [ n ] = 1 N ∑ n = 0 N − 1 X [ k ] ⋅ ( W N ) − k n

进一步，正反变换都可以看成是两个矩阵的乘法

X=W⋅x=⎡⎣⎢⎢⎢⎢⎢⎢111⋮11(WN)1(WN)2⋮(WN)N−1⋯⋯⋯⋮⋯1(WN)N−1(WN)2(N−1)⋮(WN)(N−1)2⎤⎦⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢x[0]x[1]x[2]⋮x[N−1]⎤⎦⎥⎥⎥⎥⎥⎥ X = W ⋅ x = [ 1 1 ⋯ 1 1 ( W N ) 1 ⋯ ( W N ) N − 1 1 ( W N ) 2 ⋯ ( W N ) 2 ( N − 1 ) ⋮ ⋮ ⋮ ⋮ 1 ( W N ) N − 1 ⋯ ( W N ) ( N − 1 ) 2 ] [ x [ 0 ] x [ 1 ] x [ 2 ] ⋮ x [ N − 1 ] ]

x=W−1⋅X=1N⎡⎣⎢⎢⎢⎢⎢⎢111⋮11(WN)−1(WN)−2⋮(WN)−(N−1)⋯⋯⋯⋮⋯1(WN)−(N−1)(WN)−2(N−1)⋮(WN)−(N−1)2⎤⎦⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢X[0]X[1]X[2]⋮X[N−1]⎤⎦⎥⎥⎥⎥⎥⎥ x = W − 1 ⋅ X = 1 N [ 1 1 ⋯ 1 1 ( W N ) − 1 ⋯ ( W N ) − ( N − 1 ) 1 ( W N ) − 2 ⋯ ( W N ) − 2 ( N − 1 ) ⋮ ⋮ ⋮ ⋮ 1 ( W N ) − ( N − 1 ) ⋯ ( W N ) − ( N − 1 ) 2 ] [ X [ 0 ] X [ 1 ] X [ 2 ] ⋮ X [ N − 1 ] ]

留意到 W W 矩阵求逆变换到W−1W−1只需要将每个元素求倒数，再除以 N N <script type="math/tex" id="MathJax-Element-51">N</script>即可，不需要直接求逆，所以理论上是很快的。

numpy实现

代码参考自这里。

波形显示函数

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

多个谐波叠加波形

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

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

# 傅里叶变换DFT
x = np.random.random() 
N = len(x) 
n = np.arange(N) 
k = n.reshape((N, )) 
W = np.exp(- * 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( * 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