[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