Calculating signal power from Continuous Wavelet Transform in MATLAB

Question

I would like to ask a question about the calculation of the signal power using CWT in Matlab. Assume a signal of length N points with sampling frequency $f_{s}$. Using conventional approach, the power will be calculated as follows:

Power = sum((signal).^2)/N

Once the FFT is being performed, the power can be calculated in the following way:

Power = sum(spectrum.^2)

If we use the DWT transform:

[cA,cD] = dwt(signal,'sym4');

the power can be calculated in the following way:

Power_DWT = (sum(abs(cA).^2,'all')+sum(abs(cD).^2,'all'))/(length(cA)+length(cD))

It has to be noted, that all these approaches are giving the same results. The question is how to proceed in case of CWT transformation, since my intention is to calculate power of the signal for a certain frequency range using the CWT coefficients.

Answer 1

For the DWT, the energy is preserved only when the discrete wavelet is orthogonal, since orthogonality allows $L^2$-norm isometry. When it is not, like for biorthogonal wavelets, the relationship between the energy of the signal $x$ and that of wavelet coefficients $w_x$ can be bounded:

$$ A\|x\|^2 \le \|w_x\|^2 \le B\|x\|^2$$

An orthogonal basis is a special case where $A=B=1$. In general for redundant continuous wavelet frame (basis generalization), $A\le B$. When the frame is tight, $A = B$. When your wavelet is well-chosen, and densely sampled, you can have $A\approx B$, and then you can expect pretty good energy preservation.

A read could be: A short introduction to frames, Gabor systems, and wavelet systems Christensen, Ole

Answer 2

CWT power is tricky, and is covered in detail in this post. I provide code on Github to compute energy, and power from it. TL;DR, if you're being asked to compute energy/power, then use the "energy of transform" result. Concerning conversion to physical units, see this answer.

I'll note that by default, MATLAB's CWT undertiles frequencies near Nyquist, so energy conservation won't hold in that range and energy will be under-reported. This can be accounted for by extending "Normalization" (first URL), but it's not done in the code (and that's a good thing for most purposes).

$f_s = 400\ \text{Hz}$, $T = 5\ \text{sec}$. Code output:

Between 50 and 150 Hz, DISCRETE:
994.404   -- energy     (transform)
788.931   -- energy     (signal)
0.497202  -- mean power (transform)
0.394465  -- mean power (signal)

Between 50 and 150 Hz, PHYSICAL (via Riemann integration):
2.48601   Joules -- energy     (transform)
1.97233   Joules -- energy     (signal)
0.497202  Watts  -- mean power (transform)
0.394465  Watts  -- mean power (signal)

Original signal:
1995.56          -- energy     (discrete)
0.997779         -- mean power (discrete)
4.9889    Joules -- energy     (physical)
0.997779  Watts  -- mean power (physical)