From 941a11c8ff4d4f567873a245f046c30a722ea181 Mon Sep 17 00:00:00 2001 From: Alexandre Sauve Date: Mon, 15 Apr 2019 15:31:05 +0200 Subject: [PATCH 01/13] Added suppport for cwt with fft based convolution --- pywt/_cwt.py | 59 +++++++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 52 insertions(+), 7 deletions(-) diff --git a/pywt/_cwt.py b/pywt/_cwt.py index a4e6ca536..c592d3d3d 100644 --- a/pywt/_cwt.py +++ b/pywt/_cwt.py @@ -7,7 +7,7 @@ __all__ = ["cwt"] -def cwt(data, scales, wavelet, sampling_period=1.): +def cwt(data, scales, wavelet, sampling_period=1., method='conv'): """ cwt(data, scales, wavelet) @@ -29,7 +29,17 @@ def cwt(data, scales, wavelet, sampling_period=1.): The values computed for ``coefs`` are independent of the choice of ``sampling_period`` (i.e. ``scales`` is not scaled by the sampling period). - + method : convolution method name + Can be any of + - ``conv`` uses only the ``numpyp.conv`` function + - ``fft`` uses frequency domain convolution with ``numpyp.fft.fft`` + - ``auto`` for automatic selection of the fastest convolution method + depending on the complexity at each scale. + The ``conv`` method complexity is O(len(scale)*len(data)). + The ``fft`` method is O(N*log2(N)) with N=len(scale)+len(data)-1, + it is well suited for large size signals but slower than ``conv`` on + small ones. + Returns ------- coefs : array_like @@ -74,21 +84,56 @@ def cwt(data, scales, wavelet, sampling_period=1.): wavelet = DiscreteContinuousWavelet(wavelet) if np.isscalar(scales): scales = np.array([scales]) + dt_out = None # currently keep the 1.0.2 behaviour: TODO fix in/out dtype consistency if data.ndim == 1: if wavelet.complex_cwt: - out = np.zeros((np.size(scales), data.size), dtype=complex) - else: - out = np.zeros((np.size(scales), data.size)) + dt_out = complex + out = np.zeros((np.size(scales), data.size), dtype=dt_out) precision = 10 int_psi, x = integrate_wavelet(wavelet, precision=precision) + + if method in ('auto', 'fft'): + # - to be as large as the sum of data length and and maximum wavelet + # support to avoid circular convolution effects + # - additional padding to reach a power of 2 for CPU-optimal FFT + size_pad = lambda s: 2**np.int(np.ceil(np.log2(s[0] + s[1]))) + size_scale0 = size_pad( (len(data), + np.take(scales, 0) * ((x[-1] - x[0]) + 1)) ) + fft_data = None + elif not method == 'conv': + raise ValueError("method must be in: 'conv', 'fft' or 'auto'") + for i in np.arange(np.size(scales)): step = x[1] - x[0] j = np.floor( np.arange(scales[i] * (x[-1] - x[0]) + 1) / (scales[i] * step)) if np.max(j) >= np.size(int_psi): j = np.delete(j, np.where((j >= np.size(int_psi)))[0]) - coef = - np.sqrt(scales[i]) * np.diff( - np.convolve(data, int_psi[j.astype(np.int)][::-1])) + int_psi_scale = int_psi[j.astype(np.int)][::-1] + + if method == 'conv': + conv = np.convolve(data, int_psi_scale) + else: + size_scale = size_pad( (len(data), len(int_psi_scale)) ) + if size_scale != size_scale0: + # the fft of data changes when padding size changes thus + # it has to be recomputed + fft_data = None + size_scale0 = size_scale + nops_conv = len(data) * len(int_psi_scale) + nops_fft = (2+(fft_data is None)) * size_scale * np.log2(size_scale) + if (method == 'fft') or ((method == 'auto') and (nops_fft < nops_conv)): + if fft_data is None: + fft_data = np.fft.fft(data, size_scale) + fft_wav = np.fft.fft(int_psi_scale, size_scale) + conv = np.fft.ifft(fft_wav*fft_data) + conv = conv[0:len(data)+len(int_psi_scale)-1] + else: + conv = np.convolve(data, int_psi_scale) + + coef = - np.sqrt(scales[i]) * np.diff(conv) + if not np.iscomplexobj(out): + coef = np.real(coef) d = (coef.size - data.size) / 2. if d > 0: out[i, :] = coef[int(np.floor(d)):int(-np.ceil(d))] From 75904742da0c9cfdfe4dabaa0870f3b5a6a36254 Mon Sep 17 00:00:00 2001 From: Alexandre Sauve Date: Mon, 15 Apr 2019 15:31:44 +0200 Subject: [PATCH 02/13] added test for fft based cwt --- pywt/tests/test_cwt_wavelets.py | 50 ++++++++++++++++++++++++++++++++- 1 file changed, 49 insertions(+), 1 deletion(-) diff --git a/pywt/tests/test_cwt_wavelets.py b/pywt/tests/test_cwt_wavelets.py index 4372efc6a..c3ad46c4d 100644 --- a/pywt/tests/test_cwt_wavelets.py +++ b/pywt/tests/test_cwt_wavelets.py @@ -2,7 +2,7 @@ from __future__ import division, print_function, absolute_import from numpy.testing import (assert_allclose, assert_warns, assert_almost_equal, - assert_raises) + assert_raises, assert_equal) import numpy as np import pywt @@ -371,3 +371,51 @@ def test_cwt_small_scales(): # extremely short scale factors raise a ValueError assert_raises(ValueError, pywt.cwt, data, scales=0.01, wavelet='mexh') + + +def test_cwt_method_fft(): + data = np.zeros(32, dtype=np.float32) + data[15] = 1. + scales1 = 1 + wavelet = 'cmor1.5-1.0' + + # build a reference cwt with the legacy np.conv() method + cfs_conv, _ = pywt.cwt(data, scales1, wavelet, method='conv') + + # compare with the fft based convolution + cfs_fft, _ = pywt.cwt(data, scales1, wavelet, method='fft') + assert_allclose(cfs_conv, cfs_fft, rtol=0, atol=1e-13) + + +def test_cwt_method_auto(): + np.random.seed(1) + data = np.random.randn(50) + scales = [1, 5, 25, 125] + wavelet = 'cmor1.5-1.0' + + # build a reference cwt with the legacy np.conv() method + cfs_conv, _ = pywt.cwt(data, scales, wavelet, method='conv') + + # 'fft' method switches for scale 2 with len(data)==50 + cfs_fft, _ = pywt.cwt(data, scales, wavelet, method='auto') + assert_allclose(cfs_conv, cfs_fft, rtol=0, atol=1e-13) + + +def test_cwt_dtype(): + """Currently output dtype precision is fixed in version 1.0.2""" + wavelet = 'mexh' + scales = 1 + dtype_expected = { + np.float16: np.float64, + np.float32: np.float64, + np.float64: np.float64, + np.float128: np.float64, + np.complex64: np.float64, + np.complex128: np.float64, + np.complex256: np.float64, + } + for dtype_in, dtype_out in dtype_expected.items(): + data = np.zeros(2, dtype=dtype_in) + cfs, f = pywt.cwt(data, scales, wavelet) + assert_equal(dtype_out, cfs.dtype) + From c3594b9bf13e50909bd422db202138ead84be4d8 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Tue, 16 Apr 2019 15:56:09 -0400 Subject: [PATCH 03/13] fix docstring formatting --- pywt/_cwt.py | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/pywt/_cwt.py b/pywt/_cwt.py index c592d3d3d..57cccf0c2 100644 --- a/pywt/_cwt.py +++ b/pywt/_cwt.py @@ -29,17 +29,17 @@ def cwt(data, scales, wavelet, sampling_period=1., method='conv'): The values computed for ``coefs`` are independent of the choice of ``sampling_period`` (i.e. ``scales`` is not scaled by the sampling period). - method : convolution method name - Can be any of - - ``conv`` uses only the ``numpyp.conv`` function - - ``fft`` uses frequency domain convolution with ``numpyp.fft.fft`` - - ``auto`` for automatic selection of the fastest convolution method - depending on the complexity at each scale. - The ``conv`` method complexity is O(len(scale)*len(data)). - The ``fft`` method is O(N*log2(N)) with N=len(scale)+len(data)-1, - it is well suited for large size signals but slower than ``conv`` on - small ones. - + method : {'conv', 'fft', 'auto'}, optional + The method used to compute the CWT. Can be any of: + - ``conv`` uses ``numpy.convolve``. + - ``fft`` uses frequency domain convolution via ``numpy.fft.fft``. + - ``auto`` uses automatic selection based on an estimate of the + computational complexity at each scale. + The ``conv`` method complexity is ``O(len(scale) * len(data))``. + The ``fft`` method is ``O(N * log2(N))`` with + ``N = len(scale) + len(data) - 1``. It is well suited for large size + signals but slower than ``conv`` on small ones. + Returns ------- coefs : array_like @@ -91,13 +91,13 @@ def cwt(data, scales, wavelet, sampling_period=1., method='conv'): out = np.zeros((np.size(scales), data.size), dtype=dt_out) precision = 10 int_psi, x = integrate_wavelet(wavelet, precision=precision) - + if method in ('auto', 'fft'): # - to be as large as the sum of data length and and maximum wavelet # support to avoid circular convolution effects # - additional padding to reach a power of 2 for CPU-optimal FFT size_pad = lambda s: 2**np.int(np.ceil(np.log2(s[0] + s[1]))) - size_scale0 = size_pad( (len(data), + size_scale0 = size_pad( (len(data), np.take(scales, 0) * ((x[-1] - x[0]) + 1)) ) fft_data = None elif not method == 'conv': @@ -110,7 +110,7 @@ def cwt(data, scales, wavelet, sampling_period=1., method='conv'): if np.max(j) >= np.size(int_psi): j = np.delete(j, np.where((j >= np.size(int_psi)))[0]) int_psi_scale = int_psi[j.astype(np.int)][::-1] - + if method == 'conv': conv = np.convolve(data, int_psi_scale) else: @@ -130,7 +130,7 @@ def cwt(data, scales, wavelet, sampling_period=1., method='conv'): conv = conv[0:len(data)+len(int_psi_scale)-1] else: conv = np.convolve(data, int_psi_scale) - + coef = - np.sqrt(scales[i]) * np.diff(conv) if not np.iscomplexobj(out): coef = np.real(coef) From c837c6c6956729f3a1ed3bac10bdd227b2bd5a92 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Tue, 16 Apr 2019 16:08:09 -0400 Subject: [PATCH 04/13] PEP8 fixes --- pywt/_cwt.py | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/pywt/_cwt.py b/pywt/_cwt.py index 57cccf0c2..a8f699f0b 100644 --- a/pywt/_cwt.py +++ b/pywt/_cwt.py @@ -84,7 +84,7 @@ def cwt(data, scales, wavelet, sampling_period=1., method='conv'): wavelet = DiscreteContinuousWavelet(wavelet) if np.isscalar(scales): scales = np.array([scales]) - dt_out = None # currently keep the 1.0.2 behaviour: TODO fix in/out dtype consistency + dt_out = None # TODO: fix in/out dtype consistency in a subsequent PR if data.ndim == 1: if wavelet.complex_cwt: dt_out = complex @@ -93,12 +93,12 @@ def cwt(data, scales, wavelet, sampling_period=1., method='conv'): int_psi, x = integrate_wavelet(wavelet, precision=precision) if method in ('auto', 'fft'): - # - to be as large as the sum of data length and and maximum wavelet - # support to avoid circular convolution effects + # - to be as large as the sum of data length and and maximum + # wavelet support to avoid circular convolution effects # - additional padding to reach a power of 2 for CPU-optimal FFT size_pad = lambda s: 2**np.int(np.ceil(np.log2(s[0] + s[1]))) - size_scale0 = size_pad( (len(data), - np.take(scales, 0) * ((x[-1] - x[0]) + 1)) ) + size_scale0 = size_pad((len(data), + np.take(scales, 0) * ((x[-1] - x[0]) + 1))) fft_data = None elif not method == 'conv': raise ValueError("method must be in: 'conv', 'fft' or 'auto'") @@ -114,20 +114,22 @@ def cwt(data, scales, wavelet, sampling_period=1., method='conv'): if method == 'conv': conv = np.convolve(data, int_psi_scale) else: - size_scale = size_pad( (len(data), len(int_psi_scale)) ) + size_scale = size_pad((len(data), len(int_psi_scale))) if size_scale != size_scale0: # the fft of data changes when padding size changes thus # it has to be recomputed fft_data = None size_scale0 = size_scale nops_conv = len(data) * len(int_psi_scale) - nops_fft = (2+(fft_data is None)) * size_scale * np.log2(size_scale) - if (method == 'fft') or ((method == 'auto') and (nops_fft < nops_conv)): + nops_fft = (2 + (fft_data is None)) + nops_fft *= size_scale * np.log2(size_scale) + if (method == 'fft') or ( + (method == 'auto') and (nops_fft < nops_conv)): if fft_data is None: fft_data = np.fft.fft(data, size_scale) fft_wav = np.fft.fft(int_psi_scale, size_scale) - conv = np.fft.ifft(fft_wav*fft_data) - conv = conv[0:len(data)+len(int_psi_scale)-1] + conv = np.fft.ifft(fft_wav * fft_data) + conv = conv[:data.size + int_psi_scale.size - 1] else: conv = np.convolve(data, int_psi_scale) From a767fc6a749925e1126d15d93062a58cea99e356 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Tue, 16 Apr 2019 16:56:11 -0400 Subject: [PATCH 05/13] minor refactoring of cwt --- pywt/_cwt.py | 41 ++++++++++++++++++++++------------------- 1 file changed, 22 insertions(+), 19 deletions(-) diff --git a/pywt/_cwt.py b/pywt/_cwt.py index a8f699f0b..59a6366e3 100644 --- a/pywt/_cwt.py +++ b/pywt/_cwt.py @@ -1,9 +1,12 @@ +from math import sqrt, log2, floor, ceil + import numpy as np from ._extensions._pywt import (DiscreteContinuousWavelet, ContinuousWavelet, Wavelet, _check_dtype) from ._functions import integrate_wavelet, scale2frequency + __all__ = ["cwt"] @@ -88,7 +91,7 @@ def cwt(data, scales, wavelet, sampling_period=1., method='conv'): if data.ndim == 1: if wavelet.complex_cwt: dt_out = complex - out = np.zeros((np.size(scales), data.size), dtype=dt_out) + out = np.empty((np.size(scales), data.size), dtype=dt_out) precision = 10 int_psi, x = integrate_wavelet(wavelet, precision=precision) @@ -96,33 +99,34 @@ def cwt(data, scales, wavelet, sampling_period=1., method='conv'): # - to be as large as the sum of data length and and maximum # wavelet support to avoid circular convolution effects # - additional padding to reach a power of 2 for CPU-optimal FFT - size_pad = lambda s: 2**np.int(np.ceil(np.log2(s[0] + s[1]))) - size_scale0 = size_pad((len(data), - np.take(scales, 0) * ((x[-1] - x[0]) + 1))) + size_pad = lambda s: 2**ceil(log2(s[0] + s[1])) + size_scale0 = size_pad((data.size, + scales[0] * ((x[-1] - x[0]) + 1))) fft_data = None elif not method == 'conv': raise ValueError("method must be in: 'conv', 'fft' or 'auto'") - for i in np.arange(np.size(scales)): + for i, scale in enumerate(scales): step = x[1] - x[0] - j = np.floor( - np.arange(scales[i] * (x[-1] - x[0]) + 1) / (scales[i] * step)) - if np.max(j) >= np.size(int_psi): - j = np.delete(j, np.where((j >= np.size(int_psi)))[0]) - int_psi_scale = int_psi[j.astype(np.int)][::-1] + j = np.arange(scale * (x[-1] - x[0]) + 1) / (scale * step) + j = j.astype(int) # floor + if j[-1] >= int_psi.size: + j = np.extract(j < int_psi.size, j) + int_psi_scale = int_psi[j][::-1] if method == 'conv': conv = np.convolve(data, int_psi_scale) else: - size_scale = size_pad((len(data), len(int_psi_scale))) + size_scale = size_pad((data.size, int_psi_scale.size)) if size_scale != size_scale0: # the fft of data changes when padding size changes thus # it has to be recomputed fft_data = None size_scale0 = size_scale - nops_conv = len(data) * len(int_psi_scale) - nops_fft = (2 + (fft_data is None)) - nops_fft *= size_scale * np.log2(size_scale) + if method == 'auto': + nops_conv = len(data) * len(int_psi_scale) + nops_fft = ((2 + (fft_data is None)) * + size_scale * log2(size_scale)) if (method == 'fft') or ( (method == 'auto') and (nops_fft < nops_conv)): if fft_data is None: @@ -133,22 +137,21 @@ def cwt(data, scales, wavelet, sampling_period=1., method='conv'): else: conv = np.convolve(data, int_psi_scale) - coef = - np.sqrt(scales[i]) * np.diff(conv) + coef = - sqrt(scale) * np.diff(conv) if not np.iscomplexobj(out): coef = np.real(coef) d = (coef.size - data.size) / 2. if d > 0: - out[i, :] = coef[int(np.floor(d)):int(-np.ceil(d))] + out[i, :] = coef[floor(d):-ceil(d)] elif d == 0.: out[i, :] = coef else: raise ValueError( - "Selected scale of {} too small.".format(scales[i])) + "Selected scale of {} too small.".format(scale)) frequencies = scale2frequency(wavelet, scales, precision) if np.isscalar(frequencies): frequencies = np.array([frequencies]) - for i in np.arange(len(frequencies)): - frequencies[i] /= sampling_period + frequencies /= sampling_period return out, frequencies else: raise ValueError("Only dim == 1 supported") From 0e0d41394939ddfcb7c05d9ef947bfeea0ff4270 Mon Sep 17 00:00:00 2001 From: Alexandre Sauve Date: Wed, 17 Apr 2019 10:43:39 +0200 Subject: [PATCH 06/13] Updated code (removed method 'auto'/randomstate/removed dtype checks) accordingly to #pullrequestreview-227464829 --- pywt/_cwt.py | 26 ++++++++----------- pywt/tests/test_cwt_wavelets.py | 44 +++++---------------------------- 2 files changed, 16 insertions(+), 54 deletions(-) diff --git a/pywt/_cwt.py b/pywt/_cwt.py index c592d3d3d..f7213f490 100644 --- a/pywt/_cwt.py +++ b/pywt/_cwt.py @@ -33,12 +33,10 @@ def cwt(data, scales, wavelet, sampling_period=1., method='conv'): Can be any of - ``conv`` uses only the ``numpyp.conv`` function - ``fft`` uses frequency domain convolution with ``numpyp.fft.fft`` - - ``auto`` for automatic selection of the fastest convolution method - depending on the complexity at each scale. The ``conv`` method complexity is O(len(scale)*len(data)). The ``fft`` method is O(N*log2(N)) with N=len(scale)+len(data)-1, - it is well suited for large size signals but slower than ``conv`` on - small ones. + it is well suited for large size signals but slightly slower than + ``conv`` on small ones. Returns ------- @@ -92,7 +90,8 @@ def cwt(data, scales, wavelet, sampling_period=1., method='conv'): precision = 10 int_psi, x = integrate_wavelet(wavelet, precision=precision) - if method in ('auto', 'fft'): + if method == 'fft': + # for FFT the buffer needs: # - to be as large as the sum of data length and and maximum wavelet # support to avoid circular convolution effects # - additional padding to reach a power of 2 for CPU-optimal FFT @@ -101,7 +100,7 @@ def cwt(data, scales, wavelet, sampling_period=1., method='conv'): np.take(scales, 0) * ((x[-1] - x[0]) + 1)) ) fft_data = None elif not method == 'conv': - raise ValueError("method must be in: 'conv', 'fft' or 'auto'") + raise ValueError("method must be 'conv' or 'fft'") for i in np.arange(np.size(scales)): step = x[1] - x[0] @@ -120,16 +119,11 @@ def cwt(data, scales, wavelet, sampling_period=1., method='conv'): # it has to be recomputed fft_data = None size_scale0 = size_scale - nops_conv = len(data) * len(int_psi_scale) - nops_fft = (2+(fft_data is None)) * size_scale * np.log2(size_scale) - if (method == 'fft') or ((method == 'auto') and (nops_fft < nops_conv)): - if fft_data is None: - fft_data = np.fft.fft(data, size_scale) - fft_wav = np.fft.fft(int_psi_scale, size_scale) - conv = np.fft.ifft(fft_wav*fft_data) - conv = conv[0:len(data)+len(int_psi_scale)-1] - else: - conv = np.convolve(data, int_psi_scale) + if fft_data is None: + fft_data = np.fft.fft(data, size_scale) + fft_wav = np.fft.fft(int_psi_scale, size_scale) + conv = np.fft.ifft(fft_wav*fft_data) + conv = conv[0:len(data)+len(int_psi_scale)-1] coef = - np.sqrt(scales[i]) * np.diff(conv) if not np.iscomplexobj(out): diff --git a/pywt/tests/test_cwt_wavelets.py b/pywt/tests/test_cwt_wavelets.py index c3ad46c4d..a40b8d9cb 100644 --- a/pywt/tests/test_cwt_wavelets.py +++ b/pywt/tests/test_cwt_wavelets.py @@ -2,7 +2,7 @@ from __future__ import division, print_function, absolute_import from numpy.testing import (assert_allclose, assert_warns, assert_almost_equal, - assert_raises, assert_equal) + assert_raises) import numpy as np import pywt @@ -374,48 +374,16 @@ def test_cwt_small_scales(): def test_cwt_method_fft(): - data = np.zeros(32, dtype=np.float32) + rstate = np.random.RandomState(1) + data = rstate.randn(50) data[15] = 1. - scales1 = 1 - wavelet = 'cmor1.5-1.0' - - # build a reference cwt with the legacy np.conv() method - cfs_conv, _ = pywt.cwt(data, scales1, wavelet, method='conv') - - # compare with the fft based convolution - cfs_fft, _ = pywt.cwt(data, scales1, wavelet, method='fft') - assert_allclose(cfs_conv, cfs_fft, rtol=0, atol=1e-13) - - -def test_cwt_method_auto(): - np.random.seed(1) - data = np.random.randn(50) - scales = [1, 5, 25, 125] + scales = [1, 32, 64] wavelet = 'cmor1.5-1.0' # build a reference cwt with the legacy np.conv() method cfs_conv, _ = pywt.cwt(data, scales, wavelet, method='conv') - # 'fft' method switches for scale 2 with len(data)==50 - cfs_fft, _ = pywt.cwt(data, scales, wavelet, method='auto') + # compare with the fft based convolution + cfs_fft, _ = pywt.cwt(data, scales, wavelet, method='fft') assert_allclose(cfs_conv, cfs_fft, rtol=0, atol=1e-13) - - -def test_cwt_dtype(): - """Currently output dtype precision is fixed in version 1.0.2""" - wavelet = 'mexh' - scales = 1 - dtype_expected = { - np.float16: np.float64, - np.float32: np.float64, - np.float64: np.float64, - np.float128: np.float64, - np.complex64: np.float64, - np.complex128: np.float64, - np.complex256: np.float64, - } - for dtype_in, dtype_out in dtype_expected.items(): - data = np.zeros(2, dtype=dtype_in) - cfs, f = pywt.cwt(data, scales, wavelet) - assert_equal(dtype_out, cfs.dtype) From 5cf11882a763a07842a564092c0ae93578fd4169 Mon Sep 17 00:00:00 2001 From: Alexandre Sauve Date: Wed, 17 Apr 2019 11:59:55 +0200 Subject: [PATCH 07/13] Added support for next_fast_len() to improve complexity --- pywt/_cwt.py | 77 +++++++++++++++++---------------- pywt/tests/test_cwt_wavelets.py | 7 ++- 2 files changed, 42 insertions(+), 42 deletions(-) diff --git a/pywt/_cwt.py b/pywt/_cwt.py index f7213f490..1b4ab16fc 100644 --- a/pywt/_cwt.py +++ b/pywt/_cwt.py @@ -1,9 +1,13 @@ +from math import floor, ceil +from scipy.fftpack import next_fast_len + import numpy as np from ._extensions._pywt import (DiscreteContinuousWavelet, ContinuousWavelet, Wavelet, _check_dtype) from ._functions import integrate_wavelet, scale2frequency + __all__ = ["cwt"] @@ -29,15 +33,17 @@ def cwt(data, scales, wavelet, sampling_period=1., method='conv'): The values computed for ``coefs`` are independent of the choice of ``sampling_period`` (i.e. ``scales`` is not scaled by the sampling period). - method : convolution method name - Can be any of - - ``conv`` uses only the ``numpyp.conv`` function - - ``fft`` uses frequency domain convolution with ``numpyp.fft.fft`` - The ``conv`` method complexity is O(len(scale)*len(data)). - The ``fft`` method is O(N*log2(N)) with N=len(scale)+len(data)-1, - it is well suited for large size signals but slightly slower than - ``conv`` on small ones. - + method : {'conv', 'fft'}, optional + The method used to compute the CWT. Can be any of: + - ``conv`` uses ``numpy.convolve``. + - ``fft`` uses frequency domain convolution via ``numpy.fft.fft``. + - ``auto`` uses automatic selection based on an estimate of the + computational complexity at each scale. + The ``conv`` method complexity is ``O(len(scale) * len(data))``. + The ``fft`` method is ``O(N * log2(N))`` with + ``N = len(scale) + len(data) - 1``. It is well suited for large size + signals but slightly slower than ``conv`` on small ones. + Returns ------- coefs : array_like @@ -82,65 +88,60 @@ def cwt(data, scales, wavelet, sampling_period=1., method='conv'): wavelet = DiscreteContinuousWavelet(wavelet) if np.isscalar(scales): scales = np.array([scales]) - dt_out = None # currently keep the 1.0.2 behaviour: TODO fix in/out dtype consistency + dt_out = None # TODO: fix in/out dtype consistency in a subsequent PR if data.ndim == 1: if wavelet.complex_cwt: dt_out = complex - out = np.zeros((np.size(scales), data.size), dtype=dt_out) + out = np.empty((np.size(scales), data.size), dtype=dt_out) precision = 10 int_psi, x = integrate_wavelet(wavelet, precision=precision) - + if method == 'fft': - # for FFT the buffer needs: - # - to be as large as the sum of data length and and maximum wavelet - # support to avoid circular convolution effects - # - additional padding to reach a power of 2 for CPU-optimal FFT - size_pad = lambda s: 2**np.int(np.ceil(np.log2(s[0] + s[1]))) - size_scale0 = size_pad( (len(data), - np.take(scales, 0) * ((x[-1] - x[0]) + 1)) ) + size_scale0 = -1 fft_data = None elif not method == 'conv': raise ValueError("method must be 'conv' or 'fft'") - for i in np.arange(np.size(scales)): + for i, scale in enumerate(scales): step = x[1] - x[0] - j = np.floor( - np.arange(scales[i] * (x[-1] - x[0]) + 1) / (scales[i] * step)) - if np.max(j) >= np.size(int_psi): - j = np.delete(j, np.where((j >= np.size(int_psi)))[0]) - int_psi_scale = int_psi[j.astype(np.int)][::-1] - + j = np.arange(scale * (x[-1] - x[0]) + 1) / (scale * step) + j = j.astype(int) # floor + if j[-1] >= int_psi.size: + j = np.extract(j < int_psi.size, j) + int_psi_scale = int_psi[j][::-1] + if method == 'conv': conv = np.convolve(data, int_psi_scale) else: - size_scale = size_pad( (len(data), len(int_psi_scale)) ) + # the padding is selected for + # - optimal FFT complexity + # - to be larger than the two signals length to avoid circular + # convolution + size_scale = next_fast_len(data.size + int_psi_scale.size - 1) if size_scale != size_scale0: # the fft of data changes when padding size changes thus # it has to be recomputed - fft_data = None - size_scale0 = size_scale - if fft_data is None: fft_data = np.fft.fft(data, size_scale) + size_scale0 = size_scale fft_wav = np.fft.fft(int_psi_scale, size_scale) - conv = np.fft.ifft(fft_wav*fft_data) - conv = conv[0:len(data)+len(int_psi_scale)-1] - - coef = - np.sqrt(scales[i]) * np.diff(conv) + conv = np.fft.ifft(fft_wav * fft_data) + conv = conv[:data.size + int_psi_scale.size - 1] + + coef = - np.sqrt(scale) * np.diff(conv) if not np.iscomplexobj(out): coef = np.real(coef) d = (coef.size - data.size) / 2. if d > 0: - out[i, :] = coef[int(np.floor(d)):int(-np.ceil(d))] + out[i, :] = coef[floor(d):-ceil(d)] elif d == 0.: out[i, :] = coef else: raise ValueError( - "Selected scale of {} too small.".format(scales[i])) + "Selected scale of {} too small.".format(scale)) frequencies = scale2frequency(wavelet, scales, precision) if np.isscalar(frequencies): frequencies = np.array([frequencies]) - for i in np.arange(len(frequencies)): - frequencies[i] /= sampling_period + frequencies /= sampling_period return out, frequencies else: raise ValueError("Only dim == 1 supported") diff --git a/pywt/tests/test_cwt_wavelets.py b/pywt/tests/test_cwt_wavelets.py index a40b8d9cb..e2cafcb2a 100644 --- a/pywt/tests/test_cwt_wavelets.py +++ b/pywt/tests/test_cwt_wavelets.py @@ -371,19 +371,18 @@ def test_cwt_small_scales(): # extremely short scale factors raise a ValueError assert_raises(ValueError, pywt.cwt, data, scales=0.01, wavelet='mexh') - + def test_cwt_method_fft(): rstate = np.random.RandomState(1) data = rstate.randn(50) data[15] = 1. - scales = [1, 32, 64] + scales = np.arange(1, 64) wavelet = 'cmor1.5-1.0' - + # build a reference cwt with the legacy np.conv() method cfs_conv, _ = pywt.cwt(data, scales, wavelet, method='conv') # compare with the fft based convolution cfs_fft, _ = pywt.cwt(data, scales, wavelet, method='fft') assert_allclose(cfs_conv, cfs_fft, rtol=0, atol=1e-13) - From e392672a02d9468c7cca4f0be92e64fbfbae15b9 Mon Sep 17 00:00:00 2001 From: Alexandre Sauve Date: Thu, 18 Apr 2019 20:29:18 +0200 Subject: [PATCH 08/13] Fix the scipy import by providing a fallback --- pywt/_cwt.py | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/pywt/_cwt.py b/pywt/_cwt.py index 1b4ab16fc..32bd64655 100644 --- a/pywt/_cwt.py +++ b/pywt/_cwt.py @@ -1,7 +1,4 @@ from math import floor, ceil -from scipy.fftpack import next_fast_len - -import numpy as np from ._extensions._pywt import (DiscreteContinuousWavelet, ContinuousWavelet, Wavelet, _check_dtype) @@ -11,6 +8,20 @@ __all__ = ["cwt"] +import numpy as np + +try: + from scipy.fftpack import next_fast_len +except ImportError: + # Do provide a fallback so scipy is an optional requirement + def next_fast_len(n): + """Given a number of samples `n`, returns the next power of two + following this numbe to take advantage of FFT speedup. + This fallback is less efficient as `scipy.fftpack.next_fast_len` + """ + return 2**ceil(np.log2(n)) + + def cwt(data, scales, wavelet, sampling_period=1., method='conv'): """ cwt(data, scales, wavelet) From 540d25060ba783deb13d0fdeff4d863e16f7c1e7 Mon Sep 17 00:00:00 2001 From: Alexandre Sauve Date: Thu, 18 Apr 2019 21:30:46 +0200 Subject: [PATCH 09/13] Minor: typos --- pywt/_cwt.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pywt/_cwt.py b/pywt/_cwt.py index 32bd64655..437bfbcf9 100644 --- a/pywt/_cwt.py +++ b/pywt/_cwt.py @@ -16,8 +16,8 @@ # Do provide a fallback so scipy is an optional requirement def next_fast_len(n): """Given a number of samples `n`, returns the next power of two - following this numbe to take advantage of FFT speedup. - This fallback is less efficient as `scipy.fftpack.next_fast_len` + following this number to take advantage of FFT speedup. + This fallback is less efficient than `scipy.fftpack.next_fast_len` """ return 2**ceil(np.log2(n)) From 7b24572787eca58d1fa0e06e36faccab19402f72 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Thu, 25 Jul 2019 11:54:52 -0400 Subject: [PATCH 10/13] ENH: prefer scipy.fft when available --- pywt/_cwt.py | 41 ++++++++++++++++++++++++++--------------- 1 file changed, 26 insertions(+), 15 deletions(-) diff --git a/pywt/_cwt.py b/pywt/_cwt.py index 437bfbcf9..de6c4b24d 100644 --- a/pywt/_cwt.py +++ b/pywt/_cwt.py @@ -11,15 +11,27 @@ import numpy as np try: - from scipy.fftpack import next_fast_len + # Prefer scipy.fft (new in SciPy 1.4) + import scipy.fft + fftmodule = scipy.fft + next_fast_len = fftmodule.next_fast_len except ImportError: - # Do provide a fallback so scipy is an optional requirement - def next_fast_len(n): - """Given a number of samples `n`, returns the next power of two - following this number to take advantage of FFT speedup. - This fallback is less efficient than `scipy.fftpack.next_fast_len` - """ - return 2**ceil(np.log2(n)) + try: + import scipy.fftpack + fftmodule = scipy.fftpack + next_fast_len = fftmodule.next_fast_len + except ImportError: + fftmodule = np.fft + + # provide a fallback so scipy is an optional requirement + def next_fast_len(n): + """Round up size to the nearest power of two. + + Given a number of samples `n`, returns the next power of two + following this number to take advantage of FFT speedup. + This fallback is less efficient than `scipy.fftpack.next_fast_len` + """ + return 2**ceil(np.log2(n)) def cwt(data, scales, wavelet, sampling_period=1., method='conv'): @@ -47,7 +59,7 @@ def cwt(data, scales, wavelet, sampling_period=1., method='conv'): method : {'conv', 'fft'}, optional The method used to compute the CWT. Can be any of: - ``conv`` uses ``numpy.convolve``. - - ``fft`` uses frequency domain convolution via ``numpy.fft.fft``. + - ``fft`` uses frequency domain convolution. - ``auto`` uses automatic selection based on an estimate of the computational complexity at each scale. The ``conv`` method complexity is ``O(len(scale) * len(data))``. @@ -124,18 +136,17 @@ def cwt(data, scales, wavelet, sampling_period=1., method='conv'): if method == 'conv': conv = np.convolve(data, int_psi_scale) else: - # the padding is selected for + # The padding is selected for: # - optimal FFT complexity # - to be larger than the two signals length to avoid circular # convolution size_scale = next_fast_len(data.size + int_psi_scale.size - 1) if size_scale != size_scale0: - # the fft of data changes when padding size changes thus - # it has to be recomputed - fft_data = np.fft.fft(data, size_scale) + # Must recompute fft_data when the padding size changes. + fft_data = fftmodule.fft(data, size_scale) size_scale0 = size_scale - fft_wav = np.fft.fft(int_psi_scale, size_scale) - conv = np.fft.ifft(fft_wav * fft_data) + fft_wav = fftmodule.fft(int_psi_scale, size_scale) + conv = fftmodule.ifft(fft_wav * fft_data) conv = conv[:data.size + int_psi_scale.size - 1] coef = - np.sqrt(scale) * np.diff(conv) From 823a96d8e8e1edbf2f218052728587a055069018 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Thu, 25 Jul 2019 12:21:37 -0400 Subject: [PATCH 11/13] include FFT-based cases in the CWT benchmarks --- benchmarks/benchmarks/cwt_benchmarks.py | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/benchmarks/benchmarks/cwt_benchmarks.py b/benchmarks/benchmarks/cwt_benchmarks.py index cf9cacd35..0c063b3bf 100644 --- a/benchmarks/benchmarks/cwt_benchmarks.py +++ b/benchmarks/benchmarks/cwt_benchmarks.py @@ -6,12 +6,13 @@ class CwtTimeSuiteBase(object): """ Set-up for CWT timing. """ - params = ([32, 128, 512], + params = ([32, 128, 512, 2048], ['cmor', 'cgau4', 'fbsp', 'gaus4', 'mexh', 'morl', 'shan'], - [16, 64, 256]) - param_names = ('n', 'wavelet', 'max_scale') + [16, 64, 256], + ['conv', 'fft']) + param_names = ('n', 'wavelet', 'max_scale', 'method') - def setup(self, n, wavelet, max_scale): + def setup(self, n, wavelet, max_scale, method): try: from pywt import cwt except ImportError: @@ -21,5 +22,12 @@ def setup(self, n, wavelet, max_scale): class CwtTimeSuite(CwtTimeSuiteBase): - def time_cwt(self, n, wavelet, max_scale): - pywt.cwt(self.data, self.scales, wavelet) + def time_cwt(self, n, wavelet, max_scale, method): + try: + pywt.cwt(self.data, self.scales, wavelet, method=method) + except TypeError: + # older PyWavelets does not support use of the method argument + if method == 'fft': + raise NotImplementedError( + "fft-based convolution not available.") + pywt.cwt(self.data, self.scales, wavelet) From 688e2e66dd2b060916d79b92cdbb6b6be65004e8 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Thu, 25 Jul 2019 14:26:58 -0400 Subject: [PATCH 12/13] Add a case with SciPy included on Travis-CI. change 3.7-dev to 3.7 in .travis.yml --- .travis.yml | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index 0965345f5..11df85f61 100644 --- a/.travis.yml +++ b/.travis.yml @@ -19,12 +19,15 @@ matrix: - CYTHONSPEC=cython - USE_WHEEL=1 - os: linux - python: 3.7-dev + python: 3.7 + dist: xenial # travis-ci/travis-ci/issues/9815 + sudo: true env: - NUMPYSPEC=numpy - MATPLOTLIBSPEC=matplotlib - CYTHONSPEC=cython - USE_SDIST=1 + - USE_SCIPY=1 - os: linux python: 3.5 env: @@ -59,6 +62,7 @@ before_install: - pip install pytest pytest-cov coverage codecov futures - set -o pipefail - if [ "${USE_WHEEL}" == "1" ]; then pip install wheel; fi + - if [ "${USE_SCIPY}" == "1" ]; then pip install scipy; fi - | if [ "${REFGUIDE_CHECK}" == "1" ]; then pip install sphinx numpydoc From 923a6bbc95e25b4f44328c9e2e36086b6390a90c Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Thu, 25 Jul 2019 14:34:52 -0400 Subject: [PATCH 13/13] DOC: mention optional use of scipy in the README and install docs --- README.rst | 9 ++++++--- doc/source/common_refs.rst | 3 ++- doc/source/install.rst | 5 +++-- 3 files changed, 11 insertions(+), 6 deletions(-) diff --git a/README.rst b/README.rst index 3cde6e1b6..d6d1cb2c0 100644 --- a/README.rst +++ b/README.rst @@ -65,9 +65,11 @@ For more usage examples see the `demo`_ directory in the source package. Installation ------------ -PyWavelets supports `Python`_ >=3.5, and is only dependent on `Numpy`_ +PyWavelets supports `Python`_ >=3.5, and is only dependent on `NumPy`_ (supported versions are currently ``>= 1.13.3``). To pass all of the tests, -`Matplotlib`_ is also required. +`Matplotlib`_ is also required. `SciPy`_ is also an optional dependency. When +present, FFT-based continuous wavelet transforms will use FFTs from SciPy +rather than NumPy. There are binary wheels for Intel Linux, Windows and macOS / OSX on PyPi. If you are on one of these platforms, you should get a binary (precompiled) @@ -138,7 +140,8 @@ If you wish to cite PyWavelets in a publication, you may use the following DOI. .. _Anaconda: https://www.continuum.io .. _GitHub: https://github.com/PyWavelets/pywt .. _GitHub Issues: https://github.com/PyWavelets/pywt/issues -.. _Numpy: http://www.numpy.org +.. _NumPy: https://www.numpy.org +.. _SciPy: https://www.scipy.org .. _original developer: http://en.ig.ma .. _Python: http://python.org/ .. _Python Package Index: http://pypi.python.org/pypi/PyWavelets/ diff --git a/doc/source/common_refs.rst b/doc/source/common_refs.rst index a8331b221..8fe04b957 100644 --- a/doc/source/common_refs.rst +++ b/doc/source/common_refs.rst @@ -5,7 +5,8 @@ .. _GitHub: https://github.com/PyWavelets/pywt .. _GitHub repository: https://github.com/PyWavelets/pywt .. _GitHub Issues: https://github.com/PyWavelets/pywt/issues -.. _Numpy: http://www.numpy.org +.. _NumPy: https://www.numpy.org +.. _SciPy: https://www.scipy.org .. _original developer: http://en.ig.ma .. _Python: http://python.org/ .. _Python Package Index: http://pypi.python.org/pypi/PyWavelets/ diff --git a/doc/source/install.rst b/doc/source/install.rst index 7a82d00aa..f880a8ea8 100644 --- a/doc/source/install.rst +++ b/doc/source/install.rst @@ -39,11 +39,12 @@ PyWavelets source code directory (containing ``setup.py``) and type:: The requirements needed to build from source are: - Python_ 2.7 or >=3.4 - - Numpy_ >= 1.13.3 + - NumPy_ >= 1.13.3 - Cython_ >= 0.23.5 (if installing from git, not from a PyPI source release) To run all the tests for PyWavelets, you will also need to install the -Matplotlib_ package. +Matplotlib_ package. If SciPy_ is available, FFT-based continuous wavelet +transforms will use the FFT implementation from SciPy instead of NumPy. .. seealso:: :ref:`Development guide ` section contains more information on building and installing from source code.