科学计算 | 一个怪现象引起的思考
引言
一般在研究信号时,我们不会过多研究相位谱,一方面,相位谱提供的信息比较有限,另一方面,相位谱不够直观。
但是,这学期在做一个项目时,遇到了一个奇怪的现象问题,具体如下:
设有正弦信号,系统的采样率为,根据采样定理,采样得到的连续信号能完整表征出原信号的信息。
现在有两种记录方式,由于原信号周期为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)
运行得出如下结果
根据图像看得出来,采样率为,有效谱线为,可以发现幅度谱在抽样点时有值,即对应的模拟频率为,符合实际,然而相位谱却骚得离谱,好像不属于这个世界。
非整周期截断
将上述仿真代码中37行以后改为
t_arr2 = np.arange(0, (N+1)*dt, step=dt) # 采样时间序列
sig2 = f(t_arr2) # 对信号进行采样
signal_analysis(t_arr2, sig2)
得到如下图像
信号的幅度谱发生了频谱泄露,合乎情理,但相比整周期截断而言,这个非正周期截断的相位谱却看起来有那么一丢丢规律,有那么几分合理,好像在说这个才是整周期截断的相位谱。
理论推导
这里推导整周期截断信号的。
记原连续信号为,采样率为,采样间隔,则信号序列为
根据得
计算得到(详细过程见附录)
Analysis
从理论推导可以得出,该信号的离散傅里叶变换只有在处有值。理论上相位谱也应该在两个抽样点处有值才对。可仿真事实表明并不是如此,整周期截断信号的相位谱乱七八糟😫。如果增加点数,增加采样率,事情会变得更奇怪,如下二图所示:
采样频率为40Hz 采样点数为40因此可以认为,这并不是采样率不够与采样点数少的问题(虽然从数学上看,必然不是他们的问题😁)。所以,问题究竟出现在哪?难道是Bug?
Debug
将运算得到的结果打印出来,结果如下
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,而是一个极小的数字。
可以认为,计算机在计算和存储浮点数的时候,会损失精度并增加一定的误差,这里其实和有限字长效应很相似,即计算机无法精确的表示某一些数字。
虽然不影响幅度的表示,但是由于函数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
......
运行脚本得到如下结果。
可以看到,这图像正好就是我们所要的结果。
附录
根据欧拉公式有
于是可得
当或为1时不可直接进行等比求和,又由于(n为整数),因此
于是不难得到
至于,非整周期截断的那个公式呀,我懒得推了,但为什么其表现出了一定的规律性,是因为...,算了我还是推吧,实际上就是把求和上限改为了20项而不是之前的19,于是也不能化简了
最后化简等比求和的结果会用到一个小技巧
最后结果不推了,过于丑陋(好吧是我懒😴😴)。