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

133 阅读5分钟

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

ZoroGH

引言

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

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

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

现在有两种记录方式,由于原信号周期为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)

运行得出如下结果

Figure_1.png

根据图像看得出来,采样率为20Hz20Hz,有效谱线为090\sim9,可以发现幅度谱在抽样点k=1k=1时有值,即对应的模拟频率为fs×nNn=1=1Hzf_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)

得到如下图像

Figure_2.png

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

理论推导

这里推导整周期截断信号的DFTDFT

记原连续信号为ss,采样率为fsf_s,采样间隔t0=1/fst_0=1/f_s,则信号序列为

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

根据DFTDFT

X(k)=n=019x(n)WNnk=n=019sin(2πt0n)ej2πNnk(k=0,1,,19)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)={10j,k=110j,k=190,othersX(k)= \begin{cases} -10j,\quad k=1 \\ 10j, \quad k=19 \\ 0,\quad others \end{cases}

Analysis

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

Figure_1_N40_fs40.png

采样频率为40Hz

Figure_1_N40.png

采样点数为40

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

Debug

FFTFFT运算得到的结果打印出来,结果如下

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
    ......

运行脚本得到如下结果。

Figure_3.png 可以看到,这图像正好就是我们所要的结果。

附录

X(k)=n=019x(n)WNnk=n=019sin(2πt0n)ej2πNnkX(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πt0n)=12j(ej2πt0nej2πt0n)sin(2\pi t_0n)=\frac{1}{2j}\left(e^{j2\pi t_0n}- e^{-j2\pi t_0n}\right)

于是可得

X(k)=n=019sin(2πt0n)ej2πNnk=12jn=019(ej2πt0nej2πt0n)ej2πNnk=12jn=019(ej2πn(t0kN)ej2πn(t0+kN))\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}

ej2πn(t0kN)e^{j2\pi n (t_0-\frac k N)}ej2πn(t0+kN)e^{-j2\pi n (t_0+\frac k N)}为1时不可直接进行等比求和,又由于ej2πn=1e^{j2\pi n} =1(n为整数),因此

t0kN=1k20Zk1t0+kN=1+k20Zk19t_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

于是不难得到

X(k)=12jn=0N1=19(ej2πn(t0kN)ej2πn(t0+kN))=12j(ej2πN(t0kN)1ej2πN(t0+kN)1)ej2πn=1=0(k1,19)X(1)=10jX(19)=10j\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,于是也不能化简了

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

ej2θN1ej2θ1=ejθNejθ×ejθNejθNejθejθ=ejθ(N1)sin(θN)sin(θ)\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)}

最后结果不推了,过于丑陋(好吧是我懒😴😴)。