diff --git a/doc/release/upcoming_changes/16987.improvement.rst b/doc/release/upcoming_changes/16987.improvement.rst new file mode 100644 index 00000000000..dc592a06840 --- /dev/null +++ b/doc/release/upcoming_changes/16987.improvement.rst @@ -0,0 +1,17 @@ +Arbitrary ``period`` option for `numpy.unwrap` +---------------------------------------------- +The size of the interval over which phases are unwrapped is no longer restricted to ``2 * pi``. +This is especially useful for unwrapping degrees, but can also be used for other intervals. + +.. code:: python + + >>> phase_deg = np.mod(np.linspace(0,720,19), 360) - 180 + >>> phase_deg + array([-180., -140., -100., -60., -20., 20., 60., 100., 140., + -180., -140., -100., -60., -20., 20., 60., 100., 140., + -180.]) + + >>> unwrap(phase_deg, period=360) + array([-180., -140., -100., -60., -20., 20., 60., 100., 140., + 180., 220., 260., 300., 340., 380., 420., 460., 500., + 540.]) diff --git a/numpy/lib/function_base.py b/numpy/lib/function_base.py index 42ea8c7a785..651343a2c2e 100644 --- a/numpy/lib/function_base.py +++ b/numpy/lib/function_base.py @@ -1482,26 +1482,40 @@ def angle(z, deg=False): return a -def _unwrap_dispatcher(p, discont=None, axis=None): +def _unwrap_dispatcher(p, discont=None, axis=None, *, period=None): return (p,) @array_function_dispatch(_unwrap_dispatcher) -def unwrap(p, discont=pi, axis=-1): - """ - Unwrap by changing deltas between values to 2*pi complement. +def unwrap(p, discont=None, axis=-1, *, period=2*pi): + r""" + Unwrap by taking the complement of large deltas with respect to the period. + + This unwraps a signal `p` by changing elements which have an absolute + difference from their predecessor of more than ``max(discont, period/2)`` + to their `period`-complementary values. - Unwrap radian phase `p` by changing absolute jumps greater than - `discont` to their 2*pi complement along the given axis. + For the default case where `period` is :math:`2\pi` and is `discont` is + :math:`\pi`, this unwraps a radian phase `p` such that adjacent differences + are never greater than :math:`\pi` by adding :math:`2k\pi` for some + integer :math:`k`. Parameters ---------- p : array_like Input array. discont : float, optional - Maximum discontinuity between values, default is ``pi``. + Maximum discontinuity between values, default is ``period/2``. + Values below ``period/2`` are treated as if they were ``period/2``. + To have an effect different from the default, `discont` should be + larger than ``period/2``. axis : int, optional Axis along which unwrap will operate, default is the last axis. + period: float, optional + Size of the range over which the input wraps. By default, it is + ``2 pi``. + + .. versionadded:: 1.21.0 Returns ------- @@ -1514,9 +1528,9 @@ def unwrap(p, discont=pi, axis=-1): Notes ----- - If the discontinuity in `p` is smaller than ``pi``, but larger than - `discont`, no unwrapping is done because taking the 2*pi complement - would only make the discontinuity larger. + If the discontinuity in `p` is smaller than ``period/2``, + but larger than `discont`, no unwrapping is done because taking + the complement would only make the discontinuity larger. Examples -------- @@ -1526,19 +1540,44 @@ def unwrap(p, discont=pi, axis=-1): array([ 0. , 0.78539816, 1.57079633, 5.49778714, 6.28318531]) # may vary >>> np.unwrap(phase) array([ 0. , 0.78539816, 1.57079633, -0.78539816, 0. ]) # may vary - + >>> np.unwrap([0, 1, 2, -1, 0], period=4) + array([0, 1, 2, 3, 4]) + >>> np.unwrap([ 1, 2, 3, 4, 5, 6, 1, 2, 3], period=6) + array([1, 2, 3, 4, 5, 6, 7, 8, 9]) + >>> np.unwrap([2, 3, 4, 5, 2, 3, 4, 5], period=4) + array([2, 3, 4, 5, 6, 7, 8, 9]) + >>> phase_deg = np.mod(np.linspace(0 ,720, 19), 360) - 180 + >>> np.unwrap(phase_deg, period=360) + array([-180., -140., -100., -60., -20., 20., 60., 100., 140., + 180., 220., 260., 300., 340., 380., 420., 460., 500., + 540.]) """ p = asarray(p) nd = p.ndim dd = diff(p, axis=axis) + if discont is None: + discont = period/2 slice1 = [slice(None, None)]*nd # full slices slice1[axis] = slice(1, None) slice1 = tuple(slice1) - ddmod = mod(dd + pi, 2*pi) - pi - _nx.copyto(ddmod, pi, where=(ddmod == -pi) & (dd > 0)) + dtype = np.result_type(dd, period) + if _nx.issubdtype(dtype, _nx.integer): + interval_high, rem = divmod(period, 2) + boundary_ambiguous = rem == 0 + else: + interval_high = period / 2 + boundary_ambiguous = True + interval_low = -interval_high + ddmod = mod(dd - interval_low, period) + interval_low + if boundary_ambiguous: + # for `mask = (abs(dd) == period/2)`, the above line made + # `ddmod[mask] == -period/2`. correct these such that + # `ddmod[mask] == sign(dd[mask])*period/2`. + _nx.copyto(ddmod, interval_high, + where=(ddmod == interval_low) & (dd > 0)) ph_correct = ddmod - dd _nx.copyto(ph_correct, 0, where=abs(dd) < discont) - up = array(p, copy=True, dtype='d') + up = array(p, copy=True, dtype=dtype) up[slice1] = p[slice1] + ph_correct.cumsum(axis) return up diff --git a/numpy/lib/tests/test_function_base.py b/numpy/lib/tests/test_function_base.py index eb2fc3311ac..f881ddf0079 100644 --- a/numpy/lib/tests/test_function_base.py +++ b/numpy/lib/tests/test_function_base.py @@ -1756,6 +1756,24 @@ def test_simple(self): assert_array_equal(unwrap([1, 1 + 2 * np.pi]), [1, 1]) # check that unwrap maintains continuity assert_(np.all(diff(unwrap(rand(10) * 100)) < np.pi)) + + def test_period(self): + # check that unwrap removes jumps greater that 255 + assert_array_equal(unwrap([1, 1 + 256], period=255), [1, 2]) + # check that unwrap maintains continuity + assert_(np.all(diff(unwrap(rand(10) * 1000, period=255)) < 255)) + # check simple case + simple_seq = np.array([0, 75, 150, 225, 300]) + wrap_seq = np.mod(simple_seq, 255) + assert_array_equal(unwrap(wrap_seq, period=255), simple_seq) + # check custom discont value + uneven_seq = np.array([0, 75, 150, 225, 300, 430]) + wrap_uneven = np.mod(uneven_seq, 250) + no_discont = unwrap(wrap_uneven, period=250) + assert_array_equal(no_discont, [0, 75, 150, 225, 300, 180]) + sm_discont = unwrap(wrap_uneven, period=250, discont=140) + assert_array_equal(sm_discont, [0, 75, 150, 225, 300, 430]) + assert sm_discont.dtype == wrap_uneven.dtype class TestFilterwindows: