# phase unwrapping relative to a reference

I have samples of wrapped $$(-\pi,\pi)$$ phase from a bandlimited signal $$A(t)\cos(\pi\omega t + \phi(t)+\theta_0 )$$

Frequency is constrained by very regular physics and I know a preliminary guess $$\hat{\omega}$$. I can also write the deterministic and interesting part of $$\phi(t)$$ as $$\bar{\phi}$$. So this could be written:

$$A(t)\cos((\hat{\omega} + \Delta \omega) t + \bar{\phi(t)} + \phi(t) + \theta_0)$$

The phases come from arctan formulas and $$A(t)$$ plays no further role in the question. To make some of the scales concrete, the sampling rate would be 96Hz, $$\hat{\omega}$$ would be very close to 1,2,3 or 4 Hz. The correction to the nominal frequency would be $$\Delta \omega \approx O(0.01)$$ and is a persistent over the dataset. One answer below makes the point that phase and frequency are somewhat interchangeable. I am willing to do anything efficient for the analysis, but the reporting will make frequency and phase identifiable by assuming frequency is constant.

A deterministic fluctuation in phase/frequency $$\bar{\phi(t)}$$ occurs over 5-15 cycles of the 1Hz frequency. These are of interest and are shown below. The signal also contains noise $$\phi(t)$$ polluted by noise due to instrumentation and filtering of the signal.

I would like to quantify $$\bar{phi(t)}$$ and $$\theta_0$$ and $$\Delta\omega$$. $$\bar{\phi(t)}$$ can get large episodically, but most of the time it should be in the range $$(-\pi/4,\pi/4)$$. The scalar $$\theta_0$$ is really just a datum.

In the pictures below, red is the reference signal I'd ideally recover, but it isn't typically available. The black is the corrupted version I am really dealing with. Here is a picture of a straightforward unwrapped case from a ~1Hz band with nominal frequency removed:

Here is a picture of a ~4Hz band where around $$t=48s$$ the black analyzed system misses a cycle and gets lodged around $$-2\pi$$.

• Your update is very helpful. It's not as bad as I thought and didn't fully realize those steps you described were your actual modulation. What is the magnitude and slope (freq) of those deterministic steps? Is the issue that the 2pi slip gets confused with the actual steps, and if not can't you just detect and subtract 2pi when such a slip is detected? (I'll be updating my answer in a few days once I have time again but wanted to get this comment out and let you know I'm interested). May 15, 2023 at 16:58

I first remove some possible points of confusion with wrapping,and then important assumptions since there can potentially be a lot of pitfalls with frequency and phase estimation, and then with those assumptions I will provide an optimum maximum-likelihood strategy for estimating the frequency and phase noise.

