FFT变换中的两种分辨率(物理分辨率和计算分辨率)

设时域采样序列为$x[n],0\le n<N$，则其FFT变换为
$$X[k]=\sum _{n=0}^{N-1}x[n]exp(-j\frac{2\pi }{N}kn),0\le n<N,0\le k<N$$
若采样率为$f_s$，则FFT变换后的频率轴刻度为$[-f_s/2:f_s/N:f_s/2]$。从中得到两个信息：
（1）受奈奎斯特采样定律的限制，FFT变换最大只能表征$f_s/2$频率的信号。当信号频率超过$f_s/2$，则会发生频谱混叠。设信号真实频率为$f_0$，则混叠后的频率为$\pm |n{f}_s-f_0|，n=round(f_0/f_s)$。
下图所示场景中，$f_s=10000Hz,N=2000$，包含两个正弦信号，频率分别为$1000Hz$和$6000Hz$。FFT变换后，$6000Hz$信号发生混叠，混叠后频率为4000Hz。

图1. 采样率10000Hz，两个单频信号频率分别为1000Hz和6000Hz

（2）频谱中相邻两点间的频率间隔为$\Delta f=f_s/N$，显然采样率$f_s$越小，时域序列样点数$N$越大，$\Delta f$就越小。
基于上面的结论，我们很容易想到，若在原始采样序列$x[n]$后面补上一串0，再进行FFT变换，那就通过增加时域样点数$N$达到减小FFT频谱的频率分辨率$\Delta f$的目的。那是不是意味着，单纯的补零操作就可以让FFT变换具备分辨频率差异任意小的信号的能力？显然，这是违背直觉的，因为补零本质上并没有增加任何有用信息，自然无法提升FFT对信号的分辨能力。那补零操作的本质是什么？它与FFT对信号分辨能力的关系又是什么？
其实这边涉及到两个分辨率的概念：(1)是FFT变换频谱中相邻两点之间的频率间隔，通常称为计算分辨率。补零影响的就是计算分辨率，它与信号分辨能力无关；（2）是FFT对信号的分辨能力，通常称为物理分辨率，它只与时域采样序列的有效时长有关。下面对这两种分辨率进行详细说明。
1. 计算分辨率
从上面的分析可知，计算分辨率$\Delta f=f_s/N$，因此任何降低采样率或增加时域序列点数的操作都能提高计算分辨率。补零就是最常见的操作。下面推导补零前后信号频谱的关系。
设原始时域序列为$x[n],0\le n<N$，它的FFT变换为
$$X[k]=\sum _{n=0}^{N-1}x[n]exp(-j\frac{2\pi }{N}kn),0\le n<N,0\le k<N$$
在$x[n]$后面补0使序列总点数增加至$M$，则补零后的序列为
$$x^{\prime }[n]=\{\begin{array}{c}x[n],0\le n<N\\ 0,N\le n<M\end{array}$$
补零后序列的FFT变换为
$$\begin{gathered} X^{\prime }[k]=\sum _{n=0}^{M-1}{x}^{\prime }[n]exp(-j\frac{2\pi }{M}kn)=\sum _{n=0}^{N-1}x[n]exp(-j\frac{2\pi }{M}kn) \\ =\sum _{n=0}^{N-1}[\frac{1}{N}\sum _{l=0}^{N-1}X[l]exp(j\frac{2\pi }{N}ln)]exp(-j\frac{2\pi }{M}kn) \\ =\frac{1}{N}\sum _{l=0}^{N-1}X[l]\sum _{n=0}^{N-1}exp(-j\frac{2\pi }{M}kn)exp(j\frac{2\pi }{N}ln) \\ =\frac{1}{N}\sum _{l=0}^{N-1}X[l]\Theta (k,l) \end{gathered}$$
其中
$$\begin{gathered} \Theta (k,l)=\sum _{n=0}^{N-1}exp(-j\frac{2\pi }{M}kn)exp(j\frac{2\pi }{N}ln)=\sum _{n=0}^{N-1}exp(j2\pi (\frac{l}{N}-\frac{k}{M})n) \\ =exp(j(N-1)\pi (\frac{l}{N}-\frac{k}{M}))\frac{sin(N\pi (\frac{l}{N}-\frac{k}{M}))}{sin(\pi (\frac{l}{N}-\frac{k}{M}))} \end{gathered}$$
当$\frac{l}{N}-\frac{k}{M}=0$时，$\Theta (k,l)=N$，其他情况，$\Theta (k,l)$急剧减小，所以$X^{\prime }[k]=\frac{1}{N}{\sum }_{l=0}^{N-1}X[l]\Theta (k,l)\approx X[\frac{N}{M}k]$。
基于上面的分析，可知补零后的频谱是原始频谱的插值。由于插值是基于已有信息进行操作的，所以若补零前的信号频谱不足以支撑两个信号的分辨，则补零后的频谱同样无法分辨这两个信号。
虽然补零并不能提高FFT对信号的分辨能力，但它具有以下用处：
a) 改变FFT变换的点数，提高FFT变换的计算效率。最常见的用法是将FFT变换的点数补零至$2^N$，此时可用蝶形算法保证FFT变换的效率；
b) 使信号频谱更加精细和平滑，降低栅栏效应；
c) 补零至特定点数凸显特定频率的信号。FFT变换只能表征频率为$f_s/N$整数倍的信号，若原始信号频率不在这些频率点上则可通过补零使其落在FFT可表征的频率点上。
需要注意的是：补零仅仅增加了信号长度，但没有增加信号功率，相同功率平均到更多频点上，会使得频谱的功率谱密度降低。因此若用补零后的频谱来计算某一频率范围内的信号功率时，需要进行信号功率谱密度的补偿。

