2.2 Estimating clock drift
For correcting clock drift, we need to first provide some reasonably good guesses of the true drift values.
find_clock_drift() can take care of this.
This function takes a brute force approach and uses cross correlation to see how closely sync-tag detections and emissions line up over a time series of possible clock drift values.
This works well for getting values on the scale of, say, 0.1 second resolution.
For each clock_time, we’ll check possible drift values within a 1 hour-period
int_drft <- find_clock_drift(
sync_model, clock_resolution = 0.1, clock_window = 60*60*1, quiet = T)## Warning in rlm.default(x, y, weights, method = method, wt.method = wt.method, :
## 'rlm' failed to converge in 20 steps## Receiver pair: 3 ⟷ 1
## Drift (s) xcorr
## 2021-05-08 18:00:00 0.4610390 NaN
## 2021-05-08 19:00:00 0.5000000 0.73
## 2021-05-08 20:00:00 0.4857143 NaN
## 2021-05-08 21:00:00 0.5000000 0.61
## 2021-05-08 22:00:00 0.5000000 0.70
## 2021-05-08 23:00:00 0.5000000 0.73
## 2021-05-09 00:00:00 0.5000000 0.71
## 2021-05-09 01:00:00 0.6000000 0.64
## 2021-05-09 02:00:00 0.5000000 0.70
## 2021-05-09 03:00:00 0.6000000 0.75
## 2021-05-09 04:00:00 0.6000000 0.63
## 2021-05-09 05:00:00 0.6000000 0.65
## 2021-05-09 06:00:00 0.6090909 NaN
##
## Receiver pair: 3 ⟷ 2
## Drift (s) xcorr
## 2021-05-08 18:00:00 33.3 0.78
## 2021-05-08 19:00:00 33.4 0.75
## 2021-05-08 20:00:00 33.5 0.50
## 2021-05-08 21:00:00 33.6 0.60
## 2021-05-08 22:00:00 33.7 0.70
## 2021-05-08 23:00:00 33.7 0.60
## 2021-05-09 00:00:00 33.8 0.71
## 2021-05-09 01:00:00 33.9 0.72
## 2021-05-09 02:00:00 34.0 0.78
## 2021-05-09 03:00:00 34.1 0.73
## 2021-05-09 04:00:00 34.1 0.69
## 2021-05-09 05:00:00 34.2 0.57
## 2021-05-09 06:00:00 34.3 0.41
##
## Receiver pair: 3 ⟷ 4
## Drift (s) xcorr
## 2021-05-08 18:00:00 34.7 0.72
## 2021-05-08 19:00:00 34.8 0.73
## 2021-05-08 20:00:00 34.9 0.48
## 2021-05-08 21:00:00 35.0 0.61
## 2021-05-08 22:00:00 35.1 0.65
## 2021-05-08 23:00:00 35.1 0.71
## 2021-05-09 00:00:00 35.2 0.50
## 2021-05-09 01:00:00 35.3 0.64
## 2021-05-09 02:00:00 35.4 0.74
## 2021-05-09 03:00:00 35.5 0.70
## 2021-05-09 04:00:00 35.5 0.72
## 2021-05-09 05:00:00 35.6 0.50
## 2021-05-09 06:00:00 35.7 0.29
##
## Receiver pair: 3 ⟷ 5
## Drift (s) xcorr
## 2021-05-08 18:00:00 30.2 0.78
## 2021-05-08 19:00:00 30.3 0.73
## 2021-05-08 20:00:00 30.3 0.47
## 2021-05-08 21:00:00 30.5 0.62
## 2021-05-08 22:00:00 30.5 0.65
## 2021-05-08 23:00:00 30.6 0.74
## 2021-05-09 00:00:00 30.6 0.50
## 2021-05-09 01:00:00 30.7 0.60
## 2021-05-09 02:00:00 30.8 0.74
## 2021-05-09 03:00:00 30.9 0.68
## 2021-05-09 04:00:00 31.0 0.60
## 2021-05-09 05:00:00 31.1 0.52
## 2021-05-09 06:00:00 31.1 0.33
##
## Receiver pair: 4 ⟷ 6
## Drift (s) xcorr
## 2021-05-08 18:00:00 4.5 NaN
## 2021-05-08 19:00:00 4.5 NaN
## 2021-05-08 20:00:00 4.5 0.74
## 2021-05-08 21:00:00 4.5 0.75
## 2021-05-08 22:00:00 4.5 0.69
## 2021-05-08 23:00:00 4.5 0.64
## 2021-05-09 00:00:00 4.5 0.75
## 2021-05-09 01:00:00 4.5 0.66
## 2021-05-09 02:00:00 4.5 0.80
## 2021-05-09 03:00:00 4.5 0.75
## 2021-05-09 04:00:00 4.5 0.68
## 2021-05-09 05:00:00 4.5 0.73
## 2021-05-09 06:00:00 4.5 0.81To run the cross correlation, find_clock_drift treats detections and emissions between a receiver pair as two discrete time series with a resolution of clock_resolution.
These time series are centered around each clock time, spanning a period of clock_window.
Each detection/emission is then marked as a Gaussian kernel with a standard deviation of sd.
These detection/emission time series are then cross correlated and the standardized cross correlation coefficient (ranging from 0–1) is returned for each possible drift value. The drift value which results in the best alignment is then returned along with its correlation coefficient.
For datasets with more extreme clock drift, its important that sd is set sufficiently high.
Setting a higher sd value causes sync-tag emission/detection values that are almost aligned to give higher correlation scores.
This function runs quite fast. Under the hood, we’ve implemented the cross correlations as circular convolution via Fourier transforms.
Our values look good, so we’ll go ahead and set them in our model.
# Save drift values to KaltoaClockSync model
clock_drift(sync_model) <- int_drft
# Plot sync-tag diagnostics
plot(sync_model, layout = c(4,4))
As you can see in the plot, our clocks are roughly corrected within a 0.1 second resolution. For fine-scale positioning, we need to refine this further.
fit_clock_drift will automatically refine these clock drift values.
For each receiver pair, the clock drift is treated a a Gaussian random walk along some linear trend line.
The drift values are fitted as random effects in the R package TMB and the detection error is treated as a mixture distribution (see ?ddetect).
This function can provide extremely high resolution drift values.
However, for the model to fit, it requires quite good starting values.
We suggest using find_clock_drift to get drift values with at least a 0.1 second resolution before applying fit_clock_drift.
## --- Kaltoa receiver drift model
## Detection sigma (ms) Std. Err. Drift (ms h⁻¹) Std. Err.
## 3 ⟷ 1 0.077606 NaN 33.625 NaN
## 3 ⟷ 2 0.19983 0.056821 -417.79 NaN
## 3 ⟷ 4 0.056536 NaN 58.207 NaN
## 3 ⟷ 5 0.070577 NaN 73.88 NaN
## 4 ⟷ 6 0.80613 0.056001 -2.0131 NaN
## Drift sigma (ms) Std. Err. Offset (ms) Std. Err.
## 3 ⟷ 1 695.45 NaN 0.85986 NaN
## 3 ⟷ 2 1754.9 NaN -0.81203 0.00026613
## 3 ⟷ 4 3732.8 NaN 1.7767 NaN
## 3 ⟷ 5 2.1226 NaN -0.45261 NaN
## 4 ⟷ 6 302.16 NaN -3.167 5.8179e-05
##
## Drift values
## 3 ⟷ 1 3 ⟷ 2 3 ⟷ 4 3 ⟷ 5 4 ⟷ 6
## 2021-05-08 18:00:00 0.37241 33.211 34.582 30.107 4.3171
## 2021-05-08 19:00:00 0.38304 33.291 34.667 30.181 4.3272
## 2021-05-08 20:00:00 0.39370 33.372 34.752 30.255 4.3370
## 2021-05-08 21:00:00 0.40435 33.452 34.837 30.329 4.3467
## 2021-05-08 22:00:00 0.41495 33.532 34.921 30.403 4.3575
## 2021-05-08 23:00:00 0.42554 33.613 35.006 30.477 4.3664
## 2021-05-09 00:00:00 0.43614 33.693 35.091 30.551 4.3772
## 2021-05-09 01:00:00 0.44675 33.773 35.176 30.625 4.3861
## 2021-05-09 02:00:00 0.45731 33.854 35.261 30.699 4.3962
## 2021-05-09 03:00:00 0.46796 33.934 35.346 30.773 4.4057
## 2021-05-09 04:00:00 0.47859 34.014 35.431 30.847 4.4157
## 2021-05-09 05:00:00 0.48907 34.095 35.515 30.920 4.4261
## 2021-05-09 06:00:00 0.49967 34.175 35.600 30.994 4.4368# Save drift values to KaltoaClockSync model
clock_drift(sync_model) <- fit_drft
# Plot sync-tag diagnostics
plot(sync_model, type = 'error', layout = c(4,4))
Now the sync-tag diagnostic plots show that the transmission latency is more-less constant over time. Note how some receiver pairs have a bias above or below their expected transmission latencies. These can be correct for in the next section.