# Musical Metamerism with Time-Frequency Scattering

Vincent Lostanlen, Han Han

► To cite this version:

Vincent Lostanlen, Han Han. Musical Metamerism with Time-Frequency Scattering. 2024. hal-05499572

HAL Id: hal-05499572

https://hal.science/hal-05499572v1

Preprint submitted on 8 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# MUSICAL METAMERISM WITH TIME-FREQUENCY SCATTERING

Vincent Lostanlen and Han Han

November 2024

*Abstract.* — The concept of metamerism originates from colorimetry, where it describes a sensation of visual similarity between two colored lights despite significant differences in spectral content. Likewise, we propose to call “musical metamerism” the sensation of auditory similarity which is elicited by two music fragments which differ in terms of underlying waveforms. In this technical report, we describe a method to generate musical metamers from any audio recording. Our method is based on joint time-frequency scattering in Kymatio, an open-source software in Python which enables GPU computing and automatic differentiation. The advantage of our method is that it does not require any manual preprocessing, such as transcription, beat tracking, or source separation. We provide a mathematical description of JTFS as well as some excerpts from the Kymatio source code. Lastly, we review the prior work on JTFS and draw connections with closely related algorithms, such as spectrotemporal receptive fields (STRF), modulation power spectra (MPS), and Gabor filterbank (GBFB).

## Contents

<table><tr><td>1. Introduction</td><td>2</td></tr><tr><td>2. Method</td><td>3</td></tr><tr><td>3. Implementation</td><td>6</td></tr><tr><td>4. Related work</td><td>9</td></tr><tr><td>5. Conclusion</td><td>11</td></tr><tr><td>References</td><td>12</td></tr></table>

*Key words and phrases.* — audio texture synthesis, automatic differentiation, joint time-frequency scattering (JTFS), music cognition.
# 1. Introduction

**1.1. Problem statement.** — What makes music memorable? The ability of human listeners to recognize a song “by ear” may be easy to observe, it remains difficult to explain [1]. For example, a familiar song may be identified after having undergone a certain amount of frequency transposition or time stretching. Likewise, prior research in cognitive science has demonstrated that songs may retain their familiarity after noise vocoding—a digital audio effect which locally alters spectral content while leaving the temporal envelope essentially unchanged [2].

Together, these findings hint at the idea that music familiarity is not reducible to the recognition of rigid patterns in the time-frequency domain. Instead, Henkjan Honing et al. have suggested that listeners attend to the global “contour” of a musical piece when perceiving it as familiar [3]. Here, the notion of contour is an umbrella term for various auditory attributes: i.e., melodic, harmonic, rythmic, timbral, textural, and so forth. The commonality between all these components of contour perception resides in their invariance to local alterations of the short-term power spectrum, such as those mentioned earlier: noise vocoding, time stretching, and frequency transposition. Thus, Honing’s hypothesis, which we refer to as “multicomponent contour perception”, is consistent with the available evidence on music familiarity [4].

Yet, there is a gap in scholarship concerning the relative weight of contour components in the cognition of familiar music pieces. Measuring these weights experimentally is challenging, for lack of suitable “minimal pairs”: i.e., pairs of auditory stimuli which differ from each other just enough to significantly reduce perceived familiarity while having some contour components in common [5]. The prospect of synthesizing such minimal pairs as auditory stimuli raises a difficult challenge for signal processing research, but one that is worth pursuing because of its application potential.

**1.2. Contribution.** — In this report, we present a new method for stimulus design in music cognition experiments. Our method relies on a nonlinear multiresolution approximation algorithm known as joint time-frequency scattering, or JTFS for short. JTFS decomposes the spectrogram into wavelet subbands at various temporal rates and frequential scales while retaining a fine localization of coefficients in the time-frequency domain. The Euclidean distance between two JTFS feature vectors is differentiable with respect to the underlying digital audio samples, thus lending itself to signal reconstruction via gradient descent.

Following Feather et al., we refer to the outcome of this reconstruction by “metamer” [6, 7]; i.e., a signal that is perceptually indistinguishable from the reference even so the underlying waveforms may significantly differ. Although
the synthesis of JTFS metamers has found applications in computer music [8],
the same cannot yet be said of music cognition research.

The key idea is to coarsen the spectrotemporal discretization of JTFS co-
efficients by means of Gaussian low-pass filtering over the time and frequency
variables. After this coarsening, Euclidean distances in JTFS feature space are
invariant to small shifts up to some time scale *T* and to small transposition up
to some frequency interval *F*. Such a property of invariance aligns with the
“contour hypothesis” of Honing *et al.*, i.e., the claim that music listeners primar-
ily attend to relative variations in spectrotemporal energy distribution in the
time-frequency domain as opposed to absolute time-frequency cues. There-
fore, resynthesizing a waveform from coarsened JTFS coefficients produces a
kind of “musical metamer”, that is, a signal which resembles the original in
terms of summary statistics over spectrotemporal regions of duration of the
order of *T* and pitch register of the order of *F*.

## 2. Method

This section reproduces some elements from [8], with minor adaptations.

