Linear fit to wrapped / periodic data

Question

Assume I have a number $N$ of noisy (phase) data $\varphi_i$, i.e., $\varphi_i \in [0,2\pi)$. I know that between $\varphi_{i-1}$ and $\varphi_i$ there is a constant phase shift $\Delta \varphi$ (neglecting the noise). If I want to determine $\Delta \varphi$ from the noisy data, I could do a linear fit: $$ \varphi_i = \varphi_0 + i \Delta\varphi $$ The only problem is the "wrapping" of $\varphi_i$ at $2\pi$. Actually I would have to find a least squares solution to something like: $$ \varphi_i = (\varphi_0 + i \Delta\varphi \mod 2\pi) $$ It would be possible to do an "unwrapping" of the data before by shifting $\varphi_i$ by $2n_i\pi$ such that $|\varphi_{i-1}-\varphi_i|$ gets minimal and applying a linear fit afterwards (see the first equation). However, I wonder whether there is a way to do it "directly" without the unwrapping, which might potentially affect the result (especially in case of a more complex phase shift model, which is what I am actually heading to...).

Answer 1

Let me start with the disclaimer that there seems to be a great deal of literature on the problem of phase unwrapping, and I am not familiar with any of it. So it's possible that there exist much better techniques than my answer here.

You could start by representing your phases as unit complex numbers, $z_k = e^{i\phi_k}.$ Then the linear phase model that you want to fit is just a complex sinusoid, $z_k = z_0 e^{ik\Delta\phi}$. If you take the Fourier transform of $z_k$ and pick the largest-amplitude mode, that should be very close to the optimal fit, and then you can apply a nonlinear least-squares algorithm like Levenberg-Marquardt to fine-tune the solution.

Notes:

1. The Fourier step is necessary because it's very easy for a nonlinear minimization to get stuck in a local minimum, so it's important to start somewhere close to the global optimum.

2. If the data were periodic, you could just take the largest Fourier mode and call it a day; in this case, the Fourier transform would be a spike at the true frequency and zero everywhere else (in the absence of noise, at least). Even otherwise, the mode with the largest amplitude ought to be nothing but the one closest to the true frequency, so you should get a really good initialization — I'm pretty sure this is true, but I haven't worked out the details.