PyWavelets CWT implementation

Question

I seek to understand PyWavelets' implementation of the Continuous Wavelet Transform, and how it compares to the more 'basic' version I've coded and provided here. In particular:

1. How is integrated wavelet, then differenced post-convolving, equivalent to convolving without int+diff?
2. How is resampling the integrated wavelet at increasing resolution equivalent to changing its scale?

---

Naive implementation:

def cwt(x, win_len=None, win='morlet', real=False):
    N = len(x)
    win_len = win_len or N // 8

    scales = _scales(N, nv=32)
    coef = np.zeros((len(scales), N), dtype='complex128')
    kernel = morlet_kernel(win_len)
    wl2 = win_len // 2

    for i, scale in enumerate(scales):
        coef[i, :] = np.convolve(x, kernel(scale)[::-1])[wl2:-(wl2 - 1)]
    return coef

PyWavelets (clipped to relevant parts):

def pywt_cwt(data, scales):
    out = np.empty((np.size(scales),) + data.shape)
    int_psi, x = integrate_wavelet(wavelet='morl', precision=10)
    
    for i, scale in enumerate(scales):
        step = x[1] - x[0]
        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]

        conv = np.convolve(data, int_psi_scale)
        coef = - np.sqrt(scale) * np.diff(conv, axis=-1)
        
        d = (coef.shape[-1] - data.shape[-1]) / 2.
        coef = coef[..., floor(d):-ceil(d)]
        out[i, ...] = coef
    return out

---

Comparison on $f=1, 4$ sinusoids:

Answer 1

First see "Naive Breakdown" in a below section. Onto PyWavelets: the algorithm was found on Github to stem from an old MATLAB implementation, but it provides no details on coding the wavelet kernel. So I set on to explore; various visuals and code annotations here. Key findings follow. -- All code

See "PyWavelets Breakdown" below. It is implementing the following - namely Eq 4:

---

(1) How is integrated wavelet, then differenced post-convolving, equivalent to convolving without int+diff?

