[Numpy-discussion] phase unwrapping (1d)

Neal Becker ndbecker2@gmail....
Mon Jan 14 08:39:34 CST 2013


This code should explain all:
--------------------------------
import numpy as np
arg = np.angle

def nint (x):
    return int (x + 0.5) if x >= 0 else int (x - 0.5)

def unwrap (inp, y=np.pi, init=0, cnt=0):
    o = np.empty_like (inp)
    prev_o = init
    for i in range (len (inp)):
        o[i] = cnt * 2 * y + inp[i]
        delta = o[i] - prev_o

        if delta / y > 1 or delta / y < -1:
            n = nint (delta / (2*y))
            o[i] -= 2*y*n
            cnt -= n

        prev_o = o[i]

    return o
            

u = np.linspace (0, 400, 100) * np.pi/100
v = np.cos (u) + 1j * np.sin (u)
plot (arg(v))
plot (arg(v) + arg (v))
plot (unwrap (arg (v)))
plot (unwrap (arg (v) + arg (v)))
-------------------------------

Pierre Haessig wrote:

> Hi Neal,
> 
> Le 11/01/2013 16:40, Neal Becker a écrit :
>> I wanted to be able to handle the case of
>>
>> unwrap (arg (x1) + arg (x2))
>>
>> Here, phase can change by more than 2pi.
> It's not clear to me what you mean by "change more than 2pi" ? Do you
> mean that the consecutive points of in input can increase by more than
> 2pi ? If that's the case, I feel like there is no a priori information
> in the data to detect such a "giant leap".
> 
> Also, I copy-paste here for reference the numpy.wrap code from [1] :
> 
> def unwrap(p, discont=pi, axis=-1):
>     p = asarray(p)
>     nd = len(p.shape)
>     dd = diff(p, axis=axis)
>     slice1 = [slice(None, None)]*nd # full slices
>     slice1[axis] = slice(1, None)
>     ddmod = mod(dd+pi, 2*pi)-pi
>     _nx.copyto(ddmod, pi, where=(ddmod==-pi) & (dd > 0))
>     ph_correct = ddmod - dd;
>     _nx.copyto(ph_correct, 0, where=abs(dd)<discont)
>     up = array(p, copy=True, dtype='d')
>     up[slice1] = p[slice1] + ph_correct.cumsum(axis)
>     return up
> 
> I don't know why it's too slow though. It looks well vectorized.
> 
> Coming back to your C algorithm, I'm not C guru so that I don't have a
> clear picture of what it's doing. Do you have a Python prototype ?
> 
> Best,
> Pierre
> 
> [1]
> https://github.com/numpy/numpy/blob/master/numpy/lib/function_base.py#L1117




More information about the NumPy-Discussion mailing list