From f9fb509f1f601e11bb7337809f0483b8c8ccd26f Mon Sep 17 00:00:00 2001 From: Kevin Lee Date: Wed, 4 Mar 2026 05:07:52 -0800 Subject: [PATCH 1/4] DOC: refine CWT normalization explanation --- doc/source/ref/cwt.rst | 29 +++++++++++++++++++++++++++++ pywt/_cwt.py | 10 ++++++++++ 2 files changed, 39 insertions(+) diff --git a/doc/source/ref/cwt.rst b/doc/source/ref/cwt.rst index 1910b6e1..c9d21d23 100644 --- a/doc/source/ref/cwt.rst +++ b/doc/source/ref/cwt.rst @@ -231,6 +231,35 @@ scales. The right column are the corresponding Fourier power spectra of each filter.. For scales 1 and 2 it can be seen that aliasing due to violation of the Nyquist limit occurs. + +Normalization with ``sampling_period`` +-------------------------------------- + +The ``sampling_period`` argument changes only the reported frequencies +(``scale2frequency(...)/sampling_period``). It does not rescale the CWT +coefficients themselves. + +This can be confusing when comparing against a hand-written approximation to +the continuous-time definition, + +.. math:: + + W_x(a, b) = \frac{1}{\sqrt{a}}\int x(t)\,\psi^*\left(\frac{t-b}{a}\right)\,dt. + +For sampled data ``t_n = n\,dt``, a direct Riemann-sum approximation is + +.. math:: + + W_x(a, b) \approx \frac{dt}{\sqrt{a}}\sum_n x[n]\,\psi^*\left(\frac{n\,dt-b}{a}\right). + +So a manual discrete convolution written in physical time includes the ``dt`` +factor from the integral approximation. If the scale is +parameterized in seconds (``a_sec = a*dt``), then +``dt/sqrt(a_sec) = sqrt(dt)/sqrt(a)``. Relative to a purely +sample-index-based implementation with scale ``a`` in samples, this appears as +an additional amplitude factor of ``sqrt(dt)`` (equivalently +``1/sqrt(fs)`` with ``fs = 1/dt``). + .. _Converting frequency: Converting frequency to scale for ``cwt`` diff --git a/pywt/_cwt.py b/pywt/_cwt.py index 0bea72c4..c0fd2c6d 100644 --- a/pywt/_cwt.py +++ b/pywt/_cwt.py @@ -81,6 +81,16 @@ def cwt(data, scales, wavelet, sampling_period=1., method='conv', axis=-1, Size of coefficients arrays depends on the length of the input array and the length of given scales. + The coefficients are independent of ``sampling_period``. The + ``sampling_period`` parameter only affects the values returned in + ``frequencies``. + + When comparing against a direct numerical approximation of the + continuous-time CWT integral, keep in mind that a sampled implementation in + physical time introduces a ``dt`` factor from the Riemann sum. Depending on + the normalization used in a manual calculation, this can appear as an + additional amplitude factor of ``sqrt(dt)`` (equivalently ``1/sqrt(fs)``). + Examples -------- >>> import pywt From b770b3b7a2b2c176b738b1cc2a629471a56282c9 Mon Sep 17 00:00:00 2001 From: Ralf Gommers Date: Mon, 17 Aug 2026 13:38:04 +0000 Subject: [PATCH 2/4] DOC: correct the CWT normalization explanation The previous explanation described a physical-time Riemann-sum implementation, which carries a factor sqrt(dt) *relative to* what cwt returns. The reported symptom is the opposite: a hand-written convolution that uses a scale in seconds in the 1/sqrt(a) prefactor while summing over raw samples is too large by sqrt(fs), and needs 1/sqrt(fs) to match. Both factors have the same magnitude but apply in opposite directions, so the text as written would lead readers to correct the wrong way. State instead what cwt actually computes, in units of samples: C[a, b] = 1/sqrt(a) * sum_n x[n] conj(psi((n - b)/a)) and derive the relation to the continuous-time transform from there, so that the direction of the sqrt(dt) factor is unambiguous. Name the mixed-convention mistake explicitly, since that is the one users hit. Also: - fix ``t_n = n\,dt``, which was inline literal markup and rendered the LaTeX escape verbatim - use distinct symbols for the scale in samples and in seconds; the old text switched between the two meanings of `a` within a few lines - retitle the section, since sampling_period is precisely the thing that does not affect the normalization, and drop the duplicate statement about it from the cwt docstring (it is already in the parameter docs) - move the section after the cwt_scaling_demo plot, which the text above it refers to as "the example below" - add a runnable example reproducing cwt with an explicit convolution, and document the two discretization effects that keep the agreement from being exact - add a regression test pinning the convention --- doc/source/ref/cwt.rst | 149 +++++++++++++++++++++++++------- pywt/_cwt.py | 21 +++-- pywt/tests/test_cwt_wavelets.py | 33 +++++++ 3 files changed, 165 insertions(+), 38 deletions(-) diff --git a/doc/source/ref/cwt.rst b/doc/source/ref/cwt.rst index c9d21d23..801bc3fa 100644 --- a/doc/source/ref/cwt.rst +++ b/doc/source/ref/cwt.rst @@ -231,35 +231,6 @@ scales. The right column are the corresponding Fourier power spectra of each filter.. For scales 1 and 2 it can be seen that aliasing due to violation of the Nyquist limit occurs. - -Normalization with ``sampling_period`` --------------------------------------- - -The ``sampling_period`` argument changes only the reported frequencies -(``scale2frequency(...)/sampling_period``). It does not rescale the CWT -coefficients themselves. - -This can be confusing when comparing against a hand-written approximation to -the continuous-time definition, - -.. math:: - - W_x(a, b) = \frac{1}{\sqrt{a}}\int x(t)\,\psi^*\left(\frac{t-b}{a}\right)\,dt. - -For sampled data ``t_n = n\,dt``, a direct Riemann-sum approximation is - -.. math:: - - W_x(a, b) \approx \frac{dt}{\sqrt{a}}\sum_n x[n]\,\psi^*\left(\frac{n\,dt-b}{a}\right). - -So a manual discrete convolution written in physical time includes the ``dt`` -factor from the integral approximation. If the scale is -parameterized in seconds (``a_sec = a*dt``), then -``dt/sqrt(a_sec) = sqrt(dt)/sqrt(a)``. Relative to a purely -sample-index-based implementation with scale ``a`` in samples, this appears as -an additional amplitude factor of ``sqrt(dt)`` (equivalently -``1/sqrt(fs)`` with ``fs = 1/dt``). - .. _Converting frequency: Converting frequency to scale for ``cwt`` @@ -285,3 +256,123 @@ of frequency directly. .. plot:: pyplots/cwt_scaling_demo.py + + +.. _CWT normalization: + +Normalization of the CWT coefficients +------------------------------------- + +``cwt`` works entirely in units of samples. Writing :math:`x[n]` for the input +signal and taking both the scale :math:`a` and the translation :math:`b` to be +expressed in samples, the returned coefficients are + +.. math:: + + C[a, b] = \frac{1}{\sqrt{a}}\sum_n x[n]\, + \psi^*\!\left(\frac{n - b}{a}\right). + +No sampling interval appears in this expression. In particular, the +``sampling_period`` argument of :func:`pywt.cwt` rescales only the returned +``frequencies``; the coefficients themselves do not depend on it. + +Relation to the continuous-time transform +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The continuous-time definition of the CWT, with a scale :math:`a_s` and a +translation :math:`b_s` in seconds, is + +.. math:: + + W_x(a_s, b_s) = \frac{1}{\sqrt{a_s}} + \int x(t)\,\psi^*\!\left(\frac{t - b_s}{a_s}\right)\,\mathrm{d}t. + +For data sampled at :math:`t_n = n\,\mathrm{d}t`, approximating that integral by +a Riemann sum introduces a factor :math:`\mathrm{d}t`. Substituting +:math:`a_s = a\,\mathrm{d}t` and :math:`b_s = b\,\mathrm{d}t` then gives + +.. math:: + + W_x(a_s, b_s) \approx \frac{\mathrm{d}t}{\sqrt{a\,\mathrm{d}t}} + \sum_n x[n]\,\psi^*\!\left(\frac{n - b}{a}\right) + = \sqrt{\mathrm{d}t}\, C[a, b]. + +So the coefficients returned by ``cwt`` have to be **multiplied** by +:math:`\sqrt{\mathrm{d}t}` (equivalently, divided by :math:`\sqrt{f_s}`) to be +expressed in physical-time units. PyWavelets does not apply that factor, +because doing so would make the coefficients depend on ``sampling_period``. + +A common source of confusion +~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +A frequent mistake when comparing a hand-written convolution against ``cwt`` is +to *mix* the two conventions above: to use a scale in seconds, +:math:`a_s = a\,\mathrm{d}t`, in the :math:`1/\sqrt{a_s}` prefactor, while +summing over raw samples without the :math:`\mathrm{d}t` coming from the +integral. Such a result is too large by a factor +:math:`1/\sqrt{\mathrm{d}t} = \sqrt{f_s}`, and an otherwise unexplained factor +:math:`1/\sqrt{f_s}` then has to be inserted by hand to make the amplitudes +agree with ``cwt``. Either keep everything in samples, as :math:`C[a, b]` does, +or write out the full Riemann sum including :math:`\mathrm{d}t` — but do not +combine the two. + +The example below evaluates :math:`C[a, b]` directly and compares it to +``cwt``: + +.. try_examples:: + + >>> import numpy as np + >>> import pywt + >>> fs = 3000.0 # sampling rate in Hz + >>> t = np.arange(3000) / fs + >>> x = np.sin(2 * np.pi * 40 * t) + >>> wavelet = pywt.ContinuousWavelet('cmor14-2') + >>> a = pywt.frequency2scale(wavelet, 40 / fs) # scale, in samples + >>> coefs, freqs = pywt.cwt(x, a, wavelet, sampling_period=1 / fs) + + Now the same transform, written out as a convolution over samples. ``Fb`` and + ``Fc`` are the bandwidth and center frequency of ``'cmor14-2'``: + + >>> Fb, Fc = 14.0, 2.0 + >>> lb, ub = wavelet.lower_bound, wavelet.upper_bound + >>> k = np.arange(int(a * (ub - lb)) + 1) + >>> u = lb + (k + 0.5) / a + >>> psi = np.exp(2j*np.pi*Fc*u) * np.exp(-u**2 / Fb) / np.sqrt(np.pi * Fb) + >>> conv = np.convolve(x, np.conj(psi)[::-1]) / np.sqrt(a) + >>> trim = (conv.size - x.size) // 2 + >>> manual = conv[trim:trim + x.size] + >>> rel = np.max(np.abs(manual - coefs[0])) / np.max(np.abs(coefs[0])) + >>> bool(rel < 1e-3) + True + +Note that :math:`\psi` is conjugated *and* reversed before the convolution, and +that it is evaluated at bin midpoints, ``(k + 0.5) / a``. The half-sample offset +is there because ``cwt`` convolves with the integral of :math:`\psi` and then +differences the result, which effectively averages :math:`\psi` over each +sample bin rather than sampling it pointwise. + +Limits of the agreement +~~~~~~~~~~~~~~~~~~~~~~~ + +:math:`C[a, b]` is what ``cwt`` computes in the limit of a finely sampled +wavelet. Two discretization effects prevent an analytic implementation from +reproducing it exactly: + +* :math:`\psi` is only evaluated over + ``[wavelet.lower_bound, wavelet.upper_bound]``, so its tails are truncated + (see :ref:`Choosing scales`). With the default bounds of :math:`[-8, 8]`, + ``cmor14-2`` still retains about 1% of its peak amplitude at the edges, and + that truncation dominates the residual in the example above. +* ``cwt`` resamples the precomputed integral of :math:`\psi` onto the grid for a + given scale by truncated indexing rather than by interpolation. The resulting + jitter grows with the scale and is reduced by raising ``precision``; it is the + "zipper-like" effect mentioned in the ``precision`` documentation of + :func:`pywt.cwt`. It largely averages out for narrowband signals such as the + one above, but is clearly visible for broadband input. + +Finally, note that the amplitude of the wavelets themselves does not follow a +single convention across families — ``mexh`` and ``gaus`` are normalized to unit +energy, ``morl`` carries no normalization constant at all, and ``cmor``, +``shan`` and ``fbsp`` use yet other conventions (see the formulas above). CWT +coefficient magnitudes are therefore not directly comparable between wavelet +families. diff --git a/pywt/_cwt.py b/pywt/_cwt.py index c0fd2c6d..928ef510 100644 --- a/pywt/_cwt.py +++ b/pywt/_cwt.py @@ -81,15 +81,18 @@ def cwt(data, scales, wavelet, sampling_period=1., method='conv', axis=-1, Size of coefficients arrays depends on the length of the input array and the length of given scales. - The coefficients are independent of ``sampling_period``. The - ``sampling_period`` parameter only affects the values returned in - ``frequencies``. - - When comparing against a direct numerical approximation of the - continuous-time CWT integral, keep in mind that a sampled implementation in - physical time introduces a ``dt`` factor from the Riemann sum. Depending on - the normalization used in a manual calculation, this can appear as an - additional amplitude factor of ``sqrt(dt)`` (equivalently ``1/sqrt(fs)``). + The transform is computed entirely in units of samples. With the scale + ``a`` and the translation ``b`` both expressed in samples, the returned + coefficients are:: + + C[a, b] = 1/sqrt(a) * sum_n data[n] * conj(psi((n - b)/a)) + + No sampling interval enters this expression, which is why ``coefs`` does + not depend on ``sampling_period``. To express the coefficients in + physical-time units instead, multiply them by ``sqrt(dt)`` (equivalently, + divide by ``sqrt(fs)``); that factor comes from approximating the integral + of the continuous-time transform by a Riemann sum. See the + :ref:`CWT normalization` section of the documentation for details. Examples -------- diff --git a/pywt/tests/test_cwt_wavelets.py b/pywt/tests/test_cwt_wavelets.py index 22113b33..b626fdd7 100644 --- a/pywt/tests/test_cwt_wavelets.py +++ b/pywt/tests/test_cwt_wavelets.py @@ -471,6 +471,39 @@ def test_cwt_method_fft(): assert_allclose(cfs_conv, cfs_fft, rtol=0, atol=1e-13) +def test_cwt_normalization_convention(): + # cwt works in units of samples: with the scale a and the shift b both + # given in samples, coefs[a, b] == 1/sqrt(a) * sum_n x[n] conj(psi((n-b)/a)) + # and no sampling interval enters. + fs = 3000. + Fb, Fc = 14., 2. + wavelet = pywt.ContinuousWavelet(f'cmor{Fb:g}-{Fc:g}') + x = np.sin(2 * np.pi * 40 * np.arange(3000) / fs) + scale = pywt.frequency2scale(wavelet, 40 / fs) + + cfs, freqs = pywt.cwt(x, scale, wavelet, sampling_period=1 / fs) + assert_allclose(freqs, [40.], rtol=1e-12) + + # psi is sampled at bin midpoints because cwt convolves with the integral + # of psi and then differences it, which averages psi over each sample bin. + lb, ub = wavelet.lower_bound, wavelet.upper_bound + k = np.arange(int(scale * (ub - lb)) + 1) + psi, _ = ref_cmor(lb + 0.5 / scale, ub + 0.5 / scale, k.size, Fb, Fc) + # the filter is conjugated *and* reversed; convolving with conj(psi) alone + # gives the complex conjugate of the correct result, which an abs() + # comparison would not catch + conv = np.convolve(x, np.conj(psi)[::-1]) / np.sqrt(scale) + trim = (conv.size - x.size) // 2 + manual = conv[trim:trim + x.size] + + assert_allclose(manual, cfs[0], atol=1e-3 * np.max(np.abs(cfs[0]))) + + # the coefficients themselves are unaffected by sampling_period + cfs_unit, freqs_unit = pywt.cwt(x, scale, wavelet, sampling_period=1.) + assert_allclose(cfs_unit, cfs, rtol=0, atol=0) + assert_allclose(freqs_unit * fs, freqs, rtol=1e-12) + + def test_continuous_wavelet_pickle(tmpdir): wavelet = pywt.ContinuousWavelet('cmor1.5-1.0') filename = os.path.join(tmpdir, 'cwav.pickle') From 0c4c0f111441a691a100942903b8b50f4320e6d6 Mon Sep 17 00:00:00 2001 From: Ralf Gommers Date: Mon, 17 Aug 2026 13:43:00 +0000 Subject: [PATCH 3/4] DOC: add frequency2scale to the "Other functions" reference page It is public in `__all__` and referenced from the CWT docs, but was the only one of the scale/frequency conversion helpers missing an autodoc entry. --- doc/source/ref/other-functions.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/doc/source/ref/other-functions.rst b/doc/source/ref/other-functions.rst index d9f757cc..414f1203 100644 --- a/doc/source/ref/other-functions.rst +++ b/doc/source/ref/other-functions.rst @@ -32,6 +32,8 @@ Central frequency of ``psi`` wavelet function .. autofunction:: scale2frequency +.. autofunction:: frequency2scale + Quadrature Mirror Filter ------------------------ From 445bca73c18e64475687e30e6453034f90632b0d Mon Sep 17 00:00:00 2001 From: Ralf Gommers Date: Wed, 2 Sep 2026 06:13:08 +0000 Subject: [PATCH 4/4] DOC: use the established third-level heading style in cwt.rst The new subsections used `~` underlines, but this file uses `^` for third-level headings, so docutils reported an inconsistent title style skipping from level 2 to 4. --- doc/source/ref/cwt.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/source/ref/cwt.rst b/doc/source/ref/cwt.rst index 801bc3fa..3d466adf 100644 --- a/doc/source/ref/cwt.rst +++ b/doc/source/ref/cwt.rst @@ -277,7 +277,7 @@ No sampling interval appears in this expression. In particular, the ``frequencies``; the coefficients themselves do not depend on it. Relation to the continuous-time transform -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The continuous-time definition of the CWT, with a scale :math:`a_s` and a translation :math:`b_s` in seconds, is @@ -303,7 +303,7 @@ expressed in physical-time units. PyWavelets does not apply that factor, because doing so would make the coefficients depend on ``sampling_period``. A common source of confusion -~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +^^^^^^^^^^^^^^^^^^^^^^^^^^^^ A frequent mistake when comparing a hand-written convolution against ``cwt`` is to *mix* the two conventions above: to use a scale in seconds, @@ -352,7 +352,7 @@ differences the result, which effectively averages :math:`\psi` over each sample bin rather than sampling it pointwise. Limits of the agreement -~~~~~~~~~~~~~~~~~~~~~~~ +^^^^^^^^^^^^^^^^^^^^^^^ :math:`C[a, b]` is what ``cwt`` computes in the limit of a finely sampled wavelet. Two discretization effects prevent an analytic implementation from