基于Python的信号处理（2）——离散傅里叶变换DFT

【摘要】

• 使用Python实现离散傅里叶变换
• 参考书目：《数字信号处理——原理、算法与应用》第四版
• 这是一个笔记，不一定全面
• 本文主题：DFT、IDFT、numpy

---

目录

• 1. 离散傅里叶变换(DFT)和逆离散傅里叶变换(IDFT)
• • 1.1 基本公式
  • 1.2 代码实现
  • 1.3 绘制信号时域和频域图
  • 1.4 增大 N N N，提高频谱分辨率
  • 1.5 将频谱平移至[ − π -\pi −π, π \pi π]

---

1. 离散傅里叶变换(DFT)和逆离散傅里叶变换(IDFT) 1.1 基本公式

设离散信号 x ( n ) x(n) x(n)的长度为 L L L，变换成长度为 N N N( N ≥ L N\ge L N≥L)的频谱 { X ( k ) } k = 0 N − 1 \{X(k)\}_{k=0}^{N-1} {X(k)}k=0N−1​的过程，称为离散傅里叶变换(DFT)：
X ( k ) = ∑ n = 0 N − 1 x ( n ) e − j 2 π k n / N k = 0 , ⋯   , N − 1 X(k)=\sum_{n=0}^{N-1}x(n)e^{-j2\pi kn/N}\quad k=0,\cdots,N-1 X(k)=n=0∑N−1​x(n)e−j2πkn/Nk=0,⋯,N−1
从 X ( k ) X(k) X(k)中恢复 x ( n ) x(n) x(n)的公式(IDFT)为：
x ( n ) = ∑ k = 0 N − 1 X ( k ) e j 2 π k n / N n = 0 , ⋯   , N − 1 x(n)=\sum_{k=0}^{N-1}X(k)e^{j2\pi kn/N} \quad n=0,\cdots,N-1 x(n)=k=0∑N−1​X(k)ej2πkn/Nn=0,⋯,N−1

其中 N N N越大，频谱分辨率越高。

1.2 代码实现

以 x ( n ) = [ 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 ] x(n)=[1,1,1,1,1,1,1,1,1,1] x(n)=[1,1,1,1,1,1,1,1,1,1]为原始时域信号序列，其长度为 L = 10 L=10 L=10，做 N = 100 N=100 N=100点DFT。

• 生成信号：

L, N = 10, 100
x = np.ones(L) # 时域离散信号
X = np.zeros(N) + np.zeros(N)*1j # 频域频谱

• DFT:

# DTF变换
for k in range(N):
    for n in range(L):
        X[k] = X[k] + x[n]*np.exp(-1j*n*k/N*2*np.pi)

• IDFT变换：

# IDFT变换
x_p = np.zeros(L)
for n in range(L):
    for k in range(N):
        x_p[n] = x_p[n] + 1/N*X[k]*np.exp(1j*n*k/N*2*np.pi)

1.3 绘制信号时域和频域图

fig = plt.figure(figsize=(10, 10))
a1 = plt.subplot2grid((4, 1), (0, 0))
a2 = plt.subplot2grid((4, 1), (1, 0))
a3 = plt.subplot2grid((4, 1), (2, 0))
a4 = plt.subplot2grid((4, 1), (3, 0))

a1.stem(x) #原始序列
a2.stem(np.abs(X)) #DFT频谱（幅度）
a3.stem(np.angle(X, deg=True)) # DFT频谱（相位）
a4.stem(x_p) # IDFT序列

a1.set_ylabel(r'Original sequence')

a2.set_ylabel(r'DTF:|X($\omega$)|')
a2.set_xticks([0, 25, 50,75, 99])
a2.set_xticklabels(['0', r'$\frac{\pi}{2}$',r'$\pi$',r'$\frac{3\pi}{2}$',r'\pi$'])

a3.set_ylabel(r'DFT:$\theta(\omega)$')
a3.set_xticks([0, 25, 50,75, 99])
a3.set_xticklabels(['0', r'$\frac{\pi}{2}$',r'$\pi$',r'$\frac{3\pi}{2}$',r'\pi$'])

a4.set_ylabel(r'IDFT sequence')
plt.show()

图1 做N=100点的DFT和IDFT变换

观察图1中的相位，发现其变化范围为[-180,180]。DFT变换得到的第 k k k个频谱点为X(k)=a+bi，是一个复数。在求其角度的时候可以由两种方式：
(1) 直接使用np.angle()函数。
(2) 使用np.arctan2(b/a)。注意这里使用的是arctan2，不是arctan。arctan2会更具a和b的正负情况判断角度。

1.4 增大 N N N，提高频谱分辨率

当 L < N LL<N，相当于在信号后面补了 N − L N-L N−L个零。设置 N = 1000 N=1000 N=1000，发现频谱分辨率得到了明显的提升。

图2 做N=1000点的DFT和IDFT变换

1.5 将频谱平移至[ − π -\pi −π, π \pi π]

图3 DFT后频谱进行平移（N=100）

绘图代码如下：

Y = np.roll(X, int(N/2)) # 平移频谱
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.stem(np.abs(Y))
ax.set_xticks([0, 25, 50,75, 99])
ax.set_xticklabels(['-$\pi$', r'$-\frac{\pi}{2}$',r'0',r'$\frac{\pi}{2}$',r'$\pi$'])
ax.set_ylabel(r'DTF:|X($\omega$)|')
plt.grid(True)
plt.show()