texture

Sound texture statistics and synthesis after McDermott & Simoncelli (2011). Reached as so.texture.

sonore.texture.stats

Sound texture statistics (McDermott & Simoncelli, 2011).

A texture is summarized by time-averaged statistics of an auditory model: a cochlear filterbank, compressed and downsampled subband envelopes, and modulation filterbanks applied to those envelopes.

stats = so.texture.TextureStats.measure(rain)
stats.count()          # 1515 statistics, as in the paper
stats.save("rain.npz")

Differences from the MATLAB toolbox (v1.7) are listed in one place, DIFFERENCES_FROM_TOOLBOX. This is a clean-room implementation from the paper; the toolbox was consulted for behavior only.

All filtering here is circular (pad=0): the sound is treated as one period of a periodic signal, which is what makes synthesized textures loop seamlessly.

class TextureModel(fs: float = 20000.0, rms: float = 0.01, n_bands: int = 30, f_lo: float = 20.0, f_hi: float = 10000.0, compression: float = 0.3, env_fs: float = 400.0, n_mod: int = 20, mod_lo: float = 0.5, mod_hi: float = 200.0, mod_Q: float = 2.0, n_oct: int = 7, oct_hi: float = 100.0, c1_bands: tuple[int, ...] = (1, 2, 3, 4, 5, 6), c1_offsets: tuple[int, ...] = (1, 2), corr_offsets: tuple[int, ...] = (1, 2, 3, 5, 8, 11, 16, 21))

Bases: object

Parameters of the texture model. Defaults are the paper’s.

fs: float = 20000.0

rms: float = 0.01

n_bands: int = 30

f_lo: float = 20.0

f_hi: float = 10000.0

compression: float = 0.3

env_fs: float = 400.0

n_mod: int = 20

mod_lo: float = 0.5

mod_hi: float = 200.0

mod_Q: float = 2.0

n_oct: int = 7

oct_hi: float = 100.0

c1_bands: tuple[int, ...] = (1, 2, 3, 4, 5, 6)

c1_offsets: tuple[int, ...] = (1, 2)

corr_offsets: tuple[int, ...] = (1, 2, 3, 5, 8, 11, 16, 21)

property filterbank: Filterbank

property mod_bank: ConstantQModulationFilterbank

property oct_bank: OctaveModulationFilterbank

property decimation: int

prepare(sound: Sound) → ndarray

Mono, resampled to fs, truncated to a whole number of envelope samples, RMS-normalized. Returns a 1-D array.

subbands(x: ndarray) → ndarray

Circular cochlear subbands of a prepared signal, shape (n, n_bands + 2).

envelopes(subbands: ndarray) → ndarray

