Signal Processing Stack Exchange is a question and answer site for practitioners of the art and science of signal, image and video processing. It only takes a minute to sign up.
Sign up to join this community
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).
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:
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.
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)
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)
#%%###########################################################################
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)
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).
