我有两个时间序列，我怀疑它们之间存在时间偏移，我想估计这个时间偏移。
这个问题之前已经被问过： 查找两个（非谐波）波之间的相位差并查找两个相似波形之间的时间偏移，但在我的情况下，时间偏移小于数据的分辨率。例如，数据以小时分辨率提供，时间偏移只有几分钟（见图）。
其原因是用于测量其中一个系列的数据记录器在其时间上有几分钟的变化。
有什么算法可以估计这种变化，最好不使用插值？
问问题
14751 次
6 回答
5
这是一个相当有趣的问题。这是使用傅立叶变换的部分解决方案的尝试。这依赖于适度周期性的数据。我不确定它是否适用于您的数据（端点处的导数似乎不匹配）。
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 也可以，尽管我米茫然地解释为什么会这样。
也许更熟悉数学的人可以从这里得到更好的答案......
于 2012-12-11T20:07:50.593 回答
3
您提供的链接之一有正确的想法（实际上我在这里做的几乎相同）
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
可能需要检查
请注意，相关性的 argmax() 显示了对齐的位置，它必须按长度b-a = 10-0 = 10和 N 进行缩放才能得到实际值。
检查 correlate Source的来源，从 sigtools 导入的函数的行为并不完全清楚。对于大型数据集，循环相关（通过快速傅立叶变换）比直接方法快得多。我怀疑这是在 sigtools 中实现的，但我不能确定。在我的 python2.7 文件夹中搜索该文件仅返回编译后的 C pyd 文件。
于 2012-12-11T23:21:03.053 回答
2
确实，有趣的问题，但还没有令人满意的答案。让我们试着改变它...
您说您不喜欢使用插值，但是，正如我从您的评论中了解到的那样，您真正的意思是您希望避免上采样到更高分辨率。一个基本的解决方案是利用最小二乘拟合线性插值函数，但没有上采样到更高分辨率：
import numpy as np
from scipy.interpolate import interp1d
from scipy.optimize import leastsq
def yvals(x):
return np.sin(x)+np.sin(2*x)+np.sin(3*x)
dx = .1
X = np.arange(0,2*np.pi,dx)
Y = yvals(X)
unknown_shift = np.random.random() * dx
Y_shifted = yvals(X + unknown_shift)
def err_func(p):
return interp1d(X,Y)(X[1:-1]+p[0]) - Y_shifted[1:-1]
p0 = [0,] # Inital guess of no shift
found_shift = leastsq(err_func,p0)[0][0]
print "Unknown shift: ", unknown_shift
print "Found shift: ", found_shift
样本运行给出了一个非常准确的解决方案：
Unknown shift: 0.0695701123582
Found shift: 0.0696105501967
如果在移位的 Y 中包含噪声：
Y_shifted += .1*np.random.normal(size=X.shape)
一个得到不太精确的结果：
Unknown shift: 0.0695701123582
Found shift: 0.0746643381744
当有更多数据可用时，存在噪声时的准确性会提高，例如：
X = np.arange(0,200*np.pi,dx)
一个典型的结果是：
Unknown shift: 0.0695701123582
Found shift: 0.0698527939193
于 2014-03-30T15:09:47.240 回答
2
这是一个非常有趣的问题。最初，我打算建议一个类似于 user948652 的基于互相关的解决方案。但是，根据您的问题描述，该解决方案存在两个问题：
- 数据的分辨率大于时移，并且
- 在某些日子里，预测值和测量值之间的相关性非常低
由于这两个问题，我认为直接应用互相关解决方案实际上可能会增加您的时间偏移，尤其是在预测值和测量值彼此相关性非常低的日子里。
在我上面的评论中，我问你是否有任何事件发生在两个时间序列中，你说你没有。但是，根据您的域，我认为您实际上有两个：
- 日出
- 日落
即使信号的其余部分相关性较差，日出和日落也应该有一定的相关性，因为它们将从夜间基线单调增加/减少。因此，这是一个基于这两个事件的潜在解决方案，它既应该最小化所需的插值，又不依赖于相关性差的信号的互相关。
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()具体取决于您在分析中需要多少严格性。我会使用numpy.diff，将其设置为特定值，然后取高于该值的第一个点和最后一个点。您还可以从大量文件中读取夜间数据，计算平均值和标准偏差，并查找超过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)
或者，我们可以使用三次样条对数据进行插值（参见此处）。
3. 高斯卷积
由于有一些插值，所以我们不知道实际日出和日落的预测有多精确。因此，我们可以将信号与高斯进行卷积，以表示这种不确定性。
gaussian_window = scipy.signal.gaussian(M, std)
sunrise_sunset_g = scipy.signal.convolve(sunrise_sunset, gaussian_window)
4. 互相关
使用 user948652 的答案中的互相关方法来获得时移。
这种方法有很多悬而未决的问题，需要对数据进行检查和实验才能更具体地确定，例如识别日出/日落的最佳方法是什么，高斯窗口应该有多宽等。但它是我将如何开始解决这个问题。祝你好运！
于 2012-12-13T16:34:23.080 回答
0
我已经成功使用（在 awgn 通道中）匹配滤波器方法，它在索引 n 处给出峰值能量 m[n]；然后将二次多项式 f(n) 拟合到 m[n-1]、m[n]、m[n+1] 并通过设置 f'(n)==0 找到最小值。
响应不一定是绝对线性的，特别是如果信号的自相关在 m[n-1]、m[n+1] 处没有消失。
于 2012-12-13T07:15:37.093 回答
0
优化以获得最佳解决方案
对于给定的约束，即解的相移比采样方法少一点，简单的下坡单纯形算法效果很好。我已经修改了@mgilson 的示例问题来展示如何做到这一点。请注意，此解决方案是稳健的，因为它可以处理噪声。
错误函数：可能有更多优化的东西需要优化，但效果出奇的好：
np.sqrt((X1-X2+delta_x)**2+(Y1-Y2)**2).sum()
即，仅通过调整 x 轴（相位）来最小化两条曲线之间的欧几里得距离。
import numpy as np
def yvals(x):
return np.sin(x)+np.sin(2*x)+np.sin(3*x)
dx = .1
unknown_shift = .03 * np.random.random() * dx
X1 = np.arange(0,2*np.pi,dx) #some X values
X2 = X1 + unknown_shift
Y1 = yvals(X1)
Y2 = yvals(X2) # shifted Y
Y2 += .1*np.random.normal(size=X1.shape) # now with noise
def err_func(p):
return np.sqrt((X1-X2+p[0])**2+(Y1-Y2)**2).sum()
from scipy.optimize import fmin
p0 = [0,] # Inital guess of no shift
found_shift = fmin(err_func, p0)[0]
print "Unknown shift: ", unknown_shift
print "Found shift: ", found_shift
print "Percent error: ", abs((unknown_shift-found_shift)/unknown_shift)
示例运行给出：
Optimization terminated successfully.
Current function value: 4.804268
Iterations: 6
Function evaluations: 12
Unknown shift: 0.00134765446268
Found shift: 0.001375
Percent error: -0.0202912082305
于 2012-12-13T17:19:27.677 回答