估计两个时间序列之间的小时间偏移

Question

我有两个时间序列,我怀疑他们之间有时间转换,我想估计这个时间转换.

之前已经问过这个问题: 找到两个(非谐波)波之间的相位差并找到两个相似波形之间的时间偏移,但在我的情况下,时移小于数据的分辨率.例如,数据以小时分辨率提供,时移仅为几分钟(见图).

造成这种情况的原因是,用于测量其中一个系列的数据记录器在其时间上几分钟的变化.

那里的任何算法都可以估计这种偏移,最好不使用插值？

Answer 1

这是一个相当有趣的问题。这是使用傅立叶变换的部分解决方案的尝试。这依赖于适度周期性的数据。我不确定它是否适用于您的数据（端点处的导数似乎不匹配）。

import numpy as np

X = np.linspace(0,2*np.pi,30)  #some X values

def yvals(x):
    return np.sin(x)+np.sin(2*x)+np.sin(3*x)

Y1 = yvals(X)
Y2 = yvals(X-0.1)  #shifted y values

#fourier transform both series
FT1 = np.fft.fft(Y1)
FT2 = np.fft.fft(Y2)

#You can show that analyically, a phase shift in the coefficients leads to a 
#multiplicative factor of `exp(-1.j * N * T_d)`

#can't take the 0'th element because that's a division by 0.  Analytically, 
#the division by 0 is OK by L'hopital's<sp?> rule, but computers don't know calculus :)
print np.log(FT2[1:]/FT1[1:])/(-1.j*np.arange(1,len(X)))

对打印输出的快速检查表明，功率最大的频率 (N=1,N=2) 给出了合理的估计，如果您查看绝对值 (np.absolute)，N=3 也可以，尽管我'我无法解释为什么会这样。

也许更熟悉数学的人可以从这里得到一个更好的答案......

Answer 2

您提供的链接之一有正确的想法（事实上我在这里做的事情几乎相同）

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

a,b, N = 0, 10, 1000        #Boundaries, datapoints
shift = -3                  #Shift, note 3/10 of L = b-a

x = np.linspace(a,b,N)
x1 = 1*x + shift
time = np.arange(1-N,N)     #Theoritical definition, time is centered at 0

y1 = sum([np.sin(2*np.pi*i*x/b) for i in range(1,5)])
y2 = sum([np.sin(2*np.pi*i*x1/b) for i in range(1,5)])

#Really only helps with large irregular data, try it
# y1 -= y1.mean()
# y2 -= y2.mean()
# y1 /= y1.std()
# y2 /= y2.std()

cross_correlation = correlate(y1,y2)
shift_calculated = time[cross_correlation.argmax()] *1.0* b/N
y3 = sum([np.sin(2*np.pi*i*(x1-shift_calculated)/b) for i in range(1,5)])
print "Preset shift: ", shift, "\nCalculated shift: ", shift_calculated



plt.plot(x,y1)
plt.plot(x,y2)
plt.plot(x,y3)
plt.legend(("Regular", "Shifted", "Recovered"))
plt.savefig("SO_timeshift.png")
plt.show()

这有以下输出：

Preset shift:  -3
Calculated shift:  -2.99

可能需要检查

1. Scipy 关联
2. 时滞分析

请注意，相关性的 argmax() 显示了对齐的位置，必须按b-a = 10-0 = 10和 N 的长度对其进行缩放才能获得实际值。

检查 correlate Source的来源，从 sigtools 导入的函数的行为并不完全明显。对于大型数据集，循环相关（通过快速傅立叶变换）比直接方法快得多。我怀疑这是 sigtools 中实现的，但我不能确定。在我的 python2.7 文件夹中搜索该文件仅返回编译后的 C pyd 文件。