科学计算 | 一个怪现象引起的思考

ZoroGH

引言

一般在研究信号时，我们不会过多研究相位谱，一方面，相位谱提供的信息比较有限，另一方面，相位谱不够直观。

但是，这学期在做一个项目时，遇到了一个奇怪的现象问题，具体如下：

设有 $1Hz$ 正弦信号 $sin(2\pi t)$ ，系统的采样率为 $20Hz$ ，根据采样定理，采样得到的连续信号能完整表征出原信号的信息。

现在有两种记录方式，由于原信号周期为1秒，即系统只需要采样20个点才能表征完整的一个周期的信号，第一种记录方式为采样20个点来记录该信号，第二种记录21个点来记录该信号。后续将第一种记作整周期截断，第二种为非整周期截断。现象是，整周期截断的信号其相位谱乱七八糟，而非整周期截断的相位谱却是有些许规律。这问题当时让我迷糊了一晚上，于是探究了一番。

仿真

先来看图，以免看公式会晕😵😵。

整周期截断

在Python中进行如下仿真：

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft

# 定义信号
def f(x): return np.sin(2*np.pi*x)  # 1Hz正弦波


def signal_analysis(t: np.ndarray, sig: np.ndarray) -> None:
    '''分析(数字)信号频谱并绘制图像'''
    assert len(t) == len(sig), "序列长度不匹配"
    N = len(t)
    X = fft(sig)
    A = np.abs(X)
    phi = np.angle(X)
    k_arr = np.arange(N)
    plt.figure(tight_layout=True)
    plt.subplot(211)
    plt.stem(t, sig)
    plt.xticks(t)
    plt.grid(True, alpha=0.5)
    plt.subplot(223)
    plt.stem(k_arr, A)
    plt.xticks(k_arr)
    plt.grid(True, alpha=0.5)
    plt.subplot(224)
    plt.stem(k_arr, phi)
    plt.xticks(k_arr)
    plt.grid(True, alpha=0.5)
    plt.show()


fs = 20  # 20Hz采样率
dt = 1/fs  # 采样间隔
N = 20  # 采样点数
n_arr = np.arange(N)  # 从0开始到N-1，总共N个点
t_arr1 = np.arange(0, N*dt, step=dt)  # 采样时间序列
sig1 = f(t_arr1)  # 对信号进行采样

signal_analysis(np.arange(10), sig1)

运行得出如下结果

根据图像看得出来，采样率为 $20Hz$ ，有效谱线为 $0\sim9$ ，可以发现幅度谱在抽样点 $k=1$ 时有值，即对应的模拟频率为 $f_s\times\frac{n}{N}|_{n=1}=1Hz$ ，符合实际，然而相位谱却骚得离谱，好像不属于这个世界。

非整周期截断

将上述仿真代码中37行以后改为

t_arr2 = np.arange(0, (N+1)*dt, step=dt)  # 采样时间序列
sig2 = f(t_arr2)  # 对信号进行采样
signal_analysis(t_arr2, sig2)

得到如下图像

信号的幅度谱发生了频谱泄露，合乎情理，但相比整周期截断而言，这个非正周期截断的相位谱却看起来有那么一丢丢规律，有那么几分合理，好像在说这个才是整周期截断的相位谱。

理论推导

这里推导整周期截断信号的 $DFT$ 。

记原连续信号为 $s$ ，采样率为 $f_s$ ，采样间隔 $t_0=1/f_s$ ，则信号序列为

x(n)=s(t)|_{t=n*t_0}=sin(2\pi t_0 n)\quad(n=0,1,\cdots,19)

根据 $DFT$ 得

X(k) = \sum_{n=0}^{19}x(n)W_N^{nk}=\sum_{n=0}^{19}sin(2\pi t_0 n)e^{-j\frac{2\pi}{N}nk}\quad(k=0,1,\cdots,19)

计算得到（详细过程见附录）

X(k)= \begin{cases} -10j,\quad k=1 \\ 10j, \quad k=19 \\ 0,\quad others \end{cases}

Analysis

从理论推导可以得出，该信号的离散傅里叶变换只有在 $k=1,19$ 处有值。理论上相位谱也应该在 $k=1,19$ 两个抽样点处有值才对。可仿真事实表明并不是如此，整周期截断信号的相位谱乱七八糟😫。如果增加点数，增加采样率，事情会变得更奇怪，如下二图所示：

采样频率为40Hz

采样点数为40

因此可以认为，这并不是采样率不够与采样点数少的问题（虽然从数学上看，必然不是他们的问题😁）。所以，问题究竟出现在哪？难道是Bug？

Debug

将 $FFT$ 运算得到的结果打印出来，结果如下

