2
$\begingroup$

MATLAB's icwt docs state inversion to be done by a single integral:

$$ f(t) = 2 \Re e\left\{ \frac{1}{C_{\psi, \delta}} \int_0^\infty \left< f(t), \psi(t) \right> \frac{da}{a} \tag{1} \right\} $$

Where does this expression come from, and how is integration over $b$ (translation) eliminated (as in Wiki's formula)? Does it hold if computing discretely over log-distributed $a$ (e.g. $(2^{1/32})^{[1, 2, ...]}$)?

$\endgroup$

2 Answers 2

1
$\begingroup$

The formula is premised on the wavelet being analytic, or being nonzero only over non-negative frequencies: ${\hat\psi} (\omega < 0) = 0$. (Note all wavelets also have ${\hat \psi (0)}=0$ per the admissibility criterion). From Daubechies et al,

$$ \begin{align} \int_0^\infty W_f(a, b) a^{-1} da & = \frac{1}{2\pi} \int_{-\infty}^{\infty} \int_{0}^{\infty}{\hat f}(\xi) \overline{ {\hat \psi}(a\xi)} a^{-1} da d\xi \tag{1a} \\ &= \frac{1}{2\pi} \int_0^\infty \int_0^\infty {\hat f}(\xi) \overline{ {\hat \psi} (a\xi)} e^{jb\xi} a^{-1} da d\xi \tag{1b} \\ &= \int_0^\infty \overline{ {\hat \psi}(z)} \frac{dz}{z} \cdot \frac{1}{2\pi} \int_0^\infty {\hat f}(\xi) e^{jb\xi} d\xi \tag{1c} \\ \end{align} $$

where a rather confusing change of variables was replaced with a clearer one, $z=a\xi,\ dz=\xi da$. Note that this c.o.v. is only valid if integrals share limits, which is enabled by $\psi$'s analyticity; otherwise, $z$ would integrate $-\infty $ to $\infty$, but $z$ includes $a$, which only spans $0$ to $\infty$. An identity was also utilized in splitting integrals.

Assuming $f$ is real, we have ${\hat f(\xi)} = \overline{\hat f (-\xi)}$, thus from $0$ to $\infty$,

$$ \frac{1}{4\pi} \Re e \left\{ \int_0^\infty {\hat f(\xi)} e^{jb\xi} d\xi \right\} = f(b), $$

and setting $C_\psi = \int_0^\infty \overline{ {\hat \psi (z)}} z^{-1} dz$ thusly derives

$$ \begin{align} f(b) &= \frac{1}{2} \Re e \left\{ \frac{1}{C_\psi} \int_0^\infty W_f(a, b) a^{-1} da \right\} \tag{2a} \\ &= \frac{1}{2} \Re e \left\{ \frac{1}{C_\psi} \int_0^\infty f \star \overline{\psi_{a,b}(t)} \frac{da}{a} \right\} \tag{2b} \\ &= \boxed{ \frac{1}{2} \Re e \left\{ \frac{1}{C_\psi} \int_0^\infty \left< f(t), \psi_{a,b}(t) \right> \frac{da}{a} \right\} } \tag{2c} \\ \end{align} $$

Since $W_f$ sweeps $b$ over all $t$, above equivalently holds for $f(t)$.


The paper you linked uses $a^{-3/2}$ -- Indeed, as it assumes an L2 norm on the wavelet, whereas MATLAB uses L1; further info can be found here.

What about the $\delta$? -- The $\delta$ is from MATLAB writing the expression as a special case of the more general inverse, as it shows above near it, with $\left<f, g\right>$, which is derived here.

and the $2 \Re e$? -- Good question - perhaps MATLAB's $C_{\psi, \delta}$ is a bit different, but per above definition we clearly divide by two, not multiply. If in doubt, try both.


Interpretation: it helps to observe that $b$ in all derivations can be replaced with a constant - so the relation holds pointwise; if we integrate along $a$, i.e. a column of the 2D time-frequency representation, we recover the function at that point, $f(b)$.

Notice in the double-integral iCWT we use a dual wavelet to "undo" the inner product with the CWT wavelet, analogous to $e^{+j\omega t}$ in $\mathcal{F}^{-1}$; here, we don't use any inverting kernel - the transform undoes itself when integrated over all scales. This implies some symmetry property on the wavelet $\psi$ in encoding information across scales - though I'm unsure how this exactly works out; the most important is analyticity which permits this in the first place, so one may search there for an answer.


Discretized case with exponentially distributed $a$

From same page of linked paper, the discretized linear case reads

$$ f(b) \approx \frac{1}{2} \Re e \left\{ C_\psi^{-1} \sum_k W_f (a_k, b) \frac{(\Delta a)_k}{a_k} \right\}, \tag{3} $$

with $(\Delta a)_k$ taking role of $da$. For the linear case, $(\Delta a)_k = \Delta a$, but not for the log case, as the increment size increases for greater $k$. That's the only distinction between log and linear.

In practice, however, it's better to directly use the discretized log-scaled transform - which derives as follows:

$$ \begin{align} & z = e^a, \ dz = a da \\ & \frac{da}{a^2} = \frac{d \ln (z)}{z} \tag{4} \end{align} $$

Suppose $z(a) = 2^{a/32}$. Note that the log-plot of an exponential scale is linear, so the difference $\Delta(\ln (z))$ or $d(\ln (z))$ is constant, and equal to $\ln (z(1))da$. To better see this, let's derive directly for this example, using $n_v$ for "number of voices":

$$ z(a) = 2^{a/n_v} \rightarrow dz(a) = z \frac{\ln(2)}{n_v} da = \ln(a_0)da \tag{5} $$

When discretizing, $da=1$, and we end up with

$$ \boxed{f[] \approx \frac{\ln (a_0)}{2 C_\psi} \Re e \left\{\sum_{k=1}^{K} \frac{1}{a_0^k} W_f (a_0^k, .) \right\} } \tag{6} $$

same as Mallat Eq 4.67 (note his derivation is for double-integral, but same scaling argument applies here).

$\endgroup$
0
$\begingroup$

TL;DR

Any strictly analytic or real-valued (in time) filterbank with zero phase is subject to one-integral inversion. It's as simple as the filters, in frequency domain, summing to unity, i.e. the coefficients sum to the original signal via convolution theorem. Full list of criteria:

  1. (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. 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. Full frequency tiling (i.e. non-zero); for analytic this means all positives are non-zero.

  4. Unit stride (hop_size=1)

This is applicable to any time-frequency representation (e.g. STFT), or whatever colvolves with a filterbank.


Inversion proof/derivation

The wavelet-specific derivations are useful for describing design constraints while achieving bandpassing, but are unnecessarily complicated. The definitive criterion for one-integral inversion for any filterbank is found in the convolution theorem:

$$ x \star h = \hat x \cdot \hat h $$

equivalently

$$ x \star \psi_0 + x \star \psi_1 + ... \Leftrightarrow \hat x \cdot \hat \psi_0 + \hat x \cdot \hat \psi_1 + ... = \hat x \cdot (\hat \psi_0 + \hat \psi_1 + ...) $$

Therefore, to have

$$ x \star \psi_0 + x \star \psi_1 + ... = x $$

we require

$$ \hat x \cdot (\hat \psi_0 + \hat \psi_1 + ...) = \hat x $$

i.e.

$$ \hat \psi_0 + \hat \psi_1 + ... = 1 $$

which is a more constrained version of the tight frame criterion:

$$ |\hat \phi(\omega)|^2 + \sum_{s} |\hat\psi_s(\omega)|^2 = 1, \forall \omega $$

(if $||^2$ is constant, so is $||$)) where $||$ are removed, and we see instead that the following must hold:

$$ \boxed {\hat \phi(\omega) + \sum_{s} \hat\psi_s(\omega) = 1, \forall \omega } $$

which also reveals that $\hat \psi$ must be real-valued (or at least their overlapped sum, but unsure that's possible). Reformulating fully in continuous time:

$$ \int_{-\infty}^{\infty} x \star \psi_s ds = x \Rightarrow \int_{-\infty}^{\infty} \hat x \cdot \hat \psi_s ds = \hat x \int_{-\infty}^{\infty} \hat \psi_s ds = \hat x $$ $$ \Rightarrow \boxed{ \int_{-\infty}^{\infty} \hat \psi (s \omega) ds = 1, \forall \omega } $$

(since we have infinite scale, lowpass is no longer needed)

Example / an investigation

Everything below, minus "Code", was written before I the above derivation. It's correct and provides extra intuition, but not necessary.

I've added reasoning here that shows real wavelets are a fair game for one integral reconstruction - along conditions on the entire filterbank.

Still, what exactly allows for CWT(x).sum(axis=0) == x?

First, observe that, for a wavelet symmetric about $t=0$ and $\psi(t=0) \neq 0$, like Morlet, the wavelet can be positive or negative at any other $t$ depending on scale (assuming positive-only scale for now), except at $t=0$, where it always remains some scaling of the original sign. Therefore, the sum of the wavelets at every other point can be zero, while being nonzero at $t=0$.

In convolution, we shift the wavelet such that its $t=0$ is now at some $t=t_0$ - take the product with input, and sum. Now, if all wavelets in our filterbank sum to zero at every point except $t_0$, then the sum of products of wavelets with input will sum to $C\cdot x(t_0)$, where $C = \text{sum}(\psi(t_0))$. Trivial example:

$$ \psi_0 = [-.5, 1, .5], \psi_1 =[.5, 1, -.5], x = [2, 3, 4] $$ $$ \Rightarrow \psi_0 x + \psi_1 x = (\psi_0 + \psi_1)x = [0, (1+1)\cdot 3, 0] = [0, 6, 0] $$ Repeat this for every other $t_0$, and we have the wavelet transform, and its inverse for all $t$.

I said in the other post we require a lowpass for inverse, yet we inverted fine above. That's because the $\psi$'s aren't zero-mean and thus not valid wavelets. No matter how hard we try, we cannot make $\psi$'s that have zero mean that sum to zero at everywhere but $t_0$ (as that would imply the sum of zero mean sequences is itself not zero-mean).

Illustration

I've constructed an approx tight frame ("perfect filterbank") with analytic Morlets. First begin with non-tight:

enter image description here

(n=1 means 1 sample to right of $t=0$) We see that indeed $t=0$ is dominant, and around it oscillations decay to zero. But still, they are non-negligible, and would contribute parts of input around $t_0$, preventing perfect sum to $x(t_0)$. And now, the notorious tight frame:

enter image description here

Much better. Note in particular, the (real) oscillations are all negative. This is crucial, because to end up with a perfect filterbank sum of $[0, ..., 0, 1, 0, ..., 0]$, we add the lowpass filter, which is strictly positive.

Above still isn't an exact (within float precision) tight frame because:

  1. The filterbank isn't purely per CWT, the lowest frequencies are tiled with linear center frequency spacing and constant bandwidth (i.e. like STFT); this manifests as the oscillations seen at leftmost part of LP sum
  2. Frequencies near Nyquist aren't fully tiled (see small gap in LP sum to the left of the black line, which is at Nyquist)
  3. Energy normalization used was one-sided - that is, by assming real inputs, we double LP sum to 2 for one side to conserve energy. But this neglects the tail leaks from anti-analytic wavelets that would contribute to lowest frequencies of LP sum. The effect is very small for large Q, but still worth noting.

Why sum to (<) zero at $t \neq 0$?

It follows from the concept of wavelet self-similarity, and scale:

enter image description here

That is, changing the scale is equivalent to moving along the wavelet from some reference $t = a \neq 0$ point of view ($t=0$ never changes, asm. L1 norm). If that's the case, and $\psi$ is zero-mean (zero-sum), then for $\psi$ to sum to zero, if it's positive at $t=0$, then it must sum to negative for $t \neq 0$. Thus, since traversing scales $[0, \infty)$ is, equivalently, for each $t \neq 0$, same as traversing $\psi$ for $t=(0, \infty)$, then at each such $t$ the sum of wavelets across scales will be negative.

Wrapping up

CWT is convolution with wavelets. Convolution is sum(input * kernel) at various shifts. Sum of such convolutions is equivalently the sum of the sums of products; if the filterbank sums to $0$ everywhere but at $t_0$, then the convolutions at $t_0$ will sum to $x(t_0)$.

What of the imaginary part? It's simply dropped (for real inputs), which naturally qualifies real wavelets for one-integral inversion (e.g. real part of complex Morlet).

What of wavelets what are zero at $t=0$? Unsure, but doesn't seem like a dealbreaker due to lowpass, and they can be analytic, thus satisfying the inverse derivation, but perhaps a different explanation is due.

Code

import numpy as np
from numpy.fft import ifft, ifftshift
from <redacted>.numpy import Scattering1D
from <redacted>.visuals import plot, plotscat, filterbank_scattering
# <redacted>: will update eventually

#%%###########################################################################
def viz(ts):
    filterbank_scattering(ts, lp_sum=1)

    psi_fs = [p[0] for p in ts.psi1_f]
    psis = np.array([ifftshift(ifft(p)) for p in psi_fs])
    phi = ifftshift(ifft(ts.phi_f[0])).real

    psum = psis.sum(axis=0) + phi
    N = psis.shape[-1]
    plot([], hlines=(0, {'color': 'tab:red'}))
    plotscat(psum[N//2-10:N//2+11], complex=1, show=1,
             title="wavelet sum, zoomed")

    slc = psis[:, N//2 + 1].real
    plot(slc, title="all wavelets at n=1 | sum=%.1e (real)" % slc.sum(),
         xlabel="wavelet index", show=1)

#%%###########################################################################
kw = dict(J=11, shape=2048, max_pad_factor=4)
ts = Scattering1D(Q=1, **kw)
viz(ts)
ts = Scattering1D(Q=256, r_psi=.85, **kw)
viz(ts)
$\endgroup$

Your Answer

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

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