2 つの時系列があり、それらの間に時間のずれがあると思われます。この時間のずれを推定したいと思います。
この質問は以前に尋ねられました: 2 つの (不調和な) 波の間の位相差を見つけ、2 つの類似した波形の間の時間シフトを見つけますが、私の場合、時間シフトはデータの分解能よりも小さくなります。たとえば、データは 1 時間ごとの解像度で利用でき、タイム シフトはわずか数分です (画像を参照)。
これの原因は、シリーズの 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 でも問題ありませんが、それがなぜなのかを説明するのに途方に暮れています。
たぶん、数学に詳しい人がここからそれを取って、より良い答えを出すことができます...
あなたが提供したリンクの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
確認が必要な場合があります
相関のargmax（）はアライメントの位置を示していることに注意してください。実際の値を取得するには、長さb-a = 10-0 = 10とNでスケーリングする必要があります。
相関ソースのソースを確認すると、sigtoolsからインポートされた関数がどのように動作するかが完全に明らかではありません。大規模なデータセットの場合、（高速フーリエ変換を介した）循環相関は、単純な方法よりもはるかに高速です。これがsigtoolsに実装されているものだと思いますが、はっきりとはわかりません。python2.7フォルダー内のファイルを検索すると、コンパイルされたCpydファイルのみが返されました。
これは非常に興味深い問題です。当初、user948652 のものと同様の相互相関ベースのソリューションを提案するつもりでした。ただし、問題の説明から、その解決策には 2 つの問題があります。
これら 2 つの問題の結果として、相互相関ソリューションを直接適用すると、特に予測値と測定値の相関が非常に低い日には、実際にタイム シフトが増加する可能性が高いと思います。
上記の私のコメントで、両方の時系列で発生するイベントがあるかどうかを尋ねましたが、あなたはないと答えました. ただし、ドメインに基づいて、実際には次の 2 つがあると思います。
信号の残りの部分があまり相関していなくても、日の出と日の入りは夜間のベースラインから単調に増加/減少するため、ある程度相関しているはずです。したがって、必要な補間を最小限に抑え、相関の低い信号の相互相関に依存しない、これら 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()あり、分析に必要な厳密さに応じて異なります。を使用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 次スプラインを使用してデータを補間することもできます (こちらを参照)。
3. ガウス畳み込み
補間があるため、実際の日の出と日の入りがどの程度正確に予測されたかはわかりません。したがって、信号をガウスで畳み込み、この不確実性を表すことができます。
gaussian_window = scipy.signal.gaussian(M, std)
sunrise_sunset_g = scipy.signal.convolve(sunrise_sunset, gaussian_window)
4.相互相関
user948652 の回答の相互相関法を使用して、タイム シフトを取得します。
この方法には、日の出/日の入りを特定するための最良の方法は何か、ガウスウィンドウの幅など、より具体的に特定するためにデータの調査と実験が必要な多くの未解決の質問があります。しかし、それはどのように問題に取り組み始めるか。幸運を！
確かに、興味深い問題ですが、満足のいく答えはまだありません。それを変えてみよう...
あなたは補間を使用したくないと言っていますが、あなたのコメントから私が理解しているように、あなたが本当に意味することは、より高い解像度へのアップサンプリングを避けたいということです。基本的な解決策は、線形補間関数を使用した最小二乗法を使用しますが、より高い解像度へのアップサンプリングは行いません。
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
インデックス n でピーク エネルギー m[n] を与える一致フィルター アプローチを (awgn チャネルで) 使用することに成功しました。次に、2 次多項式 f(n) を m[n-1]、m[n]、m[n+1] に当てはめ、f'(n)==0 を設定して最小値を見つけます。
応答は、特に信号の自己相関が m[n-1]、m[n+1] でゼロにならない場合、必ずしも完全に線形であるとは限りません。
与えられた制約、つまり、解がサンプリング法よりも少ない量だけ位相シフトされるという制約に対して、単純なダウンヒル シンプレックス アルゴリズムがうまく機能します。@mgilson のサンプル問題を修正して、これを行う方法を示しました。このソリューションは、ノイズを処理できるという点で堅牢であることに注意してください。
エラー関数: 最適化するより最適なものがあるかもしれませんが、これは驚くほどうまく機能します:
np.sqrt((X1-X2+delta_x)**2+(Y1-Y2)**2).sum()
つまり、x 軸 (位相) のみを調整して、2 つの曲線間のユークリッド距離を最小化します。
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