2. 物理分辨率

物理分辨率表征FFT变换对不同频率信号的分辨能力，它仅与有效信号时长有关。具体地，设信号时长为$T$，则物理分辨率为$\Delta f=1/T$。信号时长可用采样点数和采样间隔的乘积得到，即$T=N/f_s$，所以$\Delta f=f_s/N$。这个表达式与上述计算分辨率表达式的形式一致，差别在于物理分辨率表达式中$N$表示有效信号的采样点数，不包括补零长度。
下面分别从时域和频率两个角度理解物理分辨率的表达式：
(1) 时域和频域本质上是从两种不同角度去看待一个信号，在这两个域上信号所包含的信息量是相等的。因此不论是从时域还是频域角度去理解信号的分辨能力都应得到相同的结论。下面以两个单频信号为例，从时域的角度理物理分辨率。设两个单频信号的频率分别为$f_0$和$f_1$，对应的信号周期分别为$T_0=1/f_0,T_1=1/f_1$。显然，若两个信号的初始位置是对齐的，则经过$\Delta T=1/(|f_0-f_1|)$时间后，两个信号会重新对齐。如下图所示，$f_0=10Hz,f_1=12Hz$，经过0.5s后两个信号的重新回到初始状态。

图2. 两个正弦信号时域波形

将上述两个信号进行叠加，我们的目标是从叠加后的信号中判断原始两个信号的频率差。根据上面的描述，我们知道叠加信号的周期为$\Delta T=1/(|f_0-f_1|)$。所以利用周期进行判断的思路显然是可行的，但是这种方法至少需要两个周期的信号长度。下图是图2中两个正弦信号的叠加波形，时长为一个周期$\Delta T$。以$\Delta T/2$为界将信号分为两段，可见这两段信号是共轭对称的，可以利用这个性质计算$\Delta T$，进而计算两个正弦信号的频差。此时只要求信号长度为$\Delta T$即可。

图3. 图2中两个正弦信号的叠加

