如何找到信号周期（自相关与快速傅立叶变换与功率谱密度之间的关系？

Question

假设人们想找到给定正弦波信号的周期。从我在网上阅读的内容来看，似乎两种主要方法采用傅立叶分析或自相关。我正在尝试使用python自动化该过程，而我的用例是将此概念应用于类似信号，这些信号来自围绕恒星运行的模拟物体的位置（或速度或加速度）的时间序列。

对于简单示例而言，请考虑将x = sin(t)替换为0 ≤ t ≤ 10 pi。

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

## sample data
t = np.linspace(0, 10 * np.pi, 100)
x = np.sin(t)
fig, ax = plt.subplots()
ax.plot(t, x, color='b', marker='o')
ax.grid(color='k', alpha=0.3, linestyle=':')
plt.show()
plt.close(fig)

给定形式为x = a sin(b(t+c)) + d的正弦波，正弦波的周期为2 * pi / b。由于b=1（或通过目视检查），我们的正弦波的周期为2 * pi。我可以对照此基准检查从其他方法获得的结果。

尝试1：自相关

据我所知（如果我错了，请纠正我），可以使用相关性来查看一个信号是否是另一个信号的时滞副本（类似于余弦和正弦相差多少）。因此，自相关就是针对自身测试信号，以测量时滞重复该信号的时间。使用example posted here：

result = np.correlate(x, x, mode='full')

由于x和t每个都由100元素组成，result由199元素组成，所以我不确定为什么要任意选择最后一个100元素。

print("\n autocorrelation (shape={}):\n{}\n".format(result.shape, result))  

 autocorrelation (shape=(199,)):
[ 0.00000000e+00 -3.82130761e-16 -9.73648712e-02 -3.70014208e-01
 -8.59889695e-01 -1.56185995e+00 -2.41986054e+00 -3.33109112e+00
 -4.15799070e+00 -4.74662427e+00 -4.94918053e+00 -4.64762251e+00
 -3.77524157e+00 -2.33298717e+00 -3.97976240e-01  1.87752669e+00
  4.27722402e+00  6.54129270e+00  8.39434617e+00  9.57785701e+00
  9.88331103e+00  9.18204933e+00  7.44791758e+00  4.76948221e+00
  1.34963425e+00 -2.50822289e+00 -6.42666652e+00 -9.99116299e+00
 -1.27937834e+01 -1.44791297e+01 -1.47873668e+01 -1.35893098e+01
 -1.09091510e+01 -6.93157447e+00 -1.99159756e+00  3.45267493e+00
  8.86228186e+00  1.36707567e+01  1.73433176e+01  1.94357232e+01
  1.96463736e+01  1.78556800e+01  1.41478477e+01  8.81191526e+00
  2.32100171e+00 -4.70897483e+00 -1.15775811e+01 -1.75696560e+01
 -2.20296487e+01 -2.44327920e+01 -2.44454330e+01 -2.19677060e+01
 -1.71533510e+01 -1.04037163e+01 -2.33560966e+00  6.27458308e+00
  1.45655029e+01  2.16769872e+01  2.68391837e+01  2.94553896e+01
  2.91697473e+01  2.59122266e+01  1.99154591e+01  1.17007613e+01
  2.03381596e+00 -8.14633251e+00 -1.78184255e+01 -2.59814393e+01
 -3.17580589e+01 -3.44884934e+01 -3.38046447e+01 -2.96763956e+01
 -2.24244433e+01 -1.26974172e+01 -1.41464998e+00  1.03204331e+01
  2.13281784e+01  3.04712823e+01  3.67721634e+01  3.95170295e+01
  3.83356037e+01  3.32477037e+01  2.46710643e+01  1.33886439e+01
  4.77778141e-01 -1.27924775e+01 -2.50860560e+01 -3.51343866e+01
 -4.18671622e+01 -4.45258983e+01 -4.27482779e+01 -3.66140001e+01
 -2.66465884e+01 -1.37700036e+01  7.76494745e-01  1.55574483e+01
  2.90828312e+01  3.99582426e+01  4.70285203e+01  4.95000000e+01
  4.70285203e+01  3.99582426e+01  2.90828312e+01  1.55574483e+01
  7.76494745e-01 -1.37700036e+01 -2.66465884e+01 -3.66140001e+01
 -4.27482779e+01 -4.45258983e+01 -4.18671622e+01 -3.51343866e+01
 -2.50860560e+01 -1.27924775e+01  4.77778141e-01  1.33886439e+01
  2.46710643e+01  3.32477037e+01  3.83356037e+01  3.95170295e+01
  3.67721634e+01  3.04712823e+01  2.13281784e+01  1.03204331e+01
 -1.41464998e+00 -1.26974172e+01 -2.24244433e+01 -2.96763956e+01
 -3.38046447e+01 -3.44884934e+01 -3.17580589e+01 -2.59814393e+01
 -1.78184255e+01 -8.14633251e+00  2.03381596e+00  1.17007613e+01
  1.99154591e+01  2.59122266e+01  2.91697473e+01  2.94553896e+01
  2.68391837e+01  2.16769872e+01  1.45655029e+01  6.27458308e+00
 -2.33560966e+00 -1.04037163e+01 -1.71533510e+01 -2.19677060e+01
 -2.44454330e+01 -2.44327920e+01 -2.20296487e+01 -1.75696560e+01
 -1.15775811e+01 -4.70897483e+00  2.32100171e+00  8.81191526e+00
  1.41478477e+01  1.78556800e+01  1.96463736e+01  1.94357232e+01
  1.73433176e+01  1.36707567e+01  8.86228186e+00  3.45267493e+00
 -1.99159756e+00 -6.93157447e+00 -1.09091510e+01 -1.35893098e+01
 -1.47873668e+01 -1.44791297e+01 -1.27937834e+01 -9.99116299e+00
 -6.42666652e+00 -2.50822289e+00  1.34963425e+00  4.76948221e+00
  7.44791758e+00  9.18204933e+00  9.88331103e+00  9.57785701e+00
  8.39434617e+00  6.54129270e+00  4.27722402e+00  1.87752669e+00
 -3.97976240e-01 -2.33298717e+00 -3.77524157e+00 -4.64762251e+00
 -4.94918053e+00 -4.74662427e+00 -4.15799070e+00 -3.33109112e+00
 -2.41986054e+00 -1.56185995e+00 -8.59889695e-01 -3.70014208e-01
 -9.73648712e-02 -3.82130761e-16  0.00000000e+00]

尝试2：傅立叶

由于我不确定从上一次尝试到哪里去，我寻求了新的尝试。据我了解，傅立叶分析基本上是将信号从时域（x(t) vs t）移至/从频域（x(t) vs f=1/t）移出；频率空间中的信号应显示为正弦波，随时间衰减。该周期是从观察到最多的频率获得的，因为这是频率分布峰值的位置。

由于我的值都是实值，因此应用傅里叶变换应该意味着我的输出值都是复数值。除了scipy has methods for real-values的事实，我认为这不是问题。我不完全了解所有不同scipy方法之间的差异。这使我很难遵循this posted solution中提出的算法（即，如何/为什么选择阈值？）。

omega = np.fft.fft(x)
freq = np.fft.fftfreq(x.size, 1)
threshold = 0
idx = np.where(abs(omega)>threshold)[0][-1]
max_f = abs(freq[idx])
print(max_f)

此输出0.01，表示周期为1/0.01 = 100。这也没有道理。

尝试3：功率谱密度

根据scipy docs，我应该能够使用周期图（其according to wikipedia是自相关函数的傅立叶变换）来估计信号的功率谱密度（psd）。通过选择信号达到峰值的主频fmax，可以将信号的周期设为1 / fmax。

freq, pdensity = signal.periodogram(x)

fig, ax = plt.subplots()
ax.plot(freq, pdensity, color='r')
ax.grid(color='k', alpha=0.3, linestyle=':')
plt.show()
plt.close(fig)

下面显示的周期图以49.076...的频率在fmax = 0.05达到峰值。因此，period = 1/fmax = 20。这对我来说没有意义。我觉得这与采样率有关，但了解程度不足以确认或进一步发展。

我意识到我在理解这些东西的工作方式方面缺少一些根本的空白。在线上有很多资源，但是很难在大海捞针上找到这根针。有人可以帮助我了解更多吗？

Answer 1

让我们首先看您的信号（我添加了endpoint=False以使除法均匀：]

