Fixing phase unwrapping errors in Numpy

Question

I have a series of unwrapped phases, with some unwrapping errors that consist of a jump of +/- a multiple of Pi:

import numpy
a = numpy.array([0.5, 1.3, 2.4, 10.3, 10.8, 10.2, 7.6, 3.2, 2.9])

In this example there is a first jump of 2 cycles between 2.4 and 10.3, and a jump of -1 cycle between 7.6 and 3.2. I want to remove the jumps. The catch is that when you remove a jump, you need to increase or decrease the rest of the series accordingly, not just the value where the jump occurs.

Is there a cleaner way (no/less loops, faster) of doing this:

jumpsexist = 1
while jumpsexist:
    # Look for absolute differences greater than Pi
    jump = numpy.abs((numpy.roll(a,-1) -a)) > numpy.pi
    if jump[:-1].any():
        # Find the index of the first jump
        jumpind = numpy.argmax(jump) + 1
        # Calculate the number of cycles in that jump
        cycles = ((a[jumpind] - a[jumpind- 1]) / numpy.pi).astype("Int8")
        # Remove the cycles
        a[jumpind:] -= cycles * numpy.pi
    else:
        break

Note that the variable is not called "jump sexist", but "jumps exist". I'll take "s words" for 200. — Benjamin, Apr 12 '12 at 14:56
In keeping with your poorly chosen variable names, you might find the numpy function `cumsum` to be useful here. That's cumulative sum not ... nevermind. — Hooked, Apr 12 '12 at 15:00
On a side note, phase unwrapping in the presence of noise is an open research problem... You wouldn't happen to be working on InSAR data, would you? — Joe Kington, Apr 12 '12 at 15:47

score 9 · Accepted Answer · answered Apr 17 '12 at 12:37

9

NumPy offers the function numpy.unwrap() for phase unwrapping. With the default parameter values, it will correct an array of phases modulo 2π such that all jumps are less than or equal to π:

>>> a = numpy.array([0.5, 1.3, 2.4, 10.3, 10.8, 10.2, 7.6, 3.2, 2.9])
>>> numpy.unwrap(a)
array([ 0.5       ,  1.3       ,  2.4       ,  4.01681469,  4.51681469,
        3.91681469,  1.31681469,  3.2       ,  2.9       ])

answered Apr 17 '12 at 12:37

Sven Marnach

574,206
118
941
841

Changed this to accepted answer. Although not clear in my original question, this is really what I was after. – Benjamin Apr 27 '12 at 12:17
@Benjamin: I was suspecting you are – that's why I was pressing you to state your requirements clearly. :) – Sven Marnach Apr 27 '12 at 16:20

talonmies · Answer 2 · 2012-04-12T15:31:31.177

How about this:

import numpy as np 
a = np.array([0.5, 1.3, 2.4, 10.3, 10.8, 10.2, 7.6, 3.2, 2.9])
d = np.diff(a)/np.pi
b = np.empty_like(a)
b[0] = a[0]
b[1:] = a[1:]-(np.floor(np.abs(d))*np.sign(d)).cumsum()*np.pi

which gives:

In [40]: print a
[  0.5   1.3   2.4  10.3  10.8  10.2   7.6   3.2   2.9]

In [41]: print b
[ 0.5         1.3         2.4         4.01681469  4.51681469  3.91681469
  1.31681469  0.05840735 -0.24159265]

Here d holds the signed magntiude of the "jumps", and cumulative summation of the appropriately truncated "jumps" is the mulitple of pi which needs to be removed/added to each sucessive element of the series.

Is that what you meant?

Right on. For readability, `d = (numpy.diff(phases) / numpy.pi).astype("Int8")` and `b[1:] = a[1:] - numpy.cumsum(d) * numpy.pi ` might be clearer. — Benjamin, Apr 12 '12 at 17:28

Fixing phase unwrapping errors in Numpy

2 Answers2