IDFT
IDFT(Inverse Discrete Fourier Transform), 傅里叶逆变换,可以将频域信号转换到时域中, 它的公式非常简单:
x [ n ] = 1 N ∑ k = 0 N − 1 X [ k ] e j 2 π k n / N x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j2\pi kn/N} x[n]=N1k=0∑N−1X[k]ej2πkn/N
X [ k ] X[k] X[k]:离散频率下标为k时的频率大小
x [ n ] x[n] x[n]: 离散时域信号序列
N N N: 信号序列的长度,也就是采样的个数
对比我们之前讲过的DFT,两者公式类似,但是注意在DFT中指数带负号,而IDFT中不带
从矩阵的角度看IDFT
DFT的矩阵表示
讲IDFT之前,我们先复习DFT的矩阵表示形式:
[ s 0 0 s 0 1 ⋯ s 0 N − 1 ⋮ ⋮ ⋮ ⋮ s k 0 s k 1 ⋯ s k N − 1 ⋮ ⋮ ⋱ ⋮ s N − 1 0 s N − 1 1 ⋯ s N − 1 N − 1 ] [ x [ 0 ] x [ 1 ] ⋮ x [ n ] ⋮ x [ N − 1 ] ] = [ X [ 0 ] X [ 1 ] ⋮ X [ k ] ⋮ X [ N − 1 ] ] \begin{bmatrix} s_0^0 & s_0^1 & \cdots & s_0^{N-1} \\ \vdots & \vdots & \vdots & \vdots\\ s_k^0 & s_k^1 & \cdots & s_k^{N-1} \\ \vdots & \vdots & \ddots & \vdots\\ s_{N-1}^0 & s_{N-1}^1 & \cdots & s_{N-1}^{N-1} \\ \end{bmatrix} \begin{bmatrix} x[0] \\ x[1] \\ \vdots\\ x[n] \\ \vdots \\ x[N-1] \end{bmatrix} = \begin{bmatrix} X[0] \\ X[1] \\ \vdots\\ X[k] \\ \vdots \\ X[N-1] \end{bmatrix} ⎣⎢⎢⎢⎢⎢⎢⎡s00⋮sk0⋮sN−10s01⋮sk1⋮sN−11⋯⋮⋯⋱⋯s0N−1⋮skN−1⋮sN−1N−1⎦⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡x[0]x[1]⋮x[n]⋮x[N−1]⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡X[0]X[1]⋮X[k]⋮X[N−1]⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤
S S S矩阵中的每一行都是一个 S k S_k Sk向量, S k = e − j 2 π k n / N , n = 0 , 1 , ⋯   , N − 1 S_k = e^{-j2\pi kn/N}, n=0,1,\cdots,N-1 Sk=e−j2πkn/N,n=0,1,⋯,N−1,进一步简化上面的表示,得到:
[ ⋯ S 0 ⋯ ⋮ ⋯ S k ⋯ ⋮ ⋯ S N − 1 ⋯ ] [ x [ 0 ] x [ 1 ] ⋮ x [ n ] ⋮ x [ N − 1 ] ] = [ X [ 0 ] X [ 1 ] ⋮ X [ k ] ⋮ X [ N − 1 ] ] \begin{bmatrix} \cdots & S_0 & \cdots \\ & \vdots & \\ \cdots & S_k & \cdots \\ & \vdots & \\ \cdots & S_{N-1} & \cdots \\ \end{bmatrix} \begin{bmatrix} x[0] \\ x[1] \\ \vdots\\ x[n] \\ \vdots \\ x[N-1] \end{bmatrix} = \begin{bmatrix} X[0] \\ X[1] \\ \vdots\\ X[k] \\ \vdots \\ X[N-1] \end{bmatrix} ⎣⎢⎢⎢⎢⎢⎢⎡⋯⋯⋯S0⋮Sk⋮SN−1⋯⋯⋯⎦⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡x[0]x[1]⋮x[n]⋮x[N−1]⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡X[0]X[1]⋮X[k]⋮X[N−1]⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤
IDFT的矩阵表示
从IDFT的公式,可以看出,其实IDFT和DFT表示是一样的,只是对象发生了变化。具体来说,有两个变化:
- 由于指数部分不再有符号, S k S_k Sk进行了共轭操作,得到 S k ∗ S_k^* Sk∗
- 输入是频率信息X[k]
因此,矩阵表示变成了下面这样:
[ ⋯ S 0 ∗ ⋯ ⋮ ⋯ S k ∗ ⋯ ⋮ ⋯ S N − 1 ∗ ⋯ ] [ X [ 0 ] X [ 1 ] ⋮ X [ n ] ⋮ X [ N − 1 ] ] = [ x [ 0 ] x [ 1 ] ⋮ x [ k ] ⋮ x [ N − 1 ] ] \begin{bmatrix} \cdots & S_0^* & \cdots \\ & \vdots & \\ \cdots & S_k^* & \cdots \\ & \vdots & \\ \cdots & S_{N-1}^* & \cdots \\ \end{bmatrix} \begin{bmatrix} X[0] \\ X[1] \\ \vdots\\ X[n] \\ \vdots \\ X[N-1] \end{bmatrix} = \begin{bmatrix} x[0] \\ x[1] \\ \vdots\\ x[k] \\ \vdots \\ x[N-1] \end{bmatrix} ⎣⎢⎢⎢⎢⎢⎢⎡⋯⋯⋯S0∗⋮Sk∗⋮SN−1∗⋯⋯⋯⎦⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡X[0]X[1]⋮X[n]⋮X[N−1]⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡x[0]x[1]⋮x[k]⋮x[N−1]⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤
Talk is cheap, show me the code
接下来就简单多了,我们将先介绍如何使用scipy中ifft,然后自己动手实现一份ifft
导入必要的包
import numpy as np
from scipy.fftpack import fft, ifftimport matplotlib.pyplot as plt%matplotlib notebook
生成信号用于测试
def generate_sine(N, A, fs, f0, phi):'''N : number of samplesA : amplitudefs: sample ratef0: frequencyphi: initial phase'''T = 1/fsn = np.arange(N)x = A*np.cos( 2*np.pi*f0*n*T + phi )return x# generate signal
N = 501
A = 0.8
fs = 44100
f0 = 1000
phi = 0.0x = generate_sine(N, A, fs, f0, phi)plt.figure()
plt.plot(x)
plt.show()
使用scipy中的ifft
# fft the signal
N = 512 # fft size
X = fft(x, N)
mX = np.abs(X)
pX = np.angle(X)freq_axis = np.arange(N)/N * fs
plt.figure(figsize=(10, 12))
ax = plt.subplot(3,1,1)
plt.plot(freq_axis, mX)
ax.set_title('Magnitude')ax = plt.subplot(3,1,2)
plt.plot(freq_axis, pX)
ax.set_title('Phase')# ifft it
ifft_x = ifft(X)
ax = plt.subplot(3,1,3)
plt.plot(ifft_x)
ax.set_title('Synthesise')plt.show()
自己动手写ifft
只有两个地方要注意:
- 不要忘记乘上 1/N
- S k ∗ S_k^* Sk∗是 S k S_k Sk向量的共轭后的结果。反映在代码中,就是 S k ∗ S_k^* Sk∗不要共轭操作之间返回
def generate_complex_sinusoid(n, N):'''n : time index (or frequency index)N : number of sample'''k = np.arange(N)c_sin = np.exp(1j*2*np.pi*k*n/N)return c_sin# ifft loop
ifft_x = np.array([])for i in range(N):s = generate_complex_sinusoid(i, N)ifft_x = np.append(ifft_x, 1/N * np.sum(X*s))plt.figure()
plt.plot(ifft_x)
plt.show()
总结
通过自己动手,我们发现IDFT的原来和实现很简单,几乎与DFT一模一样,唯一需要注意的点就是 S k ∗ S_k^* Sk∗