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

【摘要】

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

---

---

1. 离散傅里叶变换(DFT)和逆离散傅里叶变换(IDFT)

1.1 基本公式

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

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

1.2 代码实现

以$x(n)=[1,1,1,1,1,1,1,1,1,1]$为原始时域信号序列，其长度为$L=10$，做$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'$2\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'$2\pi$'])

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

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

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

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

当$L<N$，相当于在信号后面补了$N-L$个零。设置$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()