ENH: Phase unwrapping generalized to arbitrary interval size

[scimax]

ENH: Phase unwrapping generalized to arbitrary interval size
Generalizing phase unwrapping to arbitrary interval size, such as 360 for phase in degree instead of radian. Also integer unwrapping is supported, but the return values are floats.

• closes Parameterize unwrap to allow for different ranges (not just -pi, pi) #14876
• fixes wrong unwrapping of integers in ENH: Adding parameters to unwrap #14877
• proper naming compared to fix bugs in unwrap function #16977

To-Dos:

• I have to include the tests added in ENH: Adding parameters to unwrap #14877
• Release Note
• Still to clarify: (min_val, max_val) or interval_size
• Return integers if input argument is an integer array

But I couldn't figure out yet how to unwrap integers, such that they stay integers. The only thing I came up with is checking the type of the input array at the beginning and change the output type in the return statement. So it's not really an improvement in the algorithm.

[tatsuoka999]

This change is exactly what I needed. I hope this change will be merged.
Or, for compatibility, it may be added as a new function such as unwrap2 ().

[eric-wieser]

Thanks for the patch.

This looks like a correct (or at least, as correct as before) generalization unlike the other patches, but I suspect we want to keep the API proposed in #14877

[eric-wieser]

Please delete these commented lines

[eric-wieser]

I think the API in #14877 of having separate low and high bounds makes more sense. That would also solve your problem of not returning integers, as you introduce floats by doing the /2

[scimax]

ok, in my first version I used intervalSize = max_val - min_val but then I thought how does the following make sense:

>>> unwrap([-10, -9, -8, -7, -6, -5, -11, -10, -8, -7, -6, -10, -9, -8, -7 ], min_val=2, max_val=8)
array([-10.,  -9.,  -8.,  -7.,  -6.,  -5.,  -5.,  -4.,  -2.,  -1.,   0.,
         2.,   3.,   4.,   5.])

It works as expected, just correcting the jumps. But the input is not in the given range [min_val, max_val] and, in my opinion, it just introduces unnecessary comparisons if you check whether the input data is within the bounds. But yes, if that's wanted, I can simply add the line again and change the signature back.

One concern regarding the floats: I agree that by /2 the integers become floats. But because of the modulo operation you have to shift the data by half the interval size. If the interval size is odd, the half interval size is definitely float. I will rethink this, but in my opinion the only solution to this would be case distinction. So in case of integers you could check for odd or even interval size, and then make the modulo operation work with both cases separately. But for the float case it will have to shift by half the interval size.

So (min_val, max_val) or interval_size ?

[scimax]

Regarding your requirement in #16977 I have to add that the term unwrapping is usually not used for correcting offsets. I might be wrong but from my experience in rf design I haven't seen a use of this. Is there a field where offset correction is meant?
Quoting from wikipedia on instantenous phase:

Discontinuities can then be removed by adding 2π whenever Δφ[n] ≤ −π, and subtracting 2π whenever Δφ[n] > π. That allows φ[n] to accumulate without limit and produces an unwrapped instantaneous phase.

[eric-wieser]

My understanding is that the intent of min / max is to specify the range of values to use for

• The first value
• The first value after a discontinuity, maybe?

[scimax]

What do you think about this solution?

>>> a = np.array([
    [1,  2, 3,  4,  1, 2, 3, 4],
    [17,18,19, 19, 18,17,20,19],
    [-2,-3,-4,-5, -2, -3, -4, -5]])
>>> unwrap(a, min_val=11, max_val=15)
array([[13., 14., 15., 16., 17., 18., 19., 20.],
       [13., 14., 15., 15., 14., 13., 12., 11.],
       [14., 13., 12., 11., 10.,  9.,  8.,  7.]])

It's a hybrid solution including both options. I wanted to keep the interval_size option without correcting the offset as it would break backward compatibility otherwise. I can give you an example later, in which I unwrapped over one axis first, and then over the second axis. Thereby, during the second unwrapping the first value was already outside of [-pi, pi) as it was the output of an unwrapping.

[scimax]

With backward compatibility I mean that the old function did not care whether or not the data is from the given interval. It only cared about phase differences. I will try to construct an example similar with the one below, but with integer data.

So suppose you have data which is depending on two parameters, f(x,y) = x*y, but the data wraps against both axis. In fact, instead of f(x,y) the data is z(x,y) = f(x,y) mod p. Here I choose p=11. So how do we find f(x,y)? Simple answer in terms of math, it's not possible. But assuming that we have no big jumps we can correct the data up to an unknown offset. In fact, in the example the array f_arr is almost reconstructed, up to a missing offset of 11.

>>> 
>>> x_arr = np.arange(-5,6)
>>> y_arr = np.arange(-3,4)
>>> X, Y = np.meshgrid(x_arr, y_arr)
>>> f_arr = X*Y
>>> f_arr
array([[ 15,  12,   9,   6,   3,   0,  -3,  -6,  -9, -12, -15],
       [ 10,   8,   6,   4,   2,   0,  -2,  -4,  -6,  -8, -10],
       [  5,   4,   3,   2,   1,   0,  -1,  -2,  -3,  -4,  -5],
       [  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0],
       [ -5,  -4,  -3,  -2,  -1,   0,   1,   2,   3,   4,   5],
       [-10,  -8,  -6,  -4,  -2,   0,   2,   4,   6,   8,  10],
       [-15, -12,  -9,  -6,  -3,   0,   3,   6,   9,  12,  15]])