(2) 无限长单频信号的傅里叶变换为理想的冲激函数($\delta$函数)，其频谱仅在对应频率位置上非零。然而，我们实际处理的信号是有限长度的，等效于对无限长度的信号加了如下矩形窗
$$w[t]=\{\begin{array}{c}1,0\le t<T\\ 0,others\end{array}$$

上述矩形窗函数的频谱为
$$W(f)={\int }_0^{T}exp(-j2\pi ft)dt=Texp(-j\pi fT)sinc(\pi fT)$$
显然，$|W(f)|=T|sinc(\pi fT)|$，它的示意图如下，在$1/T$的整数倍上存在零点。时域加窗等效于频域卷积，所以正弦信号的频谱表现为以其频率为中心的sinc函数形态。当两个正弦信号的频谱主瓣不发生混叠时，可在频域进行区分。

图4. sinc函数包络示意图，在1/T的整数倍上存在包络零点

3. 仿真验证

如下图所示，仿真生成两个正弦信号，观察在有效样点数不够、有效样点数足够及有效样点数不够但通过补零增加样点数目这三种情况下的频谱，各个参数设置及仿真代码见附录。下图表明，补零并不能提高FFT对信号的分辨能力，只能起到对频谱进行插值和平滑的作用，与前面的理论推导结论相符。

图5. 仿真不同样点条件下FFT的分辨能力

---

【附录1】：图1代码

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']  # 使用微软雅黑字体
plt.rcParams['axes.unicode_minus'] = False  # 处理负号显示异常

if __name__ == '__main__':
    f0 = 1000
    f1 = 24000
    fs = 10000
    N = 2000
    t = np.arange(0, N, 1) / fs
    sig = np.sin(2 * np.pi * f0 * t) + np.sin(2 * np.pi * f1 * t)
    plt.figure()
    plt.plot(np.linspace(-fs / 2, fs / 2, N), 20 * np.log10(np.fft.fftshift(np.abs(np.fft.fft(sig)))), 'b.-')
    plt.show()

【附录2】：仿真验证代码

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']  # 使用微软雅黑字体
plt.rcParams['axes.unicode_minus'] = False  # 处理负号显示异常

if __name__ == '__main__':
    f0 = 200
    f1 = 250
    fs = 100000
    df = f1 - f0
    data_time = 1 / df
    N = int(data_time * fs)
    # ##############样点数不够的情况###############
    t1 = np.arange(0, N // 2, 1) / fs
    sig_1 = np.sin(2 * np.pi * f0 * t1) + np.sin(2 * np.pi * f1 * t1)

    # ###############样点数够的情况###############
    t2 = np.arange(0, 2 * N - 1, 1) / fs
    sig_2 = np.sin(2 * np.pi * f0 * t2) + np.sin(2 * np.pi * f1 * t2)

    # ##############对样点数不够的情况进行补零###############
    sig_3 = np.zeros(10 * N)
    sig_3[0:len(sig_1)] = sig_1

    plt.figure()
    plt.subplot(3, 1, 1)
    plt.plot(np.linspace(-fs / 2, fs / 2, len(sig_1)), 20 * np.log10(np.fft.fftshift(np.abs(np.fft.fft(sig_1)))), 'b.-')
    plt.xlim([-500, 500])
    plt.title('样点数不够')
    plt.subplot(3, 1, 2)
    plt.plot(np.linspace(-fs / 2, fs / 2, len(sig_2)), 20 * np.log10(np.fft.fftshift(np.abs(np.fft.fft(sig_2)))), 'b.-')
    plt.xlim([-500, 500])
    plt.title('样点数足够')
    plt.subplot(3, 1, 3)
    plt.plot(np.linspace(-fs / 2, fs / 2, len(sig_3)), 20 * np.log10(np.fft.fftshift(np.abs(np.fft.fft(sig_3)))), 'b.-')
    plt.xlim([-500, 500])
    plt.title('样点数不够，补零')
    plt.xlabel('Frequency(Hz)')
    plt.show()