**2.1. Joint time-frequency scattering.** — Time-frequency scattering results from the cascade of two stages: an auditory filterbank and the extraction of spectrotemporal modulations with wavelets in time and log-frequency. First, we define Morlet wavelets of center frequency $\lambda > 0$ and quality factor $Q$ as

$$
(1) \qquad \psi_{\lambda}(t) = \lambda \exp\left(-\frac{\lambda^2 t^2}{2Q^2}\right) \times (\exp(2\pi i \lambda t) - \kappa),
$$

where the corrective term $\kappa$ ensures that each $\psi_{\lambda}(t)$ has one vanishing moment, i.e., a null average. Within a discrete setting, acoustic frequencies $\lambda$ are typically of the form $2^{n/Q}$ where $n$ is integer, thus covering the hearing range. For $x(t)$ a finite-energy signal, we define the CQT of $x$ as the matrix $U_1 x(t, \lambda) = |x * \psi_{\lambda}|(t)$; which stacks convolutions with all wavelets $\psi_{\lambda}(t)$ followed by the complex modulus nonlinearity.

Secondly, we define Morlet wavelets of respective center frequencies $\alpha > 0$ and $\beta \in \mathbb{R}$ with quality factor $Q = 1$. With a slight abuse of notation, we denote these wavelets by $\psi_{\alpha}(t)$ and $\psi_{\beta}(\log \lambda)$ even though they do not necessarily have the same shape as the wavelets $\psi_{\lambda}(t)$ of Equation 1. Frequencies $\alpha$, hereafter called amplitude modulation rates, are measured in Hertz (Hz) and discretized as $2^n$ with integer $n$. Frequencies $\beta$, hereafter called frequency modulation scales, are measured in cycles per octave (c/o) and discretized as $\pm 2^n$ with integer $n$. The edge case $\alpha = 0$ corresponds to $\psi_{\alpha}(t)$ being a Gaussian low-pass filter $\phi_T(t)$ of temporal width $T^{-1}$, measured in Hertz. Likewise, the edge case $\beta = 0$ corresponds to $\psi_{\beta}(\log \lambda)$ being a Gaussian low-pass filter
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Interference pattern between wavelets $\psi_{\alpha}(t)$ and $\psi_{\beta}(\log_2 \lambda)$ in the time-frequency domain $(t, \log_2 \lambda)$ for different combinations of amplitude modulation rate $\alpha$ and frequency modulation scale $\beta$. Darker shades of red (resp. blue) indicate higher positive (resp. lower negative) values of the real part.</figcaption></figure>

$\phi_F(\log \lambda)$ of bandwidth $F^{-1}$, measured in cycles per octave. These modulation scales $\beta$ play the same role as the quefrencies in a mel-frequency cepstrum. We define the fourth-order tensor

$$ (2) \quad \mathbf{U}_2\mathbf{x}(t, \lambda, \alpha, \beta) = |\mathbf{U}_1\mathbf{x}^* \psi_{\alpha}^{\log_2 \lambda} \psi_{\beta}|(t, \lambda) \\ = \left| \iint \mathbf{U}_1\mathbf{x}(\tau, s) \psi_{\alpha}(t-\tau) \psi_{\beta}(\log_2 \lambda - s) \mathrm{d}\tau \mathrm{d}s \right|, $$

which stacks convolutions in time and log-frequency with all wavelets $\psi_{\alpha}(t)$ and $\psi_{\beta}(\log_2 \lambda)$ followed by the complex modulus nonlinearity. Figure 1 shows the interference pattern of the product $\psi_{\alpha}(t-\tau)\psi_{\beta}(\log_2 \lambda - s)$ for different combinations of time $t$, frequency $\lambda$, rate $\alpha$, and scale $\beta$. We denote the multiindices $(\lambda, \alpha, \beta)$ resulting from such combinations as scattering paths [9].

**2.2. Averaging.** — Because it is a convolutional operator in the time-frequency domain, the STRF is equivariant to temporal translation $t \mapsto t + \tau$ as well as frequency transposition $\lambda \mapsto 2^s\lambda$. To produce musical metamers, we guarantee invariance to temporal translation up to some time lag $T$ and up to some frequency interval $F$. To this aim, we define time-frequency scattering as the result of a local averaging of both $\mathbf{U}_1\mathbf{x}(t, \lambda)$ and $\mathbf{U}_2\mathbf{x}(t, \lambda, \alpha, \beta)$ by a Gaussian low-pass filter $\phi_T$ of cutoff frequency equal to $T^{-1}$, yielding

$$ (3) \quad \mathbf{S}_1\mathbf{x}(t, \lambda) = (\mathbf{U}_1\mathbf{x}^* \phi_T^{\log \lambda} \phi_F)(t, \lambda) $$

and

$$ (4) \quad \mathbf{S}_2\mathbf{x}(t, \lambda, \alpha, \beta) = (\mathbf{U}_2\mathbf{x}^* \phi_T^{\log_2 \lambda} \phi_F)(t, \lambda, \alpha, \beta), $$
respectively.

