# One integral inverse CWT

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, ...]}$$)?

Revisiting this question for a more definitive/intuitive answer. 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:

(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:

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:

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.

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)

### Code

import numpy as np
from numpy.fft import ifft, ifftshift
from kymatio.numpy import Scattering1D
from kymatio.visuals import plot, plotscat, filterbank_scattering

#%%###########################################################################
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)

#%%###########################################################################
ts = Scattering1D(Q=1, **kw)
viz(ts)
ts = Scattering1D(Q=256, r_psi=.85, **kw)
viz(ts)


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$$, 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$$

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