The OP has specifically stated that the samples received were wrapped at the 0 and $$2\pi$$ boundaries. This is different from the typical wrapping assumed by unwrap in the common tools (Python's Numpy, MATLAB, Octave). In those tools the assumed boundary is -$$\pi$$ to $$\pi$$ by default (and can be set to any other value). Changing the boundary is identical to changing the assumed phase offset $$\theta_0$$. Regardless of where the wrapping is set to occur, the phase offset term $$\theta_0$$ has no effect on reducing the amount of wrapping that was done in the measurement other than trivially the first sample in the equation given by the OP: wrapping is done on the difference between successive samples, and any such offset is removed in the difference operation. Generally consider arbitrary phase samples as a phase changing versus sample $$n$$ given as $$\phi[n]$$ and a phase offset given as $$\theta_0$$:

$$\theta[n] = \phi[n] + \theta_0$$ $$\theta[n]-\theta[n-1] = \phi[n]-\phi[n-1]$$

Given the wrapping is done on phase differences, we also see that if $$\omega$$ is known or can be easily approximated, then if the frequency could have been removed in the original phase measurement prior to the wrapping that was done, there will be significantly less discontinuities (knowing these kind of measurements, I assume that was actually done and the OP's $$\omega$$ is actually a small frequency offset from the frequency in the original signal, and the $$\delta \omega$$ term is a residual static frequency error to be determined.

I am also making the following simplifying assumptions which if not valid for the OP would modify the maximum likelihood estimation strategy that I will suggest:

• The samples are uniformly spaced in time by $$T = 1/f_s$$ where $$f_s$$ is the sampling rate.

• $$A(t)$$ has been properly removed prior to phase measurement and has had no effect on the measured results.

• $$\phi(t)$$ is "small angle" phase noise; there is not a deterministic phase modulation that is being determined, and once the linear phase slope (static frequency) is removed the rms phase error is $$\phi_{rms} < 1$$ radian.

• The capture duration is sufficiently short or the signal generation and signal measurement systems are locked to a common reference. The most important point here is that we are not yet affected by frequency drift that will inevitably occur when the reference of time or frequency for the signal generation system is different than the reference for the signal measurement system.

• The noise process is assumed to be "White FM": Phase Noise (the OP's $$\phi(t)$$) and frequency noise given as $$\omega(t)$$ are related by $$f(t) = d\phi(t)/dt$$. Given this relationship, the OP's expression is given entirely in terms of phase noise with a static frequency but that doesn't mean we do not have both frequency and phase noise (they are the same noise, it is just deciding which way we want to describe it). That said, the Power Spectral Density due to Frequency Fluctuations (commonly abbreviated $$S_f(f)$$ for one-sided spectrums, and $$\mathscr{L}_f(f)$$ for two-sided spectrums) will be constant if the noise process is White FM, and the Power Spectral Density due to Phase Fluctuations (commonly abbreviated $$S_\phi(f)$$ for one-sided spectrums, and $$\mathscr{L}_\phi(f)$$ for two-sided spectrums) will go down -20 dB/decade. $$\mathscr{L}_\phi(f)$$ is most commonly used to represent spectrums of phase noise, so given that, this whole paragraph reduces to: confirm the phase noise is going with a slope of -20 dB/decade (or less). For longer captures of unlocked source and measurement, the slope of the noise will inevitably get steeper (for example from -20 dB/decade to -30 dB /decade to -40 dB/decade as the noise process migrates from white FM noise, to flicker freq noise, to random walk frequency noise), which in simpler terms will mean our estimation techniques will diverge. If this turned out to be the case, it means our estimates would improve if we shortened the capture duration used so that this does apply.

Also to avoid confusion and as motivated by this recent post, since I will be using both discrete and continuous time references to frequency, I will use $$\Omega$$ to refer to continuous time frequency as given in units of radians/sec ($$\Omega = 2\pi F$$ where $$F$$ is the frequency in Hz) and $$\omega$$ as discrete time "normalized" frequency as given in units of radians/sample, with $$\omega = \Omega/f_s$$. (And $$\omega = 2\pi f$$ where $$f$$ is normalized frequency in units of cycles/sample).

Thus we start with the OP's phase measurements as samples of the following phase versus time function:

$$\theta(t) = (\hat{\Omega}+\Delta{\Omega})t + \phi(t) + \theta_0 \tag{1}\label{1}$$

I will remove $$\theta_0$$ from consideration since it has no effect on phase wrapping as previously explained, as well as no effect on the estimate the OP seeks for static frequency offset $$\Delta \omega$$ and phase noise $$\phi(t)$$.

Thus when sampled at $$t= nT$$, we have the samples given as:

$$\theta(nT) = (\hat{\Omega}+\Delta{\Omega})nT + \phi(nT) \tag{2}\label{2}$$

The sampling rate $$f_s = 1/T$$, and thus the samples in units of normalized frequency as a function of time given simply in samples is given as:

$$\theta[n] = (\hat{\omega}+\Delta{\omega})n + \phi[n] \tag{3}\label{3}$$

With that, my recommended algorithm to determine $$\Delta \omega$$, $$\phi[n]$$ is:

(OP has mentioned in the comments that there is a deterministic phase modulation involved, so this answer will likely be modified based on the details of that)

• B Yes I am tracking this. Most of the assumptions seem fine except the "no deterministic" part. $\phi(t)$ is a modulation of interest. I am trying to apportion phase between global frequency and smallish modulation $\phi(t)$ in the presence of noise. The $\phi(t)$ I want does seem to stay within the rms bounds you describe. May 13, 2023 at 16:10
• The part about "frequency noise and phase noise are the same" resonates. $\phi$ wanders within bounds but then wanders off on a linear path towards a new modolo-two-$\pi$ value, at which point it stays there. That looks like an episode of wrong frequency. When I look at the wrapped phase, I can see it has the same number of periods as the reference phase. May 13, 2023 at 16:13
• Re first comment: Ah ok- could you please update your question with more details about your specific phase modulation of interest? Even specific quantities such as sampling rate etc can help with context and narrow down the breadth of the response. Re second comment: It would help if you can Include a time plot of that in your post and any details as to what the system is that is providing that and its function. There's a lot of reasons why you might see what you describe, but at the moment it's a very large dart board where I fear it wouldn't yet help if I started suggesting any first. May 14, 2023 at 3:58