**2.3. Reconstruction.** — By recursion over layers in the composition of wavelet modulus operators, the invertibility of wavelet modulus implies the invertibility of scattering transforms of infinite depth [10], up to a constant time shift of at most T. However, the time-frequency scattering network presented here has a finite number of layers (i.e., two layers) and is therefore not exactly invertible. Thanks to the theorem of exponential decay of scattering coefficients [11], the residual energy that is present in deeper layers can be neglected on the condition that T is small enough.

Starting from a colored Gaussian noise $y_0(t)$ whose power spectral density matches $S_1x(t, \lambda)$, we refine it by additive updates of the form $y_{n+1}(t) = y_n(t) + u_n(t)$, where the term $u_n(t)$ is defined recursively as $u_n(t) = m \times u_n(t) + \mu_n \nabla E(y_n)(t)$. In subsequent experiments, the momentum term is fixed at $m = 0.9$ while the learning rate is initialized at $\mu_0 = 0.1$ and modified at every step according to a “bold driver” heuristic [12].

Like deep neural networks, scattering networks consist of the composition of linear operators (wavelet transforms) and pointwise nonlinearities (complex modulus). Consequently, the gradient $\nabla E(y_n)$ can be obtained by composing the Hermitian adjoints of these operator in the reverse order as in the direct scattering transform – a method known as backpropagation.

We rely on the PyTorch backend of Kymatio to perform gradient backpropagation via reverse-mode automatic differentiation [13]. We outline the resulting automatic differentiation algorithm below. First, Kymatio backpropagates the gradient of Euclidean loss for second-order scattering coefficients:

$$ (5) \qquad \nabla U_2 y(t, \lambda, \alpha, \beta) = ((S_2 x - S_2 y) * \phi)(t, \lambda, \alpha, \beta). $$

Secondly, Kymatio backpropagates the second layer onto the first:

$$ \begin{aligned} \nabla U_1 y(t, \lambda) = & ((S_1 x - S_1 y) * \phi)(t, \lambda, \alpha, \beta) \\ & + \sum_{\alpha, \beta} \Re \left( \left[ \frac{(U_1 y * \bar{\psi}_{\alpha} * \bar{\psi}_{\beta})^{\log \lambda}}{|U_1 y * \bar{\psi}_{\alpha}|^{\log \lambda} |\bar{\psi}_{\beta}|} \times \nabla U_2 y \right] * \psi_{\alpha} * \psi_{\beta} \right) (t, \lambda, \alpha, \beta), \end{aligned} $$

where the symbol $\Re(z)$ denotes the real part of the complex number $z$. Lastly, Kymatio backpropagates the first layer into the waveform domain:

$$ (7) \qquad \nabla E(y)(t) = \sum_{\lambda} \Re \left( \left[ \frac{y * \psi_{\lambda}}{|y * \psi_{\lambda}|} \times \nabla U_1 y \right] * \psi_{\lambda} \right) (t) $$

Observe that the signal above is null if $S_1x = S_1y$ and $S_2x = S_2y$; i.e., if $x$ and $y$ are musical metamers for the JTFS representation.
### 3. Implementation

**3.1. Filterbanks.** — Let $(\psi_0^{\text{tm},1}, \dots, \psi_{N_1-1}^{\text{tm},1})$ be a finite sequence of complex-valued band-pass filters with center frequencies $1 > \xi_0^{\text{tm},1} > \dots > \xi_{N_1-1}^{\text{tm},1} > 0$ and bandwidths $\frac{1}{2} > \sigma_0^{\text{tm},1} > \dots > \sigma_{N_1-1}^{\text{tm},1} > 0$, where the sample rate is set equal to 1 by convention. We refer to [14] for the rationale behind the values of $\xi_0^{\text{tm},1}$ and $\sigma_0^{\text{tm},1}$.

Given a scale $J$ and a dimensionless hyperparameter $Q^{\text{tm},1}$, center frequencies and bandwidths decrease according to a geometric progression, while maintaining a quality factor of $Q^{\text{tm},1}$; until reaching an “elbow”, defined as a lower bound on the bandwidth. After the elbow, the progression of center frequencies is arithmetic, while bandwidths remain at the lower bound. Thus, filters have a constant-bandwidth property in the lower-frequency range and a constant-$Q$ property in the upper-frequency range. This construction was first proposed in [15] and is reminiscent of the construction of the mel scale in HTK$^{(1)}$. We reproduce an excerpt of the Kymatio code below.

```
1 sigma0 = 0.1
2 r_psi = math sqrt(0.5)
3
4 xi = max(1. / (1. + math.pow(2., 3. / Q)), 0.35)
5 factor = 1. / math.pow(2, 1. / Q)
6 sigma = (1 - factor) / (1 + factor) * xi / math sqrt(2 * math.log(1. / r))
7 sigma_min = sigma0 / 2**J
8
9 if sigma <= sigma_min:
10 xi = sigma
11 else:
12 yield xi, sigma
13 # High-frequency (constant-Q) region
14 while sigma > (sigma_min * math.pow(2, 1 / Q)):
15 xi /= math.pow(2, 1 / Q)
16 sigma /= math.pow(2, 1 / Q)
17 yield xi, sigma
18
19 # Low-frequency (constant-bandwidth) region
20 elbow_xi = xi
21 for q in range(Q-1):
22 xi -= 1/Q * elbow_xi
23 yield xi, sigma_min
```

