Jak mogę wyrównać / zsynchronizować dwa sygnały?

21

Robię badania, ale utknąłem na etapie analizy (powinienem był poświęcić więcej uwagi wykładom statystyk).

Zebrałem dwa jednoczesne sygnały: zintegrowane natężenie przepływu dla objętości i zmiany w rozszerzaniu klatki piersiowej. Chciałbym porównać sygnały i ostatecznie mam nadzieję uzyskać głośność na podstawie sygnału rozszerzania klatki piersiowej. Ale najpierw muszę wyrównać / zsynchronizować moje dane.

Ponieważ nagrywanie nie rozpoczyna się dokładnie w tym samym czasie, a ekspansja klatki piersiowej jest rejestrowana przez dłuższy czas, muszę znaleźć dane, które odpowiadają moim danym objętości w zestawie danych ekspansji klatki piersiowej, i mam miarę, w jakim stopniu są one wyrównane. Nie jestem do końca pewien, jak to zrobić, jeśli dwa sygnały nie zaczynają się dokładnie w tym samym czasie lub między danymi w różnych skalach i przy różnych rozdzielczościach.

Dołączyłem przykład dwóch sygnałów ( https://docs.google.com/spreadsheet/ccc?key=0As4oZTKp4RZ3dFRKaktYWEhZLXlFbFVKNmllbGVXNHc ), proszę dać mi znać, jeśli jest coś jeszcze, co mogę podać.

osoba157
źródło
Nie znam tego wystarczająco dobrze, aby udzielić odpowiedzi, i nie jestem pewien, czy to rozwiązuje pytanie, ale jedno podejście synchronizujące sygnały nazywa się „rejestracją”, która jest podzbiorem analizy danych funkcjonalnych. Temat ten omówiono w książce FDA Ramseya i Silvermana. Podstawową ideą jest to, że zaobserwowane sygnały mogą być „wypaczone” (np. Gdybyśmy byli zainteresowani mechaniką sposobu żucia ludzi, ale mamy dane o ludziach żujących przy różnych prędkościach - w tym przypadku oś czasu jest „wypaczona”) i rejestracja próbuje zdefiniować podstawowy sygnał na wspólnej, „nie wypaczonej” skali.
Makro
1
Czy zebrałeś już wszystkie swoje dane? Czy ten pilot jest przedmiotem? Jeśli dopiero zaczynasz, chciałbym oddzielić sygnał od paska i użyć go jako wyzwalacza (a nawet po prostu oznaczyć znacznik czasu) podczas rejestrowania przepływu. Zwykle systemy akwizycji mają tę zdolność z portem pomocniczym lub wyzwalającym. Jestem pewien, że istnieją sposoby na odróżnienie go po prostu przy użyciu danych, jak sugerował Macro, ale dodanie tego dodatkowego kroku zabierze wiele zgadywania.
jonsca
1
Myślę, że chcesz oszacować tylko ustalone opóźnienie. Możesz użyć korelacji krzyżowej, jak przedstawiono tutaj: stats.stackexchange.com/questions/16121/…
thias
1
Możesz zadać to pytanie na dsp.SE, gdzie myślą także o synchronizacji sygnałów.
Dilip Sarwate
1
@Thias jest poprawny, ale wygląda na to, że pierwsza seria powinna być ponownie próbkowana, aby miały wspólne przedziały.
whuber

Odpowiedzi:

24

Pytanie brzmi, jak znaleźć wielkość, o jaką jeden szereg czasowy („rozszerzenie”) jest opóźniony względem drugiego („objętość”), gdy serie są próbkowane w regularnych, ale różnych odstępach czasu.

W tym przypadku obie serie zachowują się w sposób dość ciągły, jak pokazują liczby. Oznacza to (1), że początkowe wygładzenie może być niewielkie lub wcale (2) ponowne próbkowanie może być tak proste, jak interpolacja liniowa lub kwadratowa. Kwadrat może być nieco lepszy ze względu na gładkość. Po ponownym próbkowaniu opóźnienie znajduje się poprzez maksymalizację korelacji krzyżowej , jak pokazano w wątku: W przypadku dwóch serii danych próbkowanych z przesunięciem, jakie jest najlepsze oszacowanie przesunięcia między nimi? .


Aby to zilustrować , możemy użyć danych dostarczonych w pytaniu, wykorzystując Rpseudokod. Zacznijmy od podstawowej funkcjonalności, korelacji krzyżowej i ponownego próbkowania:

cor.cross <- function(x0, y0, i=0) {
  #
  # Sample autocorrelation at (integral) lag `i`:
  # Positive `i` compares future values of `x` to present values of `y`';
  # negative `i` compares past values of `x` to present values of `y`.
  #
  if (i < 0) {x<-y0; y<-x0; i<- -i}
  else {x<-x0; y<-y0}
  n <- length(x)
  cor(x[(i+1):n], y[1:(n-i)], use="complete.obs")
}

Jest to prosty algorytm: obliczenia oparte na FFT byłyby szybsze. Ale dla tych danych (obejmujących około 4000 wartości) jest wystarczająco dobry.

resample <- function(x,t) {
  #
  # Resample time series `x`, assumed to have unit time intervals, at time `t`.
  # Uses quadratic interpolation.
  #
  n <- length(x)
  if (n < 3) stop("First argument to resample is too short; need 3 elements.")
  i <- median(c(2, floor(t+1/2), n-1)) # Clamp `i` to the range 2..n-1
  u <- t-i
  x[i-1]*u*(u-1)/2 - x[i]*(u+1)*(u-1) + x[i+1]*u*(u+1)/2
}

Pobrałem dane jako plik CSV oddzielony przecinkami i usunąłem jego nagłówek. (Nagłówek spowodował pewne problemy z R, których nie chciałem zdiagnozować).

data <- read.table("f:/temp/a.csv", header=FALSE, sep=",", 
                    col.names=c("Sample","Time32Hz","Expansion","Time100Hz","Volume"))

Uwaga: w tym rozwiązaniu założono, że każda seria danych jest uporządkowana w czasie, bez żadnych przerw w żadnej z nich. Pozwala to na wykorzystanie indeksów do wartości jako przybliżeń czasowych i skalowanie tych indeksów według częstotliwości próbkowania czasowego w celu konwersji ich na czasy.

Okazuje się, że jeden lub oba z tych instrumentów nieco z czasem dryfuje. Przed kontynuowaniem dobrze jest usunąć takie trendy. Ponadto, ponieważ na końcu występuje zwężający się sygnał głośności, powinniśmy go wyciąć.

n.clip <- 350      # Number of terminal volume values to eliminate
n <- length(data$Volume) - n.clip
indexes <- 1:n
v <- residuals(lm(data$Volume[indexes] ~ indexes))
expansion <- residuals(lm(data$Expansion[indexes] ~ indexes)

Ponownie próbkuję rzadziej występującą serię, aby uzyskać jak największą precyzję wyniku.

e.frequency <- 32  # Herz
v.frequency <- 100 # Herz
e <- sapply(1:length(v), function(t) resample(expansion, e.frequency*t/v.frequency))

Teraz można obliczyć korelację krzyżową - dla wydajności szukamy tylko rozsądnego okna opóźnień - i można określić opóźnienie, w którym można znaleźć maksymalną wartość.

lag.max <- 5       # Seconds
lag.min <- -2      # Seconds (use 0 if expansion must lag volume)
time.range <- (lag.min*v.frequency):(lag.max*v.frequency)
data.cor <- sapply(time.range, function(i) cor.cross(e, v, i))
i <- time.range[which.max(data.cor)]
print(paste("Expansion lags volume by", i / v.frequency, "seconds."))

Dane wyjściowe mówią nam, że ekspansja opóźnia głośność o 1,85 sekundy. (Jeśli ostatnie 3,5 sekundy danych nie zostały przycięte, wynik wyniósłby 1,84 sekundy.)

Dobrze jest sprawdzić wszystko na kilka sposobów, najlepiej wizualnie. Po pierwsze, funkcja korelacji krzyżowej :

plot(time.range * (1/v.frequency), data.cor, type="l", lwd=2,
     xlab="Lag (seconds)", ylab="Correlation")
points(i * (1/v.frequency), max(data.cor), col="Red", cex=2.5)

wykres korelacji krzyżowej

Następnie zarejestrujmy dwie serie w czasie i narysujmy je razem na tych samych osiach .

normalize <- function(x) {
  #
  # Normalize vector `x` to the range 0..1.
  #
  x.max <- max(x); x.min <- min(x); dx <- x.max - x.min
  if (dx==0) dx <- 1
  (x-x.min) / dx
}
times <- (1:(n-i))* (1/v.frequency)
plot(times, normalize(e)[(i+1):n], type="l", lwd=2, 
     xlab="Time of volume measurement, seconds", ylab="Normalized values (volume is red)")
lines(times, normalize(v)[1:(n-i)], col="Red", lwd=2)

Zarejestrowane działki

Wygląda całkiem nieźle! Możemy jednak lepiej zrozumieć jakość rejestracji za pomocą wykresu rozrzutu . Zmieniam kolory według czasu, aby pokazać postęp.

colors <- hsv(1:(n-i)/(n-i+1), .8, .8)
plot(e[(i+1):n], v[1:(n-i)], col=colors, cex = 0.7,
     xlab="Expansion (lagged)", ylab="Volume")

Wykres punktowy

Szukamy punktów do śledzenia tam iz powrotem wzdłuż linii: odmiany, które odzwierciedlają nieliniowości opóźnionej reakcji wzrostu na objętość. Chociaż istnieją pewne odmiany, są one dość małe. Jednak to, jak zmiany te zmieniają się w czasie, może budzić pewne zainteresowanie fizjologiczne. Cudowną cechą statystyki, zwłaszcza jej eksploracji i aspektu wizualnego, jest to, że tworzy dobre pytania i pomysły wraz z przydatnymi odpowiedziami.

Whuber
źródło
1
Cholera jasna, jesteś niesamowity. Korelacja krzyżowa jest dokładnie tym, co sobie wyobrażałem (wiedziałem, że musi istnieć jej nazwa), ale twoje odpowiedzi / wyjaśnienia wykraczały poza to. Dziękuję bardzo!
person157
Nie mam teraz czasu na pełne wyjaśnienie, ale świetne książki pojawiają się w książkach „Przepisy numeryczne”. Na przykład, spójrz na rozdziale 13.2 "Korelacja i autokorelacja Korzystanie FFT" w Numerical Recipes in C . Możesz także zajrzeć do acffunkcji R.
whuber
Nowy w 'r', proszę bądź uprzejmy: funkcja „normalizuj” zastosowana w wykresie połączonym (2. wykres ostatnio) nie działa dla mnie, czy jest aktualizacja tej funkcji od czasu opublikowania tej odpowiedzi?
CmKndy
1
@CmKndy Również byłem nowy R, kiedy opublikowałem tę odpowiedź i zapomniałem podać definicję tej funkcji. Oto oryginał:normalize <- function(x) { x.max <- max(x); x.min <- min(x); dx <- x.max - x.min; if (dx==0) dx <- 1; (x-x.min) / dx }
whuber
Idealnie, dziękuję @whuber. Jeśli możesz
napisać