t = np.linspace(0, 10*np.pi, 100, endpoint=False)
x = np.sin(t)

让弧度除掉（基本上是取t /= 2*np.pi并通过与频率相关来创建相同的信号：

fs = 20 # Sampling rate of 100/5 = 20 (e.g. Hz)
f = 1 # Signal frequency of 1 (e.g. Hz)
t = np.linspace(0, 5, 5*fs, endpoint=False)
x = np.sin(2*np.pi*f*t)

[这使f/fs == 1/20 == 0.05更加突出（即信号的周期恰好是20个采样）。正如您已经猜到的，数字信号的频率始终与其采样率有关。请注意，无论f和fs的值是什么，只要它们的比率相同，实际信号就完全相同：

fs = 1 # Natural units
f = 0.05
t = np.linspace(0, 100, 100*fs, endpoint=False)
x = np.sin(2*np.pi*f*t)

以下，我将使用这些自然单位（fs = 1）。唯一的不同是在t中，因此是所产生的频率轴。

自相关

您对自相关函数的功能的理解是正确的。它检测信号与自身时滞形式的相关性。它通过在其自身上滑动信号来实现此目的，如此处右列所示（来自Wikipedia）：

请注意，由于相关函数的两个输入都相同，结果信号必定是对称的。这就是np.correlate的输出通常从中间切成薄片的原因：

acf = np.correlate(x, x, 'full')[-len(x):]

现在索引0对应于信号的两个副本之间的0延迟。

接下来，您将要找到呈现最大相关性的索引或延迟。由于重叠部分的缩小，默认情况下它也将是索引0，因此以下内容将不起作用：

acf.argmax() # Always returns 0

相反，我建议找到最大的peak，其中将峰定义为比两个直接邻居都大的索引值。

inflection = np.diff(np.sign(np.diff(acf))) # Find the second-order differences
peaks = (inflection < 0).nonzero()[0] + 1 # Find where they are negative
delay = peaks[acf[peaks].argmax()] # Of those, find the index with the maximum value

现在delay == 20，它告诉您信号的采样率是1/20的频率：

signal_freq = fs/delay # Gives 0.05

傅立叶变换

您使用以下方法计算FFT：

omega = np.fft.fft(x)
freq = np.fft.fftfreq(x.size, 1)

这些功能已重新设计用于复数值信号。它们将适用于实值信号，但是您将获得对称输出，因为负频率分量将与正频率分量相同。 NumPy为实值信号提供单独的功能：

ft = np.fft.rfft(x)
freqs = np.fft.rfftfreq(len(x), t[1]-t[0]) # Get frequency axis from the time axis
mags = abs(ft) # We don't care about the phase information here

让我们看看：

plt.plot(freqs, mags)
plt.show()

注意两点：峰值在频率0.05处，并且轴上的最大频率为0.5（Nyquist frequency，恰好是采样率的一半）。如果我们选择了fs = 20，则为10。

现在让我们找到最大值。您尝试过的阈值方法可以工作，但是盲目选择了目标频率仓，因此在存在其他信号的情况下，该方法会受到影响。我们可以选择最大值：signal_freq = freqs[mags.argmax()] # Gives 0.05

但是，例如，如果我们有较大的DC偏移（因此索引0中的分量较大），则此操作将失败。在这种情况下，我们可以再次选择最高峰，以使其更加稳健：

inflection = np.diff(np.sign(np.diff(mags)))
peaks = (inflection < 0).nonzero()[0] + 1
peak = peaks[mags[peaks].argmax()]
signal_freq = freqs[peak] # Gives 0.05

如果我们选择了fs = 20，则由于生成频率轴的时间轴不同，因此本来会给出signal_freq == 1.0。

周期图

这里的方法本质上是相同的。 x的自相关函数具有与x相同的时间轴和周期，因此我们可以使用上述FFT来找到信号频率：

pdg = np.fft.rfft(acf)
freqs = np.fft.rfftfreq(len(x), t[1]-t[0])

plt.plot(freqs, abs(pdg))
plt.show()

[该曲线显然与x上的直接FFT的特性略有不同，但是主要要点是相同的：频率轴的范围从0到0.5*fs，并且我们发现一个峰值与之前：freqs[abs(pdg).argmax()] == 0.05。