Another filterbank is built for the second layer: $(\psi_0^{\text{tm},2}, \dots, \psi_{N_2-1}^{\text{tm},2})$, with its own hyperparameters $Q^{\text{tm},2}$. Lastly, we build a filterbank for frequential scattering: $(\psi_0^{\text{fr},1}, \dots, \psi_{N_{\text{fr}}-1}^{\text{fr},1})$, with its own hyperparameters $Q^{\text{fr},1}$ and $J^{\text{fr},1}$.

**3.2. First layer.** — Given an audio signal $x$, we extract its time-frequency representation by convolution with each filters $\psi_{n_1}^{\text{tm},1}$ where $0 \le n_1 < N_1$. For reasons of computational efficiency, we map the signal to the Fourier domain, via a real-input direct fast Fourier transform (RFFT). The filtering itself is

<footnote>(1) https://librosa.org</footnote>
an element-wise product between complex-valued vectors, implemented with the CDGMM routine. Furthermore, the filtering is immediately followed by a subsampling operation, implemented in the Fourier domain by reshaping and summation. The subsampling factor $2^{k1}$, depends on the resolution of the cutoff frequency of the filter $\psi_{n1}^{\text{tm},1}$, i.e., $j1$; and on the choice of stride, i.e., $T$. We reproduce an excerpt of the Kymatio code below.

```python
1 U_0Hat = backend.rfft(U_0)
2 for n1 in psi1:
3     j1 = psi1[n1][j']
4     k1 = min(j1, log2 stride)
5     U_1c = backend.cdgmm(U_0 hat, psi1[n1][levels][0])
6     U_1 hat = backend.subsample fourier(U_1c, 2 ** k1)
7     U_1c = backend.ifft(U_1 hat)
8     U_1m = backend.modulus(U_1c)
9     U_1 hat = backend.rfft(U_1 m)
10 U_1 hats.append({'coef': U_1 hat, 'j': (j1,), 'n': (n1,)})
```

To reduce the influence of boundary effects, this operation is preceded by padding on both sides, and followed by unpadding on both sides, not shown in the code above.

**3.3. Second layer, time variable.** — Contrary to Kymatio's time scattering module, which relies on a depth-first search over scattering paths for memory efficiency, time scattering in Kymatio's JTFS module is done according to a width-first search. In other words, the outer loop corresponds to the deeper path index ($n_2$) whereas the inner loop corresponds to the shallower path index ($n_1$). This is because all the first-order paths corresponding to a given second-order paths must be available in memory when performing scattering over the log-frequency variable (see next subsection). We reproduce an excerpt of the Kymatio code below.

```python
1 psi2 = filters[2]
2 for n2 in range(len(psi2)):
3     j2 = psi2[n2][j']
4     Y_2_list = []
5     for U_1 hat in U_1 hats:
6         j1 = U_1 hat[j'][0]
7         k1 = min(j1, log2 stride) if average_local else j1
8     if j2 > j1:
9         k2 = min(j2-k1, log2 stride) if average_local else (j2-k1)
10        U_2_c = backend.cdgmm(U_1 hat['coef'], psi2[n2][levels][k1])
11        U_2 hat = backend.subsample fourier(U_2_c, 2 ** k2)
12        U_2_c = backend.ifft(U_2 hat)
13        Y_2_list.append(U_2_c)
14     if len(Y_2_list) > 0:
15         # Stack over the penultimate dimension with shared n2.
16         # Y_2 is a complex-valued 3D array indexed by (batch, n1, time)
17         Y_2 = backendstack(Y_2_list, 2)
18         # Y_2 is a stack of multiple n1 paths so we put (-1) as placeholder.
19         # n1 ranges between 0 (included) and n1_max (excluded), which we store
20         # separately for the sake of meta() and padding/unpadding downstream.
21         yield {'coef': Y_2, 'j': (-1, j2), 'n': (-1, n2), 'n1_max': len(Y_2_list)}
```
**3.4. Second layer, frequency variable.** — Since band-pass filters $\psi_0^{\text{tm},1}$, ..., $\psi_{N_1-1}^{\text{tm},1}$ are complex-valued, scattered signals ($|x * \psi^{\text{tm},1}| * \psi^{\text{tm},2}$) are also complex-valued. Thus, it is necessary to include negative frequencies as well as positive frequencies. For this purpose, we flip the sign of all center frequencies $\xi_n^{\text{fr},1}$ while maintaining the same bandwidths $\sigma_n^{\text{fr},1}$; an operation known as wavelet *spinning*. We reproduce an excerpt of the Kymatio code below.

```python
def spin(filterbank_fn, filterbankKeywords):
    def spinned_fn(J, Q, **kwargs):
        yield from filterbank_fn(J, Q, **kwargs)
        for xi, sigma in filterbank_fn(J, Q, **kwargs):
            yield -xi, sigma
    return spinned_fn, filterbankKeywords
```

Importantly, spinning is only necessary if $\alpha \neq 0$: indeed, the case $\alpha = 0$ corresponds to a temporal low-pass filtering of the scalogram, which produces a real-valued input. We accommodate both subcases in Kymatio via the `spinned` Boolean argument:

