2
$\begingroup$

Mallat gives analytic wavelet time & frequency widths/uncertainties as

$$ \begin{align} \sigma_{ts}^2 &= \int_{-\infty}^{\infty} (t - u)^2 |\psi_{u, s}(t)|^2 dt = s^2 \sigma_t^2 \tag{4.51} \\ \sigma_{\omega s}^2 &= \frac{1}{2\pi}\int_{-\infty}^{\infty} \left(\omega - \frac{\eta}{s} \right)^2 |\hat\psi_{u, s}(\omega)|^2 d\omega = \sigma_\omega^2 / s^2 \tag{4.54} \end{align} $$

where

$$ \begin{align} \psi_{u,s}(t) = \frac{1}{\sqrt{s}} \psi^* \left( \frac{t-u}{s}\right),\ \ \hat\psi_{u,s}(\omega) &= \sqrt{s} \hat\psi(s\omega) e^{-i\omega u} \tag{4.53} \\ \eta = \frac{1}{2\pi} \int_0^\infty \omega |\hat\psi (\omega)|^2 d\omega \tag{4.52} \\ \end{align} $$ $$ \sigma_t^2 = \int_{-\infty}^{\infty}t^2 |\psi(t)|^2dt,\ \ \sigma_\omega^2 = \frac{1}{2\pi}\int_0^\infty (\omega - \eta)^2 |\hat\psi(\omega)|^2 d\omega $$

But how would one actually code this? Continuous math can't just be ported directly; the following considerations arise:

  1. How to treat $\psi(t)$ for wavelets defined in frequency domain? We ifft, which now mixes the discrete and discretized. Any normalization considerations?
  2. How to define $t$ and $\omega$, particularly in light of (1)? What are their ranges? For $\omega$ one can presume $0$ to $\pi$ - but what of $t$? Even if we go to infinity, one must still define the interval over which the wavelet decays.
  3. This answer discusses center frequency computed for wavelets sampled directly in time-domain; is the approach the same here?
  4. How to tell whether the implementation(s) is correct? One can always code integrals to spew numbers.
$\endgroup$
2
$\begingroup$

The key is units, and understanding wavelet behavior in context of application (in this case CWT). Full implementations for all discussed here is available at squeezepy.

This answer assumes analytic wavelets as per the question ($\hat \psi (\omega \leq 0) = 0$), but all ideas apply to general time-frequency atoms (e.g. windows for STFT) with some code adjustments.


Time Resolution

Firstly, note we're computing a "width", which is inversely-related to "resolution". The width is the standard deviation of the time-domain wavelet. Standard deviation is well-defined for normally-distributed data, and has desirable properties if data is a probability distribution.

By squaring modulus of $\psi(t)$ and normalizing it to integrate to 1, we interpret the result as a probability distribution, and define the time variance:

$$ \sigma_t^2 = \frac{\int_{-\infty}^{\infty} t^2 |\psi(t)|^2 dt}{\int_{-\infty}^{\infty} |\psi(t)|^2 dt}, \tag{1} $$

where $\psi(t)$ is assumed zero-mean. We thus interpret $\sigma_t$ as as the time-span of ~68% of wavelet's energy.

The purpose of the denominator is to enforce validity as probability distribution, i.e. normalizing $\psi$ if we haven't already. This is valuable for coding the integral, as the wavelet at different scales may not integrate to $1$ automatically. It also reminds us of the probabilistic context rooting $\sigma_t$. Mind that $\sigma_t$ is half the said span, as we compute the uni-lateral ('radius') variance.

The question is, how to define $t$? We can't take infinity, and even if we could, what finite width would span the non-decayed wavelet? Answer to latter is rooted in how DFT works, discussed in "Are we right?" below, but there's an extra concern: should t be zero-mean? No; it should be defined such that it aligns with the wavelet as if it were inputs to generate the wavelet.

Most wavelets will peak at $t=0$, for which t=0 should align with wavelet's peak (i.e. same array index), negative to left, positive to right. Such t will be nearly zero-mean, but not exactly, and requires different handling for even- and odd-length t.


Frequency resolution

Much of what's said for $\sigma_t$ also holds here; we define the frequency variance as

$$ \sigma_\omega^2 = \frac{\int_0^{\infty} (\omega - \tilde\omega_\psi)^2 |\Psi(\omega)|^2 d\omega } {\int_0^{\infty} |\Psi(\omega)|^2 d\omega}, \tag{2} $$

