Can finite Multi-Resolution Analysis satisfy the Littlewood-Paley Criterion (unity partition)?

Question

Suppose we are working with a multi resolution analysis (MRA) of $L^2(\mathbb{R})$ and let $\phi$ be the corresponding scaling function and $\psi$ the derived wavelets. Using standard notation with $\hat{f}$ denoting a Fourier transform of a function $f\in L^2(\mathbb{R})$ is it always true for all integer scales $J$ that

$$|\hat{\phi}(2^J\omega)|^2+\sum_{j>-J}|\hat{\psi}(2^{-j}\omega)|^2 = 1$$

(notation adapted from this paper, equation 8). Is it always true that this is the case, that is, for any wavelets derived from MRA, or do the wavelets have to be specific?

Answer 1

The Littlewood-Paley criterion is satisfied by an appropriate "generator rule" that relates a given filter to its successor, $\psi_{i + 1} = \text{rule}(\psi_i)$. The neat thing is, this holds for all $\psi$ (or let's say $\in L^2$ / "non-weirdos"), regardless of shape - but there's practical caveats.

I'm unfamiliar with 'official' justifications as I've not 'formally' studied frame theory or LP theory - below are my investigations and (sort of) proofs.

Summary

Recommended to read last unless just looking for TL;DR.

The LP criterion is satisfied for all $\psi$, but not with finite $J$ - there, the unsatisfied $\omega$ interval shrinks with growing $J$. In practice, we must additionally ensure the wavelet tiling is CQT: Q = (center frequency) / bandwidth = const., along other caveats. The LP criterion is also satisfied for STFT, with much fewer problems.

This is because, the LP sum for CWT is equivalently the convolution of $|\hat \psi|^2$ with a constant, in log space, which is a constant in log space, hence a constant in linear space. It's a convolution because of CQT: Q = freq / bandwidth = constant, and ratios are uniform in log space, which undoes the domain transformation and keeps $\hat\psi$'s shape fixed for all frequential translations.

LP: STFT

To see how it works, begin with STFT (specifically the "mod" version here): it's convolutions with windowed complex sinusoids, so all filters have the same shape in frequency domain, but are shifted. Now the magic: the LP sum is simply the sum of all $|\hat\psi|^2$, which is equivalently convolution of $\hat \psi$ with a flat line, which itself is a flat line! Here's a proof in discrete case with $\psi$ that's complex white noise:

import matplotlib.pyplot as plt
import numpy as np
from numpy.fft import fft, ifft

# completely arbitrary kernel
xf = np.random.randn(256) + 1j*np.random.randn(256)
# add up all of its shifts
sm0 = np.sum([np.roll(xf, i) for i in range(len(xf))], axis=0)
# now do it as convolution
sm1 = ifft(fft(xf) * fft(np.ones(len(xf))))

# assert equality & plot
assert np.allclose(sm0, sm1)
plt.plot(sm0.real); plt.plot(sm1.imag)

Of course if it holds for sum of $\hat\psi$, it holds for sum of $|\hat\psi|^2$.

The "generator rule" here is $\hat\psi_{i + 1}[k] = \hat\psi_i[k - 1]$ (favoring discrete notation), i.e. frequency shift without dilation, i.e. without changing bandwidth; this makes the resulting tiling uniform, i.e. the overlapped sum is the same everywhere.

LP: CWT

The argument translates elegantly: if filters sum to a constant in log-frequency, they also sum to a constant in frequency: if $f(t) = 1$, then $f(g(t)) = 1$ for all (non-weirdo) $g$, including $g(t) = \log(t)$. And in log-frequency, different wavelets have the same bandwidth - so it's same as STFT! Or... do they?

The "generator rule" in CWT is $\psi_{a}(t) = \psi(t/a)$, i.e. wavelet at scale $a$ is the mother wavelet dilated by $a$, which implies $\psi_{a_1}(t) = \psi_{a_0}(t \cdot (a_0 / a_1))$, i.e. "wavelet at scale $6$ is wavelet at scale $2$ dilated by $3$". But that's just stuff in time domain; what must be shown is, that the filterbank tiling is uniform in log-frequency.

So, what's uniform in log space? Ratios: $8$ and $4$ are as far from each other as $4$ and $2$. Incrementing in log space means multiplying. So a wavelet at frequency $8$ must be twice as wide as at $4$. Put differently, the ratio of center frequency to bandwidth must be constant: CQT (Constant-Q Transform, Q = xi/sigma).

In time, we indeed have a ratio relation. In frequency, conveniently, we have

$$ x(t/a) \Leftrightarrow \hat x (a \omega) |a| $$

and that's the key: dilation in time is contraction in frequency, for all $x$. It means that uniformity is inherent to the generator rule. In discrete, we have $\hat\psi_{i + 1} = \hat\psi_{i}[k/s]$; in continuous, $s$ can be anything, but here it must be the log equivalent of $k - 1$ in linear case, which depends on the log base, or "Q" in CQT. But as we'll see later, that's not enough either.

For a visual, here's a CQT-tiled Morlet filterbank:

LP: general

In "LP: CWT" we wrote,

$$ f(t) = 1 \Leftrightarrow f(g(t)) = 1, \forall g $$

Put in words, if a function is constant in some arbitrary transformation of its domain, then it's also constant in the original domain. Our "function" of interest is the LP sum, $\sum |\hat\psi(\omega)|^2$.

We saw that in log space, the "generator rule" must effectively undo the log transformation. It's same in the general space: a translation in $g$-space cannot change the shape of $\hat\psi$ in $g$-space. While this is achievable, the effect on time domain isn't as predictable, if at all linear, as we only have two Fourier properties pertaining to translation: $\hat x(\omega - a)$ and $\hat x(\omega / a)$!

LP: Scattering

A layer of scattering (as given in the equation) is CWT up to $J$, complemented by lowpass. Problem is... I don't see it working out.

I've simulated Morlet, alongside Gaussian. Note, "$\psi$ for $j > -J$" and "$\phi$ at $J$" means (with base $2$) $\phi$ is at twice the scale of $\psi$. So here's $\hat\psi$ at $s$, and $\hat\phi$ at $2s$, with their squared sum on right:

Uh, yeah. I've thought on this for a while and came up with nothing. Worse, I've concluded, if it's doable, it's inconsistent:

1. A big deal in scattering is that the actual, measured scale of $\phi$ does not exceed that of any $\psi$. This means we cannot define $\phi$ arbitrarily to force a complement.
2. Following "LP: General", we also see that the shape of $|\hat\phi|$ must match the shape of $|\hat\psi|$ (at same scale). Yet, combining this with $j > -J$ and "phi at $J$", we have one counter-example of things not working, meaning there's infinite counterexamples.

But don't just take my word for it:

Unfortunately this "almost" occurs over the worst possible interval - low frequencies, which are the most dominant for all but first order, per complex modulus. Now, I've read much but not all of the paper, and I think most if not all theorems are developed for $J = \infty$, for which it is all $\omega$. The cases for finite $J$ might be hand-waved, but I've not confirmed. Off-topic but, what I can confirm is, there's no such thing as "scale of invariance" for time-shifts, as is advertised everywhere, for anything but $J = \infty$ - which I'll write about later.

LP: Scattering-practical

Practical implementations complement the CQT filterbank with linearly-spaced wavelets of fixed bandwidth, i.e. STFT (except zero-mean enforcement). Separately, the linear and log portions can perfectly satisfy LP - but I suspect, the transition interval (i.e. $\omega$ between CQT and non-CQT wavelets) is incapable of this.

More importantly, I can't think of any use for $=1$, and $=1$ is more detrimental than helpful for most purposes (per excess redundancy, output size, compute, and other reasons). $\approx 1$ is totally acceptable, but it's important to stay $\leq 1$.

Practical caveats

Some info repeats:

1. CWT must be tiled according to CQT, Constant-Q Transform, Q = (center freq) / bandwidth. It's what's almost always done in literature, but not in implementation. Moreover, in the discrete case, scales must be distributed exponentially, i.e. next_freq / prev_freq = constant.

2. Wavelets must be L1 normed: sum(abs(psi)) == 1. This undoes the $\cdot |a|$ in frequency, so wavelets are only stretched and shifted. Functionally it's $a \psi (t/a)$, which is equivalent per dilation property of integration.

3. L1 norm's not enough, we must rescale by a constant factor depending on how many filters there are, to achieve unity.

4. The non-CQT tiling has a subtle problem: keeping bandwidth same while shifting frequency requires re-parametrizing the mother wavelet, meaning changing its shape, which makes perfectly satisfying LP impossible, even if neglecting the "transition" problem. An option is to relax the zero-mean criterion of wavelets, but this does more harm than good, and strongly so in scattering where DC is dominant in higher orders.

5. Discrete perfect LP is possible with STFT but not CWT, because we cannot attain perfect sampling uniformity in log space over all intervals. There are again workarounds, that again do more harm than good.

6. Nyquist requires careful handling - both the region near it, and itself.

Useful facts

1. The more redundant the discrete CWT filterbank, the better LP is satisfied - which improves one-integral inversion, hence performance of synchrosqueezing.

2. Someone's gone through the grind to account for all of 1-6 caveats. Me! Result below, with library to be released "soon" - "Follow" this post or "Watch" ssqueezepy to be notified.