从 FFT 中找出信号的周期

Question

我有一个周期信号，我想找到周期。

由于有边界效应，我先剪掉边界，看第一个和最后一个最小值，保留N个周期。

然后，我计算 FFT。

代码：

import numpy as np

from matplotlib import pyplot as plt

# The list of a periodic something
L = [2.762, 2.762, 1.508, 2.758, 2.765, 2.765, 2.761, 1.507, 2.757, 2.757, 2.764, 2.764, 1.512, 2.76, 2.766, 2.766, 2.763, 1.51, 2.759, 2.759, 2.765, 2.765, 1.514, 2.761, 2.758, 2.758, 2.764, 1.513, 2.76, 2.76, 2.757, 2.757, 1.508, 2.763, 2.759, 2.759, 2.766, 1.517, 4.012]
# Round because there is a slight variation around actually equals values: 2.762, 2.761 or 1.508, 1.507
L = [round(elt, 1) for elt in L]

minima = min(L)
min_id = L.index(minima)

start = L.index(minima)
stop = L[::-1].index(minima)

L = L[start:len(L)-stop]

fft = np.fft.fft(np.asarray(L))/len(L)
fft = fft[range(int(len(L)/2))]

plt.plot(abs(fft))

我知道我在列表的 2 个点之间有多少时间（即采样频率，在本例中为 190 Hz）。我认为 fft 应该给我一个与周期中点数相对应的值的尖峰，从而给我点数和周期数。然而，这根本不是我观察到的输出：

我目前的猜测是 0 处的尖峰对应于我的信号的平均值，并且这个 7 左右的小尖峰应该是我的时期（尽管重复模式只包括 5 个点）。

我究竟做错了什么？谢谢！

Answer 1

您的数据是正确的，只是您没有正确对其进行预处理：

1. 第一个巨峰是信号的 DC/平均值。如果你在 DFT 之前减去它，它就会消失
2. 在采用 DFT 之前不加窗信号会在 DFT 频谱中产生振铃，降低峰值并提高“非峰值”。

如果包含这两个步骤，结果应该更符合您的预期：

import numpy as np
import scipy.signal

from matplotlib import pyplot as plt

L = np.array([2.762, 2.762, 1.508, 2.758, 2.765, 2.765, 2.761, 1.507, 2.757, 2.757, 2.764, 2.764, 1.512, 2.76, 2.766, 2.766, 2.763, 1.51, 2.759, 2.759, 2.765, 2.765, 1.514, 2.761, 2.758, 2.758, 2.764, 1.513, 2.76, 2.76, 2.757, 2.757, 1.508, 2.763, 2.759, 2.759, 2.766, 1.517, 4.012])
L = np.round(L, 1)
# Remove DC component
L -= np.mean(L)
# Window signal
L *= scipy.signal.windows.hann(len(L))

fft = np.fft.rfft(L, norm="ortho")

plt.plot(L)
plt.figure()
plt.plot(abs(fft))

您会注意到，您会在 附近看到一个峰值8，而在两倍处看到另一个峰值，16。这也是意料之中的：周期信号在n*period采样后始终是周期性的，其中 n 是任何自然数。在你的情况下：n*8。

Answer 2

一旦信号的直流部分被去除，该函数就可以与其自身进行卷积以捕获周期。事实上，卷积将在每个周期倍数处出现峰值。FFT 可用于计算卷积。

fft = np.fft.rfft(L, norm="ortho")

def abs2(x):
    return x.real**2 + x.imag**2

selfconvol=np.fft.irfft(abs2(fft), norm="ortho")

第一个输出不太好，因为图像的大小不是周期的倍数。

正如 Nils Werner 所指出的，可以应用窗口来限制频谱泄漏的影响。作为替代方案，可以使用周期的第一个粗略估计来中继信号，并且可以按照我在如何缩放基于 FFT 的互相关使其峰值等于 Pearson's rho 中的回答重复该过程。

从这里开始，获取周期归结为找到第一个最大值。这是一种可以完成的方法：

import numpy as np
import scipy.signal

from matplotlib import pyplot as plt

L = np.array([2.762, 2.762, 1.508, 2.758, 2.765, 2.765, 2.761, 1.507, 2.757, 2.757, 2.764, 2.764, 1.512, 2.76, 2.766, 2.766, 2.763, 1.51, 2.759, 2.759, 2.765, 2.765, 1.514, 2.761, 2.758, 2.758, 2.764, 1.513, 2.76, 2.76, 2.757, 2.757, 1.508, 2.763, 2.759, 2.759, 2.766, 1.517, 4.012])
L = np.round(L, 1)
# Remove DC component, as proposed by Nils Werner
L -= np.mean(L)
# Window signal
#L *= scipy.signal.windows.hann(len(L))

fft = np.fft.rfft(L, norm="ortho")

def abs2(x):
    return x.real**2 + x.imag**2

selfconvol=np.fft.irfft(abs2(fft), norm="ortho")
selfconvol=selfconvol/selfconvol[0]

plt.figure()
plt.plot(selfconvol)
plt.savefig('first.jpg')
plt.show()


# let's get a max, assuming a least 4 periods...
multipleofperiod=np.argmax(selfconvol[1:len(L)/4])
Ltrunk=L[0:(len(L)//multipleofperiod)*multipleofperiod]

fft = np.fft.rfft(Ltrunk, norm="ortho")
selfconvol=np.fft.irfft(abs2(fft), norm="ortho")
selfconvol=selfconvol/selfconvol[0]

plt.figure()
plt.plot(selfconvol)
plt.savefig('second.jpg')
plt.show()


#get ranges for first min, second max
fmax=np.max(selfconvol[1:len(Ltrunk)/4])
fmin=np.min(selfconvol[1:len(Ltrunk)/4])
xstartmin=1
while selfconvol[xstartmin]>fmin+0.2*(fmax-fmin) and xstartmin< len(Ltrunk)//4:
    xstartmin=xstartmin+1

xstartmax=xstartmin
while selfconvol[xstartmax]<fmin+0.7*(fmax-fmin) and xstartmax< len(Ltrunk)//4:
    xstartmax=xstartmax+1

xstartmin=xstartmax
while selfconvol[xstartmin]>fmin+0.2*(fmax-fmin) and xstartmin< len(Ltrunk)//4:
    xstartmin=xstartmin+1

period=np.argmax(selfconvol[xstartmax:xstartmin])+xstartmax

print "The period is ",period