[WIP] Introduce xi_min

[OverLordGoldDragon]

Introduces explicit guard against bandpasses ending up as pure sines.

See images. Attaining a proper wavelet can be very costly on padding, but at the same time we want to tile the entire frequency axis; I thus introduced max_pad_factor in latest JTFS, so user can choose faster compute at expense of boundary effects in small fraction of coefficients.

Rationale for xi_min = 2/N is, minimal possible wavelet is made with 3 samples in frequency domain (arguably 2), so peak cannot lie on bin 1, thus we put it at bin 2. As shown in "this PR only", this alone doesn't guarantee a good wavelet - it does, however, spare much padding that'd be required to make psi1_f[-1] a proper wavelet.

Note that other PRs in progress will allow a bandpass to become exactly a pure sine, exacerbating the problem.

Current behavior

After this PR + others

After this PR only

Plots code

show

# run from different branches
import numpy as np
from kymatio.numpy import Scattering1D
from ssqueezepy import ifft
from ssqueezepy.visuals import plot, scat

N = 2048 * 4
J = int(np.log2(N) - 1)
Q = 16
scattering = Scattering1D(J, N, Q, max_pad_factor=None)

p_fr = scattering.psi1_f[-1][0]
p_t = ifft(p_fr)

scat(p_fr[:30], show=1, title="psi1_f[-1][:30]")
plot(p_t, complex=1)
plot(p_t, abs=1, linestyle='--', color='k', title="ifft(psi1_f[-1])")

[OverLordGoldDragon]

Problem is worse than I thought, it's not just "last few filters" but last few octaves; for large Q dozens of filters never decay to zero. With J, Q = 11, 16, here's psi1_f[-30], meaning 29 more wavelets have worse decay than this:

This is remedied by padding more, rather than via xi_min (though it helps by requiring less padding).

[OverLordGoldDragon]

This can merge at any point (preferably after #673).

[lostanlen]

so this PR effectively reduces the number of filters in the filterbank while leaving everything else the same?

i see your line of thinking for this bounds check but why allow the user to override it? wouldn't that add unnecessary complexity to filter design?

[OverLordGoldDragon]

It's passed internally by scattering_filter_factory, not user-exposed.

[OverLordGoldDragon]

Post inspired by #800.

Design concerns A) feature quality, and B) preserving of information. As we pad more, we lower the lowest possible frequency, and hence add more room for low-frequency wavelets. Likewise, more padding means populating more lower frequencies of the input signal, fft(x_pad). Incomplete tiling == inability to invert, hence loss of information... yet at the same time, padding == adding no information. Who "wins"?

First, the wavelets: top row setting the lowest frequency at input's minimum (DFT bin 1), bottom row setting it at padded's bin 1:

As a reminder, at this scale (J = log2(N)) we require padding, else we get a pure sine for a wavelet: (note, while the two psi_t look similar, only one suffers from circular conv.)

Taking x = randn(N) and comparing with padded:

Padding is by a factor of 16, hence unpadded bin 1 moves to padded bin 16. Yet, bins 1 to 15 are filled, and if we don't tile these, information is lost. However,

has nothing in 1-15, hence the answer: what's lost is what's needed to recover the padded portion. The "extra" frequencies in zero or whichever padding are required to encode the fact that our transform is to invert to zeros outside the original input's interval.

• Caveat, bit of a lapse in reasoning: the padded portion serves not only as an "optional" segment to be inverted, but also as a contributing interval to a large kernel convolution. Since xi=1/N is all that's needed to fully tile the unpadded signal, and since the trimmed temporal waveforms are close in decay and identical in modulation... in zero-pad case we have N slidings of the wavelet with same carrier but different windowings, so all that'll differ for recovery is the un-windowing.

This covers B). For A, we ask whether a wavelet that completes less than one cycle over the input is important. This wavelet correlates to zero with periodic padding, as evident from freq-domain, but not with other padding. While both 1/N and 1/N_pad for xi have the same support (since non-CQT portion) of which the original input comprises only a small fraction, a lower frequency requires a longer observation interval to ascertain the measurement; in this sense 1/N_pad is worse.

Favoring 1/N also helps keep filterbank design independent of padding settings, convenient for a number of reasons, particularly in JTFS. Last point concerns: 1/N or 2/N? Latter helps guarantee a wavelet as opposed to a pure tone, as a wavelet is minimally three (arguably two) points in frequency domain, and we can't plant at DC so the peak must be at 2. More importantly, insufficient padding with J=log2(N) yields a near-pure sine regardless, which risks not tiling bin 1 at all. Hence, 1/N.

[OverLordGoldDragon]

"WIP" just to change to 1/N and see if tests bork

[OverLordGoldDragon]

Missed edge case; 'reflect' can meaningfully increase variability:

code

import numpy as np
from numpy.fft import rfft
import matplotlib.pyplot as plt

t = np.linspace(0, 1, 129, 1)
g = np.cos(2*np.pi * .5 * t)
gp = np.pad(g, [64, 63], mode='reflect')
gf, gpf = rfft(g), rfft(gp)

tkw = dict(weight='bold', fontsize=17, loc='left')
fig, axes = plt.subplots(1, 2, sharey=True)
axes[0].plot(g)
axes[1].plot(gp)
axes[0].set_title("x", **tkw)
axes[1].set_title("x_pad", **tkw)
fig.subplots_adjust(left=0, right=1, wspace=.02)
fig.set_size_inches(11, 4.5)
plt.show()

fig, axes = plt.subplots(1, 2, sharey=True)
axes[0].plot(np.abs(gf))
axes[1].plot(np.abs(gpf))
axes[1].scatter(np.arange(len(gpf))[:3], np.abs(gpf)[:3])
axes[0].set_title("|rfft(x)|", **tkw)
axes[1].set_title("|rfft(x_pad)|", **tkw)
fig.subplots_adjust(left=0, right=1, wspace=.02)
fig.set_size_inches(11, 4.5)
plt.show()

Half cycle -> full cycle, and the padded wavelet whose frequency is now lesser than it was possible without padding will capture this frequency more accurately than the fully invertible filterbank of unpadded, winning in "feature quality" and tying in invertibility.

Whether this is plausible or not is a secondary matter; the only important downside is making the filterbank depend upon padding settings. But this is meaningful, since the padded part cannot simply be "detached" from the original (despite me trying to do this). It's also arguable whether "feature quality" is improved since the half cycle persisting for a full cycle is a prior - but by retaining full tiling of padded we assure everything's captured for any padding choice. The added wavelets are max j, so again minimal speed and no memory overhead.

Hence, xi_min = 1/2**J_pad. This PR does 2/2**J_pad, but also in JTFS I made it so that xi_min is computed inside of compute_params_filterbank. Might update later, or just keep this as a "concept" PR and merge the real thing with JTFS.