```python
def frequency_scattering(X, backend, filters_fr, log2Stride_fr,
                          average_local_fr, spinned):
    # Unpack filters_fr list
    phi, psis = filters_fr
    log2_F = phi['j']

    # Swap time and frequency axis
    X_T = backend swp_time_frequency(X['coef'])

    # Zero-pad frequency domain
    pad_right = phi['N'] - X['n1_max']
    X_pad = backend.pad_frequency(X_T, pad_right)

    # Spinned case switch
    if spinned:
        # Complex-input FFT
        X_hat = backend.cfft(X_pad)
    else:
        # Real-input FFT
        X_hat = backend.rfft(X_pad)
        # Restrict to nonnegative spins
        psis = filter(lambda psi: psi['xi']>=0, psis)
    for n_fr, psi in enumerate(psis):
        j_fr = psi['j']
        spin = np.sign(psi['xi'])
        k_fr = min(j_fr, log2Stride_fr) if average_local_fr else j_fr
        Y_fr_hat = backend.cdgmm(X_hat, psi['levels'][0])
        Y_fr_sub = backend.subsample_fourier(Y_fr_hat, 2 ** k_fr)
        Y_fr = backend.ifft(Y_fr_sub)
        Y_fr = backend swp_time_frequency(Y_fr)
        # If not spinned, X['n'] = S1['n'] = (-1,) is a 1-tuple.
        # If spinned, X['n'] = Y2['n'] = (-1,2) is a 2-tuple.
        # In either case, we elide the "-1" placeholder and define the new 'n'
        # as (X['n'][1:] + (n_fr,)), i.e., n=(n_fr,) is not spinned
        # and n=(2, n_fr) if spinned. This 'n' tuple is unique.
        yield {**X, 'coef': Y_fr, 'n': (X['n'][1:] + (n_fr,)),
                'j_fr': (j_fr,), 'n_fr': (n_fr,), 'n1 stride': (2**j_fr),
                'spin': spin}
```
**3.5. Main algorithm.** — The main JTFS algorithm has only six lines of code, thanks to the `yield` coroutine and the `yield from` expression, added in Python 3.3.

```python
def joint_timefrequency_scattering(U_0, backend, filters, log2Stride,
                                   average_local, filters_fr, log2Stride_fr, average_local_fr):
    """
    4 Yields
    5
    6 # Zeroth order
    7 if average_local:
    8     * S_0 indexed by (batch, time[log2_T])
    9 else:
   10     * U_0 indexed by (batch, time)
   11
   12 # First order
   13 for n_fr < len filters_fr):
   14     * Y_1_fr indexed by (batch, n1[n_fr], time[log2_T]), complex-valued,
   15     where n1 has been zero-padded to size N_fr before convolution
   16
   17 # Second order
   18 for n2 < len(psi2):
   19     for n_fr < len filters_fr):
   20         * Y_2_fr indexed by (batch, n1[n_fr], time[n2]), complex-valued,
   21         where n1 has been zero-padded to size N_fr before convolution
   22
   23 time_gen = time_scattering_widthfirst(
   24     U_0, backend, filters, log2Stride, average_local)
   25 yield next(time_gen)
   26
   27 # First order: S1(n1, t) = (|x*psi_{n1}|*phi)(t[log2_T])
   28 S_1 = next(time_gen)
   29
   30 # Y_1_fr_{n_fr}(n1, t[log2_T]) = (|x*psi_{n1}|*phi*psi_{n_fr})(t[log2_T])
   31 yield from frequency_scattering(S_1, backend,
   32                               filters_fr, log2Stride_fr, average_local_fr, spinned=False)
   33
   34 # Second order: Y_2_{n2}(t[log2_T], n1[j_fr])
   35 # = (|x*psi_{n1}|*psi_{n2})(t[j2], n1[j_fr])
   36 for Y_2 in time_gen:
   37
   38 # Y_2_fr_{n2,n_fr}(t[j2], n1[j_fr])
   39 # = (|x*psi_{n1}|*psi_{n2}*psi_{n_fr})(t[j2], n1[j_fr])
   40 yield from frequency_scattering(Y_2, backend, filters_fr,
   41                               log2Stride_fr, average_local_fr, spinned=True)
```

## 4. Related work

**4.1. Applications to audio classification.** — The first implementation of JTFS was made in 2014 by Joakim Andén and was written in MATLAB<sup>(2)</sup>. In 2015, Vincent Lostanlen wrote another MATLAB implementation enabling custom routines for gradient backpropagation <sup>(3)</sup>. The first JTFS article relied on the first implementation for speech classification and on the second implementation for texture synthesis of bird sounds [16]. A mathematical link between JTFS and deep convolutional network was outlined by [17] and

<footnote><sup>(2)</sup> https://github.com/scatnet/scatnet</footnote>

