我有兩個時間序列，我懷疑他們之間存在時間偏移，我想估計這個時間偏移。估計兩個時間序列之間的小時間偏移
此問題已被問： Find phase difference between two (inharmonic) waves和find time shift between two similar waveforms但在我的情況下，時移小於數據的分辨率。例如數據以小時分辨率提供，時間偏移僅爲幾分鐘（見圖）。
原因在於用於測量其中一個系列的數據記錄器在其時間內有幾分鐘的變化。
任何可以估計這種偏移的算法，最好不使用插值？
這是一個非常有趣的問題。這是一個嘗試使用傅里葉變換的部分解決方案。這依賴於中等週期性的數據。我不確定它是否可以與您的數據一起工作（在端點的衍生物似乎不匹配）。
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），得到合理的估計，N = 3確實行太多，如果你看 絕對值（np.absolute），儘管我無法解釋爲什麼會這樣。
也許有人更熟悉的數學可以把它從這裏給一個更好的答案...
一個你所提供的鏈接的看法是正確的（其實我在這裏做幾乎同樣的事情）
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得到實際值。
檢查關聯源Source從sigtools導入的函數的行爲並不完全清楚。對於大數據集，循環相關（通過快速傅里葉變換）比直接法更快。我懷疑這是在sigtools中實現的，但我無法確定。在我的python2.7文件夾中搜索文件只返回編譯後的C pyd文件。
我已經成功地使用了（在awgn通道中）匹配濾波器方法，它在索引n處給出峯值能量m [n];然後將二次多項式f（n）擬合到m [n-1]，m [n]，m [n + 1]，並通過設置f'（n）== 0來找到最小值。
響應不一定是絕對線性的，特別是如果信號的自相關不會在m [n-1]，m [n + 1]處消失。
這是一個非常有趣的問題。最初，我打算建議一個類似於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)
可替代地，我們可以使用三次樣條內插的數據（見here）。
3高斯卷積
由於有一定的插值，那麼我們不知道實際的日出和日落如何準確地進行了預測。所以，我們可以用高斯信號來卷積信號來表示這種不確定性。
gaussian_window = scipy.signal.gaussian(M, std)
sunrise_sunset_g = scipy.signal.convolve(sunrise_sunset, gaussian_window)
4.互相關
使用在user948652的回答獲得時移的互相關方法。
在這種方法中有很多未解答的問題需要對數據進行檢查和實驗以更具體地確定，比如什麼是識別日出/日落的最佳方法，高斯窗應該有多寬，等等。但是我就是這樣開始攻擊這個問題的。 祝你好運！
對於給定的約束，即該解決方案是相移少量小於採樣方法，簡單的單純形算法工作得很好。我修改了@mgilson的示例問題以顯示如何執行此操作。請注意，該解決方案非常強大，因爲它可以處理噪音。
Error函數：有可能是更理想的東西，以優化過，但這種工作得非常好：
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
爲什麼不簡單地執行X2-X1？沒有迭代需要和完美的結果！不，嚴重的是，X2是未知的，所以當你在你的err_func中使用它時，你實際上是在作弊！雖然我必須承認你激勵了我的答案...... – kadee
事實上，有趣的問題，但還沒有令人滿意的答案。讓我們嘗試改變...
你說你不喜歡使用插值，但正如我從你的評論中瞭解的那樣，你真正的意思是你想避免向更高分辨率上採樣。鹼性溶液利用最小二乘擬合與線性插值函數的，但沒有上採樣到更高的分辨率：
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
（+1）尼斯問題。出於興趣，你爲什麼禁止使用插值？ – NPE
我只是想，如果你想估計轉換到高精度，那麼你需要內插到一個非常高的分辨率。由於我有很多數據，我想避免這種情況。 – omar
在我看來，如果您的數據大致是週期性的，那麼fourier系列可能會有所幫助...... – mgilson