where $\tilde\omega_\psi$ is the energy center frequency (see "Center Frequency"). We thus interpret it as frequency-span of ~68% of wavelet's spectral energy. In CWT filter bank terms, a small $\sigma_\omega$ for a given scale means that most of input signal's frequencies will be band-passed about $\tilde\omega_\psi$, and the rest attenuated/rejected.


Center Frequency

The cited answer's approach won't work here - at least not directly. Either way it's suboptimal; since discretization favors us in frequency domain, integrating is best. Ref [1] defines three meaningful frequencies one can associate with a wavelet's scales, and corresponding center frequencies; two of those for center frequency are the energy center frequency and the peak center frequency, latter defined as frequency at which $\Psi(\omega)$ peaks, and former as:

$$ \tilde\omega_\psi = \frac{\int_0^{\infty}\omega |\Psi(\omega)|^2 d\omega} {\int_0^{\infty}|\Psi(\omega)|^2d\omega} \tag{3} $$

Here we don't assume $\Psi(\omega)$ is zero-mean (in fact it can't be). If $\Psi(\omega)$ is even-symmetric, then energy center frequency == peak center frequency. For finite and coded $\Psi(\omega)$, this distinction becomes significant and crucial to any wavelet regardless of shape ("Edge Cases").


Are we right?

Three sanity checks:

  1. Can results be interpreted meaningfully?
  2. Is $\sigma_t \sigma_\omega \approx 0.5$?
  3. Does $\sigma_t \sigma_\omega$ scale correctly with varying wavelet parameters ($\mu$, scale, etc)?

#2 is our main check; if the product is much greater or lesser than 0.5, we're off.

To cut it short, I had to 'cheat' to find the proper range for $t$. I started assuming $\omega_\psi$ to be correct since it made sense and correctly gave center frequency as the exact location of the mode=peak. For $\sigma_t$ I ranged $t\in [-.5, .5)$, as that's interpreted as fraction of frame's length. $\sigma_t \sigma_\omega$, however, did not multiply to $0.5$.

Working backwards, $\sigma_t = 0.5 / \sigma_\omega$, I found the correct $t$ to range $t \in [-N/2, N/2)$! And then it all clicked: units.

The units of $\sigma_\omega$ are samples per cycle. Then, to cancel in Heisenberg's relation, $\sigma_t$ must have the reciprocal units: cycles per sample. This is unified by below relation for DFT:

$$ \begin{align} f_p \left[ \frac{\text{cycles}}{\text{second}} \right] & = \left( f_{\text{DFT}} \left[ \frac{\text{cycles}}{\text{samples}} \right] \right) \cdot \left( f_s \left[ \frac{\text{samples}}{\text{second}} \right] \right) \end{align} \tag{4} $$

which can also be used to obtain resolutions in physical units (Hertz, seconds, etc).

Above is also a bit simplified; the full units are:

$$ \begin{align} \left[\left[ \sigma_\omega \right]\right] &= \frac{\text{cycles}\cdot \text{radians}}{\text{samples}} \tag{5a} \\ \left[\left[ \sigma_t \right]\right] &= \frac{\text{samples}}{\text{cycles}\cdot \text{radians}} \tag{5b} \\ \end{align} $$

Radians in $(5b)$ are mandated by those in $(5a)$ to cancel in product to yield unitless Heisenberg area. They're also sensibly interpreted once considering the nondimensional variants.


Edge cases

Discretizing problems don't end at finite precision and scaling conversions; consider what happens with scale pushed to extremes.

High scale

Freq-domain not doing too poorly - mainly finite precision imperfections. But time-domain is trouble: it doesn't fit the frame! Two problems here:

  1. How do we handle $\psi (t)$? Should we extend the frame by adding samples until decay?
  2. How does our choice in 1 affect CWT?

#2 is key to #1. Recall in CWT with freq-domain wavelets, we do a form of discretized convolution theorem, getting sort of the exact continuous-time result, then discretizing. An imperfect time-domain wavelet is characterized by not fitting in the frame - but what if the frequency-domain wavelet fits just fine?

I've yet to meet a time-domain wavelet that doesn't decay sufficiently if the frequency-domain wavelet does; problems arise when the frequency-domain wavelet doesn't. At the high scale extreme, the freq-dom wavelet is defined by a single sample, which is a time-domain pure sine with $\sigma_t=\infty$. Proper decay demands even symmetry of the freq-domain wavelet, which is bare minimally attainable by two adjacent bins at equal values. How exactly various scenarios are handled can be found in squeezepy's tests/props_test.py.

