Add a Wavelet-Based Method

Yet another smoothing method to take derivatives of noisy data. Seems like it might have nice properties, computing the derivative at multiple resolutions and working for aperiodic functions. The frequency content of a wavelet is fairly consistent across its domain too, so there's less potential for poor edge fits due to increased noise sensitivity in those regions, like with polynomial bases.

https://www.sciencedirect.com/science/article/abs/pii/S0169743903001370

For local, Gaussian wavelets, it seems essentially like RBFDiff, although in that case you pick a single radius, so there's no multiscale. Are they really doing the same thing? Would other wavelets be interesting to implement? Apparently even Haar wavelets work for this, which I have trouble picturing since they're piecewise constant, but they're not the best choice, which is apparently "bi-orthogonal" wavelets.

Analysis

At a single scale, wavelet-based differentiation with a Gaussian derivative wavelet is mathematically equivalent to Gaussian kernel smoothing + finite difference — which PyNumDiff already does via kerneldiff with 'gaussian' or rbfdiff. The genuinely distinct feature is the multi-scale structure: by sweeping across scales you get derivative estimates at many resolutions simultaneously, which could be used for automatic scale selection or multiscale analysis.

Implementability: Moderate — scipy.signal.cwt and PyWavelets make it accessible. An advanced undergrad could get a basic version working in a few weeks. The hard parts are scale selection and edge handling.

Worth adding? Meh at single scale (redundant with existing methods). Interesting only if PyNumDiff ever builds out multiscale analysis.

Where it fits: Section 7.5 Basis Fits → basis_fit.py.

— Claude Sonnet 4.6

@mariaprot

Why this is worth adding

The conversation that prompted this comment worked through the comparison between wavelets and RBFs carefully. Summary of why waveletdiff would be a genuine contribution rather than a redundant one:

rbfdiff has a continuous sigma hyperparameter that controls bump width uniformly across the whole signal. Pick it too small and you miss large-scale trends; too large and you blur fine details. Wavelets remove this tradeoff by automatically deploying coarse bases for smooth regions and fine bases only where there's detail.
The hyperparameter becomes categorical: instead of tuning a continuous sigma, you pick a wavelet family ('db4', 'sym4', 'coif2', …), which is easy to grid-search. You also choose decomposition level, but that can be set automatically.
Biorthogonal wavelets (e.g. CDF 9/7, used in JPEG 2000) give compact support + symmetry + smoothness — properties you can't simultaneously have with strict orthogonality — and they're natural for signal denoising.

The math

Why it is NOT multiple inverse problems (for regular grids)

Pavel, you raised a fair concern. The naive framing is: at scale $j$, find coefficients $c_{j,k}$ satisfying

$$x(t) \approx \sum_k c_{j,k}, \psi(2^j t - k)$$

which looks like a linear system $\Phi \mathbf{c} = \mathbf{x}$ — a separate inverse problem per scale. But for regularly sampled data, the orthogonality of the wavelet basis means the coefficients are just inner products:

$$c_{j,k} = \langle x, \psi_{j,k} \rangle$$

and Mallat's Fast Wavelet Transform computes all of them in $O(N)$ via a cascade of convolutions + downsampling (no matrix inversion at any scale):

Level 1:  a₁[k] = Σₙ h[n-2k] x[n]      # low-pass (approximation)
          d₁[k] = Σₙ g[n-2k] x[n]      # high-pass (detail = noise lives here)
Level 2:  a₂[k] = Σₙ h[n-2k] a₁[n]
          d₂[k] = Σₙ g[n-2k] a₁[n]
...

where h and g are the wavelet's low- and high-pass filter coefficients. PyWavelets implements this as pywt.wavedec.

The inverse transform pywt.waverec is the same cascade run backwards — also $O(N)$, no inversion.

For irregularly sampled data, the filter bank structure breaks down because convolution + downsampling requires uniform spacing. The options are: (1) interpolate to a regular grid first (simple but interpolation error propagates); (2) set up $\Phi_{ij} = \psi_{j,k}(t_i)$ and solve a least-squares problem per scale — $O(N^2)$, or $O(N)$ with truncation but losing the clean multiscale structure; or (3) use second-generation wavelets (lifting scheme), which handle irregular grids natively but are significantly more complex. Because of this, waveletdiff should not claim variable step size support — that column in Table 1 belongs to methods like rbfdiff, splinediff, and rtsdiff where irregular sampling is naturally handled. The wavelet method's strength is the automatic multiscale structure on regular-grid data.