<footnote><sup>(3)</sup> https://github.com/lostanlen/scattering.m</footnote>
formally built by [18], with new applications to the supervised classification of music instruments and urban sounds. Since then, JTFS has been applied to other domains: playing techniques from Chinese bamboo flute [19], bioacoustics [20], and medical acoustics [21]. Another publication has shown that Euclidean distances in JTFS feature space predict auditory judgments of timbre similarity within a large vocabulary of instrumental playing techniques, as collected from a group of 31 professional composers and non-expert music listeners [22]. An article for the French academy of sciences, written on the anniversary of Joseph Fourier's birth, discusses the historical background of JTFS [23].

**4.2. Use in computer music.** — The earliest use of JTFS in computer music was *FAVN*, an opera which was commissioned to composer Florian Hecker by the Alte Oper in Frankfurt in 2016, followed by an exhibition at the Kunsthalle in Vienna in 2017. The catalog of this exhibition discusses the underpinnings of JTFS in terms of programming language theory, as well as connections with auditory neurophysiology [24]. This discussion was expanded in [8] on the occasion of *The Shape of RemiXXXes to Come*: i.e., a Lorenzo Senni record, released by Warp Records, containing a JTFS-based remix by Florian Hecker. JTFS was also used in other pieces, such as *Synopsis Seriation* [25], released by Urbanomic, and *Resynthese FAVN*, released as a 10-CD box set by Blank Forms. The use of JTFS in music has been discussed by art historian Ina Blom [26] and by philosopher Amy Ireland [27].

**4.3. Kymatio implementation.** — JTFS was re-implemented in the Python language by Mathieu Andreux in 2018 [28], with an application to generative audio models. The next year, the Kymatio consortium was officially created, with the goal of providing a stable, portable, and efficient implementation of scattering transforms in Python<sup>(4)</sup>. The development of time scattering in Kymatio was led by Andén, Andreux, and Lostanlen, with important contributions by John Muradeli, Changhong Wang, and Philip Warrick [13]. A first implementation of GPU-enabled JTFS from Kymatio's Pytorch routines was written by John Muradeli, under the name of wavespin<sup>(5)</sup>, with new applications to manifold learning and music instrument classification [29]. The current implementation of JTFS in Kymatio was written by Vincent Lostanlen, Cyrus Asfa, Han Han, Changhong Wang, Muawiz Chaudhary, and Joakim Andén<sup>(6)</sup>. Since then, Chris Mitcheltree has reimplemented JTFS in PyTorch for pedagogical purposes<sup>(7)</sup>.

<footnote>(4) https://kymat.io</footnote>

<footnote>(5) https://github.com/OverLordGoldDragon/wavespin/tree/dafx2022-jtfs</footnote>

<footnote>(6) https://github.com/kymatio/kymatio</footnote>

<footnote>(7) https://github.com/christhetree/jtfst_implementation</footnote>
Reverse-mode automatic differentiation in Kymatio has been successfully applied to perceptual sound matching [30]. The combination of JTFS, differentiable physical modeling synthesis, and Riemannian geometry has led to an original technique for training neural networks to solve inverse problems, known as perceptual-neural-physical sound matching [31, 32]. This is a kind of model-based deep learning for music information research—we refer to [33] for an introduction.

**4.4. Spectrotemporal receptive fields.** — Besides JTFS, there are other algorithms which extract spectrotemporal modulations at various scales and rates. As early as 1998, [34] have implemented spectro-temporal modulation transfer functions so as to fit neurophysiological responses to ripple noise stimuli, as previously measured in the auditory cortex of ferrets [35]. We refer to [36] for a complete mathematical description of this implementation, known as “full cortical model” [37], and to [38] for a review of applications. A MATLAB implementation of STRF, including the synthesis of sounds from time-averaged (but not frequency-averaged) coefficients, was achieved by Sam Norman-Haigneré [39] (8). We refer to [40] for an introduction to STRF in the interdisciplinary context of music cognition and music information retrieval (MIR). Recently, STRF has been applied to other tasks, such as bioacoustics [41] and industrial acoustics [42].

Another implementation of STRF was proposed in 2012 by [43], under the name of Gabor filterbanks (GBFB), with applications to automatic keyword spotting in noisy environments. In music cognition, STRF is sometimes known under the name of “modulation power spectrum”, with applications to timbre perception research [44].

**5. Conclusion**

We have described a method to generate musical metamers from any audio recording. Our method is based on joint time-frequency scattering in Kymatio, an open-source software in Python which enables GPU computing and automatic differentiation. The advantage of our method is that it does not require any manual preprocessing, such as transcription, beat tracking, or source separation. Hence, it is suitable to any kind of audio signal, without need to postulate the presence of “Western classical” cues such as tonality, counterpoint, meter, or theme. Lastly, because the resynthesis operates by iterative refinement from a random initial guess, it is compatible with the methodological guidelines of reproducible research: by running the algorithm anew, one

<footnote>(8) https://github.com/snormanhaignere/spectrotemporal-synthesis-v2</footnote>
may generate an infinite number of *independent* metamers without compromising the practitioner's ability to conduct counterfactual reasoning. In future work, we plan to use the generated metamers in the context of music cognition to test Honing's hypothesis of multicomponent contour perception at various spectrotemporal scales.

## References

