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