2nd (actually nth) order clock offset correction?

[dmedine]

Has anyone ever evaluated the accuracy of load_xdf's drift correction algorithm? It is a non-standard algorithm (i.e. not part of Matlab) for fitting a first order curve to correct for clock drift between sending and receiving devices

https://github.com/xdf-modules/xdf-Matlab/blob/master/load_xdf.m#L505-L516

xdf-Matlab/load_xdf.m

Lines 773 to 815 in 0cdf054

	function x = robust_fit(A,y,rho,iters)
	% Perform a robust linear regression using the Huber loss function.
	% x = robust_fit(A,y,rho,iters)
	%
	% Input:
	% A : design matrix
	% y : target variable
	% rho : augmented Lagrangian variable (default: 1)
	% iters : number of iterations to perform (default: 1000)
	%
	% Output:
	% x : solution for x
	%
	% Notes:
	% solves the following problem via ADMM for x:
	% minimize 1/2*sum(huber(A*x - y))
	%
	% Based on the ADMM Matlab codes also found at:
	% https://web.stanford.edu/~boyd/papers/admm_distr_stats.html
	%
	% Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
	% 2013-03-04
	if ~exist('rho','var')
	rho = 1; end
	if ~exist('iters','var')
	iters = 1000; end
	x_offset = min(A(:, 2));
	A(:, 2) = A(:, 2) - x_offset;
	Aty = A'*y;
	L = chol( A'*A, 'lower' );
	L = sparse(L);
	U = sparse(L');
	z = zeros(size(y)); u = z;
	for k = 1:iters
	x = U \ (L \ (Aty + A'*(z - u)));
	d = A*x - y + u;
	z = rho/(1+rho)*d + 1/(1+rho)*max(0,(1-(1+1/rho)./abs(d))).*d;
	u = d - z;
	end
	x(1) = x(1) - x(2)*x_offset;
	end

The reason I ask is that it has come to my attention that 1st order curves may not be appropriate for typical drift. For example, Here is the difference between LSL Markers and a screen flip signal recorded via an EEG device.

The total marker/EEG drift changes from about 27 to 23 ms over the course of several hours. This is definitely enough to muck up ERP/ERSP analysis. Such second order curves are pretty typical but an nth order curve is probably what is best.

If I were to rewrite this, I would use polyfit and do away with the ADMM method using the Huber loss function as minimization criteria but maybe there is a compelling reason for using such a robust tool. I confess that the main reason I would choose polyfit is that I don't fully understand how to adapt the code to return an nth order set of mapping coefficients. scipy has analogs to polyfit but I don't know quite what to do in C++, but there must a be a good lib out there for that stuff.

[cboulay]

Please make sure this GitHub issue in pyxdf in on your radar when thinking about this. As Matthew mentioned in that thread, he prefers a rolling window linear fit.

I kind of like the rolling window linear fit idea because it's closer to what liblsl does if someone enables online clock sync with dejittering.

Either method is a major improvement, and there's no reason the xdf importer can't do both:

• Have a parameter for rolling window duration (0 is the whole file).
• Have a separate parameter for the fit order.

27 to 23 ms over the course of several hours. This is definitely enough to muck up ERP/ERSP analysis.

That's interesting. I wouldn't think a 4 msec shift can affect ERSP. What's the application? I understand how this would affect gamma and phase alignment but that's not normally in the "ERSP" bucket. I agree the current implementation should be improved; I'm just looking to learn about other high-temporal-precision applications.

[dmedine]

Cool, I wasn't aware of that other issue. I may have some time/incentive to implement this option in Matlab. I will aim for the rolling window if that is the preferred interpolation.

It could be that 4ms isn't a show stopper with ERSPs. I am not an expert in EEG analysis. I do know that even a few ms will attenuate your ERPs. It depends, of course, on how many trials you have etc. Also, setup time and experiment time. At SCCN, a tolerance of 1ms was always the mantra so I just assumed that it was important for ERSP as well---but maybe it was just an arbitrary figure, like 10000 steps per day. Then again, if you are timelocking ERSPs to things like people jogging on a treadmill or turning their head, which will probably require time warping to do the averaging, it seems to me that you want to be as accurate as possible.

[chkothe]

Hey, I agree that offering a sliding window (or adaptive RLS) fit would be a somewhat neater solution for a non-linear fit than trying to do a polynomial fit, mainly since it'd be the same that LSL does online (which we'd indeed want to have as a choice for offline simulation in any case), and more alternatives beyond that would probably mainly serve to create confusion (there'd be a number of ways to do a polynomial fit, e.g., Legendre polynomials or a series of band-limited sinusoids etc, and if the fits are supposed to be robust to outliers that'd add even more complexity and uncertainty as to whether the fit actually worked or got derailed somehow).

A note re the 4ms above -- that should be without the first-order correction (but correct me if I'm wrong, because otherwise there shouldn't be a linear / slope component in there), so the actual residual error between linear and a hypothetical non-linear fit should be more like 1ms at the point of max deviation.

As for the RLS time constant, the one that LSL uses had been tuned to be a more or less optimal tradeoff between the achievable precision (sub-ms), which requires a longer integration window, and the ability to track more or less rapid rate changes (which requires a shorter window, and which we reasoned are mostly driven by device/ambient temperature changes that play out over multi-minute time scales). I recall having written something on that but not sure if it was on github or in some code comments.

[cboulay]

I'm working with a user who happens to have recordings where the outlet was from a PC with a pretty bad clock (drifts 1.5 s / hr and highly variable). They are using pyxdf and the _robust_fit linear regression is doing its best but I think they'd be better off with an adaptive fit.

I agree that offering a sliding window (or adaptive RLS) fit would be a somewhat neater solution for a non-linear fit than trying to do a polynomial fit, mainly since it'd be the same that LSL does online (which we'd indeed want to have as a choice for offline simulation in any case)

@chkothe , I couldn't find a sliding window for clock offset in the liblsl source code. All that I found was that online clock correction uses the "best" offset from the last probe burst (see result_aggregation_scheduled), which could be locally bad if there's a particularly bad burst.

The dejitterer updates smoothly, and for continuous streams this probably compensates for outlier offset measurements, but only if enabled and it doesn't apply to irregular streams.

Should the online clock offset estimation also have some smoothing?

[chkothe]

Interesting idea, yeah. One could have an option in the config for doing online smoothing for particularly bad clocks. For recording of course one would want to have raw offsets.

Theoretically one could use the same RLS code for that, but getting the time constant right can be a bit crafty, especially considering that those offsets would normally be at much lower sampling rate (0.2 hz). One could up the offset sampling rate maybe to 1hz and do the math to see what would be a good default time constant, or alternatively look into sliding window estimation.