让我们首先看看你的信号(我已经添加了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
and 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
以及生成的频率轴。
自相关
您对自相关函数作用的理解是正确的。它检测信号与其自身的时滞版本的相关性。它通过将信号滑过自身来实现这一点,如右列所示(来自维基百科):
请注意,由于相关函数的两个输入相同,所得信号必然是对称的。这就是为什么输出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
Now 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(奈奎斯特频率,恰好是采样率的一半)。如果我们选择了fs = 20
,这将是 10。
现在让我们找出最大值。您尝试过的阈值方法can工作,但目标频率仓是盲目选择的,因此该方法在存在其他信号的情况下会受到影响。我们可以只选择最大值:
signal_freq = freqs[mags.argmax()] # Gives 0.05
然而,如果例如我们有一个大的直流偏移(因此索引 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()
这条曲线显然与直接 FFT 的特性略有不同x
,但主要要点是相同的:频率轴的范围是0
to 0.5*fs
,我们在与之前相同的信号频率处找到一个峰值:freqs[abs(pdg).argmax()] == 0.05
.
Edit:
测量实际的周期性np.sin
,我们可以只使用我们传递给的“角度轴”np.sin
生成频率轴时代替时间轴:
freqs = np.fft.rfftfreq(len(x), 2*np.pi*f*(t[1]-t[0]))
rad_period = 1/freqs[mags.argmax()] # 6.283185307179586
虽然这看起来毫无意义,对吧?我们传入2*np.pi
我们得到2*np.pi
。然而,我们可以对任何常规时间轴做同样的事情,而不需要预先假设pi
在任何一点:
fs = 10
t = np.arange(1000)/fs
x = np.sin(t)
rad_period = 1/np.fft.rfftfreq(len(x), 1/fs)[abs(np.fft.rfft(x)).argmax()] # 6.25
当然,现在真正的价值位于两个箱子之间。这就是插值的用武之地以及选择合适的窗口函数的相关需求。