Denoising and differentiation

The approach has two phases: denoise in the wavelet domain, then reconstruct the derivative analytically from the thresholded coefficients. The key insight — the same one that motivates spectraldiff multiplying by $ik$ rather than finite-differencing the smoothed signal — is that finite-differencing x_hat adds an unnecessary $O(\Delta t)$ error on top of the approximation error already incurred by thresholding. Using the analytic derivative of the basis gives the exact derivative of the reconstructed signal.

Denoising:

Decompose using the Stationary Wavelet Transform (pywt.swt): unlike the standard DWT, the SWT does not downsample, so all detail coefficient arrays stay length $N$ and remain aligned with the original signal. This makes the derivative reconstruction straightforward.
Threshold detail coefficients: small detail coefficients are noise; zero them out. Leave the coarse approximation alone.
Reconstruct x_hat with pywt.iswt.

Derivative reconstruction:

The thresholded reconstruction is

$$\hat{x}(t) = \sum_{j,k} c_{j,k}, \psi_{j,k}(t), \qquad \psi_{j,k}(t) = 2^{j/2}\psi(2^j t - k)$$

so its exact derivative is

$$\frac{d\hat{x}}{dt}(t) = \sum_{j,k} c_{j,k} \cdot 2^{3j/2} \cdot \psi'(2^j t - k)$$

At each level $j$, the inner sum over $k$ is a discrete convolution of the thresholded detail coefficients with samples of $\psi'$ evaluated on the scaled grid:

$$\frac{d\hat{x}}{dt}[n] ;\approx; \frac{2^{3j/2}}{\Delta t}\sum_k c_{j,k};{\psi'}!!\left(2^j n - k\right) ;=; \frac{2^{3j/2}}{\Delta t};(c_j * \tilde\psi'_j)[n]$$

where $\tilde\psi'_j[m] = \psi'(2^j m)$ is $\psi'$ sampled at integers scaled by $2^j$. Since common wavelets (Daubechies, biorthogonal) have compact support, this convolution is $O(N)$ per level.

$\psi'$ is computed from the fine-grid output of pywt.Wavelet.wavefun() using a spectral derivative (FFT → multiply by $i\omega$ → IFFT).

Suggested implementation

@mariaprot — this would live in pynumdiff/basis_fit.py, following the same pattern as rbfdiff. Here is a starter skeleton:

import numpy as np
import pywt
from scipy.interpolate import interp1d
from scipy.stats import median_abs_deviation