[1] H. Honing, “Unravelling the origins of musicality: Beyond music as an epiphenomenon of language,” *Behavioral and Brain Sciences*, vol. 44, 2021.

[2] J. Li, D. Baker, J. Burgoyne, and H. Honing, “Is pitch information indispensable for music recognition? A pilot study based on a musical matching pairs game,” in *Proceedings of the International Conference on Music Perception and Cognition (ICMPC)*, 2023.

[3] H. Honing, “Musicality as an upbeat to music: Introduction and research agenda,” in *The Origins of Musicality*, H. Honing, Ed. MIT Press, 2018, pp. 3–20.

[4] H. Honing, C. ten Cate, I. Peretz, and S. E. Trehub, “Without it no music: Cognition, biology and evolution of musicality,” p. 20140088, 2015.

[5] H. Honing, “Evidence for tempo-specific timing in music using a web-based experimental setup.” *Journal of Experimental Psychology: Human Perception and Performance*, vol. 32, no. 3, p. 780, 2006.

[6] J. Feather, A. Durango, R. Gonzalez, and J. McDermott, “Metamers of neural networks reveal divergence from human perceptual systems,” *Advances in Neural Information Processing Systems (NeurIPS)*, vol. 32, 2019.

[7] J. Feather, G. Leclerc, A. Mądry, and J. H. McDermott, “Model metamers reveal divergent invariances between biological and artificial neural networks,” *Nature Neuroscience*, vol. 26, no. 11, pp. 2017–2034, 2023.

[8] V. Lostanlen and F. Hecker, “The Shape of RemiXXXes to Come: Audio Texture Synthesis with Time–Frequency Scattering,” in *Proceedings of the International Conference on Digital Audio Effects (DAFx)*, 2019.

[9] S. Mallat, “Group invariant scattering,” *Communications in Pure and Applied Mathematics*, vol. 65, no. 10, pp. 1331–1398, 2012.

[10] I. Waldspurger, “Wavelet transform modulus: Phase retrieval and scattering,” Ph.D. dissertation, École normale supérieure, 2015.

[11] ——, “Exponential decay of scattering coefficients,” in *Proceedings of the International Conference on Sampling Theory and Applications (SampTA)*. IEEE, 2017, pp. 143–146.

[12] I. Sutskever, J. Martens, G. Dahl, and G. Hinton, “On the importance of initialization and momentum in deep learning,” in *Proceedings of the International Conference on Machine Learning (ICML)*, 2013, pp. 1139–1147.
[13] M. Andreux, T. Angles, G. Exarchakis, R. Leonarduzzi, G. Rochette, L. Thiry, J. Zarka, S. Mallat, J. Andén, E. Belilovsky et al., “Kymatio: Scattering transforms in python,” *Journal of Machine Learning Research*, vol. 21, no. 60, pp. 1–6, 2020.

[14] V. Lostanlen, G. Lafay, J. Andén, and M. Lagrange, “Relevance-based quantization of scattering features for unsupervised mining of environmental audio,” *EURASIP Journal on Audio, Speech, and Music Processing*, vol. 2018, no. 1, pp. 1–10, 2018.

[15] J. Andén and S. Mallat, “Multiscale scattering for audio classification.” in *Proceedings of the International Society for Music Information Retrieval Society (ISMIR) Conference*. Miami, Florida, 2011, pp. 657–662.

[16] J. Andén, V. Lostanlen, and S. Mallat, “Joint time-frequency scattering for audio classification,” in *Proceedings of the IEEE International Workshop on Machine Learning for Signal Processing (MLSP)*. IEEE, 2015, pp. 1–6.

[17] S. Mallat, “Understanding deep convolutional networks,” *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences*, vol. 374, no. 2065, p. 20150203, 2016.

[18] J. Andén, V. Lostanlen, and S. Mallat, “Joint time-frequency scattering,” *IEEE Transactions on Signal Processing*, vol. 67, no. 14, pp. 3704–3718, 2019.

[19] C. Wang, V. Lostanlen, E. Benetos, and E. Chew, “Playing technique recognition by joint time-frequency scattering,” in *Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)*. IEEE, 2020, pp. 881–885.

[20] C. Wang, E. Benetos, S. Wang, and E. Versace, “Joint scattering for automatic chick call recognition,” in *Proceedings of the European Signal Processing Conference (EUSIPCO)*. IEEE, 2022, pp. 195–199.

[21] G. Sharma, K. Umapathy, and S. Krishnan, “Time-frequency scattergrams for biomedical audio signal representation and classification,” *IEEE/ACM Transactions on Audio, Speech, and Language Processing*, vol. 32, pp. 564–576, 2023.

[22] V. Lostanlen, C. El-Hajj, M. Rossignol, G. Lafay, J. Andén, and M. Lagrange, “Time-frequency scattering accurately models auditory similarities between instrumental playing techniques,” *EURASIP Journal on Audio, Speech, and Music Processing*, vol. 2021, no. 1, p. 3, 2021.

[23] V. Lostanlen, J. Andén, and M. Lagrange, “Fourier at the heart of computer music: From harmonic sounds to texture,” *Comptes Rendus Physique*, vol. 20, no. 5, pp. 461–473, 2019.

