1
$\begingroup$

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?

$\endgroup$
0

1 Answer 1

8
$\begingroup$

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:

  1. Analyticity: strict vs pseudo
  2. Aliasing: insufficient sampling rate or duration
  3. Zero-mean: required for wavelets and why so
  4. Zero-phase: purely real in Fourier, and why it's desired
  5. Frequency coverage: whether we capture every frequency, whether we do so excessively or insufficiently, and
  6. Frequency-bandwidth tiling: whether frequencies are tiled according to CQT
  7. Redundancy: whether filters overlap excessively or duplicate
  8. Decay: whether filters decay sufficiently and smoothly in time and in frequency; whether filters are pure sines
  9. Temporal peaks: whether the peak is at $t=0$, and there's only one peak (if applicable)

Baseline: Generic filterbank

  • Energy conservation: fail. An analytic filterbank tiles only the positive frequency axis, which is half of input's energy.
  • Analyticity: fail. The two highest frequency wavelets leak past Nyquist.
== 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:

  • Energy conservation: same, except now the lowest frequencies have a drastic spike, despite the filterbank plot not looking any different.
== 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: multiple wavelets have the same peak frequency! Moreover, they overlap excessively, making them essentially duplicates.
  • Temporal decay: these "wavelets" are pure sines! They aren't even wavelets. A wavelet requires, super minimum, two samples in frequency - the wavelet scale is too great, creating a "sampling artifact".
  • Aliasing: consistent with previous point, these wavelets are actually aliased in frequency domain.
== 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:

Sufficient padding

Much better. Time-domain:

Hence, padding isn't only about boundary effects, but wavelet behavior.

Energy norm: simple $\sqrt{2}$

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]

L2 norm

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.

Incomplete tiling

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 ...]

Incorrect frequency-bandwidth tiling

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 ...

Redundancy

Earlier filterbanks had a sort of "critical" redundancy. Below is high redundancy:

Whether this is beneficial is application-dependent:

  • Classification of individual 1D CWT vectors: bad news, as our features are highly correlated.
  • Synchrosqueezing: good news. Synchrosqueezing is very redundancy-hungry. Note the completely flat LP sum as far as eye can tell: with a constant rescaling, and disregarding the region near Nyquist, we have a tight frame (ideal), which perfects one-integral inversion (which synchrosqueezing uses in reassignment).
  • Classification (other): greater redundancy comes with greater time localization and time-warp stability.
== 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 ...

Non-zero phase, non-zero mean

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: matters wherever complex outputs do. This includes one-integral inversion, and instantaneous frequency and amplitude localization (e.g. synchrosqueezing). Depending on the kind of non-zero phase, it could also completely misalign the output in time.
  • Non-zero mean: mixes modulations with offsets (analogously, "rates and offsets", or "weights and biases"). A distinct advantage of CWT over STFT is, not doing this. "Frequency" means frequency, period.

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 ...

Bonus: PyWavelets

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).

Temporal peak

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

Decay

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[]

Aliasing

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]

Analyticity

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.

Lowpass filter

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.

Proper normalization

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.

Perfect filterbank?

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.

Noise = filterbank?

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.

Real-valued wavelets?

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.

Anti-analytic pair?

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:

No direct convolution?

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.

Closed source implementation?

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.

Code

To be released soon.

$\endgroup$
1
  • 1
    $\begingroup$ Well, I'll be coming back to this for the foreseeable future, so thank you. $\endgroup$
    – Jdip
    Jan 6, 2023 at 13:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.