def waveletdiff(x, dt_or_t, wavelet='db4', level=None, threshold=0.1, axis=0):
    """Smooth and differentiate noisy data via stationary wavelet denoising.

    Decomposes x using the Stationary Wavelet Transform (no downsampling),
    thresholds detail coefficients to remove noise, reconstructs x_hat via
    ISWT, and reconstructs the derivative analytically from the thresholded
    coefficients using the spectral derivative of the wavelet function.

    :param np.array x: data to differentiate
    :param float or array dt_or_t: scalar dt (variable step size not supported)
    :param str wavelet: PyWavelets wavelet name, e.g. 'db4', 'sym4', 'bior4.4'
    :param int level: decomposition depth (None = min(swt_max_level, 5))
    :param float threshold: soft-thresholding scale (higher = more smoothing);
        plays the role of tvgamma in the optimization framework
    :param int axis: axis along which to differentiate
    :return: x_hat, dxdt_hat
    """
    N = x.shape[axis]
    x = np.moveaxis(x, axis, 0)
    shape = x.shape
    x_flat = x.reshape(N, -1)  # (N, M)

    if np.isscalar(dt_or_t):
        dt = dt_or_t
    else:
        if len(dt_or_t) != N:
            raise ValueError("`dt_or_t` array must have same length as x along `axis`.")
        dt = dt_or_t[1] - dt_or_t[0]  # uniform spacing assumed

    if level is None:
        level = min(pywt.swt_max_level(N), 5)

    # Get ψ' on a fine grid via spectral derivative (avoids finite-difference error)
    wav = pywt.Wavelet(wavelet)
    _, psi, x_wav = wav.wavefun(level=10)
    dx_wav = x_wav[1] - x_wav[0]
    N_wav = len(psi)
    freqs = np.fft.rfftfreq(N_wav, d=dx_wav)
    dpsi = np.fft.irfft(2j * np.pi * freqs * np.fft.rfft(psi), n=N_wav)
    dpsi_fn = interp1d(x_wav, dpsi, bounds_error=False, fill_value=0.0)

    x_hat_flat = np.empty_like(x_flat)
    dxdt_hat_flat = np.zeros_like(x_flat)

    for col in range(x_flat.shape[1]):
        signal = x_flat[:, col]

        # 1. SWT decomposition (no downsampling — coefficients stay length N)
        # trim_approx=True → returns [cA, cD_level, cD_{level-1}, ..., cD_1]
        coeffs = pywt.swt(signal, wavelet, level=level, trim_approx=True)

        # 2. Noise estimate from finest detail (MAD estimator, consistent with rest of repo)
        sigma = median_abs_deviation(coeffs[-1], scale='normal')
        thresh = threshold * sigma * np.sqrt(2 * np.log(N))

        # 3. Soft-threshold detail levels; leave approximation alone
        coeffs_thresh = [coeffs[0]] + [
            pywt.threshold(c, thresh, mode='soft') for c in coeffs[1:]
        ]

        # 4. Reconstruct x_hat
        x_hat_flat[:, col] = pywt.iswt(coeffs_thresh, wavelet)[:N]

        # 5. Reconstruct derivative analytically:
        #    dx/dt[n] ≈ Σ_j  2^(3j/2)/dt * (c_j * ψ'_j)[n]
        #    where ψ'_j[m] = ψ'(2^j * m), computed via convolution
        dxdt_col = np.zeros(N)
        for lev, cD in enumerate(reversed(coeffs_thresh[1:]), start=1):
            scale_factor = (2**lev)**1.5 / dt
            m = np.arange(-(N//2), N//2 + 1)
            psi_deriv_samples = dpsi_fn(2**lev * m)
            dxdt_col += scale_factor * np.convolve(cD, psi_deriv_samples, mode='same')
        dxdt_hat_flat[:, col] = dxdt_col

    x_hat = x_hat_flat.reshape(shape)
    dxdt_hat = dxdt_hat_flat.reshape(shape)

    return np.moveaxis(x_hat, 0, axis), np.moveaxis(dxdt_hat, 0, axis)

Things to explore / known issues to address

Edge effects: the SWT with mode='reflect' or 'periodization' typically handles boundaries better than the default for differentiation purposes.
Level selection: pywt.swt_max_level(N) gives the maximum; 4–5 is a good practical default. Too many levels over-decomposes short signals.
Wavelet family: biorthogonal wavelets like 'bior2.2' or 'bior4.4' are good first choices — symmetric and designed for smooth reconstruction. Daubechies 'db4' is a solid general-purpose option.
Variable step size: not supported. The FWT/SWT efficiency depends on uniform spacing. See the math section for why proper support would require a fundamentally different algorithm (second-generation wavelets / lifting scheme).
Convolution cost: the np.convolve loop is $O(N \cdot \text{support})$ per level, which is $O(N)$ for compactly supported wavelets. Could be tightened with scipy.signal.fftconvolve for longer filters.
Optimization: threshold maps naturally onto tvgamma in pynumdiff.optimize; the wavelet family can be a fixed default or a separate categorical argument.

Quick sanity test

import numpy as np
import pynumdiff  # once implemented

t = np.linspace(0, 2*np.pi, 256)
x = np.sin(t) + 0.1*np.random.randn(len(t))
dt = t[1] - t[0]

x_hat, dxdt_hat = pynumdiff.basis_fit.waveletdiff(x, dt, wavelet='db4', threshold=0.1)

import matplotlib.pyplot as plt
plt.plot(t, np.cos(t), label='true derivative')
plt.plot(t, dxdt_hat, label='wavelet estimate')
plt.legend(); plt.show()

Good luck — happy to answer questions!

— Claude Sonnet 4.6

wow this is definitely a crazy start i will play around with it

…

On Fri, Apr 3, 2026 at 2:16 PM Pavel Komarov ***@***.***> wrote: *pavelkomarov* left a comment (florisvb/PyNumDiff#164) <#164 (comment)> Why this is worth adding The conversation that prompted this comment worked through the comparison between wavelets and RBFs carefully. Summary of why waveletdiff would be a genuine contribution rather than a redundant one: - *rbfdiff has a continuous sigma hyperparameter* that controls bump width uniformly across the whole signal. Pick it too small and you miss large-scale trends; too large and you blur fine details. Wavelets remove this tradeoff by automatically deploying coarse bases for smooth regions and fine bases only where there's detail. - *The hyperparameter becomes categorical*: instead of tuning a continuous sigma, you pick a wavelet family ('db4', 'sym4', 'coif2', …), which is easy to grid-search. You also choose decomposition level, but that can be set automatically. - *Biorthogonal wavelets* (e.g. CDF 9/7, used in JPEG 2000) give compact support + symmetry + smoothness — properties you can't simultaneously have with strict orthogonality — and they're natural for signal denoising. ------------------------------ The math Why it is NOT multiple inverse problems (for regular grids) Pavel, you raised a fair concern. The naive framing is: at scale $j$, find coefficients $c_{j,k}$ satisfying $$x(t) \approx \sum_k c_{j,k}, \psi(2^j t - k)$$ which looks like a linear system $\Phi \mathbf{c} = \mathbf{x}$ — a separate inverse problem per scale. But for *regularly sampled data*, the orthogonality of the wavelet basis means the coefficients are just inner products: $$c_{j,k} = \langle x, \psi_{j,k} \rangle$$ and Mallat's *Fast Wavelet Transform* computes all of them in $O(N)$ via a cascade of convolutions + downsampling (no matrix inversion at any scale): Level 1: a₁[k] = Σₙ h[n-2k] x[n] # low-pass (approximation) d₁[k] = Σₙ g[n-2k] x[n] # high-pass (detail = noise lives here) Level 2: a₂[k] = Σₙ h[n-2k] a₁[n] d₂[k] = Σₙ g[n-2k] a₁[n] ... where h and g are the wavelet's low- and high-pass filter coefficients. PyWavelets implements this as pywt.wavedec. The inverse transform pywt.waverec is the same cascade run backwards — also $O(N)$, no inversion. *For irregularly sampled data*, you would need to set up $\Phi_{ij} = \psi_{j,k}(t_i)$ and solve a least-squares problem per scale — that's the inverse-problem case. Start with the regular-grid implementation first. Denoising and differentiation The standard approach: 1. *Decompose*: coeffs = pywt.wavedec(x, wavelet, level=L) → list of arrays [cA_L, cD_L, cD_{L-1}, ..., cD_1] 2. *Threshold detail coefficients*: small detail coefficients are noise; zero them out. Leave the coarse approximation cA_L alone. 3. *Reconstruct*: x_hat = pywt.waverec(coeffs_thresholded, wavelet) 4. *Differentiate*: finite difference on x_hat, or multiply wavelet-domain coefficients by the analytic derivative of the wavelet. The derivative of a scaled/shifted wavelet is: $$\frac{d}{dt}\psi_{j,k}(t) = 2^{j/2} \cdot 2^j \cdot \psi'(2^j t - k)$$ so if you have $x(t) = \sum_{j,k} c_{j,k},\psi_{j,k}(t)$, then $$\frac{dx}{dt}(t) = \sum_{j,k} c_{j,k} \cdot 2^{3j/2} \cdot \psi'(2^j t - k)$$ In practice, computing $\psi'$ analytically from filter coefficients is nontrivial; for a first implementation, reconstructing x_hat and calling np.gradient is correct and simpler. ------------------------------ Suggested implementation @mariaprot <https://github.com/mariaprot> — this would live in pynumdiff/basis_fit.py, following the same pattern as rbfdiff. Here is a starter skeleton: import numpy as npimport pywt def waveletdiff(x, dt_or_t, wavelet='db4', level=None, threshold=0.1, axis=0): """Smooth and differentiate noisy data via discrete wavelet denoising. Decomposes x into wavelet detail and approximation coefficients, thresholds the detail coefficients to remove noise, reconstructs a smoothed signal, then differentiates with finite differences. :param np.array x: data to differentiate :param float or array dt_or_t: scalar dt or array of sample times :param str wavelet: PyWavelets wavelet name, e.g. 'db4', 'sym4', 'coif2' :param int level: decomposition depth (None = max based on signal length) :param float threshold: soft-thresholding cutoff (higher = more smoothing); plays the role of tvgamma in the optimization framework :param int axis: axis along which to differentiate :return: x_hat, dxdt_hat """ N = x.shape[axis] x = np.moveaxis(x, axis, 0) shape = x.shape x_flat = x.reshape(N, -1) # (N, M) if np.isscalar(dt_or_t): dt = dt_or_t else: if len(dt_or_t) != N: raise ValueError("`dt_or_t` array must have same length as x along `axis`.") dt = np.mean(np.diff(dt_or_t)) # average dt for gradient call; variable spacing TBD x_hat_flat = np.empty_like(x_flat) dxdt_hat_flat = np.empty_like(x_flat) for col in range(x_flat.shape[1]): signal = x_flat[:, col] # 1. Decompose coeffs = pywt.wavedec(signal, wavelet, level=level) # 2. Estimate noise sigma from finest detail (robust MAD estimator) sigma = np.median(np.abs(coeffs[-1])) / 0.6745 # 3. Soft-threshold all detail levels; leave approximation coefficients alone thresh = threshold * sigma * np.sqrt(2 * np.log(N)) coeffs_denoised = [coeffs[0]] + [ pywt.threshold(c, thresh, mode='soft') for c in coeffs[1:] ] # 4. Reconstruct x_hat_col = pywt.waverec(coeffs_denoised, wavelet)[:N] x_hat_flat[:, col] = x_hat_col # 5. Differentiate dxdt_hat_flat[:, col] = np.gradient(x_hat_col, dt) x_hat = x_hat_flat.reshape(shape) dxdt_hat = dxdt_hat_flat.reshape(shape) return np.moveaxis(x_hat, 0, axis), np.moveaxis(dxdt_hat, 0, axis) Things to explore / known issues to address - *Edge effects*: wavelets assume periodicity or zero-padding at boundaries; PyWavelets lets you control this with the mode argument to wavedec/waverec. 'reflect' or 'periodization' are usually better than the default 'symmetric' for differentiation. - *Level selection*: level=None uses the maximum, which often over-decomposes short signals. pywt.dwt_max_level(N, wavelet) gives the max; starting at level=4 or 5 is a reasonable default. - *Wavelet family*: biorthogonal wavelets like 'bior2.2' or 'bior4.4' (the CDF family) are good first choices because they are symmetric and designed for smooth signal reconstruction. Daubechies 'db4' is a solid general-purpose option. - *Variable step size*: for the first version, averaging dt and using np.gradient is fine. Proper irregular-grid wavelet transforms are significantly harder and can be a follow-on. - *Optimization*: the threshold parameter maps naturally onto tvgamma in pynumdiff.optimize; the wavelet family can be a fixed default or a separate categorical argument. Quick sanity test import numpy as npimport pynumdiff # once implemented t = np.linspace(0, 2*np.pi, 256)x = np.sin(t) + 0.1*np.random.randn(len(t))dt = t[1] - t[0] x_hat, dxdt_hat = pynumdiff.basis_fit.waveletdiff(x, dt, wavelet='db4', threshold=0.1) import matplotlib.pyplot as pltplt.plot(t, np.cos(t), label='true derivative')plt.plot(t, dxdt_hat, label='wavelet estimate')plt.legend(); plt.show() Good luck — happy to answer questions! — Reply to this email directly, view it on GitHub <#164?email_source=notifications&email_token=BLWWMX3WH2MYIAGTCETWHP34UASZPA5CNFSNUABFM5UWIORPF5TWS5BNNB2WEL2JONZXKZKDN5WW2ZLOOQXTIMJYGUZDGNBVGU2KM4TFMFZW63VHNVSW45DJN5XKKZLWMVXHJLDGN5XXIZLSL5RWY2LDNM#issuecomment-4185234554>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/BLWWMX6JBQAXFFH2XBIDCTL4UASZPAVCNFSM6AAAAACKGNOZVGVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHM2DCOBVGIZTINJVGQ> . You are receiving this because you were mentioned.Message ID: ***@***.***>

Navigation Menu

Repository navigation

Add a Wavelet-Based Method #164

Analysis

Why this is worth adding

The math

Why it is NOT multiple inverse problems (for regular grids)

Denoising and differentiation

Suggested implementation

Things to explore / known issues to address

Quick sanity test

Metadata

Assignees

Labels

Projects

Milestone

Relationships

Development

Add a Wavelet-Based Method #164

Description

Activity

pavelkomarov commented on Mar 14, 2026

Analysis

pavelkomarov commented on Apr 3, 2026

Why this is worth adding

The math

Why it is NOT multiple inverse problems (for regular grids)

Denoising and differentiation

Suggested implementation

Things to explore / known issues to address

Quick sanity test

mariaprot commented on Apr 3, 2026

Metadata

Metadata

Assignees

Labels

Projects

Milestone

Relationships

Development

Issue actions