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:
- Naive higher: wavelet dilates too much, exiting the sampling frame far before its zero tails.
- 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).
- Naive lower: all's good.
- pywt lower: wavelet has small resolution.
- 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:
- 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.
pywt
integrates real Morlet via np.cumsum(psi) * step
, accounting for the differential step size
- The integrated wavelet,
int_psi
, is reused for all scales
- 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
:
- The resulting wavelet for increasing
scale
, side-by-side with naive (for later):
- The convolution result first
diff
'd,
- negated,
- multiplied by
sqrt(scale)
,
- 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:
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.