DSP: Noisy Modulo Unwrapping Technique

Overview

This page will discuss a simple algorithm to unwrap a noisy discrete real waveform transformed under a modulo operation. The algorithm can be implemented as a real-time algorithm on a waveform stream.

Setup

Let \vec v = (a_1,a_2, ..., a_N) be a discrete real waveform sampled from a continuous waveform and for a fixed real constant M>0, consider its modulo-M waveform \vec w = (b_1,b_2, ..., b_N) such that b_k = a_k\:mod\:M. The algorithm seeks to recover the unwrapped waveform \vec v from the wrapped waveform \vec w up to a constant offset.

Algorithm – (Noiseless Version)

For any given integer L, define LM_s = [LM+\frac{sM}{3}, LM+\frac{(s+1)M}{3}). Note that LM_s forms a disjoint partition of the entire real line. This means that for any given term a_k, there exist a unique modulo level L(a_k) and the subdivision level s(a_k) such that x \in L(a_k)M_{s(a_k)}. Furthermore, if L(a_k) and b_k are known, the original a_k can be recovered via the equation a_k=L(a_k)M + b_k.

In the absence of noise, the modulo level can either change by one level, exactly when the subdivision levels corresponding to two contiguous terms change from between 0 and 2 across two modulo levels. And the continuity of the source waveform prevents the subdivision level from jumping up or down by two steps at once.

This observation allows the correct assignment of the modulo level for any given term up to a constant offset using the following recursive rule:

L_{0} = 0

L_{k+1}=L_k+1, if s_{k}=2 and s_{k+1}=0

L_{k+1}=L_k-1, if s_{k}=0 and s_{k+1}=2

L_{k+1}=L_k, otherwise

The rule can be compactly rewritten as:

L_{0} = 0

L_{k+1}=L_k + \lfloor \frac{s_{k} - s_{k+1}}{2} \rfloor

In turn, this would allow the correct unwrapping of the modulo waveform up to a constant offset via the relationship a_k=L(a_k)M + b_k.

Algorithm – (Noisy Version)

The level-recovery method presented above may be unsuitable in the presence of a high noise level. If the noise corresponding to a specific term is large enough to make the term jump multiple subdivision levels, the level-recovery method will fail by creating a discontinuity point in the recovered waveform. The modified algorithm described below can safeguard against this type of shortcoming.

Pick an integer-valued scanning window size r>0 and initialize the states as follows:

W_0 = (b_1,b_2,....,b_r)

L_{W_0} = 0

and define s_{W_k} to be the majority subdivision level of W_k.

Then, recover the level sequences of the terms recursively as follows until the window is shifted all the way to the end (k \leq N-r-1):

  1. L_{k}=L_{W_k} + \lfloor \frac{s_{W_k} - s_{k}}{2} \rfloor
  2. Define W_{k+1} by removing the first element of W_{k} and appending b_{k+r+1}.
  3. L_{W_{k+1}}=L_{W_k} + \lfloor \frac{s_{W_k} - s_{W_{k+1}}}{2} \rfloor

Once the window is shifted all the way to the right, N-r \leq k \leq N, recover the level sequences of the terms by using the last window states:

L_{k}=L_{W_{N-r-1}} + \lfloor \frac{s_{W_{N-r-1}} - s_{k}}{2} \rfloor

The window size should be carefully chosen by weighing the following trade-off:

  1. It should be small enough that as the window is shifted from the beginning of the waveform (with the left edge of the window aligned with the first term) to the end of the waveform (with the right edge of the window aligned with the last term), the window correctly traces the general trajectory of the original waveform. If the window size is too wide, it may miss subdivision jumps of the original waveform.
  2. It should be chosen large enough that the windowed mean subdivision level prevents noisy terms from introducing discontinuity.

Leave a comment