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

Question

我有两个时间序列，我怀疑它们之间有一个时间偏移，我想估计一下这个时间偏移。

这个问题以前在： Find phase difference between two (inharmonic) waves和find time shift between two similar waveforms但是在我的例子中，时间偏移小于数据的分辨率。例如，数据的分辨率为每小时一次，而时间偏移仅为几分钟（见图）。

原因是用于测量其中一个系列的数据记录器的时间有几分钟的偏移。

有没有哪种算法可以估计这种偏移，最好不用插值？

Answer 1

你提供的其中一个链接有正确的想法（事实上我在这里做的几乎是相同的事情）

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 Correlate
2. Time Delay Analaysis

注意，相关性的argmax（）表示对齐的位置，它必须按b-a = 10-0 = 10和N的长度缩放才能得到实际值。

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

Answer 2

这是一个非常有趣的问题。最初，我打算提出一个类似于用户948652的基于互相关的解决方案。但是，从您的问题描述来看，该解决方案有两个问题：

1. 数据的分辨率大于时移，并且
2. 有时，预测值和测量值之间的相关性很低

由于这两个问题，我认为直接应用互相关解决方案可能实际上会增加您的时间偏移，特别是在预测值和测量值之间相关性非常低的日子。

在我上面的评论中，我问你在这两个时间序列中是否发生过任何事件，你说你没有。但是，根据你的领域，我认为你实际上有两个：

1. 日出
2. 日落

即使信号的其余部分相关性很差，日出和日落也应该具有一定的相关性，因为它们从/减少到夜间基线将单调地增加。所以这里有一个潜在的解决方案，基于这两个事件，这两个都应该最小化所需的插值，而不是依赖于低相关信号的互相关。

1。寻找大致的日出/日落

这应该很简单，只需取第一个和最后一个数据点，这些数据点高于夜间平面线，并将它们标记为日出和日落的近似值。然后，我将重点关注这些数据，以及任何一边的点，即：

width=1
sunrise_index = get_sunrise()
sunset_index = get_sunset()

# set the data to zero, except for the sunrise/sunset events.
bitmap = zeros(data.shape)
bitmap[sunrise_index - width : sunrise_index + width] = 1
bitmap[sunset_index - width : sunset_index + width] = 1
sunrise_sunset = data * bitmap 

有几种方法可以实现get_sunrise()和get_sunset()，这取决于您的分析需要多严格。我将使用^{}，以特定的值设置阈值，并取该值之上的第一个和最后一个点。您还可以从大量文件中读取夜间数据，计算平均值和标准偏差，并查找超过夜间数据的第一个和最后一个数据点，例如0.5 * st_dev。您还可以进行某种基于集群的模板匹配，特别是如果一天中的不同类别（即晴天与部分阴天与非常阴天）具有非常典型的日出/日落事件。

2。重新采样数据

我认为没有插值就没有办法解决这个问题。我将使用重采样的数据，以更高的采样率比移位。如果轮班时间在分钟范围内，则向上取样至1分或30秒。

num_samples = new_sample_rate * sunrise_sunset.shape[0]
sunrise_sunset = scipy.signal.resample(sunrise_sunset, num_samples)

或者，我们可以使用三次样条插值数据（参见here）。

3。高斯卷积

因为有一些插值，所以我们不知道实际的日出和日落是如何精确预测的。所以，我们可以用高斯函数来卷积信号，来表示这种不确定性。

gaussian_window = scipy.signal.gaussian(M, std)
sunrise_sunset_g = scipy.signal.convolve(sunrise_sunset, gaussian_window)

4。互相关

使用用户948652答案中的互相关方法获得时间偏移。

在这种方法中，有很多未回答的问题需要对数据进行检查和实验，以便更具体地确定，例如识别日出/日落的最佳方法是什么、高斯窗口应该有多宽等等，但这是我将如何开始解决问题的。 祝你好运！

Answer 3

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

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），尽管我无法解释为什么会这样。

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