I've seen many different implementations of CWT, especially in generating the wavelets. Often these implementations produce significantly different outputs - e.g. scipy vs PyWavelets vs ssqueezepy.
How can we tell which implementation is "better", or if an implementation is valid in the first place? Are there any implementable mathematical criteria or visualizations?
Wavelets isn't just sampling with scales from some min to max - but it is what many implementations do, including scipy and PyWavelets, which yield serious flaws.
Along understanding what's possible, it helps to automate the process - to this end I've written a function, validate_filterbank, that prints a comprehensive report for the input filterbank. Below each section is supplemented with a relevant excerpt of this report.
This post covers:
== ANALYTICITY =====================================================================
Found not strictly analytic filter(s); threshold for ratio of `spectral sum`
to `spectral sum past Nyquist` is 1000.0 - got (less is worse):
psi_fs[0]: 4.2
psi_fs[1]: 283.0
psi_fs[63]: 58.5
== LP-SUM (with phi) ===============================================================
LP sum falls below threshold of 1.0 (for real inputs) by at most 0.530 (more is worse;
~1.0 implies ~zero capturing of the frequency!) at following frequency bin indices
(0 to 2048, shown skipping every 1 values):
[ 1 3 5 7 9 11 14 16 2012 2014 2016 2018 2020 2022
2024 2026 2028 2030 2032 2034 2036 2038 2040 2042 2044 2046 2048]
Above is rather generous for simplicity. Actually generic filterbank:
== LP-SUM (with phi) ===============================================================
LP sum exceeds threshold of 2 (for real inputs) by at most 0.871 (more is worse)
at following frequency bin indices (0 to 2048):
[2 4 5]
The lowest frequencies are worth dedicated inspection:
== REDUNDANCY ======================================================================
Found filters with redundancy exceeding 0.4 (energy overlap relative to sum of
individual energies) -- This isn't necessarily bad. Showing up to 20 filters:
psi_fs[63] & psi_fs[64]: 0.789
psi_fs[66] & psi_fs[67]: 0.804
psi_fs[68] & psi_fs[69]: 0.924 ...
Found filters with duplicate peak frequencies! Showing up to 20 filters:
psi_fs[63], peak_idx=8
psi_fs[64], peak_idx=8
psi_fs[66], peak_idx=6 ...
== DECAY (check for pure sines) ====================================================
Found filter(s) that are pure sines! Threshold for ratio of Fourier peak to
next-highest value is 1000.0 - got (more is worse):
psi_fs[69]: 1.04e+03
psi_fs[71]: 1.03e+03
psi_fs[74]: 1.16e+03 ...
== DECAY (boundary effects) ========================================================
Found filter(s) with incomplete decay (will incur boundary effects), with
following ratios of amplitude max to edge (less is worse; threshold is 1000.0):
psi_fs[52]: 580.6
psi_fs[53]: 107.7
psi_fs[54]: 51.9 ...
Lowpass filter has incomplete decay (will incur boundary effects), with
following ratio of amplitude max to edge: 1.0 > 1000.0
== ALIASING ========================================================================
Found Fourier peaks that are spaced neither exponentially nor linearly,
suggesting possible aliasing.
psi_fs[n], n=[20 18 17]
As another perspective on temporal decay, we plot the time-domain envelopes for all wavelets:
Much better. Time-domain:
Hence, padding isn't only about boundary effects, but wavelet behavior.
Since all wavelets peak at 1, we might reason to make the energy overlap peak at 2, just multiply them by $\sqrt{2}$:
but due to the overlaps, they'll exceed 2. However, the non-CQT region
lucks out, as it will generally have less overlap. Depending on design, it's also easy to under-tile it:
LP sum falls below threshold of 1.0 (for real inputs) by at most 0.699
(more is worse; ~1.0 implies ~zero capturing of the frequency!) at following
frequency bin indices (0 to 16384):
[1 2 3 7 8 11 12 13 16 17 21 22 16384]
Commonly used, makes all wavelets' energies the same, which makes intuitive sense on surface. Yet, and as MATLAB agrees ("L1 Norm"), it's a bad idea:
Ouch! Note, we can directly invert the (unsubsampled) CWT via sum(CWT, axis=0), i.e. combining all rows (one intergral inverse). Hence, the overlap-added wavelets in frequency domain are CWT's transfer function, and the LP sum is the energy transfer function, so we've dramatically amplified and attenuated various frequencies. L2's inverse formula also has a variable rescaling term, while L1's doesn't, hence L1's is "closer to inversion". Thus, L1 more faithfully captures the input.
But why's it happen, especially if energy norm is supposed to "equalize" the wavelets? Briefly, L1 norm is actually the "information norm". Feel free to open a question.
Despite Morlets being non-compact, above is not invertible, as Gauss decays below float epsilon very quickly - the untiled regions are "float zero". Even if they were invertible, it defeats the purpose of the transform to attenuate time-frequency geometry into negligibility.
Common cause is mis-specifying the scales parameter to cwt functions.
== LP-SUM (with phi) ===============================================================
LP sum falls below threshold of 1.0 (for real inputs) by at most 1.000
(more is worse; ~1.0 implies ~zero capturing of the frequency!) at following ...
[1 548 1095 1642 2190 2737 3282 3827 4377 4924 5471 6012 ...]
The Constant-Q Transform, or CQT, is defined as having constant $Q=\xi / \sigma = $ (center frequency) / bandwidth. It's what's automatically satisfied by exponentially spaced scales plugged into standard Morlet formulae, but what can be easy to miss or realize: CWT is not just exponential in frequency, or exponential in both frequency and bandwidth, but exponential in both in a fixed ratio.
== FREQUENCY-BANDWIDTH TILING ======================================================
Found non-CQT wavelets in upper quarters of frequencies - i.e.,
`(center freq) / bandwidth` isn't constant:
psi_fs[32], Q=1.9328358208955223
psi_fs[34], Q=2.2686567164179103
psi_fs[35], Q=2.462686567164179 ...
Earlier filterbanks had a sort of "critical" redundancy. Below is high redundancy:
Whether this is beneficial is application-dependent:
== REDUNDANCY ======================================================================
Found filters with redundancy exceeding 0.4 (energy overlap relative to sum of
individual energies) -- This isn't necessarily bad. Showing up to 20 filters:
psi_fs[0] & psi_fs[1]: 0.679
psi_fs[1] & psi_fs[2]: 0.679
psi_fs[2] & psi_fs[3]: 0.679 ...
Note: showing full frequency axis, unlike in most plots above (where right half is appropriately zero)
This is scipy's filterbank, with scales picked favorably for scipy (it has no auto-scales, and if we exceed bounds things get even worse).
Non-zero phase illustration:
Operating on scipy's complex outputs will be bad news.
Non-zero mean illustration:
Look similar? They shouldn't:
x1's lower freqs should be completely empty.
scipy's problem is sampling in time domain - described further here. Moreover, it centers wavelets in time as appropriate for FFT convolution after ifftshift, but this isn't how np.convolve(, mode='same') computes for even-length kernels, so the wavelet is never even centered at input's index.
== ZERO-MEAN =======================================================================
Found non-zero mean filter(s)!:
psi_fs[97][0] == 2.62e+00+0.00e+00j
psi_fs[98][0] == 1.30e+00+0.00e+00j
psi_fs[99][0] == -8.23e-01+0.00e+00j ...
== PHASE ===========================================================================
Found filters with non-zero phase, with following absolute mean imaginary values:
psi_fs[0]: 2.9e-01
psi_fs[1]: 4.2e-01
psi_fs[2]: 3.3e-01 ...
== TEMPORAL PEAK ===================================================================
Found filters with temporal peak not at t=0!, with following peak locations:
psi_fs[0]: 1
psi_fs[1]: 1
psi_fs[2]: 1 ...
Note: showing full frequency axis, unlike in most plots above (where right half is appropriately zero)
The reader is left with this one as an exercise.
The problem is, not only do we sample in time, but we don't bother to keep doing so, instead resampling from the same sampled sequence. Except for very specific choices of scales, this results in nonuniform sampling and other issues, and can reproduce every problem described in this post (except redundancy). np.diff drastically amplifies existing flaws, which were already great. Note pywt doesn't provide L1 norm as an option, only L2; it's done here manually for clearer plotting. scales was also fixed to be logscaled, while docs examples are linear (which do reproduce redundancy, minus any benefits).
We now inspect individual wavelets (which can occur in a filterbank).
A common price to pay for zero-mean Morlet that's not of excessive scale. Two peaks, or a non-centered peak, adversely affect the mapping to center frequency.
== TEMPORAL PEAK ===================================================================
Found filters with temporal peak not at t=0!, with following peak locations:
psi_fs[]: 178
Found filters with multiple temporal peaks (or incomplete/non-smooth decay)!
(more precisely, >1 inflection points) with following number of inflection points:
psi_fs[]: 3
A crucial property of wavelets is smooth decay in both domains. Lack thereof can be detected with inflection points, or second-order derivatives (or second finite differences). If our wavelets are unimodal, we expect at most one inflection point, where the derivative (instantaneous slope) goes from positive to negative when traversing a wavelet left to right:
== TEMPORAL PEAK ===================================================================
Found filters with multiple temporal peaks (or incomplete/non-smooth decay)!
(more precisely, >1 inflection points) with following number of inflection points:
psi_fs[]: 101
It's also possible to end up with something like
and if our wavelets are supposed to be unimodal (one mode/peak), it's bad news. "Frequency decay" also includes never sufficiently decaying in the first place, but that's easy to imagine so we don't show it.
== DECAY (frequency) ===============================================================
Found filter(s) that decay then rise again in frequency:
psi_fs[]
Unlike with STFT, subsampling is generally trouble with CWT, and should rarely if ever be done without modulus after the transform. No modulus means analytic becomes anti-analytic, and high frequency becomes low frequency:
Aliasing generally will manifest as misordered peaks, so if we intend to order them higher to lower, then it won't be such any longer. It can also break exponential peak spacing.
Even with modulus, there will be aliasing (but as the post shows, we can measure and account for it).
== ALIASING ========================================================================
Found Fourier peaks that are spaced neither exponentially nor linearly,
suggesting possible aliasing.
psi_fs[n], n=[64]
Morlet is pseudo-analytic. If the goal is extracting amplitudes and instantaneous frequencies, strict analyticity is critical. Example for amplitude, CWT of a high frequency pure sine:
It can be achieved by simply zeroing the negative frequencies, where we pay the price in temporal localization, but still worth it. Better yet, use Generalized Morse Wavelets, which suffers less from temporal penalties per having more vanishing moments:
Note, in line with the Hilbert transform, the Nyquist bin should be halved - else, among other things, we get bad time decay: time resolution L2 integral:
zooming on tail:
== ANALYTICITY =====================================================================
Nyquist bin isn't halved for strictly analytic wavelet; yields improper analyticity
with bad time decay.
The DC bin, or input's mean, is at infinity for CWT, so it must be processed separately (and not concatenated with CWT in a 2D representation, except for 1D purposes e.g. 1D convs). If the filterbank is properly designed, the lowpass is redundant in capturing low frequencies, and its only role becomes capturing the DC, also some phase information (as we don't usually pass it through modulus even if we do |CWT|).
For scattering, it's acceptable to omit the lowpass entirely from LP sum. We note, if log2_T < J, then lowpass's energy will overflow the filterbank, and accounting for this overflow will attenuate low frequencies needlessly. For these reasons, validate_filterbank reports LP sum with phi and without phi.
With plenty of examples on what not to do, a "good" filterbank is shown below:
It's achieved with careful consideration of minimum and maximum center frequencies and bandwidths, and a highly heuristic energy normalization algorithm.
What's not been mentioned so far is, all shown filterbanks (except scipy/pywt's) have STFT-like tiling of lowest frequencies: fixed bandwidth, linear spaced frequencies. This is desired, else we're onslaught with redundancy problems. Minus the fact that we enforce zero-mean, this tiling is STFT.
Can we build a filterbank whose LP sum is a flat line - a tight frame?
Very close, but trouble near Nyquist and DC: whether it can all be flat, I'm unsure.
Note, this is "perfection" mainly in a "showy" sense; the practical use for 256 high-redundancy wavelets per octave isn't great.
All of the above shows us why our filterbank can't be just made up of anything - what really makes wavelets wavelets, and if generated in right relation to each other, what makes a proper wavelet filterbank.
All described, except analyticity, also applies to real-valued wavelets. While validate_filterbank only supports analytic, one can simply zero the negatives, which are mirrored duplicates of positives, and still obtain an accurate report - but to be sure, one should verify whether said report relies on the time-domain wavelet being complex.
The complete analytic filterbank includes the anti-analytic pair, but that's only needed for complex inputs. In that case, we pass in each part separately, as validate_filterbank already treats the analytic-only input as having an anti-analytic pair - but with the assumption that this pair is simply a mirrored duplicate in negative frequencies, which isn't universally true (so the function doesn't support this case).
Example from JTFS, with DC centered:
PyWavelets's compute logic isn't directly that of convolution, but it's indeed equivalent to convolution. If a CWT implementation doesn't directly convolve wavelets with input, then we must modify our wavelets in frequency domain to exactly represent the additional operations. Regardless the implementation, CWT is always convolution - or it's not CWT in the first place.
If code's hidden, simply pass in unit impulse - x=zeros(2048); x[1024]=1 - and the (non-modulus) output will be the wavelets themselves, each row being an individual wavelet in time domain (see visual just above "Sufficient padding").
If the implementation pads, and there are large scale wavelets whose support is greater than the input's - then those wavelets are trimmed, but we can still recover up to double the input's support by passing x[0]=1 and x[-1]=1, separately, which give us the wavelet's two halves separately. In fact we can repeat this trick indefinitely if we control scales or they scale with input length. Also ensure if there's padding, it's zero-padding.
To be released soon.