# Joint Time-Frequency Scattering structure & implementation?

How does JTFS differ from wavelet time scattering in its computation graph, and how does FDTS discriminability work, at a lower level? How is it implemented in practice, and how can one visualize involved ops? Can the coefficients be normalized (for e.g. ML applications)?

JTFS overview provided in this post.

### Computational structure

JTFS breaks the tree structure of time scattering by convolving along frequency, exploiting the joint time-frequency geometry:

\begin{align} S_{(J, J_{fr})}^{{(0)}} x(t) &= x \star \phi_T(t), \\ S_{(J, J_{fr})}^{{(1)}} x(t, \lambda) &= |x \star \psi_\lambda^{{(1)}}| \star \phi_T, \\ S_{(J, J_{fr})}^{{(2)}} x(t, \lambda, \mu, l, s) &= ||x \star \psi_\lambda^{{(1)}}| \star \Psi_{{\mu, l, s}}| \star \Phi_{(T, F)}. \end{align}

$$\Psi_{{\mu, l, s}}$$ comprises of five kinds of joint wavelets:

\begin{align} \Psi_{{\mu, l, +1}}(t, \lambda) &= \psi_\mu^{{(2)}}(t) \psi_{{l, s}}(+\lambda) && \text{spin up bandpass} \\ \Psi_{{\mu, l, -1}}(t, \lambda) &= \psi_\mu^{{(2)}}(t) \psi_{{l, s}}(-\lambda) && \text{spin down bandpass} \\ \Psi_{{\mu, -\infty, 0}}(t, \lambda) &= \psi_\mu^{{(2)}}(t) \phi_F(\lambda) && \text{tm bandpass, fr lowpass} \\ \Psi_{{-\infty, l, 0}}(t, \lambda) &= \phi_T(t) \psi_{{l, s}}(\lambda) && \text{tm lowpass, fr bandpass} \\ \Psi_{{-\infty, -\infty, 0}}(t, \lambda) &= \phi_T(t) \phi_F(\lambda) && \text{joint lowpass} \\ \end{align}

and $$\Phi_{(T, F)}$$ optionally does temporal and/or frequential averaging:

$$\Phi_{(T, F)}(t, \lambda) = \phi_T(t) \phi_F(\lambda)$$

where

• $$t =$$ time
• $$\lambda =$$ log-frequency (wavelet frequency along $$t$$ is $$2^\lambda$$). Rate of sinusoidal oscillation ($$=$$ frequency) of the signal.
• $$\mu =$$ log-frequency (wavelet frequency along $$t$$ is $$2^\mu$$). Frequency of amplitude envelopes of sinusoidal oscillations of the signal (frequency of AM), logscaled.
• $$l =$$ log-quefrency (wavelet quefrency along $$\lambda$$ is $$2^l$$). Frequency of FM bands - rate of variation of frequential bands (independent components/modes) of varying widths, decay factors, and temporal evolutions.

Each $$\Psi$$ captures a distinct structure:

• Spinned: FDTS
• Time lowpass: vertical time-freq geometry (e.g. spikes)
• Freq lowpass: horizontal time-freq geometry (e.g. pure sines)
• Joint lowpass: time-freq mean (dc component)

Higher orders are computed by iterating on unaveraged slices (see "Energy analysis").

Note: $$S^{(2)}$$ with $$\Phi$$ is the implementation variant; one consistent with time scattering's breakdown omits it. $$S^{(1)}$$ also has an alternate breakdown. See "Energy analysis".

Units & quefrency:

• frequency has physical units $$\text{cycles/second} = \text{Hz}$$, discrete units $$\text{cycles/sample}$$
• "quefrency" has units $$\text{cycles/octave}$$. It's a dimensionless, log unit that encodes the physical units in offset of the log axis (1mm vs 1m is a constant shift in log space), and encodes the base as 2 by definition (octave == base 2). It's the rate of variation along the log-frequency axis.

### Pseudocode

Unoptimized implementation is fairly straightforward; showing only up to spin up, unaveraged:

# first order time scattering
U1 = []
for p1t in psi1_t:
U1.append(abs(conv(x, p1t)))

# joint scattering
U2 = []
for p2t in psi2_t:
# second order time scattering
U2_tm = []
for u1 in U1:
U2_tm.append(conv(u1, p2t))
U2_tm = array(U2_tm)  # shape == (freq, time)

# frequential scattering
for p1f in psi1_fr_up:
U2.append(abs(conv(U2_tm, p1f, axis=-2)))


### vs Wavelet Time Scattering

High quefrency coefficients approximate second order time scattering; hence, little discriminative power is lost (without freq averaging). The catch is much greater output size for same time averaging configuration.

### Separable convolutions