[24] V. Lostanlen, “On time-frequency scattering and computer music,” in *Florian Hecker: Halluzination, Perspektive, Synthese*, V. J. Müller, Ed. Sternberg Press, 2018.

[25] V. Lostanlen and F. Hecker, “Synopsis Seriation: A computer music piece made with time-frequency scattering and information geometry,” in *Journées d’Informatique Musicale 2021*, 2021.

[26] I. Blom, “Sound Effect: Ina Blom on the art of Florian Hecker,” *Artforum International*, vol. 59, no. 6, pp. 100–107, 2021.
[27] A. Ireland and P. Gómez, "Al lado de una estatua fría," *Revista de Filosofía de la Universidad de Costa Rica*, vol. 63, no. 166, pp. 193-201, 2024.

[28] M. Andreux and S. Mallat, "Music generation and transformation with moment matching-scattering inverse networks," in *Proceedings of the International Conference on Music Information Retrieval (ISMIR)*, 2018, pp. 327-333.

[29] J. Muradeli, C. Vahidi, C. Wang, H. Han, V. Lostanlen, M. Lagrange, and G. Fazekas, "Differentiable time-frequency scattering on GPU," in *Proceedings of the International Conference on Digital Audio Effects (DAFx)*, 2022.

[30] C. Vahidi, H. Han, C. Wang, M. Lagrange, G. Fazekas, and V. Lostanlen, "Mesostructures: Beyond spectrogram loss in differentiable time-frequency analysis," *Journal of the Audio Engineering Society*, vol. 71, no. 9, pp. 577-585, 2023.

[31] H. Han, V. Lostanlen, and M. Lagrange, "Perceptual-neural-physical sound matching," in *Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)*. IEEE, 2023, pp. 1-5.

[32] —, "Learning to solve inverse problems for perceptual sound matching," *IEEE/ACM Transactions on Audio, Speech, and Language Processing*, 2024.

[33] G. Richard, V. Lostanlen, Y.-H. Yang, and M. Müller, "Model-based deep learning for music information research," *arXiv preprint arXiv:2406.11540*, 2024.

[34] T. Chi, Y. Gao, M. C. Guyton, P. Ru, and S. Shamma, "Spectro-temporal modulation transfer functions and speech intelligibility," *The Journal of the Acoustical Society of America*, vol. 106, no. 5, pp. 2719-2732, 1999.

[35] N. Kowalski, D. A. Depireux, and S. A. Shamma, "Analysis of dynamic spectra in ferret primary auditory cortex. i. characteristics of single-unit responses to moving ripple spectra," *Journal of Neurophysiology*, vol. 76, no. 5, pp. 3503-3523, 1996.

[36] T. Chi, P. Ru, and S. A. Shamma, "Multiresolution spectrotemporal analysis of complex sounds," *The Journal of the Acoustical Society of America*, vol. 118, no. 2, pp. 887-906, 2005.

[37] K. Patil, D. Pressnitzer, S. Shamma, and M. Elhilali, "Music in our ears: the biological bases of musical timbre perception," *PLOS Computational Biology*, vol. 8, no. 11, p. e1002759, 2012.

[38] T. Lindeberg and A. Friberg, "Idealized computational models for auditory receptive fields," *PLOS ONE*, vol. 10, no. 3, 2015.

[39] S. V. Norman-Haignere and J. H. McDermott, "Neural responses to natural and model-matched stimuli reveal distinct computations in primary and nonprimary auditory cortex," *PLOS Biology*, vol. 16, no. 12, p. e2005127, 2018.

[40] K. Siedenburg, I. Fujinaga, and S. McAdams, "A comparison of approaches to timbre descriptors in music information retrieval and music psychology," *Journal of New Music Research*, vol. 45, no. 1, pp. 27-41, 2016.

[41] D. Min, H. Nam, and Y.-H. Park, "Few-shot bioacoustic event detection utilizing spectro-temporal receptive field," in *Proceedings of the INTER-NOISE and NOISE-CON Congress and Conference*, vol. 270, no. 10. Institute of Noise Control Engineering, 2024, pp. 1694-1702.
[42] K. Li, K. Zaman, X. Li, M. Akagi, and M. Unoki, “Machine anomalous sound detection using spectral-temporal modulation representations derived from machine-specific filterbanks,” *arXiv preprint arXiv:2409.05319*, 2024.

[43] M. R. Schädler, B. T. Meyer, and B. Kollmeier, “Spectro-temporal modulation subspace-spanning filter bank features for robust automatic speech recognition,” *The Journal of the Acoustical Society of America*, vol. 131, no. 5, pp. 4134–4151, 2012.

[44] E. Thoret, P. Depalle, and S. McAdams, “Perceptually salient regions of the modulation power spectrum for musical instrument identification,” *Frontiers in Psychology*, vol. 8, p. 587, 2017.

<footnote>VINCENT LOSTANLEN AND HAN HAN, NOVEMBER 2024, Nantes Université, École Centrale Nantes, CNRS, LS2N, UMR 6004, F-44000 Nantes, France
*E-mail : vincent.lostanlen@cnrs.fr • Url : https://audio.ls2n.fr*</footnote>