3.2 Estimating clock drift: Permutation
For correcting clock drift, we need to first provide some reasonably good guesses of the true drift values.
permute_clock_drift() can take care of this using a brute-force approach.
For each clock time and receiver pair combination, a series of plausible drift values will be permuted over. For these drift values, an alignment score will be calculated indicating how well the emission and detection events line up between the receivers when this drift has been taken into account. The drift value with the highest alignment score will be returned for that receiver pair and clock time.
There are two available methods for measuring how well the emissions and detections line up and both have their particular use case.
- xcorr: Emission and detection events are represented as Gaussian kernels placed along two time-series (one for each receiver). Cross correlation is applied and the most likely drift value is the one which results in the most overlap between these emission and detection kernels (i.e. the highest cross-correlation score).
- likelihood: Emissions are arranged in a time-series and a mixed detection error distribution is placed over each emission event. Detections are then sampled from this and the joint likelihood that all detection events resulted from any emission event is returned for each drift value. The drift value resulting in the maximum likelihood is chosen.
The above methods can also be used for time-difference-of-arrival detections. That is, sync-tag detections on the receiver pair which are emitted from another receiver in the array.
The xcorr method is preferred if some emission events may be missing from the dataset. However, if there are more reflections than direct detections between the receiver pair, this method risks using the late-arriving reflected sync-tag transmissions to select the drift values.
Conversely, the likelihood method directly accounts for the presence of reflections by treating the detection error as a mixed distribution (see ?ddetect()).
This method assumes that all detection events have a corresponding emission event.
If a few emission events are missing from the dataset, this can cause low likelihood scores.
Both methods are coded by treating the underlying math as convolution operations in the Fourier domain, hence they run extremely fast for the amount of data they process. These work well for getting values on the scale of, say, 0.1 second resolution.
In the following examples, we’ll mostly stick with the default parameters for the xcorr and likelihood methods. These should work fine for most use cases. For more information on how to make use of these optional parameters, see Campbell et al. (in prep) for a detailed description of these methods.
Below, we’ll use the xcorr method.
For each clock_time, we’ll check possible drift values using sync-tag emissions/detections from within a 1 hour-period (the clock_window) centered around each clock_time.
int_drft_xcor <- permute_clock_drift(
sync_model,
method = "xcorr",# Evaluate emission/detection alignment with cross-correlation
clock_resolution = 0.1,# (s) Resolution of the emission/detection time-series
clock_window = 60*60*1,# (s) Length of window to estimate drift values from
quiet = T)
print(int_drft_xcor)## Clock drift permutation values
## Receiver pair: 3 ⟷ 1
## Drift (s) xcorr
## 2021-05-08 18:00:00 0.4 0.95
## 2021-05-08 19:00:00 0.4 0.96
## 2021-05-08 20:00:00 0.4 0.91
## 2021-05-08 21:00:00 0.4 0.95
## 2021-05-08 22:00:00 0.4 0.91
## 2021-05-08 23:00:00 0.4 0.92
## 2021-05-09 00:00:00 0.4 0.94
## 2021-05-09 01:00:00 0.4 0.88
## 2021-05-09 02:00:00 0.5 0.89
## 2021-05-09 03:00:00 0.5 0.88
## 2021-05-09 04:00:00 0.5 0.89
## 2021-05-09 05:00:00 0.5 0.91
## 2021-05-09 06:00:00 0.5 0.91
##
## Receiver pair: 3 ⟷ 2
## Drift (s) xcorr
## 2021-05-08 18:00:00 33.2 0.88
## 2021-05-08 19:00:00 33.3 0.93
## 2021-05-08 20:00:00 33.4 0.85
## 2021-05-08 21:00:00 33.5 0.93
## 2021-05-08 22:00:00 33.5 0.85
## 2021-05-08 23:00:00 33.6 0.88
## 2021-05-09 00:00:00 33.7 0.91
## 2021-05-09 01:00:00 33.8 0.85
## 2021-05-09 02:00:00 33.9 0.84
## 2021-05-09 03:00:00 33.9 0.88
## 2021-05-09 04:00:00 34.0 0.84
## 2021-05-09 05:00:00 34.1 0.82
## 2021-05-09 06:00:00 34.1 0.89
##
## Receiver pair: 3 ⟷ 4
## Drift (s) xcorr
## 2021-05-08 18:00:00 34.6 0.93
## 2021-05-08 19:00:00 34.7 0.95
## 2021-05-08 20:00:00 34.8 0.89
## 2021-05-08 21:00:00 34.8 0.91
## 2021-05-08 22:00:00 34.9 0.92
## 2021-05-08 23:00:00 35.0 0.92
## 2021-05-09 00:00:00 35.1 0.88
## 2021-05-09 01:00:00 35.2 0.87
## 2021-05-09 02:00:00 35.3 0.91
## 2021-05-09 03:00:00 35.3 0.84
## 2021-05-09 04:00:00 35.4 0.79
## 2021-05-09 05:00:00 35.5 0.82
## 2021-05-09 06:00:00 35.6 0.83
##
## Receiver pair: 3 ⟷ 5
## Drift (s) xcorr
## 2021-05-08 18:00:00 30.1 0.92
## 2021-05-08 19:00:00 30.2 0.92
## 2021-05-08 20:00:00 30.3 0.88
## 2021-05-08 21:00:00 30.3 0.92
## 2021-05-08 22:00:00 30.4 0.94
## 2021-05-08 23:00:00 30.5 0.93
## 2021-05-09 00:00:00 30.5 0.91
## 2021-05-09 01:00:00 30.6 0.85
## 2021-05-09 02:00:00 30.7 0.88
## 2021-05-09 03:00:00 30.8 0.83
## 2021-05-09 04:00:00 30.9 0.90
## 2021-05-09 05:00:00 30.9 0.87
## 2021-05-09 06:00:00 31.0 0.89
##
## Receiver pair: 4 ⟷ 6
## Drift (s) xcorr
## 2021-05-08 18:00:00 4.3 0.85
## 2021-05-08 19:00:00 4.3 0.81
## 2021-05-08 20:00:00 4.3 0.73
## 2021-05-08 21:00:00 4.4 0.83
## 2021-05-08 22:00:00 4.3 0.87
## 2021-05-08 23:00:00 4.4 0.79
## 2021-05-09 00:00:00 4.4 0.80
## 2021-05-09 01:00:00 4.4 0.79
## 2021-05-09 02:00:00 4.4 0.78
## 2021-05-09 03:00:00 4.4 0.78
## 2021-05-09 04:00:00 4.4 0.78
## 2021-05-09 05:00:00 4.4 0.78
## 2021-05-09 06:00:00 4.4 0.62Note that larger clock_windows will result in more sync-tag events being used, giving more accurate drift values.
Additionally, lowering the clock_resolution will yield more precise drift values.
The trade-off is that doubling clock_window (or halving clock_resolution) will result in a roughly doubled processing time.
Next we can plot out our estimated drift values.
Let’s also do an example of using the likelihood method.
phi is more-less the maximum possible detection error that a late-arriving reflection can result in.
For example, if you expect a tag transmission can travel no father than 1km, then a conservative value for phi can be set with phi = 1000 / 1500.
In other words, we expect reflected transmissions will travel less than one additional kilometer farther than any direct detection.
int_drft_lik <- permute_clock_drift(
sync_model,
method = "like",# Evaluate emission/detection alignment with cross-correlation
clock_resolution = 0.1,# (s) Resolution of the emission/detection time-series
clock_window = 60*60*1,# (s) Length of window to estimate drift values from
phi = 1000 / 1500,# (s) Maximum additional latency from reflected detections
quiet = T)
print(int_drft_lik)## Clock drift permutation values
## Receiver pair: 3 ⟷ 1
## Drift (s) LogLik.
## 2021-05-08 18:00:00 0.3 1.63
## 2021-05-08 19:00:00 0.3 10.47
## 2021-05-08 20:00:00 0.3 5.25
## 2021-05-08 21:00:00 0.4 11.57
## 2021-05-08 22:00:00 0.4 -1.82
## 2021-05-08 23:00:00 0.4 5.93
## 2021-05-09 00:00:00 0.4 9.89
## 2021-05-09 01:00:00 0.4 3.43
## 2021-05-09 02:00:00 0.4 4.63
## 2021-05-09 03:00:00 0.4 2.19
## 2021-05-09 04:00:00 0.4 -0.37
## 2021-05-09 05:00:00 0.5 0.88
## 2021-05-09 06:00:00 0.5 2.39
##
## Receiver pair: 3 ⟷ 2
## Drift (s) LogLik.
## 2021-05-08 18:00:00 33.2 2.90
## 2021-05-08 19:00:00 33.3 6.42
## 2021-05-08 20:00:00 33.3 3.21
## 2021-05-08 21:00:00 33.4 9.70
## 2021-05-08 22:00:00 33.5 4.72
## 2021-05-08 23:00:00 33.6 1.57
## 2021-05-09 00:00:00 33.7 10.89
## 2021-05-09 01:00:00 33.7 -2.74
## 2021-05-09 02:00:00 33.8 -0.99
## 2021-05-09 03:00:00 33.9 1.07
## 2021-05-09 04:00:00 34.0 -12.37
## 2021-05-09 05:00:00 34.0 -13.35
## 2021-05-09 06:00:00 34.1 3.12
##
## Receiver pair: 3 ⟷ 4
## Drift (s) LogLik.
## 2021-05-08 18:00:00 34.6 5.42
## 2021-05-08 19:00:00 34.6 14.16
## 2021-05-08 20:00:00 34.7 9.41
## 2021-05-08 21:00:00 34.8 9.21
## 2021-05-08 22:00:00 34.9 6.21
## 2021-05-08 23:00:00 35.0 6.65
## 2021-05-09 00:00:00 35.1 4.31
## 2021-05-09 01:00:00 35.1 4.97
## 2021-05-09 02:00:00 35.2 6.70
## 2021-05-09 03:00:00 35.3 4.74
## 2021-05-09 04:00:00 35.4 -3.48
## 2021-05-09 05:00:00 35.5 2.35
## 2021-05-09 06:00:00 35.5 3.78
##
## Receiver pair: 3 ⟷ 5
## Drift (s) LogLik.
## 2021-05-08 18:00:00 30.1 4.86
## 2021-05-08 19:00:00 30.1 11.72
## 2021-05-08 20:00:00 30.2 2.79
## 2021-05-08 21:00:00 30.3 11.86
## 2021-05-08 22:00:00 30.4 7.86
## 2021-05-08 23:00:00 30.5 10.58
## 2021-05-09 00:00:00 30.5 8.58
## 2021-05-09 01:00:00 30.6 6.30
## 2021-05-09 02:00:00 30.7 5.52
## 2021-05-09 03:00:00 30.7 -4.45
## 2021-05-09 04:00:00 30.8 4.65
## 2021-05-09 05:00:00 30.9 10.71
## 2021-05-09 06:00:00 30.9 1.11
##
## Receiver pair: 4 ⟷ 6
## Drift (s) LogLik.
## 2021-05-08 18:00:00 4.3 1.11
## 2021-05-08 19:00:00 4.3 8.38
## 2021-05-08 20:00:00 4.3 -4.62
## 2021-05-08 21:00:00 4.3 6.54
## 2021-05-08 22:00:00 4.3 2.09
## 2021-05-08 23:00:00 4.3 9.59
## 2021-05-09 00:00:00 4.3 -5.13
## 2021-05-09 01:00:00 4.3 -3.83
## 2021-05-09 02:00:00 4.4 -0.47
## 2021-05-09 03:00:00 4.4 -1.30
## 2021-05-09 04:00:00 4.4 -5.92
## 2021-05-09 05:00:00 4.4 -5.21
## 2021-05-09 06:00:00 4.4 -21.49
Here, the quality of the drift estimates are given by a log likelihood score. Clock drift values that reduce the detection error will result in relatively higher log likelihood scores.
Now we’ll store our initial drift values in our KaltoaClockSync object.
# Save drift values to KaltoaClockSync model
clock_drift(sync_model) <- int_drft_lik
# 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’ll need to refine this further using fit_clock_drift.