Oszacowanie małego przesunięcia czasowego między dwiema szeregami czasowymi

Question

Mam dwie serie czasowe i podejrzewam, że istnieje między nimi przesunięcie czasowe i chcę oszacować przesunięcie czasowe.

Na to pytanie zadano wcześniej: Find phase difference between two (inharmonic) waves i find time shift between two similar waveforms, ale w moim przypadku przesunięcie czasowe jest mniejsze niż rozdzielczość danych. na przykład dane są dostępne w rozdzielczości godzinowej, a przesunięcie czasowe to tylko kilka minut (patrz zdjęcie).

Powodem tego jest to, że rejestrator używany do pomiaru jednej z serii ma kilka minut zmiany w swoim czasie.

Jakieś algorytmy, które mogą oszacować to przesunięcie, najlepiej bez użycia interpolacji?

Answer 1

Jest to dość ciekawy problem. Oto próba częściowego rozwiązania przy użyciu transformacji Fouriera. Opiera się to na danych umiarkowanie okresowych. Nie jestem pewien, czy to będzie działać z twoimi danymi (gdzie derywaty na punktach końcowych nie pasują do siebie).

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))) 

Szybka inspekcja wydruku pokazuje, że częstotliwości z najbardziej mocy (n = 1, n = 2) dają uzasadnione szacunki, N = 3 robi OK zbyt jeśli spojrzeć na wartości bezwzględnej na (np.absolute), chociaż nie jestem w stanie wyjaśnić, dlaczego tak się stało.

Może ktoś bardziej zaznajomieni z matematyki może go stąd dać lepszą odpowiedź ...

Answer 2

Jeden z linków dostarczonych ma prawo pomysł (w rzeczywistości robię prawie to samo tutaj)

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() 

ten ma następujące dane wyjściowe:

Preset shift: -3 
Calculated shift: -2.99 

Może być konieczne sprawdzenie

1. Scipy Correlate
2. Time Delay Analaysis

Uwaga że The argmax() korelacji pokazuje położenie wyrównania, to musi być skalowane przez długość z b-a = 10-0 = 10 i N, aby uzyskać aktualną wartość.

Sprawdzanie źródła korelacji Source nie jest całkowicie oczywiste, jaka jest zachowana funkcja z sigtools. W przypadku dużych zbiorów danych korelacja kołowa (za pomocą szybkich transformacji Fouriera) jest znacznie szybsza niż metoda bezpośrednia. Podejrzewam, że jest to zaimplementowane w sigtools, ale nie jestem pewien. Poszukiwanie pliku w moim folderze python2.7 tylko zwróciło skompilowany plik C pyd.

Answer 3

Z powodzeniem zastosowałem (w kanale ANN) dopasowane podejście filtrowe, które daje energię szczytową m [n] w indeksie n; następnie dopasowanie wielomianu 2 stopnia f (n) do m [n-1], m [n], m [n + 1] i znalezienie minimum przez ustawienie f '(n) == 0.

Odpowiedź niekoniecznie jest absolutnie liniowa, szczególnie jeśli autokorelacja sygnału nie znika w m [n-1], m [n + 1].

Answer 4

To bardzo interesujący problem. Pierwotnie chciałem zaproponować rozwiązanie oparte na korelacji krzyżowej podobne do user948652. Jednak z opisem problemu, istnieją dwa problemy z tego rozwiązania:

1. Rozdzielczość danych jest większy niż przesunięciem czasowym, a
2. W niektóre dni, przewidywana wartość i wartości pomiarowe mają bardzo niska korelacja siebie

w wyniku tych dwóch kwestii, myślę, że bezpośrednio stosując rozwiązanie cross-korelacji może rzeczywiście zwiększyć przesunięcie czasu, zwłaszcza w dniach, w których przewidywane i zmierzone wartości mają bardzo niska korelacja względem siebie.

W powyższym komentarzu zapytałem, czy miałeś jakieś wydarzenia, które występują w obu szeregach czasowych, i powiedziałeś, że nie. Jednak na podstawie domeny, myślę, że rzeczywiście mają dwa:

1. Sunrise
2. Sunset

Nawet jeśli reszta sygnału jest słabo skorelowane, wschody i zachody słońca powinien być nieco skorelowane, ponieważ monotonicznie wzrosną one od/poniżej poziomu podstawowego w nocy. Oto potencjalne rozwiązanie, oparte na tych dwóch zdarzeniach, które powinno zarówno zminimalizować potrzebną interpolację, jak i nie być zależne od korelacji krzyżowej słabo skorelowanych sygnałów.