Compressed, downsampled envelopes, shape (n // decimation, n_bands + 2): Hilbert magnitude, raised to compression at the full rate, resampled (FFT, circular) to env_fs, clipped at 0.

class TextureStats(model: TextureModel, env_mean: ndarray, env_var: ndarray, env_skew: ndarray, env_kurt: ndarray, env_corr: ndarray, mod_power: ndarray, c1: ndarray, c2: ndarray, subband_var: ndarray, duration: float = 0.0)

Bases: View

The statistics of one texture. Arrays are indexed by cochlear channel first (n_bands + 2 channels, including the lowpass and highpass edges, as in the paper’s count of 1515 statistics).

env_mean, env_var, env_skew, env_kurt

Weighted moments of the compressed envelopes; env_var is var/mean^2.

Type:

(B,)

env_corr

env_corr[j, i] = correlation of channels j and j + corr_offsets[i].

Type:

(B, len(corr_offsets))

mod_power

Modulation power in each constant-Q band, relative to envelope variance.

Type:

(B, n_mod)

c1

Correlation of channels j and j + d within an octave modulation band.

Type:

(B, len(c1_bands), len(c1_offsets))

c2

Correlation of each octave modulation band (frequency-doubled) with the next one up in the same channel; real and imaginary parts are the in-phase and quadrature components.

Type:

(B, n_oct - 1), complex

subband_var

Variance of each cochlear subband (used to set levels in synthesis).

Type:

(B,)

discards = "TextureStats keeps only time-averaged statistics of the envelopes and their modulations: it discards the phase, the fine structure and every event's timing."

back_to_sound = 'sonore.texture.synth.synthesize draws a new sound with these statistics, which is not the measured one.'

no_plot = 'TextureStats has no plot: it holds several kinds of statistic with no single picture, so plot the arrays you need (env_mean, mod_power, c1, ...).'

model: TextureModel

env_mean: ndarray

env_var: ndarray

env_skew: ndarray

env_kurt: ndarray

env_corr: ndarray

mod_power: ndarray

c1: ndarray

c2: ndarray

subband_var: ndarray

duration: float = 0.0

classmethod measure(sound: Sound, model: TextureModel | None = None, window: str = 'ramped') → TextureStats

Measure the statistics of sound. window="ramped" (the default, for recorded originals) downweights the ends with measurement_window(); "uniform" weights all samples equally (as for circular synthetic signals).

classmethod from_subbands(sb: ndarray, model: TextureModel, window: str = 'ramped') → TextureStats

classmethod from_envelopes(env: ndarray, subband_var: ndarray, model: TextureModel, window: str = 'uniform') → TextureStats

Statistics of compressed, downsampled envelopes (n_env, B) (subband variances are passed through). Used during synthesis.

get(name: str) → ndarray

count(classes=('env_mean', 'env_var', 'env_skew', 'env_kurt', 'env_corr', 'mod_power', 'c1', 'c2')) → int

Number of statistics (defined entries; a complex C2 value counts as one, as in the paper). The default counts the paper’s classes: 1515.

replace(**changes) → TextureStats

A copy with some classes replaced, e.g. for hybrid textures: a.replace(mod_power=b.mod_power).

channel_mask(range_db: float = 30.0) → ndarray

Channels whose subband variance is within range_db of the loudest. Quieter channels are ignored by snr() (as in the toolbox): their statistics are dominated by noise and inaudible.

snr(other: TextureStats, classes=('env_mean', 'env_var', 'env_skew', 'env_kurt', 'env_corr', 'mod_power', 'c1', 'c2'), range_db: float = 30.0) → dict[str, float]

How well other matches these (target) statistics, per class: 10*log10(sum |target|**2 / sum |target - other|**2) in dB, over channels within range_db of the loudest. Pairwise classes (C, C1) count a pair only if both channels qualify.

The SNR is relative to the target’s own magnitude, so classes whose targets are near zero (C, C1, C2 and skew of noise-like textures) score low even between two samples of the same texture: sampling fluctuation dominates the ratio.

save(path) → None

classmethod load(path) → TextureStats

measurement_window(n: int, n_seconds: int) → ndarray

The toolbox’s window for measuring an original: flat, with a raised-cosine ramp of n // (n_seconds + 1) samples at each end (so about duration / (seconds + 1)). Normalized to sum to 1.

sonore.texture.grad

Texture statistics of one envelope channel, with analytic gradients.

Synthesis (McDermott & Simoncelli, 2011) adjusts one cochlear channel’s compressed, downsampled envelope s at a time, holding the others fixed. Each function here computes one class of statistics as a function of s and returns (value, vjp): the statistics, exactly as sonore.texture.TextureStats.measure() defines them, and a vector-Jacobian product vjp(g) = J(s).T @ g, so the gradient of 0.5 * ||value - target||**2 is vjp(value - target).

Every gradient is checked against finite differences in the tests.

Notation: w is the measurement window (sums to 1), d = s - w@s, m_p = w @ d**p. The modulation filters are real, zero-phase and circular, so they are symmetric operators (their own adjoints). The analytic filter A (the octave bank with negative frequencies removed) is Hermitian, so the adjoint of s -> Re(A s) is g -> Re(A g) and of s -> Im(A s) is g -> -Im(A g).

class ChannelContext(model: TextureModel, w: ndarray, H_mod: ndarray, A_oct: ndarray)

Bases: object

Everything that depends only on the model and the envelope length: the window and the filter responses on the FFT grid.

model: TextureModel

w: ndarray

H_mod: ndarray

A_oct: ndarray

classmethod build(model: TextureModel, n: int, w: ndarray | None = None) → ChannelContext

property n: int

mod_filter(x: ndarray) → ndarray

(n,) -> (n, n_mod); symmetric, so also its own adjoint.

mod_adjoint(G: ndarray) → ndarray

Adjoint of mod_filter(): (n, n_mod) -> (n,).

analytic(x: ndarray, bands) → ndarray

(n,) -> (n, len(bands)) complex analytic octave bands.

analytic_many(X: ndarray, bands) → ndarray

(n, k) -> (n, len(bands), k): analytic() of each column.

analytic_adjoint(g_re: ndarray, g_im: ndarray, bands) → ndarray

Adjoint of x -> (Re A x, Im A x) for the given bands, summed over bands.

Since A is Hermitian, this is Re(A g_re) - Im(A g_im) = Re(A (g_re + i g_im)): one complex FFT, and the sum over bands is taken before the inverse transform.

env_moments(s: ndarray, context: ChannelContext)

[mean, var/mean**2, skew, kurtosis] of the envelope.

mod_power(s: ndarray, context: ChannelContext)

Modulation power in each constant-Q band, relative to envelope variance.

env_corr(s: ndarray, others: ndarray, context: ChannelContext)

Correlation of s with each column of others (n, k) (fixed envelopes of other channels).

c1(s: ndarray, others: ndarray, context: ChannelContext, bands=None)

C1: correlation of s’s octave modulation bands with the same bands of each column of others (n, k). Shape (len(bands), k). No mean subtraction (paper Eq. 6).

c2(s: ndarray, context: ChannelContext)

C2 within one channel: correlation of each octave band, frequency-doubled, with the next band up. Complex (n_oct - 1,): real part against the band’s real part, imaginary part against its imaginary (quadrature) part.

mod_power_core(B: ndarray, d: ndarray, w: ndarray)

mod_power() from the filtered envelope B (n, M) and the centered envelope d. vjp(g) returns (G, direct): the cotangent of B (to be passed through the filter adjoint) and the gradient term that reaches s directly through the variance.

c1_core(R: ndarray, Ro: ndarray, w: ndarray)

c1() from the octave bands of s (R, (n, K)) and of the neighbors (Ro, (n, K, k)). vjp(g) returns the cotangent of R.

c2_core(A: ndarray, w: ndarray)

c2() from all analytic octave bands A (n, K). vjp(g) returns the cotangents (g_re, g_im) of A.real and A.imag.

Each band pair is a lower band (bands 0..K-2) and the band above it.

sonore.texture.synth

Imposing texture statistics on envelopes (McDermott & Simoncelli, 2011).

Milestone 3 of the synthesis: adjust one channel’s compressed envelope so its statistics approach the targets, holding the other channels fixed. The full loop (channel ordering, fine structure, subband rescaling, iterations) builds on impose_channel().

The objective for channel k is the unweighted sum of squared errors of the statistics that depend on k:

• its envelope moments, modulation power, and C2;

• its envelope correlation C and C1 with channels already adjusted in this pass (adjusted), at the model’s offsets in both directions.

It is minimized by nonlinear conjugate gradient (SciPy’s Polak-Ribiere CG) with a fixed number of iterations, then the envelope is clipped at 0.

The loss is nearly flat along two directions: a uniform shift of the envelope and a scaling of its deviations from the mean. Only the mean and variance depend on them, so a few CG iterations barely correct those two. Since every other statistic is invariant to shifting and scaling, they are set exactly after CG (see affine in impose_channel()); this is a deviation from the toolbox, listed in sonore.texture.DIFFERENCES_FROM_TOOLBOX.

class ChannelObjective(context: ChannelContext, target: TextureStats, k: int, env: ndarray, adjusted: ndarray, classes: tuple[str, ...] = ('env_mean', 'env_var', 'env_skew', 'env_kurt', 'env_corr', 'mod_power', 'c1', 'c2'))

Bases: object

The squared-error objective for channel k given the other channels’ envelopes env (n, B) and which of them count as adjusted.

context: ChannelContext

target: TextureStats

k: int

env: ndarray

adjusted: ndarray

classes: tuple[str, ...] = ('env_mean', 'env_var', 'env_skew', 'env_kurt', 'env_corr', 'mod_power', 'c1', 'c2')

terms()

[(name, f(s) -> (value, vjp), target)] for this channel.

reference(s: ndarray, per_term: bool = False)

The objective computed class by class from terms() (the straightforward, separately tested path; __call__() must agree).

impose_channel(target: TextureStats, env: ndarray, k: int, adjusted: ndarray | None = None, n_iter: int = 5, context: ChannelContext | None = None, classes=('env_mean', 'env_var', 'env_skew', 'env_kurt', 'env_corr', 'mod_power', 'c1', 'c2'), affine: bool = True) → tuple[ndarray, dict]

Adjust channel k of env (n, B) toward target.

adjusted marks channels whose envelopes are already final in this pass (default: none). Runs n_iter conjugate-gradient iterations and clips the result at 0. With affine (default), CG is followed by an exact affine correction s -> a + b*(s - mean) that sets the envelope mean and variance to their targets. Every other statistic is invariant to that map (skew, kurtosis and all correlations are scale- and shift-invariant; modulation power is normalized by the variance and its filters reject DC), so the correction costs nothing elsewhere, apart from clipping. It is needed because those two directions are nearly flat in the objective. Returns the new envelope (n,) and a report with the objective before and after (total and per term).

synthesize(target: TextureStats, duration: float = 5.0, classes=('env_mean', 'env_var', 'env_skew', 'env_kurt', 'env_corr', 'mod_power', 'c1', 'c2'), rng=None, max_iter: int = 60, stop_db: float = 30.0, converged_db: float = 20.0, init: Sound | None = None, callback=None, progress: bool = False) → tuple[Sound, dict]

Synthesize a sound whose statistics match target.

Starts from Gaussian noise (or init) and iterates: decompose into subbands; take compressed envelopes (downsampled, plus the high-rate residual the downsampling removes) and fine structure; impose the statistics channel by channel (impose_channel(), in channel_order(), each channel’s correlations measured against the channels already adjusted); restore the full-rate envelopes, decompress, recombine with the fine structure; rescale each subband to its target variance; resynthesize. All analysis is circular, so the result loops seamlessly.

Stops when every class in classes is at least stop_db SNR, or after max_iter iterations. The synthesis counts as converged if the average SNR over classes is at least converged_db.

Run time. Cost is linear in duration and in the number of iterations: about 0.4 s per iteration per second of sound on one core (2 s per iteration for 5 s), so the default 60 iterations of 5 s take about 2 minutes. Most textures are close to their final quality by 20-30 iterations; use max_iter to trade quality for time, progress=True to print one line per iteration, or callback(iteration, sound, snr) to watch or stop from your own code.

Returns the sound (at target.model.fs, RMS target.model.rms) and a report: snr (per-iteration dicts), converged, iterations, best_iteration. The returned sound is the iterate with the best average SNR.

channel_order(env_mean: ndarray) → list[int]

Start at the channel with the largest target envelope mean and alternate outward: k, k+1, k-1, k+2, k-2, ....