PyWavelets CWT: normalization? Vs Scipy?

Question

Related. The equation being implemented normalizes by sqrt(1 / scale):

$$ C_{a, b} = \frac{1}{\sqrt{a}} \sum_k s(k)\left( \int_{-\infty}^{k+1} \overline{\psi \left(\frac{t-b}{a} \right) }dt - \int_{-\infty}^{k} \overline{\psi \left(\frac{t-b}{a} \right) }dt \right) $$

so why does pywt._cwt do * sqrt(scale)? Also, how is it normalizing (L1, L2, etc)?

---

Edit: After reviewing Scipy's cwt, I realized it's implementing exactly what I described as "ideal" in the answer, thus it's same as comparing the two (except a few obvious details). I don't know in which scenarios pywt's algorithm is better, but generally - easy win for Scipy. I've opened a PR to address pywt's key flaw.

Edit2: Found an apparent edge of pywt's over scipy's in the least expected place: low scales. Addressed here.

Edit 3: after fixing the low precision in pywt, I actually found pywt superior to scipy in all of a variety of synthetic signals, but squeezepy's always beat both. Tests here, with more coming.

Answer 1

Let's follow the math from incubation to delivery. It begins with psi, a rescaled morlet2 (as shown previously) at a scale $a=64$, and $\sigma=5$:

$$ \psi = \psi_{\sigma}(t/a) = e^{j\sigma (t/a)} e^{-(t/a)^2/2} \tag{2} $$

gets integrated and L1-normalized: -- (see caveat2 below)

$$ \psi = \psi_{\text{int}}(t) = \frac{1}{a} \int \psi_{\sigma}(t/a)\ dt \tag{3} $$

The norm is off by a scalar factor, but what matters is it's same for different $a$; you can verify it's so (for large $a$, increase integration limits). Note that it's the wavelet that's L1-normalized, not the integrated wavelet $\psi_\text{int}(t)$; the distinction is essential, as we'll see.

$\psi$ is then convolved with the signal $s$, and finite-differenced and negated; from here, this is:

$$ -\text{DIFF}(\text{conv}) = [s \cdot (\psi_{\text{int}}(t - 0) - \psi_{\text{int}}(t - 1)),\ ...] \tag{4} $$

Inspecting closely, this is equivalent to undoing the former integration, with certain numeric differences (addressed later), so we're left with

$$ \psi = \frac{1}{a} \psi_{\sigma}(t/a) \tag{5} $$

Lastly, * sqrt(scale):

$$ \boxed{\psi = \frac{1}{\sqrt{a}} \psi_{\sigma}(t/a)} \tag{6} $$

which is thus ultimately L2-normalized. We can verify this on a test signal; L2 is biased in amplitude/power representation toward lower frequencies, suppressing high frequencies. If above is right, not only can we see this, but apply a fix by removing * sqrt(scale) from code:

Implementation-specific caveats discussed next.

---

How reliable is the normalization across scales? -- Not great (-- code2):

On right-hand side we see the best case using limited number of samples, on left we see the resampling implementation. Clearly due to inherently poor resolution of small number of samples, our norms don't do well for low scales - but the recomputed variant at least maintains a mean norm approximately constant across them.

More interestingly, the resampled one suddenly appears quadratic after $\log_2(a)=6$; what's the deal? Recall "It begins ... at scale $a=64$"; this is when we resample with more samples than in original (int_psi):

whose diff is

• In 1LR we see int_psi at scale=256 is resampling the same values from int_psi multiple times, which manifests as zeros in diff (2LR).
• Next we square the values to get an idea of the L2 norm, and compare against the L2-normalized values in 3LR; we find normed values blowing unnormed values out of water. Is this excessive, or appropriate?
• Compare against recomputed at same scale; in 4LR, we at least see recomputed's shape - it's no longer negligible. This suggests the resampling to be overestimating wavelet amplitude, and since in norm computation this is significantly offset by half the values being zero due to diff, then the overestimation must be great.

How great? We can compare peak wavelet amplitudes across scales of the unnormed resampled diff(int_psi_256) and recomputed w (morlet2/scale); norming changes both by the same factor (*sqrt(scale)), won't make a difference.

Surely enough, again at $a=64$ we see critical behavior; peak amplitude of resampled saturates at a minimum, whereas the recomputed continues to decrease. The difference at $a=256$ is more than great enough (~12db=4x, if squared -> 16x) to offset the zeros in diff, hence why the L2 norm grows.

But notice the L1 norm remains well-behaved at higher scales. To find out why, we must precisely compare the proportion of zeros with the overestimation factor. We already know latter to be 4x; for former, I took a non-flat region (otherwise it's zeros anyway) and found ratio of zeros to number of samples: 0.75. Which means, 1 in 4 values are non-zero, and each are 4x greater ... the exact proportion needed to retain the L1 norm.

That zeros proportion isn't very clear from above plots - so instead, dots:

---

If diff undoes cumsum, why do either at all?

Good question. In fact, by ridding of both, I found the quality degradation (zipper) effect disappear:

But something else also disappears - the high frequency scales (left), and the sign flips (I undid this by removing minus at coef = -). I partly offset this by replacing * sqrt(scale) by 1 / sqrt(scale) (right), but that gets us to the L2 normed case (1 / scale is overkill), and is a baseless modification in light of everything found so far. So can it be done? Sure, one can implement without the int-diff dance, but I won't be exploring it here.

---

Caveat: effective norm

Recall, conv gets trimmed to len(x), so the entire wavelet may not participate in computation, whereas above computations assume otherwise. The assumption is truer the longer x is - but here's what happens if len(x) == 512, with same wavelet precision=10 (2^10=1024 int_psi), compared against recomputed wavelet:

A mess! Both distort significantly, possibly implying it's inherent to great enough scales pinned against short signals. A simple workaround is zero-padding the signal, but I wouldn't vouch for it without understanding the caveats (which I don't for wavelets, but do for DFT).

Furthermore, note that the norm isn't only per-scale, but also per-shift ($b$); less of the wavelet is in the signal when it's just entering as opposed to when its center is at signal's center.

But keeping it per-scale, I confess I cheated in the earlier example with $f=4$; initially I picked $f=1$, but the L1 norm plot didn't look nice, and I didn't understand why. Since now I do, here are both L1 and L2 per scale for $f=1,64$:

Compare against earlier plot; note how the left modes (low scales) of both norms are unaffected, only the high scale (low frequency) ones are now much greater (L2) and lower (L1), just as expected.

---

Caveat 2: * step doesn't do the L1-normalizing

Instead, it results from diffing the resampled integrated wavelet; the * step only sets the correct initial scale, and can be omitted without losing the norm. This does not happen with recomputing the wavelet, or resampling the unintegrated wavelet.

(Note that recomputed scipy.signal.morlet2 is L2-normalized, but only because it's explicitly normalized; the diff-resample-int is "self-normalized").