diff --git a/doc/source/conf.py b/doc/source/conf.py index 5254eb8e8..b830351c5 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -36,6 +36,7 @@ extensions = ['sphinx.ext.doctest', 'sphinx.ext.autodoc', 'sphinx.ext.todo', 'sphinx.ext.extlinks', 'sphinx.ext.mathjax', 'sphinx.ext.autosummary', 'numpydoc', + 'sphinx.ext.intersphinx', 'matplotlib.sphinxext.plot_directive'] # Add any paths that contain templates here, relative to this directory. @@ -224,3 +225,9 @@ plot_formats = [('png', 96), 'pdf'] plot_html_show_formats = False plot_html_show_source_link = False + +# -- Options for intersphinx extension --------------------------------------- + +# Intersphinx to get Numpy and other targets +intersphinx_mapping = { + 'numpy': ('https://docs.scipy.org/doc/numpy/', None)} diff --git a/doc/source/pyplots/plot_boundary_modes.py b/doc/source/pyplots/plot_boundary_modes.py index cbb5f6eaa..940822f76 100644 --- a/doc/source/pyplots/plot_boundary_modes.py +++ b/doc/source/pyplots/plot_boundary_modes.py @@ -5,7 +5,7 @@ In practice, which signal extension mode is beneficial will depend on the signal characteristics. For this particular signal, some modes such as -"periodic", "antisymmetric" and "zeros" result in large discontinuities that +"periodic", "antisymmetric" and "zero" result in large discontinuities that would lead to large amplitude boundary coefficients in the detail coefficients of a discrete wavelet transform. """ @@ -28,5 +28,5 @@ boundary_mode_subplot(x, 'periodization', axes[5], symw=False) boundary_mode_subplot(x, 'smooth', axes[6], symw=False) boundary_mode_subplot(x, 'constant', axes[7], symw=False) -boundary_mode_subplot(x, 'zeros', axes[8], symw=False) +boundary_mode_subplot(x, 'zero', axes[8], symw=False) plt.show() diff --git a/doc/source/ref/signal-extension-modes.rst b/doc/source/ref/signal-extension-modes.rst index f8e466bb4..e05f2cce1 100644 --- a/doc/source/ref/signal-extension-modes.rst +++ b/doc/source/ref/signal-extension-modes.rst @@ -136,3 +136,15 @@ periodization per N/A antisymmetric asym, asymh N/A antireflect asymw reflect, reflect_type='odd' ================== ============= =========================== + +Padding using PyWavelets Signal Extension Modes - ``pad`` +--------------------------------------------------------- + +.. autofunction:: pad + +Pywavelets provides a function, :func:`pad`, that operate like +:func:`numpy.pad`, but supporting the PyWavelets signal extension modes +discussed above. For efficiency, the DWT routines in PyWavelets do not +expclitly create padded signals using this function. It can be used to manually +prepad signals to reduce boundary effects in functions such as :func:`cwt` and +:func:`swt` that do not currently support all of these signal extension modes. diff --git a/pywt/_doc_utils.py b/pywt/_doc_utils.py index 20c3fafbb..ee906aeab 100644 --- a/pywt/_doc_utils.py +++ b/pywt/_doc_utils.py @@ -4,8 +4,10 @@ import numpy as np from matplotlib import pyplot as plt +from ._dwt import pad + __all__ = ['wavedec_keys', 'wavedec2_keys', 'draw_2d_wp_basis', - 'draw_2d_fswavedecn_basis', 'pad', 'boundary_mode_subplot'] + 'draw_2d_fswavedecn_basis', 'boundary_mode_subplot'] def wavedec_keys(level): @@ -149,63 +151,6 @@ def draw_2d_fswavedecn_basis(shape, levels, fmt='k', plot_kwargs={}, ax=None, return fig, ax -def pad(x, pad_widths, mode): - """Extend a 1D signal using a given boundary mode. - - Like numpy.pad but supports all PyWavelets boundary modes. - """ - if np.isscalar(pad_widths): - pad_widths = (pad_widths, pad_widths) - - if x.ndim > 1: - raise ValueError("This padding function is only for 1D signals.") - - if mode in ['symmetric', 'reflect']: - xp = np.pad(x, pad_widths, mode=mode) - elif mode in ['periodic', 'periodization']: - if mode == 'periodization' and x.size % 2 == 1: - raise ValueError("periodization expects an even length signal.") - xp = np.pad(x, pad_widths, mode='wrap') - elif mode == 'zeros': - xp = np.pad(x, pad_widths, mode='constant', constant_values=0) - elif mode == 'constant': - xp = np.pad(x, pad_widths, mode='edge') - elif mode == 'smooth': - xp = np.pad(x, pad_widths, mode='linear_ramp', - end_values=(x[0] + pad_widths[0] * (x[0] - x[1]), - x[-1] + pad_widths[1] * (x[-1] - x[-2]))) - elif mode == 'antisymmetric': - # implement by flipping portions symmetric padding - npad_l, npad_r = pad_widths - xp = np.pad(x, pad_widths, mode='symmetric') - r_edge = npad_l + x.size - 1 - l_edge = npad_l - # width of each reflected segment - seg_width = x.size - # flip reflected segments on the right of the original signal - n = 1 - while r_edge <= xp.size: - segment_slice = slice(r_edge + 1, - min(r_edge + 1 + seg_width, xp.size)) - if n % 2: - xp[segment_slice] *= -1 - r_edge += seg_width - n += 1 - - # flip reflected segments on the left of the original signal - n = 1 - while l_edge >= 0: - segment_slice = slice(max(0, l_edge - seg_width), l_edge) - if n % 2: - xp[segment_slice] *= -1 - l_edge -= seg_width - n += 1 - elif mode == 'antireflect': - npad_l, npad_r = pad_widths - xp = np.pad(x, pad_widths, mode='reflect', reflect_type='odd') - return xp - - def boundary_mode_subplot(x, mode, ax, symw=True): """Plot an illustration of the boundary mode in a subplot axis.""" @@ -236,7 +181,7 @@ def boundary_mode_subplot(x, mode, ax, symw=True): left -= 0.5 step = len(x) rng = range(-2, 4) - if mode in ['smooth', 'constant', 'zeros']: + if mode in ['smooth', 'constant', 'zero']: rng = range(0, 2) for rep in rng: ax.plot((left + rep * step) * o2, [xp.min() - .5, xp.max() + .5], 'k-') diff --git a/pywt/_dwt.py b/pywt/_dwt.py index ea2d05047..bf2a3bbb2 100644 --- a/pywt/_dwt.py +++ b/pywt/_dwt.py @@ -12,7 +12,7 @@ __all__ = ["dwt", "idwt", "downcoef", "upcoef", "dwt_max_level", - "dwt_coeff_len"] + "dwt_coeff_len", "pad"] def dwt_max_level(data_len, filter_len): @@ -135,7 +135,6 @@ def dwt(data, wavelet, mode='symmetric', axis=-1): Axis over which to compute the DWT. If not given, the last axis is used. - Returns ------- (cA, cD) : tuple @@ -211,7 +210,6 @@ def idwt(cA, cD, wavelet, mode='symmetric', axis=-1): Axis over which to compute the inverse DWT. If not given, the last axis is used. - Returns ------- rec: array_like @@ -401,3 +399,119 @@ def upcoef(part, coeffs, wavelet, level=1, take=0): if part not in 'ad': raise ValueError("Argument 1 must be 'a' or 'd', not '%s'." % part) return np.asarray(_upcoef(part == 'a', coeffs, wavelet, level, take)) + + +def pad(x, pad_widths, mode): + """Extend a 1D signal using a given boundary mode. + + This function operates like :func:`numpy.pad` but supports all signal + extension modes that can be used by PyWavelets discrete wavelet transforms. + + Parameters + ---------- + x : ndarray + The array to pad + pad_widths : {sequence, array_like, int} + Number of values padded to the edges of each axis. + ``((before_1, after_1), … (before_N, after_N))`` unique pad widths for + each axis. ``((before, after),)`` yields same before and after pad for + each axis. ``(pad,)`` or int is a shortcut for + ``before = after = pad width`` for all axes. + mode : str, optional + Signal extension mode, see :ref:`Modes `. + + Returns + ------- + pad : ndarray + Padded array of rank equal to array with shape increased according to + ``pad_widths``. + + Notes + ----- + The performance of padding in dimensions > 1 may be substantially slower + for modes ``'smooth'`` and ``'antisymmetric'`` as these modes are not + supported efficiently by the underlying :func:`numpy.pad` function. + + Note that the behavior of the ``'constant'`` mode here follows the + PyWavelets convention which is different from NumPy (it is equivalent to + ``mode='edge'`` in :func:`numpy.pad`). + """ + x = np.asanyarray(x) + + # process pad_widths exactly as in numpy.pad + pad_widths = np.array(pad_widths) + pad_widths = np.round(pad_widths).astype(np.intp, copy=False) + if pad_widths.min() < 0: + raise ValueError("pad_widths must be > 0") + pad_widths = np.broadcast_to(pad_widths, (x.ndim, 2)).tolist() + + if mode in ['symmetric', 'reflect']: + xp = np.pad(x, pad_widths, mode=mode) + elif mode in ['periodic', 'periodization']: + if mode == 'periodization': + # Promote odd-sized dimensions to even length by duplicating the + # last value. + edge_pad_widths = [(0, x.shape[ax] % 2) + for ax in range(x.ndim)] + x = np.pad(x, edge_pad_widths, mode='edge') + xp = np.pad(x, pad_widths, mode='wrap') + elif mode == 'zero': + xp = np.pad(x, pad_widths, mode='constant', constant_values=0) + elif mode == 'constant': + xp = np.pad(x, pad_widths, mode='edge') + elif mode == 'smooth': + def pad_smooth(vector, pad_width, iaxis, kwargs): + # smooth extension to left + left = vector[pad_width[0]] + slope_left = (left - vector[pad_width[0] + 1]) + vector[:pad_width[0]] = \ + left + np.arange(pad_width[0], 0, -1) * slope_left + + # smooth extension to right + right = vector[-pad_width[1] - 1] + slope_right = (right - vector[-pad_width[1] - 2]) + vector[-pad_width[1]:] = \ + right + np.arange(1, pad_width[1] + 1) * slope_right + return vector + xp = np.pad(x, pad_widths, pad_smooth) + elif mode == 'antisymmetric': + def pad_antisymmetric(vector, pad_width, iaxis, kwargs): + # smooth extension to left + # implement by flipping portions symmetric padding + npad_l, npad_r = pad_width + vsize_nonpad = vector.size - npad_l - npad_r + # Note: must modify vector in-place + vector[:] = np.pad(vector[pad_width[0]:-pad_width[-1]], + pad_width, mode='symmetric') + vp = vector + r_edge = npad_l + vsize_nonpad - 1 + l_edge = npad_l + # width of each reflected segment + seg_width = vsize_nonpad + # flip reflected segments on the right of the original signal + n = 1 + while r_edge <= vp.size: + segment_slice = slice(r_edge + 1, + min(r_edge + 1 + seg_width, vp.size)) + if n % 2: + vp[segment_slice] *= -1 + r_edge += seg_width + n += 1 + + # flip reflected segments on the left of the original signal + n = 1 + while l_edge >= 0: + segment_slice = slice(max(0, l_edge - seg_width), l_edge) + if n % 2: + vp[segment_slice] *= -1 + l_edge -= seg_width + n += 1 + return vector + xp = np.pad(x, pad_widths, pad_antisymmetric) + elif mode == 'antireflect': + xp = np.pad(x, pad_widths, mode='reflect', reflect_type='odd') + else: + raise ValueError( + ("unsupported mode: {}. The supported modes are {}").format( + mode, Modes.modes)) + return xp diff --git a/pywt/_utils.py b/pywt/_utils.py index 48f814e21..291e85070 100644 --- a/pywt/_utils.py +++ b/pywt/_utils.py @@ -2,6 +2,7 @@ # # See COPYING for license details. import inspect +import numpy as np import sys from collections.abc import Iterable @@ -17,7 +18,7 @@ def _as_wavelet(wavelet): - """Convert wavelet name to a Wavelet object""" + """Convert wavelet name to a Wavelet object.""" if not isinstance(wavelet, (ContinuousWavelet, Wavelet)): wavelet = DiscreteContinuousWavelet(wavelet) if isinstance(wavelet, ContinuousWavelet): diff --git a/pywt/tests/test_dwt_idwt.py b/pywt/tests/test_dwt_idwt.py index 1779f8093..7e9e206bf 100644 --- a/pywt/tests/test_dwt_idwt.py +++ b/pywt/tests/test_dwt_idwt.py @@ -2,8 +2,8 @@ from __future__ import division, print_function, absolute_import import numpy as np -from numpy.testing import assert_allclose, assert_, assert_raises - +from numpy.testing import (assert_allclose, assert_, assert_raises, + assert_array_equal) import pywt # Check that float32, float64, complex64, complex128 are preserved. @@ -228,8 +228,72 @@ def test_error_on_continuous_wavelet(): def test_dwt_zero_size_axes(): # raise on empty input array assert_raises(ValueError, pywt.dwt, [], 'db2') - + # >1D case uses a different code path so check there as well x = np.ones((1, 4))[0:0, :] # 2D with a size zero axis assert_raises(ValueError, pywt.dwt, x, 'db2', axis=0) + +def test_pad_1d(): + x = [1, 2, 3] + assert_array_equal(pywt.pad(x, (4, 6), 'periodization'), + [1, 2, 3, 3, 1, 2, 3, 3, 1, 2, 3, 3, 1, 2]) + assert_array_equal(pywt.pad(x, (4, 6), 'periodic'), + [3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3]) + assert_array_equal(pywt.pad(x, (4, 6), 'constant'), + [1, 1, 1, 1, 1, 2, 3, 3, 3, 3, 3, 3, 3]) + assert_array_equal(pywt.pad(x, (4, 6), 'zero'), + [0, 0, 0, 0, 1, 2, 3, 0, 0, 0, 0, 0, 0]) + assert_array_equal(pywt.pad(x, (4, 6), 'smooth'), + [-3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9]) + assert_array_equal(pywt.pad(x, (4, 6), 'symmetric'), + [3, 3, 2, 1, 1, 2, 3, 3, 2, 1, 1, 2, 3]) + assert_array_equal(pywt.pad(x, (4, 6), 'antisymmetric'), + [3, -3, -2, -1, 1, 2, 3, -3, -2, -1, 1, 2, 3]) + assert_array_equal(pywt.pad(x, (4, 6), 'reflect'), + [1, 2, 3, 2, 1, 2, 3, 2, 1, 2, 3, 2, 1]) + assert_array_equal(pywt.pad(x, (4, 6), 'antireflect'), + [-3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9]) + + # equivalence of various pad_width formats + assert_array_equal(pywt.pad(x, 4, 'periodic'), + pywt.pad(x, (4, 4), 'periodic')) + + assert_array_equal(pywt.pad(x, (4, ), 'periodic'), + pywt.pad(x, (4, 4), 'periodic')) + + assert_array_equal(pywt.pad(x, [(4, 4)], 'periodic'), + pywt.pad(x, (4, 4), 'periodic')) + + +def test_pad_errors(): + # negative pad width + x = [1, 2, 3] + assert_raises(ValueError, pywt.pad, x, -2, 'periodic') + + # wrong length pad width + assert_raises(ValueError, pywt.pad, x, (1, 1, 1), 'periodic') + + # invalid mode name + assert_raises(ValueError, pywt.pad, x, 2, 'bad_mode') + + +def test_pad_nd(): + for ndim in [2, 3]: + x = np.arange(4**ndim).reshape((4, ) * ndim) + if ndim == 2: + pad_widths = [(2, 1), (2, 3)] + else: + pad_widths = [(2, 1), ] * ndim + for mode in pywt.Modes.modes: + xp = pywt.pad(x, pad_widths, mode) + + # expected result is the same as applying along axes separably + xp_expected = x.copy() + for ax in range(ndim): + xp_expected = np.apply_along_axis(pywt.pad, + ax, + xp_expected, + pad_widths=[pad_widths[ax]], + mode=mode) + assert_array_equal(xp, xp_expected)