1. Znajdź przybliżona Sunrise/Sunset

ta powinna być na tyle łatwe, po prostu podjąć pierwsze i ostatnie punkty danych, które są wyższe niż płaskiej linii nocnej i oznaczyć te przybliżoną wschód i zachód słońca. Następnie, chciałbym skupić się na tych danych, a także natychmiast punkty po obu stronach, tj .:

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 

Istnieje kilka sposobów, aby wdrożyć get_sunrise() i get_sunset() w zależności od tego, ile trzeba rygoru w swojej analizie. Chciałbym użyć numpy.diff, ustawić go na określoną wartość i przyjąć pierwszy i ostatni punkt powyżej tej wartości. Można również odczytać dane nocne z dużej liczby plików, obliczyć średnie odchylenie standardowe i odszukać pierwsze i ostatnie punkty danych, które przekraczają, powiedzmy, 0.5 * st_dev danych nocnych. Można również wykonać dopasowanie oparte na klastrze, w szczególności jeśli różne klasy dnia (np. Słoneczne lub częściowo zachmurzone lub bardzo słabe) mają wysoce stereotypowe zdarzenia wschodu i zachodu słońca.

2. Resample danych

Nie sądzę, że istnieje jakiś sposób, aby rozwiązać ten problem bez jakiejś interpolacji. Chciałbym użyć ponownie próbki danych do wyższej częstotliwości próbkowania niż przesunięcie. Jeśli przesunięcie jest w skali minut, a następnie wzrośnie do 1 minuty lub 30 sekund.

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

Alternatywnie, możemy użyć sześcienny splajn do interpolacji danych (patrz here).

3. Gaussa Splot

Ponieważ istnieją pewne interpolacja, to nie wiemy, jak dokładnie faktyczny wschodu i zachodu słońca były przewidywane. Tak więc możemy skonfigwać sygnał za pomocą gaussa, aby przedstawić tę niepewność.

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

4. Cross-Korelacja

pomocą metody korelacji krzyżowej w odpowiedzi user948652 do uzyskania przesunięcia czasowego.

Istnieje wiele pytań bez odpowiedzi w tej metodzie, które wymagałyby zbadania i eksperymentowania z danymi w celu dokładniejszego ich umiejscowienia, np. Jaka jest najlepsza metoda określania wschodu/zachodu słońca, jak szerokie powinno być okno gaussowskie, itp. Ale w ten sposób zacznę atakować problem. Powodzenia!

Answer 5

Optymalizacja najlepszego rozwiązania

Dla ograniczenia podane, a mianowicie, że rozwiązanie jest przesunięty w fazie o małej ilości mniejszej niż metody pobierania próbek, prosty algorytm górki simplex działa dobrze. Zmodyfikowałem przykładowy problem @mgilsona, aby pokazać, jak to zrobić. Należy pamiętać, że to rozwiązanie jest solidne, ponieważ może poradzić sobie z hałasem.

funkcja Error: Nie może być bardziej optymalne rzeczy do optymalizacji skończyła, ale to działa zaskakująco dobrze:

np.sqrt((X1-X2+delta_x)**2+(Y1-Y2)**2).sum() 

Oznacza to, zminimalizować odległość euklidesową między dwiema krzywymi tylko przez regulację osi x (faza).

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) 

Trasa próbka daje:

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 

Answer 6

Rzeczywiście, ciekawy problem, ale jeszcze nie satysfakcjonująca odpowiedź. Spróbujmy to zmienić ...

Mówisz, że wolisz nie używać interpolacji, ale, jak rozumiem z twojego komentarza, naprawdę masz na myśli to, że chciałbyś uniknąć upsamplingu do wyższej rozdzielczości.Podstawowym rozwiązaniem korzysta z najmniejszych kwadratów pasuje do funkcji liniowej interpolacji, ale bez upsampling do wyższej rozdzielczości:

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 

Trasa próbka daje dość dokładne rozwiązanie:

Unknown shift: 0.0695701123582 
Found shift: 0.0696105501967 

Jeśli zawiera hałasu w przesuniętej Y:

Y_shifted += .1*np.random.normal(size=X.shape) 

Jeden dostaje nieco mniej dokładne wyniki:

Unknown shift: 0.0695701123582 
Found shift: 0.0746643381744 

Dokładność w obecności szumu poprawia się, gdy dostępnych jest więcej danych, np. z:

X = np.arange(0,200*np.pi,dx) 

Typowym wynikiem jest:

Unknown shift: 0.0695701123582 
Found shift: 0.0698527939193