array([-1.66533454e-16-0.00000000e+00j,  5.55111512e-17-1.00000000e+01j, 
       -6.65500331e-16+6.24143231e-16j, -9.29002411e-17-1.02444732e-15j, 
        5.55111512e-17+0.00000000e+00j, -1.11022302e-16-1.05471187e-15j, 
        9.43689571e-16+0.00000000e+00j,  4.89045368e-16+4.23316153e-16j, 
        4.91399809e-16+2.95549861e-16j, -8.32667268e-16+8.88178420e-16j, 
       -4.99600361e-16-0.00000000e+00j, -8.32667268e-16-8.88178420e-16j, 
        4.91399809e-16-2.95549861e-16j,  4.89045368e-16-4.23316153e-16j, 
        9.43689571e-16-0.00000000e+00j, -1.11022302e-16+1.05471187e-15j, 
        5.55111512e-17-0.00000000e+00j, -9.29002411e-17+1.02444732e-15j, 
       -6.65500331e-16-6.24143231e-16j,  5.55111512e-17+1.00000000e+01j])

所以啊，就是这么个玩意让图像变丑。

从理论推导来看，除了k=1,19处应该有值外，其他地方应该均是0，但是在计算机运算时，结果不严格等于0，而是一个极小的数字。

可以认为，计算机在计算和存储浮点数的时候，会损失精度并增加一定的误差，这里其实和有限字长效应很相似，即计算机无法精确的表示某一些数字。

误差流图.png

虽然不影响幅度的表示，但是由于函数np.angle是求复数的幅角的，它不会在乎这个复数的模多小，只会在乎虚部与实部之比，因此，即是计算得到的复数模值很小，但虚部实部相当时，误差通过角度函数便会放大，以至于影响后续判断。另外，从数学角度来看，一个复数的模值为0时，其幅角是无意义的。

Solution

因此，为了解决由于计算机精度而产生的问题，需要后续数据进行修正。于是不妨取一个极小数，对计算结果进行阈值过滤（借用图像处理的术语是：阈值取零），然后再进行图像的绘制。

修正后的代码如下（只修改signal_analysis函数）

def signal_analysis(t: np.ndarray, sig: np.ndarray,epsilon = 1e-4) -> None:
    ......
    X = fft(sig)
    A = np.abs(X)
    mask = A>epsilon
    A = A*mask
    phi = np.angle(X) * mask
    ......

运行脚本得到如下结果。

可以看到，这图像正好就是我们所要的结果。

附录

X(k) = \sum_{n=0}^{19}x(n)W_N^{nk}=\sum_{n=0}^{19}sin(2\pi t_0 n)e^{-j\frac{2\pi}{N}nk}

根据欧拉公式有

sin(2\pi t_0n)=\frac{1}{2j}\left(e^{j2\pi t_0n}- e^{-j2\pi t_0n}\right)

于是可得

\begin{aligned} X(k) &= \sum_{n=0}^{19}sin(2\pi t_0 n)e^{-j\frac{2\pi}{N}nk} \\ &=\frac{1}{2j}\sum_{n=0}^{19}\left(e^{j2\pi t_0n}- e^{-j2\pi t_0n}\right)e^{-j\frac{2\pi}{N}nk}\\ &=\frac{1}{2j}\sum_{n=0}^{19}\left(e^{j2\pi n (t_0-\frac k N)}- e^{-j2\pi n (t_0+\frac k N)}\right) \end{aligned}

当 $e^{j2\pi n (t_0-\frac k N)}$ 或 $e^{-j2\pi n (t_0+\frac k N)}$ 为1时不可直接进行等比求和，又由于 $e^{j2\pi n} =1$ (n为整数)，因此

t_0-\frac k N=\frac{1-k}{20} \notin Z \rightarrow k\ne 1\\ t_0+\frac k N=\frac{1+k}{20} \notin Z \rightarrow k\ne 19

于是不难得到

\begin{aligned} X(k)&=\frac{1}{2j}\sum_{n=0}^{N-1=19}\left(e^{j2\pi n (t_0-\frac k N)}- e^{-j2\pi n (t_0+\frac k N)}\right)\\ &=\frac{1}{2j}\left( \frac{e^{j2\pi N (t_0-\frac k N)}-1}{\cdots} -\frac{e^{-j2\pi N (t_0+\frac k N)}-1}{\cdots} \right)\\ &e^{j2\pi n}=1\\ &= 0\quad(k\ne 1,19)\\ X(1)&=-10j\\ X(19)&= 10j \end{aligned}

至于，非整周期截断的那个公式呀，我懒得推了，但为什么其表现出了一定的规律性，是因为...，算了我还是推吧，实际上就是把求和上限改为了20项而不是之前的19，于是也不能化简了

最后化简等比求和的结果会用到一个小技巧

\frac{e^{j2\theta N}-1}{e^{j2\theta}-1}=\frac{e^{j\theta N}}{e^{j\theta}}\times\frac{e^{j\theta N}-e^{-j\theta N}}{e^{j\theta}-e^{-j\theta}}=e^{j\theta(N-1)}\frac{\sin(\theta N)}{\sin(\theta)}

最后结果不推了，过于丑陋（好吧是我懒😴😴）。