Low scale

Now the time-domain wavelet fares fine, while freq-domain is clearly troubled: half (or more) of it is entirely missing! We have the "trimmed bell" problem, contrary to "insufficiently sampled bell" for high scale case. The center frequency distinction becomes critical here; peak will always equal $\pi$, while energy will lie to its left. Investigating the behavior of the $\sigma_\omega$ integral, I found the energy variant to yield stabler and at times more sensible results, but depending on application, peak can also be used (with due re-interpreting of $\sigma_\omega$).

Nondimensional measures

Doing everything right above, one can take the Morlet wavelet and vary its mu only to notice that time and frequency resolutions don't change! What's the deal?

In short, we need a nondimensional measure; this is obtained by normalizing above as follows (thus completing formulae in ref[2]): $\sigma_t \cdot= \omega_\psi,\ \sigma_\omega /= \omega_\psi$ in relevant equations above, where $\omega_\psi$ is peak center frequency. As to how this is interpreted or why it works - it's already a lot for one question, so briefly: try making sense of $\sigma_t$ relative to center frequency's time-domain period, and $\sigma_\omega$ relative to spread around center frequency in log space (or simply ratio-wise).

Non-dimensional measures change with wavelet parameters, but not with scales, so each has its own use and interpretation.

Where's the $2\pi$?

Mallat's equations divide by $2\pi$, we don't, but there's no contradiction: this is per Parseval-Plancherel's theorem relating time-domain and frequency-domain energies. Since we explicitly divide by said energies in $(1)$ and $(2)$, the $2\pi$ is automatically accounted for.

Further reading

Not everything was explained here, nor can be in a single post. ssqueezepy contains a complete implementation, with accompanying tests that provide further explanation - particularly on handling scale extrema, which is a whole topic of its own.

References

  1. Higher-Order Properties of Analytic Wavelets - J. M. Lilly, S. C. Olhede
  2. Wavelet Tour, Ch1-4 - S. Mallat
$\endgroup$
1
$\begingroup$

We can visualize said resolutions with ssqueezepy. Let us compare the Morlet wavelet for mu=5 vs mu=20 parameters. For reference, I also recommend this Desmos. All plots below (unless otherwise shown) via:

wavelet = ssqueezepy.wavelets.Wavelet(('morlet', {'mu': 5}))  
wavelet.viz()  # repeat w/ {'mu': 20}

Waveforms across scales

enter image description here

  • Heatmaps show top-down view of time- & frequency-domain wavelets across different scales
  • mu=20 is lot more "wiggly" than mu=5 in time-domain, as expected for higher center frequency
  • mu=20 is much "thinner" than mu=5 in freq-domain (and wider in time-domain), showing a difference in resolutions (see other answer)
  • The two span different range of scales; this is per ssqueezepy 'intelligently' picking minimum and maximum, similar to (but not same as) MATLAB's cwtfreqbounds

Time & frequency joint resolution

  • Textbook-like plots, showing joint time & frequency widths (widths drawn to scale, but not amplitudes)
  • Dimensional widths are same for scale=100, but not non-dimensional
  • Box areas not accurately represented; bounds precision is within 1 pixel

Time & frequency widths vs scales

enter image description here

  • Plots compress a lot of information; see respective methods in ssqueezepy.wavelets to see how values are computed
  • Notice "well-behaved" and "poorly-behaved" ranges of scales; mu=20 is harder to tame
  • Near extrema (low/high scales), distinct discretization problems begin to show (see other answer)
  • Heisenberg product less than 0.5 is misleading here; the design decision was such that std_t and std_w by default more accurately reflect their own widths, but not joint. Defaults can be changed to ensure area >=0.5 (not recommended)
  • Further comments on edge cases in tests/props_test.py

Joint resolution across scales

  • wavelet.viz('anim:time-frequency', N=2048) (omitted mu=5; see here)
  • Adaptive resolution nature of the wavelet is clearly seen
  • Scale extrema problems are partially observed; can you guess what happens if we push them even further (higher & lower)?

How scales determine waveforms

enter image description here

  • wavelet.viz('waveforms', scale=1.5, N=1024) (repeat for scale=400)
  • Changing scales amounts to changing the spacing of the sampling grid in frequency domain (red dots)
  • Partial insight into extreme behavior is obtained by sweeping scale near extrema and observing waveform behavior
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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