The formulation is in the right, deriving from the left, equations; Instead of continuously integrating over all $\mathcal{R}$, the integral is split into segments ssumming over $k$. The input is assumed constant over $[k, k+1]$ (though shouldn't it be $)$?), as it's discrete, whereas the wavelet is kept continuous. This enables to move the input outside the integral; visually:

So how do we code this? It's clear from Eq 4 that we need a diff somewhere. pywt does np.diff(conv), which is basically np.diff(coef); why does it seem to take difference of the would-be coefficient in Naive?

In Eq 4 we have the integrated wavelet up-to k + 1, minus up-to k, multiplied with the wavelet. This is equivalent to product with the wavelet at b - 1 and b, respectively (plug in t=k+1 and t=k, compare with plugging in b=b-1 and b=b instead`).

The answer's then in in conv; assume it was already trimmed so that it's only when the wavelet is at least halfway inside the signal. Then, every point in conv is the product of the wavelet with a portion of the input, at a different tau (or rather b). So two adjacent points in conv are products with wavelets at b and b + 1:

$$ \begin{align} \text{conv} &= [s \cdot \psi_{\text{int}}(t - 0),\ s \cdot \psi_{\text{int}}(t - 1),\ ...] \\ \Rightarrow \text{DIFF}(\text{conv}) &= [s \cdot (\psi_{\text{int}}(t - 1) - \psi_{\text{int}}(t - 0)),\ ...] \end{align} $$

However, we seek minus between t=k+1 and t, or a greater t and lower t, which is the negative of above, thus explaining step 7 in PW Breakdown.

---

(2) How is resampling the integrated wavelet at increasing resolution equivalent to changing its scale?

See figure under (5) in PW Breakdown. pywt appears to "scale" via the number of samples that define the wavelet, rather than recomputing the wavelet with a dilated argument. ... though, visually, this resampling clearly does dilate the wavelet, stretching it relative to the input.

So it's indeed a form of scaling; but how does it compare to naive? Naive is recomputed more intuitively per equations, using fixed wavelet length, and passing a scaled argument. Then, which is better? The opposite question is clearer to answer; consider what happens at higher and lower scales:

1. Naive higher: wavelet dilates too much, exiting the sampling frame far before its zero tails.
2. pywt higher: wavelet is just fine; the problem's with conv. Wavelet length is fixed at 1024, so if the input is any shorter, then higher scale wavelets can never fully multiply the signal. The greater the disparity, the more the wavelet is "seen" similar to "Naive higher" by the signal; this can be seen in the question's heatmaps differing by vertical shifts.
   • Maximum scale of the wavelet is also unlimited; len(j) grows proportionally with scale (but there's a normalization problem, see here).
3. Naive lower: all's good.
4. pywt lower: wavelet has small resolution.
5. Fixed vs variable window length: higher scales must correspond to greater dilation, or being-nonzero, of the wavelet relative to the input. Both accomplish this, except Naive discontinues to zero abruptly at higher scales.

The higher scale problems in both are solved by limiting the maximum scale; if both are in their completely safe zones, the Naive implementation seems preferable due to higher resolution. However, this "safe zone" for Naive is prohibitively small, if we are to consider "safe" as not jumping to zero at tails.

To contrary, pywt's problem is more easily solved by increasing maximum wavelet length, which was suggested to rid of the zipper-like artifacts (which can be seen in the question's plots). That'd, however, worsen higher-scale performance for short inputs - but in practice it's unlikely to be shorter than 1024.

---

PyWavelets Breakdown:

1. Wavelet, prior to integration, matches exactly with the shown code blob, which is an approximation of the complete real Morlet (used by Naive) assuming $\sigma > 5$ in the Wiki.
2. pywt integrates real Morlet via np.cumsum(psi) * step, accounting for the differential step size
3. The integrated wavelet, int_psi, is reused for all scales
4. For each scale, the same int_psi is resampled at increasing resolution, with j as an integer varying linearly between same min and max (almost) for all scales:

5. The resulting wavelet for increasing scale, side-by-side with naive (for later):

6. The convolution result first diff'd,
7. negated,
8. multiplied by sqrt(scale),
9. then trimmed on both sides to either (a) only include parts where the wavelet is at least halfway "inside the signal", or (b) force its length to len(x).

---

Naive Breakdown: if unfamiliar with CWT/STFT, I strongly recommend parts 1-3 here.

We need (1) wavelet type; (2) wavelet length; (3) wavelet increment/overlap. We'll use real Morlet, 96 samples, and increment by 20.

win_len = 96  # "win" == window == wavelet
win_inc = 20
n_wins = (len(x) - win_len) // win_inc + 1

coef = np.zeros((n_wins, len(scales)), dtype='complex128')
kernel = morlet_kernel(win_len)

for tau in range(n_wins):
    start = tau * win_inc
    end   = start + win_len
    coef[tau, :] = _transform(x[start:end], kernel, scales)

def _transform(x, kernel, scales):
    coef = np.zeros(len(scales), dtype='complex128')
    for i, scale in enumerate(scales):
        psi = np.conj(kernel(scale))
        coef[i] = np.sum(x * psi / np.sqrt(scale))
    return coef

Visualized for scale=2:

The difference between this and naive implementation in the question is, latter uses win_inc=1, equivalent to convolving with a flipped wavelet and trimming the result to only include parts where the wavelet was fully "inside the signal". Also, here we have tau ("increment") in the outer loop, there scale, but the two are equivalent.

Lastly, shapes along timeshift dim differ; question's implementation trims so that only parts of convolution where wavelet is at least halfway "inside the input" are included; this one trims further to where the full wavelet is included.

---

UNSOLVED:

1. Why * sqrt(scale) (step 8 pywt)? Eq 4 shows 1 / sqrt(scale). I considered * scale per chain rule, except we've integrated the wavelet numerically... Makes sense to have per-scale treatment; do we get away with not integrating a per-scale wavelet because the per-scale wavelet is... the same resampled wavelet? Resolved.