Power/Energy from Continuous Wavelet Transform

Question

How can power or energy be computed from Continuous Wavelet Transform? Is it just $\sum |\text{CWT}(x)|^2$, or are there other considerations, particularly if interested in a subset of frequencies? Do the results interpret differently from what's computed from DFT?

Answer 1

CWT power is tricky. We must distinguish between energy of transform and energy of signal, and the two can be very different. As energy is the more fundamental and less conceptually intricate quantity, I develop this answer for energy, then relate to power. Answer applies to CWT and STFT, with differences described near bottom.

CWT here is defined "semi-discrete" (continuous domain, finite number of wavelets) for notation, but all generalizes to fully continuous and fully discrete where applicable.

The two compared, briefly

Each function above is used to compute energy. Again, very different.

TL;DR?

Read "Which to use?" and TL;DR in "Conditions on $\psi$", and adjust "Code example" to your use, but there's no shortcuts to avoiding mistakes.

Energy of transform

There is a "Parseval's equivalent" for each; for transform, we have

$$ E(\text{CWT}_\psi(x)) = \sum_i \|x \psi_i\|^2, \tag{1} $$

which is the sum of energies of coefficients. Note this is simply Parseval-Plancherel's theorem over the entire transform; for any one coefficient from wavelet e.g. $\psi_4$, the convolution in frequency domain is $\hat x(\omega) \hat \psi_4(\omega)$, its energy is $\|\hat x(\omega) \hat \psi_4(\omega)\|^2$, which via theorem is $\|x(t) \psi_4(t)\|^2$, and since $|AB| = |A||B|$, that's $\int |x(t)|^2 |\psi_4(t)|^2$ (can't factor the aggregation too, i.e. $\|x\|^2 \|\psi\|^2$).

As the expression is factorable,

$$ \int_{-\infty}^{\infty} |x|^2 \sum_i |\psi_i|^2, \tag{2} $$

we have the "energy transfer function" for the transform:

$$ E_{\psi, \text{CWT}}(\omega) = \sum_i |\hat\psi_i(\omega)|^2, \tag{3} $$

the sum of energies of wavelets, pointwise, known as the Littlewood-Paley sum. For a subset of frequencies, just have $i$ correspondingly range from $\omega_\min$ to $\omega_\max$.

Energy of signal

CWT is a decomposition: the total is the combination of its parts. Conveniently, the "combination" here is simply a sum of its coefficients, for each timestep - the one-integral inverse, x == sum(CWT(x), axis=0). For a frequency subset of interest, then, "signal energy" is energy of that subset's signal, e.g. $\|x * \psi_4 + x * \psi_5\|^2$ which is $\|\hat x \hat\psi_4 + \hat x \hat\psi_5\|^2$. Repeating above logic, the expression for the entire signal is now

$$ E(x_\text{inv}) = \left\| \sum_i \hat x(\omega) \hat \psi_i(\omega) \right\|^2, \tag{4} $$

which factors as

$$ \int_{-\infty}^{\infty} |x|^2 \left| \sum_i \psi_i \right|^2, \tag{5} $$

so the energy transfer function is

$$ E_{\psi, x_\text{inv}}(\omega) = \left| \sum_i \hat\psi_i(\omega) \right|^2, \tag{6} $$

the energy of sum of wavelets, pointwise. Compared, we have $\sum \|\psi \|^2$ vs $\|\sum \psi \|^2$: not same! (Consider $|2|^2 + |3|^2$ vs $|2 + 3|^2$). For a subset of frequencies, just have $i$ correspondingly range from $\omega_\min$ to $\omega_\max$.

Which to use?

TL;DR If you're asked to find energy/power, then almost certainly ET. Else, depends:

• Energy of transform (ET) is the energy of the coefficients of the transform.
• Energy of signal (ES) is the energy of the signal that's obtained by inverting the transform.
• Each can be applied over a part of the transform, or the whole.

Which to use depends on application. Use ET if the transform modulus (scalogram) is of interest. Use ES if the signal underlying the transform is of interest. Examples:

• Feeding scalogram to neural net: ET
• Synthesizing signal after modifying CWT: ES

There's certainly more I'm forgetting. Important notes, ES and ET have different properties:

1. ET is non-decreasing: including more coefficients never decreases energy, unlike ES
2. ET is separable: energies of different bands can be analyzed independently. ES can't do this per 1, which also means we can't plot "Energy vs Frequency" for ES.
3. ES lacks linearity (w.r.t. input freqs), which just restates 2: $\text{ES(5Hz to 50Hz) + ES(50Hz to 100Hz)}\ \neq \text{ES(5Hz to 100Hz)}$.

Again there's more I'm missing, just compare the respective functions. Lastly, small note on ES - I made it up. Because, it makes sense. Maybe it's already a thing, but case in point, most sources mean ET when they say "energy".

Conditions on $\psi$

TL;DR use filters that are real-valued in frequency, are L1 normalized, are either real-valued in time or are strictly analytic (negative frequencies = 0), use hop_size=1. Also "more wavelets" is better (see "Normalization").

The following one-integral inverse conditions apply (strict analyticity, unit stride) apply:

1. ET & ES: (A) strictly analytic (no - freqs) -- (B) anti-analytic (no + freqs) -- (C) real-valued in time

   • Real-valued $x$: (A) OR (B) OR (C)
   • Complex-valued $x$: (A) AND (B), OR (C)

2. ET & ES: Real-valued in frequency (zero-phase). More precisely, the sum of all filters in freq domain must not have an imag part, but it's unlikely if individual filters have imag part.

3. ES: Unit stride (hop_size=1)

4. ES: filters are L1 normalized. (note this isn't in the linked post)

Ignoring these conditions:

• ET: still yields the energy of the transform, but then the interpretation changes as the transform itself mirepresents the signal.
• ES: no longer yields the energy of the inverted signal, though it may be close.

Lastly, for calculations to match exactly, no unpadding. This is necessary, as the padded filters' and signal's DFT aren't same as unpadded or non-padded. However, we should still pad, and results should be quite close. Subject of another post.

Normalization

As ET and ES functions can differ significantly in shape and norm, it's usually impossible to design filterbanks that satisfy both perfectly. The most important trait to optimize for is flatness in frequency domain - flat in ET is at least approximately flat in ES, and vice versa, and makes ET and ES within scalar multiples of each other. For this, more wavelets = better.

Whether these adjustments are "valid" depends on use.

• If a result (e.g. convolving along frequency of modulus of CWT) expects the CWT to have certain energy, then our adjustments of said energy must also be done on the coefficients.
• If we're simply inspecting, or aiming for CWT that's a sampled version of the continuous result, then the coefficients shouldn't be adjusted, and filters shouldn't be rescaled to satisfy energy measures.
• If we're trying to compare different signals, or compare ET to ES in a scaling-agnostic manner, or obtain certain normalized measures, then it makes sense to adjust both ET and ES.
• If we seek ET & ES that's comparable to the signal (e.g. physical quantities), then it's mandatory.
• For operating on the scalogram, I doubt it's ever valid to adjust L2-normed CWT except for scalar scaling.

A scalar multiple adjustment is done in code examples below so our transform-based results don't contradict the original signal (i.e. exceeding its energy/power).

Energy conservation

ET provides a simple two-number measure of energy conservation; if

$$ A \leq \sum_i |\hat\psi_i(\omega)|^2 \leq B $$

then multiplying by $|\hat x(\omega)|^2$ and applying Parseval-Plancherel's theorem

$$ A \|x\|^2 \leq \| \text{CWT}(x) \|^2 \leq B \|x\|^2 $$

bounds worst-case energy losses/gains by $A, B$.

• $A = B = 1$ is a tight frame. This might imply ET=ES, I'm unsure (it is the case for Morlets).
• $A = 0$ means the transform is lossy (non-invertible).
• $B > 1$ means the transform is expansive. $B < 1$ makes it contractive.
• The DC bin is excluded from compute as it makes the measure useless (all wavelets are zero-mean; finite CWT discards signal offset).

STFT?

STFT is convolution with windowed complex sinusoids. With hop_size=1, all logic here translates smoothly. With hop_size > 1, this answer is already too long, but one should mind that the spectrogram is always aliased (to different extents, sometimes negligibly), so ET won't yield consistent results, and ES is either done differently or not at all.

ET vs ES, when should I really care?

When ET and ES plots disagree in flatness. In provided codes below, plot ET_tfn and ES_tfn (done at top of this post).

Physical/instantaneous energy/power?

See this answer. In general case, |cwt(x)|^2 is already instantaneous power by definition, but if physical units make it instantaneous energy instead - apply answer to differentiate the ET or ES that lacks the temporal aggregation step, along accounting for sampling period and duration, and normalization (this answer).

Note, the complex modulus heavily shifts energy toward lower frequencies, meaning differentiating CWT/STFT modulus rows suffers minimal aliasing in most cases. This yields excellent estimation, at expense of loss in temporal resolution - so the result shouldn't quite be treated as "instantaneous".

This also summarizes the case vs DFT - time-frequency results are localized, in both time and frequency. However, with simple-spread wavelets like Morlet and high time frequency resolution, the difference along frequency may be inconsequential; here, synchrosqueezing can offer significant enhancement.

Code example

Available on Github, in Python and MATLAB. $f_s = 400\ \text{Hz}$, $T = 5\ \text{sec}$. Output (Python):

Between 50 and 150 Hz, DISCRETE:
970.433   -- energy     (transform)
883.466   -- energy     (signal)
0.485217  -- mean power (transform)
0.441733  -- mean power (signal)

Between 50 and 150 Hz, PHYSICAL (via Riemann integration):
2.42608   Joules -- energy     (transform)
2.20867   Joules -- energy     (signal)
0.485217  Watts  -- mean power (transform)
0.441733  Watts  -- mean power (signal)

Original signal:
1913.8           -- energy      (discrete)
0.956901         -- mean power  (discrete)
4.7845    Joules -- energy      (physical)
0.956901  Watts  -- mean power  (physical)

and as a sanity check, include full transform:

Between -1 and 201 Hz, DISCRETE:
1895.69   -- energy     (transform)
1878.75   -- energy     (signal)
0.947844  -- mean power (transform)
0.939374  -- mean power (signal)

Between -1 and 201 Hz, PHYSICAL (via Riemann integration):
4.73922   Joules -- energy     (transform)
4.69687   Joules -- energy     (signal)
0.947844  Watts  -- mean power (transform)
0.939374  Watts  -- mean power (signal)

Code validation

Demonstrating that ET and ES work as described. Available on Github in Python.

Appendum: Shortcomings: unpadding, stride

Not discussed is the need to account for padding and non-unity stride. That's its own topic, but major points on unpadding:

• ET does not conserve energy! With padding, energy is conserved while padded, but then contracted upon unpadding. The effect is the most severe with zero-padding, and hard to quantify with other schemes, but I found 'reflect' to work best for aiding conservation. For worst case, consider the case x = zeros(N); x[0] = 1; half of ET is gone.
• ES conserves energy! The sum of a tight frame, zero-phase filterbank, is an all-pass filter. It's actually the unit impulse. That gives it perfect temporal localization, hence unpadded energy exactly matches the signal's energy.
  • Caveat: for a subset of the transform, the subset must fully capture a "component" (AM/FM, or intended signal) for results to be meaningful. Example, unit impulse: energy should be zero everywhere but at one time instant, but this won't happen unless all frequencies are included.
• Unpadding is aliasing, making exact correction for ET impossible.
• Not padding circumvents the problem, at expense of feature quality.

On stride - with large strides, the exact unpad index may be fractional. This means we either under-unpad, or over-unpad. Either can be accounted for as follows:

• Suppose input length is 64, and stride is 8. Then, len(out) = 8. Suppose input length is 65, and stride is 8. Then, len(out) = 9, yet the exact length is 65/8 = 8.125, which is unachievable. Hence, because of just one sample, we have an up to 9/8 energy difference relative to the first case - or, to be exact, 9/8.125. Make up for it with coef *= sqrt(8.125/9).
• This assumes a uniform energy extension of the missing continuum, i.e. that it's accurately predicted by the mean of the rest. Not flawless but pretty good in unaliased convolution settings.
• Besides inexact unpadding, the basic adjustment is coeffs *= sqrt(hop_size), via Parseval's, assuming energy-unaliased convolutions. "Energy-unaliased" here means the result doesn't fold on itself (subsampling in time $\Leftrightarrow$ folding in Fourier) over non-zero regions (changes shape of Fourier modulus), so aliasing that's just a frequency shift is fine.

Full derivations

Note, again, $|AB| = |A||B|$ and $|AB|^2 = (|AB|)^2 = (|A||B|)^2 = |A|^2|B|^2$. We derive "transfer functions" in sense of

$$ \texttt{OUT} = \int_{-\infty}^{\infty} \texttt{IN}(x) \cdot \texttt{TRANSFER}\ (x)\ dx $$

Energy of transform:

$$ \begin{align} E(\text{CWT}_\psi(x)) &= \sum_i \|x \psi_i\|^2 \\ &= \sum_i \int_{-\infty}^{\infty} |x(t) \psi_i(t)|^2 dt \\ &= \sum_i \int_{-\infty}^{\infty} |x(t)|^2 |\psi_i(t)|^2 dt \\ &= \int_{-\infty}^{\infty} |x(t)|^2 |\psi_0(t)|^2 dt + \int_{-\infty}^{\infty} |x(t)|^2 |\psi_1(t)|^2 dt\ +\ ... \\ &= \int_{-\infty}^{\infty} |\hat x(\omega)|^2 |\hat \psi_0(\omega)|^2 d\omega + \int_{-\infty}^{\infty} |\hat x(\omega)|^2 |\hat \psi_1(\omega)|^2 d\omega\ +\ ... \\ &= \int_{-\infty}^{\infty} |\hat x(\omega)|^2 \left(|\hat \psi_0(\omega)|^2 + |\hat \psi_1(\omega)|^2 +\ ...\right) d\omega \\ &= \int_{-\infty}^{\infty} |\hat x(\omega)|^2 \sum_i |\hat \psi_i(\omega)|^2 d\omega \\ &= \int_{-\infty}^{\infty} |\hat x(\omega)|^2 E_{\psi, \text{CWT}}(\omega) d\omega \\ \end{align} $$

Energy of signal:

$$ \begin{align} E(x_\text{inv}) &= \left\| \sum_i \hat x(\omega) \hat \psi_i(\omega) \right\|^2 \\ &= \int_{-\infty}^{\infty} \left| \sum_i \hat x(\omega) \hat \psi_i(\omega) \right|^2 d\omega \\ &= \int_{-\infty}^{\infty} \left| \hat x(\omega) \sum_i \hat \psi_i(\omega) \right|^2 d\omega \\ &= \int_{-\infty}^{\infty} \left| \hat x(\omega) \right|^2 \left| \sum_i \hat \psi_i(\omega) \right|^2 d\omega \\ &= \int_{-\infty}^{\infty} \left| \hat x(\omega) \right|^2 E_{\psi, x_\text{inv}}(\omega) d\omega \\ \end{align} $$

Notation

• $\| f(t) \| = \sqrt{\int_{-\infty}^{\infty} |f(t)|^2 dt}$, the L2 norm
• $E(f(t)) = \text{Energy}(f(t)) = \|f(t)\|^2$
• $\hat f(\omega) = \mathcal{F}(f(t))$, Fourier transform