>>> z_arr = np.mod(f_arr, 11)
>>> z_arr
array([[1, 2, 3, 4, 0, 1, 2],
       [2, 4, 1, 3, 0, 2, 4],
       [3, 1, 4, 2, 0, 3, 1],
       [4, 3, 2, 1, 0, 4, 3],
       [0, 0, 0, 0, 0, 0, 0]], dtype=int32)
>>> # Unwrap against x axis
>>> interval_size=11
>>> dd = np.diff(z_arr, axis=1)
>>> ddmod = np.mod(dd + interval_size/2, interval_size) - interval_size/2
>>> np.copyto(ddmod, interval_size/2, where=(ddmod == -interval_size/2) & (dd > 0))
>>> ph_correct = ddmod - dd
>>> np.copyto(ph_correct, 0, where=abs(dd) < interval_size/2)
>>> up = np.array(z_arr, copy=True, dtype='d')
>>> up[:, 1:] = z_arr[:, 1:] + ph_correct.cumsum(1) 
>>> z_arr_2 = up
>>> z_arr_2
array([[  4.,   1.,  -2.,  -5.,  -8., -11., -14., -17., -20., -23., -26.],
       [ 10.,   8.,   6.,   4.,   2.,   0.,  -2.,  -4.,  -6.,  -8., -10.],
       [  5.,   4.,   3.,   2.,   1.,   0.,  -1.,  -2.,  -3.,  -4.,  -5.],
       [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
       [  6.,   7.,   8.,   9.,  10.,  11.,  12.,  13.,  14.,  15.,  16.],
       [  1.,   3.,   5.,   7.,   9.,  11.,  13.,  15.,  17.,  19.,  21.],
       [  7.,  10.,  13.,  16.,  19.,  22.,  25.,  28.,  31.,  34.,  37.]])
>>> # Unwrap against y axis (axis 0 of the array)
>>> dd = np.diff(z_arr_2, axis=0)
>>> ddmod = np.mod(dd + interval_size/2, interval_size) - interval_size/2
>>> np.copyto(ddmod, interval_size/2, where=(ddmod == -interval_size/2) & (dd > 0))
>>> ph_correct = ddmod - dd
>>> np.copyto(ph_correct, 0, where=abs(dd) < interval_size/2)
>>> up = np.array(z_arr_2, copy=True, dtype='d')
>>> up[1:, :] = z_arr_2[1:, :] + ph_correct.cumsum(0) 
>>> up
array([[  4.,   1.,  -2.,  -5.,  -8., -11., -14., -17., -20., -23., -26.],
       [ -1.,  -3.,  -5.,  -7.,  -9., -11., -13., -15., -17., -19., -21.],
       [ -6.,  -7.,  -8.,  -9., -10., -11., -12., -13., -14., -15., -16.],
       [-11., -11., -11., -11., -11., -11., -11., -11., -11., -11., -11.],
       [-16., -15., -14., -13., -12., -11., -10.,  -9.,  -8.,  -7.,  -6.],
       [-21., -19., -17., -15., -13., -11.,  -9.,  -7.,  -5.,  -3.,  -1.],
       [-26., -23., -20., -17., -14., -11.,  -8.,  -5.,  -2.,   1.,   4.]])

>>> 

Now if you check up + 11 == f_arr all entries will be True. So the unwraping will give us f(x,y) up to an unknown offset. As you can see in the second unwraping the input is not at all from the interval [0, 11) anymore. But this is not important for the unwraping.

I can really recommend taking a look at the different intermediate steps. I tried to make a shorter version. But unwraping is quite tricky. So this is in my opinion the smallest 2D unwraping example where you can actually see what is going on. If the input arrays are smaller it's not clear anymore.

Finally, this is what I mean by backwards compatibility. The current version is doing exactly this. Apart from the fact that it's only working with an interval size of 2*pi, not 11.

I'm sorry for the late response. I didn't expect this to consume so much time.

[scimax]

as it would break backward compatibility otherwise.

Backwards compatibility with what? My understanding was that unwrap(..., interval_size=i) would be equivalent to unwrap(..., min_val=i/2, max_val=i/2) (

Adding a short notice to unwrap(..., interval_size=i) being equivalent to unwrap(..., min_val=i/2, max_val=i/2). This is not really true since it's not really clear to me, what you would expect from this:

>>> a = np.array([1,  2, 3,  4,  1, 2, 3, 4])
>>> unwrap(a, min_val=11, max_val=15)

What I have implemented now (see this comment) would map the first value to the interval [min_val, max_val). Regarding the fact that the current unwrap function is not map the first value to the interval [-pi, pi), I would say this breaks backward compatibility. My preferred way would be to get rid of min_val and max_val as I mentioned at the beginning. Unwrapping does not care about offsets. But I'm definitely open for other suggestions.

[eric-wieser]

Regarding the fact that the current unwrap function is not map the first value to the interval [-pi, pi),

Oh, it doesn't? I assumed that today np.unwrap([7*pi]) would give [pi]. I'm not at a terminal so can't check...

[scimax]

>>> np.unwrap([7*np.pi])
array([21.99114858])
>>> 

which is 7*pi. (Currently I'm using version 1.19.1, but unwrap hasn't changed since then.)

[eric-wieser]

Thanks - I think that invalidates much of my point then.

[charris]

Needs a release note and tests.

[scimax]

Example with quasi 2d unwrapping where the input data is not in the interval [-pi,pi) anymore (sorry, code is not very pythonic):

from scipy.constants import c
import matplotlib.pyplot as plt

freqs = np.linspace(220, 330, 111)
l= np.linspace(10,11, 11)

phi = -2* 2*np.pi*np.outer(l, freqs*1e6/c)
phi = np.mod(phi + np.pi, 2*np.pi) - np.pi

plt.figure()
plt.plot(freqs, phi[0,:], label=r"l={:2.2f}".format(l[0]))
plt.plot(freqs, phi[l.shape[0]//2,:], label=r"l={:2.2f}".format(l[l.shape[0]//2]))
plt.plot(freqs, phi[-1,:], label=r"l={:2.2f}".format(l[-1]))
plt.grid()
plt.xlabel("f [MHz]")
plt.ylabel(r"$\varphi$ [rad]")
plt.legend()
plt.ylim((-4,6))

After unwrapping over frequency you can see that the phase still wraps along the parameter l:

phi_un = np.unwrap(phi, axis=1)

plt.figure()
for i in range(9):
    plt.plot(freqs, phi_un[i,:], label=r"l={:2.2f}".format(l[i]))

plt.grid()
plt.xlabel("f [MHz]")
plt.ylabel(r"$\varphi$ [rad]")
plt.legend()
plt.ylim((-10,3))
plt.xlim((220, 240))

Now phases at around 240 MHz are clearly outside of [-pi,pi). But since unwrapping only cares about differences, not about offsets, it can still be unwrapped:

phi_un_2d = np.unwrap(phi_un, axis=0)

plt.figure()
for i in range(9):
    plt.plot(freqs, phi_un_2d[i,:], label=r"l={:2.2f}".format(l[i]))
plt.grid()
plt.xlabel("f [MHz]")
plt.ylabel(r"$\varphi$ [rad]")
plt.legend()
plt.ylim((-10,3))
plt.xlim((220, 240))

[eric-wieser]

I need to think a bit more about min/max vs interval_size. I don't fully understand your example yet, anything you can do to simplify it would greatly help. Avoiding pi and using rational wrapping bounds would help a lot. Admittedly I haven't looked at it yet.

In the mean-time, some small comments

[eric-wieser]

Thoughts on using the shorter name period?

[scimax]

Yeah, sounds also good! I just picked the first name which came to my mind! Shall I rename it?

[eric-wieser]

Yeah, I think so - we use period already in np.interp, and as far as I can tell the meaning is the same.

[eric-wieser]

We should fix this to return integers before merging, will think about the best way to handle it.

[eric-wieser]

These doctest outputs now need updating to integers

[scimax]

Just this comment left to address: #16987 (comment).

I cleaned up the whitespace myself

@eric-wieser, I'm sorry, I got distracted from this PR a little bit. Do you remember to which comment you wanted to refer? The link does not jump to a certain comment.

[eric-wieser]

Wow, that's frustrating. I'll have another look at that on desktop tomorrow, but I also can't work out where that permalink was supposed to go

[scimax]

Wow, that's frustrating. I'll have another look at that on desktop tomorrow, but I also can't work out where that permalink was supposed to go

Thanks a lot! I will try to go through the comments as well tonight, and see, if something is still missing.

[eric-wieser]

I worked it out, I somehow trimmed a digit off the URL. I've corrected that link now.

[eric-wieser]

I've tweaked the docs slightly, but I think this is ready to go in once CI is green.

[scimax]

I've tweaked the docs slightly, but I think this is ready to go in once CI is green.

Hi @eric-wieser , I just came across and was wondering what the status ist. Do you have an idea, why the azure-pipelines fail and the circleCI tests? I can't take a look at the logs.

[charris]

close/reopen

[charris]

There are 7 lines that are too long, hit details for the lint test. They need to be < 80 characters long.

[charris]

Thanks @scimax and @eric-wieser .

[dpesios]

Hello,

It seems to me that unwrap is returning degrees (or equivalent) instead of radians. Could please someone confirm this ? Maybe, @scimax ?

I'm not really familiar with the algorithm you implemented. The impression I have stems from using the Euler formula to obtain the real and imaginary component of a phasor. Judging by the result, in order to obtain the correct one, I have to convert into radians its output. I'm using floats as input.

Thanks in advance.

[scimax]

Hi @dpesios, I'm not sure, what you mean. Could you give a minimal example? By default, it unwraps radians. The goal of this PR was to add the option the unwrapping of degrees, but it's not the default.