The joint filters are defined separably, i.e. as product of two 1D functions. This enables inheriting time scattering's properties along log-frequency, and FDTS discriminability via independent scaling along time and log-frequency. 2D Morlets, as in image wavelet scattering, aren't appropriate, as they introduce rotation, and a rotation does not preserve the relationship between time and frequency (a rotated scalogram generally isn't a scalogram of some other signal).

Further, this yields significantly faster convolutions, and a separable interpretation. To understand FDTS sensitivity, consider echirp after second order time scattering: Each row is formed by convolving psi2_t with a U1 coefficient: Since echirp is a straight line in a scalogram, Gaussian-scaled along time, the result is as shown in the first image.

Frequential scattering will then convolve with columns of U1 * psi2_t; taking the central slice: It's a Morlet! Morlet convolved with Morlet is another Morlet with twice the time and half the frequency support (Gauss * Gauss in freq). Note that the opposite spinned wavelet convolves to zero: conv(analytic, anti_analytic) == 0.

This Morlet, or more importantly the odd symmetry of the imaginary part along frequency, is formed only by a relative misalignment of rows in the same direction - i.e. FDTS. Spin resonance is strongest when the time-log-frequency geometry is a straight line, but other forms of FDTS also work (e.g. linear chirp).

Different FDTS rates are sloped differently in time-frequency, which will resonate with joint wavelets of different relative scalings along time and frequency. This is akin to rotation, but isn't same, as "orientation" is generated by relative stretching of the wavelet.

### Invertibility

Scalogram is invertible within a global phase shift $$e^{j\phi}$$ (see "Invertibility").

Time scattering is invertible within a global phase shift, and time shift $$T$$.

JTFS is invertible within a global phase shift, time shift, and log-frequency shift $$F$$.

Loss of information within a log-frequency shift is much more severe than within a phase or time shift. For a moderately large $$F$$, two orders (which is where most implementations stop) are insufficient to recover information lost to averaging, and time-frequency ridges on scales $$ are blended together. Further, a log frequency shift is more nonlinear than a constant shift, which can significantly recharacterize a signal.

In most cases, such loss is undesired, and unaveraged JTFS is preferred; appropriate amount of averaging, or any further manipulation, can instead be learned.

### Energy analysis

Recall for standard scattering energy decomposition follows a lowpass-bandpass split. For JTFS, it's the same - just in 2D, and a mixture: Right is pure bandpass, left is mix of lowpass and bandpass that gets filed under "lowpass", with pure lowpass at very center. The breakdown is hence:

• S0: x * phi_t
• U1: |x * psi_t|
• S1: U1 * phi_t * phi_f, |U1 * phi_t * psi_f|, |U1 * psi_t * phi_f|
• U2: |U1 * psi_t * psi_f|
• S2: U2 * phi_t * phi_f, |U2 * phi_t * psi_f|, |U2 * psi_t * phi_f|
• ...

An alternative breakdown is

• S1: U1 * phi_t, |U1 * psi_t * phi_f|

replacing phi_t * phi_f with phi_t * psi_f. Note the combined energies and frequency tilings are identical (but first inherits JTFS's properties), as it's phi_t -> phi_t * phi_f, phi_t * psi_f, where * phi_f and * psi_f combined equate to * 1 in frequency response. Likewise * psi_t * psi_f and * psi_t * phi_f are * psi_t's (second order time scattering) energy equivalents.

Lastly the LP sum of the joint filterbank:  This is with simple L1 norm, no energy renormalization, which exceeds the intended upper bound of 1. Deviations from 1 (tight frame) amplify quintically for JTFS, a square up from time scattering.

### Up vs down

Wavelet transform can interpret two ways:

• convolution: response of a physical system to a stimulus, with impulse response being the wavelet filterbank
• cross-correlation: similarities of wavelets with the input (inner product with conjugate of wavelet)

Implementation follows the former per FFT convolution, but latter is the true mathematical formulation and favored for understanding coefficients. Since the wavelets are purely real in Fourier, the two can be shown to be equivalent. As such, the time-domain wavelet is what directly correlates (conj + inner product) with its input (the scalogram).

Observing them in time domain, spin down is crest up (2D maxima are sloped diagonally up). Hence, down resonates with up, as in rising chirp with crest up.

Consistent with separable treatment, we pad convolutions along time and frequency, treating frequential slices same as a signal input to time scattering. JTFS padding, along time or frequency, is tricky, it can break FDTS discriminability; the standard is zero padding for both, but I made reflect work.

### Normalization

Ref 2 proposes a channel log normalization:

$$\widetilde{S}x_n(p) = \log{\left( 1 + \frac{\int_\mathbb{R}Sx_n(t, p) dt}{\epsilon \mu(p)} \right)}$$

where

$$\mu(p) = \underset{1\leq n \leq N}{\text{median}} \int_{\mathbb{R}} Sx_n (t, p) dt$$

is median of path/channel $$p$$ aggregated over $$N$$ samples. This has the effect of "Gaussianizing" the histogram of coefficients, and works well in the proposed ML framework, with coefficients fed as frequency vectors and time integrated out.

Time integration can be omitted in the general case, leaving

$$\widetilde{S}x_n(t, p) = \log{\left( 1 + \frac{Sx_n(t, p)}{\epsilon \mu(p)} \right)}$$

and $$\epsilon$$ is a user-chosen constant which controls the extent of log scaling (bringing different orders of magnitude closer together).

This approach bears two problems, however:

1. Relative scaling of paths changes. If dataset lacks content in high frequencies, then high frequency features will be uninformative or be noise, while channel norm will bring their energy to same level as informative ones, lowering SNR.
2. Breaks spatial coherence necessary for higher-dimensional convs, like scale-shifting every individual timestep for 1D convs Blue = before, orange = after. Path 1 (max amplitude 20) is clearly dominant to path 0 (max 1), but post-norm they're of ~equal intensity.

A solution is sample normalization, but it's not a win-all. Log also has drawbacks, which PCEN improves on -- see further discussion.

### Example: pure sine Joint slices: Why are all spinned coefficients zero?

Wavelet convolved with a constant is zero, since wavelets are zero mean and a constant == DC. |conv(psi1, sine)| = U1 = const (first order yields constants along rows), and conv(psi2, U1) = 0 (second order convolves over those constants). Only the phi_t * pairs are sensitive to constant horizontal geometries, as they're non-zero mean along time.

JTFS is doper