Audio Signal Analysis, Indexing and Transformations

Contents

• Introduction and overview
• 0. Review: signals, Fourier analysis, filters and random signals
• 1. Time–frequency representations and audio features
• 2. Audio indexing I: classifying genres and instruments
• 3. Audio indexing II: onsets, tempo, beats and downbeats
• 4. Digital audio effects
• 5. Time-scale, pitch and timbre modifications
• 6. Spatial audio: binaural, transaural and wave field synthesis
• 7. Source separation I: linear and convolutive mixtures
• 8. Source separation II: underdetermined and sparse models, DUET
• 9. High-resolution analysis I: sinusoidal models beyond Fourier
• 10. High-resolution analysis II: MUSIC, ESPRIT and model order
• 11. Complements for self-study
• 12. Research directions and open problems
• Exercises
• Projects and assessment
• References

Introduction and overview

Audience: a student with a good M1 in mathematics, computer science or physics, who knows the Fourier transform and linear filters, and wants to analyse, index, separate and transform real audio signals.

Empirical first: two microphones record three people talking at once. No linear demixing matrix can undo a mixture with more sources than sensors, yet the DUET algorithm of ch. 8 raises the signal-to-distortion ratios of the three voices from −4.6, −3.1 and −1.6 dB to 14.7, 10.7 and 14.9 dB.

In the same spirit, two tones closer than the Fourier resolution of a 64-sample frame give a periodogram with two peaks at the wrong frequencies, while ESPRIT (ch. 10) finds both within 0.0003 cycles per sample.

This course covers four families of audio problems. Indexing, or machine listening, extracts descriptors from a recording to recognise its genre, its instruments, its tempo and its bars. Transformations change a sound: digital effects, reverberation, time stretching, pitch shifting, and 3D rendering over headphones or loudspeaker arrays.

Source separation recovers individual sources from a mixture. High-resolution spectral analysis estimates the frequencies of a sum of sinusoids beyond the Fourier limit. The course is taught by Gaël Richard and Roland Badeau (Télécom Paris, per master-mva.com).

It shares tools with Audio Signal Processing: Time–Frequency Analysis, Algorithms for Speech and Natural Language Processing and Machine Learning for Time Series.

What you will be able to do

• Compute and invert a short-time Fourier transform, and choose its window, length and hop for a given audio task.
• Extract MFCC, chroma and spectral descriptors, and build a genre or instrument classifier evaluated without session or artist leakage.
• Detect onsets, estimate the tempo, track beats by dynamic programming and locate downbeats.
• Design flangers, phasers, distortions and artificial reverberators, and predict their frequency responses and decay times.
• Stretch a sound in time and shift its pitch with a phase vocoder, while preserving its spectral envelope.
• Explain how HRTFs, crosstalk cancellation and wave field synthesis render a sound in space, and compute their limits.
• Separate determined linear mixtures by independent component analysis, and explain the permutation problem of convolutive mixtures.
• Separate more sources than microphones with sparse time–frequency models and DUET.
• Estimate the frequencies, damping factors and amplitudes of sinusoids with Prony, MUSIC and ESPRIT, and select the model order.
• Read a research paper of the field, reimplement its method and present it critically, as the official assessment asks.

Official objectives and syllabus (2026–2027)

The course description on master-mva.com states that the course aims to span several domains of audio signal analysis. It names audio indexing (or machine listening), high-resolution audio spectral analysis, audio source separation, and audio transformations: 3D sound rendering, sound effects and sound modifications.

It announces 13.5 h of theoretical lectures and 10.5 h of practical sessions in Matlab, or Python if preferred; this page uses Python throughout. Assessment is by papers reading and analysis, with a written report and an oral presentation. The Université Paris-Saclay course detail gives only the assessment type, and no programme of its own.

The table has one row per item of the “Thèmes abordés” list, split at its listed sub-topics, and one row per announced organisation item. The list ends with an item “Les intervenants” (the teachers), a heading rather than a topic, which is why it has no row.

Every official topic has a section in these notes. The extensions beyond the syllabus are the review of ch. 0, the self-study complements of ch. 11 (EM, feedback delay networks, linear prediction, STFT inversion, amplitude panning and Ambisonics, NMF, harmonic–percussive separation, the Cramér–Rao bound, neural models) and the research directions of ch. 12.

Where this course fits in the programme

Official prerequisite: basics of signal processing. In practice the course also uses linear algebra (eigendecompositions, least squares), elementary probability (Gaussian models, maximum likelihood) and some Python; ch. 0 reviews the signal-processing part.

• Refreshers — Linear Algebra, Analysis and Optimization, Probability and Statistics and Scientific Python.
• Semester 1 courses that prepare it — Deep Learning and Signal Processing and Machine Learning for Time Series (spectral analysis, filtering); Introduction to Probabilistic Graphical Models (EM, mixtures, hidden Markov models); Convex Optimization (the ℓ1 problems of ch. 8).
• Related courses of semester 2 — Audio Signal Processing: Time–Frequency Analysis, Speech and Natural Language Processing, Kernel Methods (the SVM of ch. 2), Large Random Matrices (eigenvalues of sample covariances, ch. 9–10), Detection Theory and Brain Imaging and BCI, which also uses ICA.
• Afterwards — the Research Internship, for which the paper-study format of the assessment is a rehearsal.

Course outline

• 0. Review — DFT, convolution, filters, windows and leakage, stationary random signals (self-study before session 1).
• 1. Time–frequency representations and audio features — STFT and its inversion, mel and constant-Q scales, MFCC, chroma, spectral descriptors.
• 2. Audio indexing I — the bag-of-frames pipeline, Gaussian mixtures and SVMs, genre and instrument recognition, honest evaluation.
• 3. Audio indexing II — onset detection, tempo, beat tracking by dynamic programming, downbeats.
• 4. Digital audio effects — comb filters, flanger, phaser, distortion, artificial reverberation.
• 5. Time-scale, pitch and timbre modifications — overlap-add, the phase vocoder, pitch shifting, spectral envelopes.
• 6. Spatial audio — interaural cues, HRTFs, binaural and transaural rendering, wave field synthesis.
• 7. Source separation I — mixing models, independent component analysis, convolutive mixtures, evaluation.
• 8. Source separation II — time–frequency sparsity, ℓ1 recovery, DUET, Wiener masks.
• 9. High-resolution analysis I — sinusoidal models, the Fourier limit, Prony, the structure of Hankel and covariance matrices.
• 10. High-resolution analysis II — MUSIC, ESPRIT, the ESTER order criterion.
• 11. Complements for self-study — EM, feedback delay networks, linear prediction, STFT inversion, panning and Ambisonics, NMF, harmonic–percussive separation, the Cramér–Rao bound, neural models.
• 12. Research directions and open problems — a self-study map of current work.

Session plan

24 h in class = 8 sessions of 3 h. The 24 h come from the Université Paris-Saclay programme table (14 h of lectures, 10 h of labs). master-mva.com announces the same total, split into 13.5 h of lectures and 10.5 h of practical sessions; the plan follows that split, with one 3 h lecture and seven sessions of 1.5 h of lecture plus 1.5 h of lab. The split into sessions is a proposal.

Level 1 = foundations, 2 = core methods, 3 = advanced and research-level material. The level never goes down from one session to the next.

Self-study, outside the 24 h (about 10 h in total): ch. 0 Review of the signal-processing prerequisites, before session 1 (≈ 3 h, with Ex. 1) · ch. 11 Complements for self-study (≈ 4 h, with Ex. 18) · ch. 12 Research directions and open problems (≈ 3 h, with Ex. 19).

How to study this module

• Work through ch. 0 before session 1: it reviews the official prerequisite (basics of signal processing), and every later chapter uses the DFT, windows and the z-transform without warning.
• After each session, do the exercises listed in its row, then run the code of the chapter and change one parameter (window length, hop, delay, number of sources, signal-to-noise ratio).
• Listen to what you compute: write the signals of the code to WAV files (scipy.io.wavfile.write) and compare them by ear.
• The code uses NumPy, SciPy and scikit-learn on synthetic signals; every block runs in a few seconds on a laptop.
• Start a project of the last section after session 4; the official assessment is a paper study with a report and a talk.

0. Review: signals, Fourier analysis, filters and random signals

Every chapter of this course computes spectra, filters signals and reasons about noise. This chapter collects the facts it uses, with the two that matter most for audio: the Fourier resolution limit and the leakage of windows.

The discrete Fourier transform

A digital audio signal is a sequence x(n) sampled at f_s hertz, typically 16 kHz for speech and 44.1 kHz for music. A normalised angular frequency \omega (radians per sample) corresponds to f = \omega f_s / 2\pi hertz.

By the sampling theorem, a continuous signal whose spectrum vanishes above f_s / 2 is determined by its samples. Frequencies above f_s / 2 fold back, or alias, which matters for the distortion effects of ch. 4.

Linear filters

Windows, leakage and resolution

A frame of N samples is the product of the signal with a window w supported on \{0, \dots, N - 1\}. By Proposition 0.2, the spectrum of the frame is the spectrum of the signal convolved with that of the window.

Two consequences organise the rest of the course.

• Resolution. Two exponentials closer than about one main-lobe half-width merge into one peak. This is the Rayleigh limit: a frequency gap of f_s / N hertz for the rectangle. Zero-padding interpolates the DTFT, but it does not resolve anything.
• Leakage. A strong component spreads energy far away through its sidelobes and can hide a weak one. Windows trade main-lobe width (resolution) for sidelobe level (dynamic range).

The code below adds a tone 50 dB below a strong one, 150 Hz away, in a 50 ms frame.

import numpy as np

fs, N = 8000, 400                         # a 50 ms analysis window, fs/N = 20 Hz
n = np.arange(N)
strong = np.cos(2 * np.pi * 1000 * n / fs)
weak = 10 ** (-50 / 20) * np.cos(2 * np.pi * 1150 * n / fs)   # 50 dB below

def level_db(x, w, f0, pad=16):
    X = np.abs(np.fft.rfft(x * w, pad * N))   # zero padding interpolates only
    f = np.fft.rfftfreq(pad * N, 1 / fs)
    X_db = 20 * np.log10(X / X.max())
    return X_db[np.argmin(np.abs(f - f0))]

for name, w in [("rectangular", np.ones(N)), ("Hann", np.hanning(N))]:
    without = level_db(strong, w, 1150)
    with_ = level_db(strong + weak, w, 1150)
    print(f"{name:12s} level at 1150 Hz: {without:6.1f} dB without the weak tone,"
          f" {with_:6.1f} dB with it")

With the rectangular window the level at 1150 Hz is −27.0 dB without the weak tone and −26.9 dB with it: the weak tone stays invisible under a sidelobe. With the Hann window it rises from −62.2 dB to −49.9 dB: the weak tone appears at its true level.

Random signals

Filtering a stationary signal by H multiplies its spectral density by |H(e^{i\omega})|^2. The covariance matrix \mathbf{R} of L consecutive samples is Hermitian, Toeplitz and positive semidefinite.

Chapters 7, 9 and 10 also use two facts from linear algebra. A Hermitian matrix has an orthonormal basis of eigenvectors with real eigenvalues. Any matrix has a singular value decomposition \mathbf{X} = \mathbf{U} \boldsymbol{\Sigma} \mathbf{V}^H, and its first K left singular vectors span the column space of its best rank-K approximation.

1. Time–frequency representations and audio features

Audio signals are not stationary: notes start and stop, voices glide. This chapter builds the short-time Fourier transform, proves when it can be inverted, and derives from it the descriptors used for indexing in ch. 2–3.

The short-time Fourier transform

The window length sets a trade-off. At 16 kHz, N = 1024 gives frames of 64 ms and bins 15.6 Hz apart: good for pitch, poor for drum attacks. Speech analysis typically uses frames of 20–40 ms and music analysis 40–100 ms.

For the periodic Hann window and R = N/4, the shifted windows sum to 2 and their squares to 3/2. So w_s = \tfrac{2}{3} w reconstructs exactly; the code below normalises by \sum_m w^2(n - mR), which is the same thing.

A modified STFT (masked in ch. 8, re-phased in ch. 5) is generally not the STFT of any signal; ch. 11 explains what overlap-add then returns.

Perceptual scales

The ear resolves low frequencies finely and high frequencies coarsely: its critical bands are about 100 Hz wide below 500 Hz and about 20 % of the centre frequency above, roughly a third of an octave. Audio descriptors imitate this with nonlinear frequency scales.

The constant-Q transform follows the same idea with geometrically spaced bins, f_k = f_{\min}\, 2^{k / b} with b bins per octave, so that every bin spans the same musical interval.

The cosine transform of a logarithm comes from the source–filter view of sound production.

Keeping the low quefrencies (liftering) keeps the envelope, which carries timbre and vowel identity, and discards pitch. MFCCs do the same on the mel scale.

Descriptors and chroma

The code computes an STFT, checks the reconstruction of Proposition 1.2, then extracts MFCCs and chroma from four notes a semitone apart, starting on A3 (220 Hz).

import numpy as np
from scipy.fft import dct

fs, N, R = 16000, 1024, 256                  # window 64 ms, hop 16 ms
w = np.hanning(N + 1)[:N]                    # periodic Hann window
print("COLA sum (hop N/4):", np.round(sum(np.roll(w, k * R) for k in range(N // R))[:3], 6))

def stft(x):
    frames = [x[i:i + N] * w for i in range(0, len(x) - N + 1, R)]
    return np.fft.rfft(np.array(frames), axis=1)            # shape (M, N/2 + 1)

def istft(X, length):
    y, norm = np.zeros(length), np.zeros(length)
    for m, frame in enumerate(np.fft.irfft(X, n=N, axis=1)):
        y[m * R:m * R + N] += frame * w                      # synthesis window
        norm[m * R:m * R + N] += w ** 2
    return y / np.maximum(norm, 1e-12)

def mel(f):  return 2595 * np.log10(1 + f / 700)
def imel(m): return 700 * (10 ** (m / 2595) - 1)

def mel_filterbank(n_mels=40, fmin=0.0, fmax=fs / 2):
    edges = imel(np.linspace(mel(fmin), mel(fmax), n_mels + 2))
    f = np.fft.rfftfreq(N, 1 / fs)
    H = np.zeros((n_mels, len(f)))
    for j in range(n_mels):
        lo, c, hi = edges[j], edges[j + 1], edges[j + 2]
        H[j] = np.clip(np.minimum((f - lo) / (c - lo), (hi - f) / (hi - c)), 0, None)
    return H

rng = np.random.default_rng(0)
t = np.arange(2 * fs) / fs
f0 = 220 * 2 ** (np.floor(t / 0.5) / 12)                     # four notes, a semitone apart
x = sum(np.sin(2 * np.pi * k * np.cumsum(f0) / fs) / k for k in range(1, 9))
x += 0.01 * rng.standard_normal(len(x))

X = stft(x)
y = istft(X, len(x))
inner = slice(N, len(x) - N)
print("reconstruction error:", f"{np.max(np.abs(x[inner] - y[inner])):.1e}")

logmel = np.log(np.abs(X) ** 2 @ mel_filterbank().T + 1e-10)    # (M, 40)
mfcc = dct(logmel, type=2, norm="ortho", axis=1)[:, :13]
print("frames:", X.shape[0], "| MFCC shape:", mfcc.shape)
f = np.fft.rfftfreq(N, 1 / fs)
pc = np.round(12 * np.log2(np.maximum(f, 1) / 440)) % 12      # pitch class, 0 = A
keep = (f > 50) & (f < 2000)
chroma = np.stack([(np.abs(X[:, keep]) ** 2)[:, pc[keep] == p].sum(1) for p in range(12)], 1)
names = "A A# B C C# D D# E F F# G G#".split()
centres = ((np.arange(4) * 0.5 + 0.25) * fs - N / 2) // R   # one frame inside each note
print("dominant pitch class per note:", [names[np.argmax(chroma[int(m)])] for m in centres])

The shifted windows sum to 2, the reconstruction error is 1.1 \times 10^{-15} (rounding only), the 2 s signal gives 122 frames of 13 MFCCs, and the chroma follows the notes A, A#, B, C.

2. Audio indexing I: classifying genres and instruments

Content-based indexing labels a recording from its signal alone: its genre, the instruments playing, the sound events it contains. This chapter builds the classical pipeline on the features of ch. 1 and shows how easily its evaluation can be fooled.

The bag-of-frames pipeline

It ignores the order of frames, hence its name: it suits timbre-driven tasks (instrument, genre, acoustic scene), not rhythm, which ch. 3 treats separately.

Gaussian mixture models

The independence assumption is false (consecutive frames overlap), which makes the posterior probabilities overconfident, but the argmax is often still good.

The class densities are fitted to the training frames of each class by the EM algorithm, whose updates and monotonicity are proved in ch. 11.

Discriminative classifiers and feature selection

Genre recognition and instrument identification

Evaluating an indexing system

The code builds three synthetic instruments recorded in 18 sessions (one instrument per session, each with its own random equaliser), describes each note by the mean and standard deviation of 13 MFCCs, and compares a random split with a session-disjoint split.

import numpy as np
from scipy.fft import dct
from sklearn.mixture import GaussianMixture
from sklearn.svm import SVC
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import StratifiedKFold, StratifiedGroupKFold

fs, N, R = 16000, 512, 256
rng = np.random.default_rng(0)
w = np.hanning(N + 1)[:N]
f = np.fft.rfftfreq(N, 1 / fs)
mel = lambda x: 2595 * np.log10(1 + x / 700)
edges = 700 * (10 ** (np.linspace(0, mel(fs / 2), 42) / 2595) - 1)
H = np.array([np.clip(np.minimum((f - a) / (b - a), (c - f) / (c - b)), 0, None)
              for a, b, c in zip(edges[:-2], edges[1:-1], edges[2:])])   # 40 mel filters

def note(kind, f0, eq):
    """Synthetic note: harmonic amplitudes and envelope depend on the instrument."""
    t = np.arange(int(0.5 * fs)) / fs
    k = np.arange(1, 21)
    amp = [np.where(k % 2 == 1, 1 / k, 0.3 / k),       # 'clarinet': weak even harmonics
           1 / k,                                       # 'bowed string': all harmonics
           1 / k ** 1.5][kind]                          # 'flute': steeper spectral slope
    amp = amp * eq(k * f0) * (k * f0 < fs / 2)          # recording-session colouring
    attack = (0.06, 0.12, 0.04)[kind] * np.exp(0.3 * rng.standard_normal())
    env = np.minimum(t / attack, 1) * np.exp(-t * (1.0, 0.5, 2.0)[kind])
    x = env * (amp[:, None] * np.sin(2 * np.pi * np.outer(k, f0 * t)
                                     + rng.uniform(0, 2 * np.pi, (20, 1)))).sum(0)
    return x + (0.002, 0.01, 0.03)[kind] * env * rng.standard_normal(len(t))

def features(x):
    frames = np.array([x[i:i + N] * w for i in range(0, len(x) - N, R)])
    P = np.abs(np.fft.rfft(frames, axis=1)) ** 2
    mfcc = dct(np.log(P @ H.T + 1e-10), norm="ortho", axis=1)[:, 1:14]
    return np.concatenate([mfcc.mean(0), mfcc.std(0)])  # temporal integration: 26 numbers

X, y, groups = [], [], []
for s in range(18):                                     # 18 sessions, one instrument each
    kind = s % 3
    a, ph = rng.standard_normal(4), rng.uniform(0, 2 * np.pi, 4)
    eq = lambda fr, a=a, ph=ph: np.exp(0.4 * sum(a[j] * np.cos((j + 1) * np.log2(fr / 100) + ph[j])
                                                 for j in range(4)))
    for _ in range(20):
        f0 = 196 * 2 ** (rng.integers(0, 19) / 12)      # G3 to C#5
        X.append(features(note(kind, f0, eq))); y.append(kind); groups.append(s)
X, y, groups = np.array(X), np.array(y), np.array(groups)

def gmm_predict(Xtr, ytr, Xte, k=2):                    # one GMM per class + Bayes rule
    sc = StandardScaler().fit(Xtr)
    L = [GaussianMixture(k, covariance_type="diag", random_state=0)
         .fit(sc.transform(Xtr[ytr == c])).score_samples(sc.transform(Xte)) for c in range(3)]
    return np.argmax(L, axis=0)                         # equal priors

def svm_predict(Xtr, ytr, Xte):
    return make_pipeline(StandardScaler(), SVC(C=10)).fit(Xtr, ytr).predict(Xte)

splits = {"random split": StratifiedKFold(3, shuffle=True, random_state=0).split(X, y),
          "session-disjoint": StratifiedGroupKFold(3, shuffle=True, random_state=0).split(X, y, groups)}
for split, folds in splits.items():
    folds = list(folds)
    for name, predict in [("GMM", gmm_predict), ("SVM", svm_predict)]:
        acc = [np.mean(predict(X[tr], y[tr], X[te]) == y[te]) for tr, te in folds]
        print(f"{split:16s} {name}: accuracy {np.mean(acc):.3f}")

The random split reports 0.828 for the GMM and 0.983 for the SVM. On sessions never seen in training, the accuracies drop to 0.758 and 0.925: the gap is the part of the score that recognised sessions rather than instruments.

3. Audio indexing II: onsets, tempo, beats and downbeats

Rhythm is the temporal side of indexing: where notes start, how fast the pulse is, where each beat falls and where each bar begins. The chain goes from an onset detection function to a tempo, then to beat times and downbeats.

Onset detection

The logarithm makes soft and loud notes comparable, and the half-wave rectification ignores decays. Bello et al. (2005) compare this function with energy-, phase- and complex-domain variants.

Tempo estimation

Hence the prior \pi. Comb-filter banks (Scheirer, 1998) and Fourier tempograms (Grosche and Müller, 2011) replace the autocorrelation in other systems and face the same ambiguity.

Beat tracking by dynamic programming

The restriction to \Delta \in [\tau/2, 2\tau] is part of the model: a slowly changing tempo is followed, but a tempo change by more than a factor 2 is not.

Downbeat estimation

Probabilistic models make the bar position a hidden state that advances with each beat, decoded by the Viterbi algorithm (Klapuri, Eronen and Astola, 2006). Current systems learn the salience with neural networks and decode beats and downbeats jointly (Böck, Krebs and Widmer, 2016; Durand et al., 2017).

Evaluation

A detected beat is correct if it lies within ±70 ms of an annotated one, and precision, recall and their harmonic mean, the F-measure, follow. Tempo accuracy is reported with a 4 % tolerance, with or without accepting factors 2 and 3. The mir_eval library (Raffel et al., 2014) implements these measures.

The code synthesises 10 bars of a drum pattern at 128 BPM with 5 ms timing jitter (kick drum on beat 1 and between beats 2 and 3, snare on beats 2 and 4, hi-hat on eighth notes), then runs the whole chain.

import numpy as np

fs, N, R = 16000, 1024, 160                       # ODF frame rate fs/R = 100 Hz
rng = np.random.default_rng(1)
bpm, n_bars = 128, 10
beat = 60 / bpm
t = np.arange(int((4 * n_bars + 1) * beat * fs)) / fs
x = 0.003 * rng.standard_normal(len(t))

def hit(t0, kind):
    i = int(t0 * fs); L = int(0.25 * fs); tt = np.arange(L) / fs
    if kind == "kick":  s = np.sin(2 * np.pi * (50 + 80 * np.exp(-tt / 0.03)) * tt) * np.exp(-tt / 0.12)
    if kind == "snare": s = 0.5 * rng.standard_normal(L) * np.exp(-tt / 0.06)
    if kind == "hat":   s = 0.15 * np.diff(rng.standard_normal(L + 1)) * np.exp(-tt / 0.02)
    x[i:i + L] += s[:len(x) - i]

beats = []
for b in range(4 * n_bars):
    tb = b * beat + 0.005 * rng.standard_normal()   # 5 ms human timing jitter
    beats.append(tb)
    if b % 4 == 0: hit(tb, "kick")                  # kick on beat 1 ...
    if b % 4 == 1: hit(tb + beat / 2, "kick")       # ... and on the "and" of 2
    if b % 2 == 1: hit(tb, "snare")                 # snare on 2 and 4
    hit(tb, "hat"); hit(tb + beat / 2, "hat")       # eighth-note hi-hats
beats = np.array(beats)

# 1. Onset detection function: half-wave rectified log-spectral flux
w = np.hanning(N)
S = np.abs(np.array([np.fft.rfft(x[i:i + N] * w) for i in range(0, len(x) - N, R)]))
L = np.log1p(100 * S)
odf = np.maximum(np.diff(L, axis=0), 0).sum(1)
odf = np.concatenate([[0], (odf - odf.mean()) / odf.std()])
times = (np.arange(len(odf)) * R + N / 2) / fs

# 2. Tempo: autocorrelation weighted by a log-Gaussian prior around 120 BPM
lags = np.arange(25, 101)                           # 60 to 240 BPM at 100 frames/s
ac = np.array([odf[:-l] @ odf[l:] for l in lags])
prior = np.exp(-0.5 * np.log2(6000 / lags / 120) ** 2)
lag = lags[np.argmax(ac * prior)]
print(f"tempo estimate: {6000 / lag:.1f} BPM (true {bpm})")

# 3. Beat tracking by dynamic programming (Ellis, 2007)
alpha, tau = 100.0, lag
C, back = odf.copy(), np.zeros(len(odf), int)
for n in range(len(odf)):
    prev = np.arange(n - 2 * tau, n - tau // 2 + 1)
    prev = prev[prev >= 0]
    if len(prev):
        score = C[prev] - alpha * np.log((n - prev) / tau) ** 2
        k = np.argmax(score); C[n] += score[k]; back[n] = prev[k]
n = np.argmax(C[-2 * tau:]) + len(C) - 2 * tau
path = [n]
while back[path[-1]] > 0: path.append(back[path[-1]])
est = times[np.array(path[::-1])]

def f_measure(est, ref, tol=0.07):
    hits = sum(np.min(np.abs(ref - e)) < tol for e in est)
    p, r = hits / len(est), hits / len(ref)
    return 2 * p * r / (p + r)
print(f"beats found: {len(est)} (true {len(beats)}), F-measure at 70 ms: {f_measure(est, beats):.3f}")

# 4. Downbeat: the bar position whose beats carry the most low-frequency onset energy
low = np.maximum(np.diff(L[:, :8], axis=0), 0).sum(1)       # 0 to 125 Hz
energy = np.interp(est, times[1:], low)
phase = np.argmax([energy[p::4].mean() for p in range(4)])
down_est, down_ref = est[phase::4], beats[::4]
print(f"downbeats: {len(down_est)} estimated, F-measure at 70 ms: {f_measure(down_est, down_ref):.3f}")

The autocorrelation peaks at a lag of 47 frames, that is 127.7 BPM. The dynamic programme finds 39 of the 40 beats (F-measure 0.987), and the low-frequency salience selects the right bar position (downbeat F-measure 0.947). The simplified code always keeps the best predecessor, a slight departure from the \max(0, \cdot) of Proposition 3.5.

4. Digital audio effects

Most studio effects are small systems whose behaviour follows from their transfer function. This chapter builds delay-based effects, the phaser, distortion and artificial reverberation from the filters of ch. 0.

Comb filters, flanger and chorus

By Proposition 4.2 the flanger’s notches are harmonically spaced and sweep together. A delay varying at the rate D'(t) also scales the instantaneous frequency by 1 - D'(t): this Doppler effect is the pitch modulation of vibrato and chorus.

Phaser

Distortion

Since \tanh is not a polynomial, its harmonics never stop, and digital distortion is computed at an oversampled rate, then low-pass filtered and decimated.

Artificial reverberation

For a 10 × 8 × 4 m room with \alpha = 0.2 everywhere, V = 320 m³ and A = 0.2 \times 304 = 60.8 m², so T_{60} \approx 0.85 \text{ s}.

Schroeder’s reverberator (1962) combines both: four parallel feedback combs with mutually prime delays of 30–45 ms, tuned by Proposition 4.9, then two short all-pass filters that multiply the echo density without colouring the spectrum.

Feedback delay networks and convolution with measured responses, the modern alternatives, are treated in ch. 11.

import numpy as np
from scipy.signal import lfilter, freqz

fs = 16000
# 1. Flanger at one instant: y = x + g x(n - D), D = 2 ms -> notches at odd multiples of 250 Hz
D, g = int(0.002 * fs), 0.9
f, Hc = freqz(np.r_[1, np.zeros(D - 1), g], worN=8192, fs=fs)
notch = f[1:-1][(np.abs(Hc[1:-1]) < np.abs(Hc[:-2])) & (np.abs(Hc[1:-1]) < np.abs(Hc[2:]))]
print("flanger notches (Hz):", np.round(notch[:4]))

# 2. Phaser: y = (x + A(z) x) / 2 with 4 first-order all-pass sections
a = -0.6
num, den = np.array([1.0]), np.array([1.0])
for _ in range(4):
    num, den = np.convolve(num, [a, 1]), np.convolve(den, [1, a])
f, A = freqz(num, den, worN=8192, fs=fs)
P = np.abs(1 + A) / 2
notch = f[1:-1][(P[1:-1] < P[:-2]) & (P[1:-1] < P[2:])]
print("all-pass |A| in [", np.abs(A).min().round(6), ",", np.abs(A).max().round(6), "]")
print("phaser notches (Hz):", np.round(notch))

# 3. Distortion: tanh soft clipping of a 500 Hz tone, harmonic levels
n = np.arange(fs)
y = np.tanh(4 * np.sin(2 * np.pi * 500 * n / fs))
Y = np.abs(np.fft.rfft(y * np.hanning(fs)))
lev = [20 * np.log10(Y[500 * k] / Y[500]) for k in range(2, 8)]
print("harmonics 2..7 (dB re fundamental):", np.round(lev, 1))

# 4. Schroeder reverberator: 4 parallel combs + 2 all-pass, target T60 = 1.5 s
T60 = 1.5
x = np.zeros(2 * fs); x[0] = 1.0
wet = np.zeros_like(x)
for ms in (29.7, 37.1, 41.1, 43.7):
    d = int(ms * fs / 1000)
    gc = 10 ** (-3 * d / (T60 * fs))          # -60 dB after T60 seconds
    wet += lfilter([1.0], np.r_[1, np.zeros(d - 1), -gc], x)
for ms, ga in ((5.0, 0.7), (1.7, 0.7)):
    d = int(ms * fs / 1000)
    wet = lfilter(np.r_[-ga, np.zeros(d - 1), 1], np.r_[1, np.zeros(d - 1), -ga], wet)
edc = 10 * np.log10(np.cumsum(wet[::-1] ** 2)[::-1] / np.sum(wet ** 2))  # Schroeder integration
i5, i35 = np.argmax(edc < -5), np.argmax(edc < -35)
print(f"measured T60 (from the -5 to -35 dB decay): {2 * (i35 - i5) / fs:.2f} s")

The flanger’s notches fall at 250, 750, 1250 and 1750 Hz (Proposition 4.2), and the four all-pass sections give 2 notches, at 525 and 2766 Hz.

The tanh distortion has odd harmonics only (−11.3, −18.8 and −25.7 dB for the 3rd, 5th and 7th; even ones below −190 dB). The reverberator decays in the target 1.50 s, measured on Schroeder’s backward integral (ch. 11).

5. Time-scale, pitch and timbre modifications

Playing a recording faster raises its pitch: duration and frequency are coupled. This chapter decouples them, first in time (time-scale modification), then in pitch, then in timbre, by working on the STFT of ch. 1 and on the spectral envelope.

Resampling couples time and pitch

A pitch shift by s semitones uses \beta = 2^{s/12}; four semitones (a major third) give \beta = 1.26. To keep the duration, the signal must first be stretched by \beta without changing its pitch.

Time-scale modification by overlap-add

On polyphonic music no single shift aligns every component; the phase vocoder works per frequency bin instead.

The phase vocoder

Each bin is thus re-synthesised at its true frequency with a phase consistent at the new hop. The method still smears transients and adds a reverberant “phasiness”, because the phases of neighbouring bins of one partial drift apart. Phase locking (Laroche and Dolson, 1999) propagates the phase of each spectral peak and keeps the bins around it locked to it.

Pitch shifting and the spectral envelope

This moves the spectral envelope along with the partials: the formants of a voice shift too, and a voice raised by an octave sounds like a cartoon character. Preserving timbre requires the source–filter view of Proposition 1.5.

import numpy as np

fs, N, Ra = 16000, 1024, 256                  # Hann window, analysis hop N/4
w = np.hanning(N + 1)[:N]
omega = 2 * np.pi * np.arange(N // 2 + 1) / N  # bin centre frequencies (rad/sample)
princarg = lambda p: (p + np.pi) % (2 * np.pi) - np.pi

def phase_vocoder(x, alpha):
    """Time-scale modification by a factor alpha (> 1 = slower), pitch unchanged."""
    Rs = int(round(alpha * Ra))
    frames = [np.fft.rfft(x[i:i + N] * w) for i in range(0, len(x) - N, Ra)]
    y = np.zeros(Rs * len(frames) + N); norm = np.zeros_like(y)
    psi = np.angle(frames[0])
    for m, X in enumerate(frames):
        if m > 0:
            dphi = np.angle(X) - np.angle(frames[m - 1]) - omega * Ra
            inst = omega + princarg(dphi) / Ra       # instantaneous frequency per bin
            psi = psi + Rs * inst                    # phase propagation at the synthesis hop
        frame = np.fft.irfft(np.abs(X) * np.exp(1j * psi), n=N) * w
        y[m * Rs:m * Rs + N] += frame; norm[m * Rs:m * Rs + N] += w ** 2
    return y / np.maximum(norm, 1e-3)

def peaks(x, k=2):
    X = np.abs(np.fft.rfft(x * np.hanning(len(x)), 8 * len(x)))
    f = np.fft.rfftfreq(8 * len(x), 1 / fs)
    loc = [i for i in range(1, len(X) - 1) if X[i] >= X[i - 1] and X[i] > X[i + 1]]
    return np.sort(f[sorted(loc, key=lambda i: -X[i])[:k]]).round(1)

t = np.arange(fs) / fs                             # 1 s of a two-note chord
x = np.sin(2 * np.pi * 440 * t) + 0.5 * np.sin(2 * np.pi * 660 * t + 1)
y = phase_vocoder(x, 1.5)
print(f"duration {len(x) / fs:.2f} s -> {len(y) / fs:.2f} s; peaks {peaks(x)} -> {peaks(y[N:-N])} Hz")

beta = 2 ** (4 / 12)                               # pitch shift by a major third
z = phase_vocoder(x, beta)
z = np.interp(np.arange(0, len(z) - 1, beta), np.arange(len(z)), z)   # resample by 1/beta
print(f"pitch shift: duration {len(z) / fs:.2f} s, peaks {peaks(z[N:-N])} Hz "
      f"(expected {np.round(440 * beta, 1)} and {np.round(660 * beta, 1)})")

The stretched chord lasts 1.48 s (1.5 s minus the edge frames) with its peaks still at 440 and 660 Hz. The pitch-shifted version lasts 1.00 s with peaks at 554.4 and 831.5 Hz, as predicted.

The second block estimates the envelope of a synthetic vowel (formants at 700 and 1200 Hz, pitch 120 Hz) both ways.

import numpy as np
from scipy.signal import lfilter, freqz
from scipy.linalg import solve_toeplitz

fs, f0 = 16000, 120
# Source-filter vowel: impulse train through two resonances (formants at 700 and 1200 Hz)
poles = []
for F, B in ((700, 80), (1200, 100)):
    r = np.exp(-np.pi * B / fs)
    poles += [r * np.exp(2j * np.pi * F / fs), r * np.exp(-2j * np.pi * F / fs)]
a_true = np.real(np.poly(poles))
e = np.zeros(fs // 2); e[::fs // f0] = 1.0
x = lfilter([1.0], a_true, e)

seg = x[2000:2000 + 1024] * np.hanning(1024)
# 1. Cepstral envelope: keep quefrencies below the pitch period (fs/f0 = 133 samples)
logS = np.log(np.abs(np.fft.fft(seg, 4096)) + 1e-9)
c = np.fft.ifft(logS).real
lifter = np.zeros(4096); lifter[:50] = 1; lifter[-49:] = 1
env_cep = np.exp(np.fft.fft(c * lifter).real[:2049])
# 2. LPC envelope of order 8 (autocorrelation method, Yule-Walker equations)
r = np.correlate(seg, seg, "full")[1023:1023 + 9]
a = solve_toeplitz(r[:8], -r[1:9])
f, Hl = freqz([1.0], np.r_[1, a], worN=2049, fs=fs)
for name, env in (("cepstrum", env_cep), ("LPC", np.abs(Hl))):
    pk = [i for i in range(1, 600) if env[i] > env[i - 1] and env[i] > env[i + 1]]
    print(f"{name:8s} envelope peaks below 2400 Hz: {np.round(f[pk])} Hz")

The cepstral envelope peaks at 726 and 1163 Hz and the LPC envelope at 711 and 1199 Hz. A shorter lifter (20 coefficients) smooths the two formants into a single peak, the usual trade-off between smoothness and detail.

6. Spatial audio: binaural, transaural and wave field synthesis

3D rendering makes a listener hear a sound at a chosen position. Perceptual approaches reproduce the signals at the two eardrums; physical approaches reproduce the sound field over a whole area. This chapter derives both, with their limits.

Localisation cues

With a = 8.75 cm the ITD reaches 0.66 ms at 90°. The ITD dominates below about 1.5 kHz, and the ILD, caused by head shadowing, above it (the duplex theory). Both cues are equal for all directions on a cone of confusion; the spectral filtering of the pinnae resolves front–back and elevation ambiguities.

Head-related transfer functions and binaural synthesis

HRTFs are measured on a grid of directions, on dummy heads such as KEMAR (Gardner and Martin, 1995) or on people (CIPIC, Algazi et al., 2001). Non-individual HRTFs cause front–back confusions and in-head localisation; head tracking reduces both.

Transaural reproduction

Transaural rendering plays binaural signals over two loudspeakers. Each ear hears both loudspeakers, so the crosstalk must be cancelled.

In free field, G_c \approx G_i\, e^{-2i\pi f \Delta / c}, where \Delta is the path difference between the far and the near loudspeaker. The inverse then blows up at low frequency and at the multiples of c / 2\Delta. The cancellation also holds only in a small “sweet spot” around the assumed head position.

Wave field synthesis

Wave field synthesis (Berkhout, de Vries and Vogel, 1993) reproduces a sound field inside a region by Huygens’ principle: a wavefront can be re-emitted by secondary sources placed on a surface.

For a plane boundary, choosing a Green’s function with zero normal derivative on the plane removes the dipoles and doubles the monopoles: this is the Rayleigh I integral. A virtual source behind the plane is rendered in the half-space in front by monopoles with driving function -2\, \partial P / \partial n, the normal now pointing into the listening region. WFS samples this layer with a line of loudspeakers.

With 15 cm between loudspeakers, waves arriving at 45° alias above 1339 Hz.

import numpy as np

c, a = 343.0, 0.0875                                  # speed of sound (m/s), head radius (m)
# 1. Woodworth's interaural time difference
for deg in (0, 30, 60, 90):
    th = np.radians(deg)
    print(f"azimuth {deg:2d} deg: ITD = {1e3 * a / c * (th + np.sin(th)):.3f} ms")

# 2. Crosstalk cancellation, free-field model: loudspeakers at +-30 deg, 2 m; ears at +-a
ears = np.array([[-a, 0.0], [a, 0.0]])
spk = 2.0 * np.array([[-np.sin(np.pi / 6), np.cos(np.pi / 6)], [np.sin(np.pi / 6), np.cos(np.pi / 6)]])
dist = np.linalg.norm(ears[:, None, :] - spk[None, :, :], axis=2)    # ear x loudspeaker
def G(f):                                            # acoustic transfer matrix (point sources)
    k = 2 * np.pi * f / c
    return np.exp(-1j * k * dist) / (4 * np.pi * dist)
delta = dist[0, 1] - dist[0, 0]                       # contralateral minus ipsilateral path
print(f"path difference {100 * delta:.2f} cm: G is singular at multiples of c/(2 delta) = {c / (2 * delta):.0f} Hz")
for f in (100, 1000, 1900, 3000):
    Gf = G(f)
    beta = 1e-3 * np.linalg.norm(Gf, 2) ** 2           # Tikhonov regularisation
    Cf = np.linalg.solve(Gf.conj().T @ Gf + beta * np.eye(2), Gf.conj().T)
    E = Gf @ Cf                                       # ear signals = E @ binaural signals
    sep = 20 * np.log10(np.abs(E[0, 0]) / np.abs(E[1, 0]))
    effort = 20 * np.log10(np.linalg.norm(Cf, 2) * np.linalg.norm(Gf, 2))
    print(f"f = {f:4d} Hz: cond(G) = {np.linalg.cond(Gf):6.1f}, separation {sep:5.1f} dB, "
          f"filter gain {effort:5.1f} dB")

# 3. Wave field synthesis of a plane wave with a line of monopoles (Rayleigh I, 2.5D)
dx = 0.15                                             # loudspeaker spacing (m)
xs = np.arange(-4, 4 + 1e-9, dx)
grid = np.array([(u, v) for u in np.linspace(-0.5, 0.5, 11) for v in np.linspace(1, 2, 11)])
r = np.linalg.norm(grid[:, None, :] - np.c_[xs, 0 * xs][None], axis=2)
for deg in (0, 45):
    th = np.radians(deg)
    print(f"plane wave at {deg} deg: aliasing above c / (dx (1 + sin)) = {c / (dx * (1 + np.sin(th))):.0f} Hz")
    for f in (500, 1000, 2000, 4000):
        k = 2 * np.pi * f / c
        drive = np.exp(-1j * k * xs * np.sin(th))           # delays of the incoming wave
        p = (drive * np.exp(-1j * k * r) / (4 * np.pi * r)).sum(1)
        target = np.exp(-1j * k * (grid[:, 0] * np.sin(th) + grid[:, 1] * np.cos(th)))
        gain = np.vdot(target, p) / np.vdot(target, target)   # best complex scale
        err = np.linalg.norm(p - gain * target) / np.linalg.norm(gain * target)
        print(f"   f = {f:4d} Hz: relative error of the synthesised field {err:.2f}")

The ITDs are 0.26, 0.49 and 0.66 ms at 30°, 60° and 90°. The crosstalk canceller separates the ears by 90 and 80 dB at 1000 and 3000 Hz, but only by 24 and 17 dB at 100 and 1900 Hz, near its singular frequencies, despite 20–23 dB of extra gain.

The array reproduces the frontal wave with a relative error of 0.12–0.18 up to 2000 Hz, and 1.60 at 4000 Hz, beyond its aliasing frequency (2287 Hz). At 45° the error jumps from 0.16 at 1000 Hz to 0.89 at 2000 Hz.

The residual of about 12 % at low frequency is the 2.5D amplitude error, not an effect of the array length. A line of point sources makes the wave decay as 1/\sqrt{r}, by a factor \sqrt2 between 1 and 2 m, and a 128 m array leaves the error at 0.11.

7. Source separation I: linear and convolutive mixtures

Source separation, or demixing, recovers individual sources from recordings of their mixture. This chapter treats mixtures with at least as many microphones as sources, which can be inverted once identified blindly from source statistics.

Mixing models

Independent component analysis

The blind assumption is that the sources are mutually independent. Second-order statistics alone do not suffice.

Maximum-likelihood ICA, which updates the whole demixing matrix at once, is presented in ch. 11.

Convolutive mixtures in the frequency domain

Frequency-domain ICA separates each bin with complex-valued ICA (Smaragdis, 1998). Proposition 7.2 now applies per bin. Each bin returns the sources in its own order (the permutation problem), which is resolved by aligning the bins on the correlation of the source envelopes across frequencies or on the directions of arrival.

Each bin also has its own scale, which the minimal distortion principle fixes: \mathbf{W}(k) \leftarrow \mathrm{diag}\big(\mathbf{W}(k)^{-1}\big)\, \mathbf{W}(k) . Each output is then the spatial image of its source on one microphone.

Measuring separation quality

The code mixes a sawtooth and a modulated Laplacian noise with a fixed 2 \times 2 matrix, compares whitening alone with FastICA, then measures the narrowband error for a 32 ms room-like filter.

import numpy as np
from sklearn.decomposition import FastICA, PCA
from scipy.signal import lfilter

rng = np.random.default_rng(0)
fs = 8000
t = np.arange(2 * fs) / fs
s1 = 2 * ((110 * t) % 1) - 1                                   # sawtooth (sub-Gaussian)
s2 = rng.laplace(size=len(t)) * (0.6 + 0.4 * np.sin(2 * np.pi * 3 * t))   # speech-like, super-Gaussian
S = np.c_[s1, s2] / np.c_[s1, s2].std(0)
A = np.array([[1.0, 0.6], [0.5, 1.0]])
X = S @ A.T                                                     # instantaneous mixture x = A s

def sir_db(Y, S):
    """Signal-to-interference ratio after resolving permutation and scale."""
    C = np.abs(np.corrcoef(Y.T, S.T)[:2, 2:])                    # |corr(y_i, s_j)|
    perm = [0, 1] if C[0, 0] + C[1, 1] > C[0, 1] + C[1, 0] else [1, 0]
    out = []
    for i, j in enumerate(perm):
        G = np.linalg.lstsq(S, Y[:, i], rcond=None)[0]         # y_i = G_j s_j + G_other s_other
        out.append(20 * np.log10(abs(G[j]) / abs(G[1 - j])))
    return np.round(out, 1)

print("mixtures      SIR (dB):", sir_db(X, S))
print("PCA whitening SIR (dB):", sir_db(PCA(2, whiten=True).fit_transform(X), S))
print("FastICA       SIR (dB):", sir_db(FastICA(2, whiten="unit-variance", random_state=0).fit_transform(X), S))

# Narrowband approximation of a convolutive mixture: X(m,k) ~ A(k) S(m,k)?
h = rng.standard_normal(256) * np.exp(-np.arange(256) / 40)   # a 32 ms room-like filter
s = rng.standard_normal(4 * fs)
xconv = lfilter(h, [1.0], s)
for N in (256, 1024, 4096):
    w = np.hanning(N)
    fr = lambda z: np.array([np.fft.rfft(z[i:i + N] * w) for i in range(N, len(z) - N, N // 4)])
    Xs, Ss = fr(xconv), fr(s)
    Hk = np.fft.rfft(h, N)
    err = np.linalg.norm(Xs - Hk * Ss) / np.linalg.norm(Xs)
    print(f"window {N:4d} samples: relative error of X = H S per bin: {err:.2f}")

The mixtures have SIRs of 4.4 and 6.0 dB, whitening alone 0.4 dB (the residual rotation of Proposition 7.3), FastICA 42.5 and 44.1 dB. The narrowband error falls from 0.36 to 0.10 and 0.02 for windows of 256, 1024 and 4096 samples.

8. Source separation II: underdetermined and sparse models, DUET

With fewer microphones than sources, the mixing matrix has no inverse and independence is not enough. The extra assumption that makes separation possible is sparsity: in the time–frequency plane, few sources are active at each point.

Sparsity in the time–frequency plane

Then each time–frequency point belongs to one source, and separation reduces to classifying points and masking. A softer version of the assumption is a sparse prior: each coefficient follows a heavy-tailed density, such as the Laplacian, whose maximum a posteriori estimate penalises the ℓ1 norm.

For complex STFT coefficients the same problem is a second-order cone programme. In practice \mathbf{A} is unknown and must first be estimated, which is what DUET does from the geometry of anechoic mixtures.

DUET

The symmetric attenuation treats the two microphones alike: swapping them changes the sign of \alpha. Steps 3 and 4 are maximum-likelihood decisions under white Gaussian noise on both channels.

import numpy as np
from scipy.signal import stft, istft
from sklearn.cluster import KMeans

fs, rng = 16000, np.random.default_rng(0)
t = np.arange(3 * fs) / fs

def voice(f0, rate):
    """Harmonic source with vibrato and random syllable-like on/off envelope."""
    f = f0 * (1 + 0.02 * np.sin(2 * np.pi * 5 * t))
    x = sum(np.sin(2 * np.pi * k * np.cumsum(f) / fs) / k for k in range(1, 16) if k * f0 < 7000)
    gates = np.repeat(rng.random(int(3 * rate) + 1) < 0.6, int(fs / rate))[:len(t)]
    env = np.convolve(gates.astype(float), np.hanning(800) / 400, "same")
    return x * env

src = [voice(f0, r) for f0, r in ((110, 4), (175, 5), (260, 6))]
a_true, d_true = (0.6, 1.0, 1.5), (-0.7, 0.0, 0.6)          # attenuations and delays (samples)
omega = 2 * np.pi * np.fft.rfftfreq(len(t))
x1 = sum(src)
x2 = sum(np.fft.irfft(a * np.fft.rfft(s) * np.exp(-1j * omega * d), len(t))
         for s, a, d in zip(src, a_true, d_true))

N = 1024
_, _, X1 = stft(x1, nperseg=N); _, _, X2 = stft(x2, nperseg=N)
w = 2 * np.pi * np.arange(X1.shape[0])[:, None] / N            # rad/sample per bin
keep = (np.abs(X1) > 1e-3 * np.abs(X1).max())
keep[0] = False                                                # DC carries no delay information
R = X2[keep] / X1[keep]
alpha = np.abs(R) - 1 / np.abs(R)                              # symmetric attenuation
delta = -np.angle(R) / np.broadcast_to(w, X1.shape)[keep]
weight = (np.abs(X1[keep] * X2[keep]))                         # DUET weighting, p = 1, q = 0
feats = np.c_[alpha, delta]
km = KMeans(3, n_init=10, random_state=0).fit(feats, sample_weight=weight)
for c in km.cluster_centers_[np.argsort(km.cluster_centers_[:, 1])]:
    a_hat = (c[0] + np.sqrt(c[0] ** 2 + 4)) / 2                # invert alpha = a - 1/a
    print(f"estimated attenuation {a_hat:.2f}, delay {c[1]:+.2f} samples")

centres = km.cluster_centers_
a_est = (centres[:, 0] + np.sqrt(centres[:, 0] ** 2 + 4)) / 2
cost = [np.abs(a_est[j] * np.exp(-1j * w * centres[j, 1]) * X1 - X2) ** 2 / (1 + a_est[j] ** 2)
        for j in range(3)]                                     # residual if only source j were active
label = np.argmin(cost, axis=0)

def sdr(est, ref):
    g = est @ ref / (ref @ ref)
    return 10 * np.log10(np.sum((g * ref) ** 2) / np.sum((est - g * ref) ** 2))
for j in np.argsort(centres[:, 1]):
    S = (X1 + a_est[j] * np.exp(1j * w * centres[j, 1]) * X2) / (1 + a_est[j] ** 2)
    est = istft(np.where(label == j, S, 0), nperseg=N)[1][:len(t)]
    jt = np.argmin([abs(a_est[j] - a) for a in a_true]); ref = src[jt]
    print(f"source {jt + 1}: SDR of the mixture {sdr(x1, ref):5.1f} dB -> after DUET {sdr(est, ref):5.1f} dB")

Three harmonic “voices” (110, 175 and 260 Hz, with vibrato and random pauses) are mixed onto two channels. The weighted clustering, used here instead of histogram peak picking, recovers the attenuations 0.61, 1.02, 1.48 and the delays −0.76, −0.02, +0.63 samples (true: 0.6, 1.0, 1.5 and −0.7, 0, 0.6).

The masks raise the SDRs from −4.6, −3.1 and −1.6 dB in the mixture to 14.7, 10.7 and 14.9 dB.

Soft masks and the Wiener filter

The ratio v_j / \sum_l v_l is a soft mask; when one variance dominates at each point it tends to the binary mask of DUET. The difficulty moves to estimating the variances v_j(m, k). Nonnegative matrix factorisation (ch. 11), the multichannel local Gaussian model (Duong, Vincent and Gribonval, 2010) and neural networks are the main tools.

9. High-resolution analysis I: sinusoidal models beyond Fourier

Many musical sounds are sums of slowly varying sinusoids, and their analysis needs frequency estimates finer than the Fourier resolution of a short frame. This chapter states the model, shows where Fourier analysis fails, and exposes the algebraic structure that high-resolution methods exploit.

Sinusoidal models

Sinusoidal analysis–synthesis (McAulay and Quatieri, 1986) tracks the peaks of successive short-time spectra as partials; spectral modelling synthesis (Serra and Smith, 1990) adds a stochastic residual. Both rely on accurate per-frame estimates of (\omega_k, \delta_k, a_k, \phi_k), which is what high-resolution methods provide.

What Fourier analysis cannot do

The limit is one of resolution, not of precision. For a single tone, the precision has a floor set by the noise:

At 10 dB of signal-to-noise ratio and N = 64, the bound on the standard deviation is 2.4 \times 10^{-4} cycles per sample, 65 times smaller than the Fourier resolution 1/N.

Linear prediction and Prony’s method

Prony’s method is exact without noise, but noise makes the equations inconsistent, and the least-squares coefficients, hence the roots, are biased. The subspace methods of ch. 10 are its robust replacement.

Hankel and covariance matrices

import numpy as np

rng = np.random.default_rng(0)
N = 64
n = np.arange(N)
f_true = np.array([0.200, 0.208])                 # cycles/sample, gap 0.008 < 1/N = 0.0156
z = np.exp(2j * np.pi * f_true)
alpha = np.array([1.0, 0.8 * np.exp(1j * np.pi / 3)])
x0 = (alpha[None, :] * z[None, :] ** n[:, None]).sum(1)

# 1. Periodogram (zero-padded): where are its peaks?
X = np.abs(np.fft.fft(x0, 8192)) ** 2
f = np.fft.fftfreq(8192)
band = (f > 0.15) & (f < 0.25)
pk = [i for i in np.where(band)[0] if X[i] > X[i - 1] and X[i] > X[i + 1]]
print("periodogram peaks above half the maximum:", np.round(f[pk][X[pk] > 0.5 * X.max()], 4))

# 2. Prony / linear prediction: x(n) + c1 x(n-1) + c2 x(n-2) = 0
def prony(x, K=2):
    Amat = np.array([x[k - K:k][::-1] for k in range(K, len(x))])
    c = np.linalg.lstsq(Amat, -x[K:], rcond=None)[0]
    return np.sort(np.angle(np.roots(np.r_[1, c])) / (2 * np.pi))   # frequencies of the roots
print("Prony, no noise      :", np.round(prony(x0), 4))
for snr_db in (60, 30):
    sigma = 10 ** (-snr_db / 20)
    x = x0 + sigma * (rng.standard_normal(N) + 1j * rng.standard_normal(N)) / np.sqrt(2)
    print(f"Prony, SNR {snr_db} dB    :", np.round(prony(x), 4))

# 3. Eigenvalues of the empirical covariance (Hankel data matrix, L = 16)
L = 16
x = x0 + 0.1 * (rng.standard_normal(N) + 1j * rng.standard_normal(N)) / np.sqrt(2)
H = np.array([x[i:i + N - L + 1] for i in range(L)])
ev = np.linalg.eigvalsh(H @ H.conj().T / (N - L + 1))[::-1]
print("largest eigenvalues:", np.round(ev[:4], 3), "| median of the others:",
      np.round(np.median(ev[2:]), 4), "| noise power 0.01")

Two tones at 0.200 and 0.208 cycles per sample, half the Fourier resolution apart, give periodogram peaks at 0.1957 and 0.2147, each wrong by more than half the gap (Proposition 9.2).

Prony is exact without noise and nearly so at 60 dB, but returns 0.0233 and 0.2018 at 30 dB. With noise power 0.01, the covariance has two eigenvalues above the noise (9.42 and 0.395) and the others near 0.01 (Proposition 9.7); the second is small because the two Vandermonde vectors are nearly parallel.

10. High-resolution analysis II: MUSIC, ESPRIT and model order

The signal subspace of ch. 9 contains the Vandermonde vectors of the poles. MUSIC finds the poles as the Vandermonde vectors orthogonal to the noise subspace; ESPRIT reads them off a shift invariance of the signal subspace. Both need the model order K, which the ESTER criterion estimates.

Estimated subspaces

With \mathbf{v}(z) = (1, z, \dots, z^{L-1})^\top the Vandermonde vector of a candidate pole, the two methods follow from Propositions 9.6–9.7.

MUSIC

MUSIC was designed for direction-of-arrival estimation with sensor arrays; for time series it assumes undamped poles, since it searches the unit circle.

ESPRIT

For a matrix \mathbf{A} with L rows, write \mathbf{A}_\downarrow for \mathbf{A} without its last row and \mathbf{A}_\uparrow for \mathbf{A} without its first row.

The cost is dominated by the SVD; fast subspace tracking (Badeau, David and Richard, 2005) updates \mathbf{W} from frame to frame at a lower cost. Unlike MUSIC, ESPRIT estimates damping factors, which audio partials have.

Choosing the model order

Information criteria such as MDL (Wax and Kailath, 1985) count the eigenvalues above the noise level, which assumes white noise. In practice audio noise is whitened first, and signals are split into subbands that each hold few components.

import numpy as np

rng = np.random.default_rng(0)
N, L = 64, 21                                     # signal length, Hankel rows (about N/3)
n = np.arange(N)
f_true = np.array([0.200, 0.208, 0.300])          # two close tones and a damped one
damp = np.array([0.0, 0.0, -0.02])
alpha = np.array([1.0, 0.8 * np.exp(1j * np.pi / 3), 0.5])
z = np.exp(damp + 2j * np.pi * f_true)
x0 = (alpha * z ** n[:, None]).sum(1)
sigma = 10 ** (-30 / 20)                          # SNR 30 dB for the first tone
x = x0 + sigma * (rng.standard_normal(N) + 1j * rng.standard_normal(N)) / np.sqrt(2)

H = np.array([x[i:i + N - L + 1] for i in range(L)])          # L x (N-L+1) Hankel matrix
U, s, _ = np.linalg.svd(H)

def esprit(K):
    W = U[:, :K]                                             # signal subspace basis
    Phi = np.linalg.lstsq(W[:-1], W[1:], rcond=None)[0]      # W_up = W_down Phi
    zk = np.linalg.eigvals(Phi)
    return zk, np.linalg.norm(W[1:] - W[:-1] @ Phi, 2)

zk, _ = esprit(3)
order = np.argsort(np.angle(zk))
print("ESPRIT frequencies:", np.round(np.angle(zk[order]) / (2 * np.pi), 4),
      "damping:", np.round(np.log(np.abs(zk[order])), 3))

V = zk[None, order] ** n[:, None]                             # amplitudes by least squares
print("amplitudes |alpha|:", np.round(np.abs(np.linalg.lstsq(V, x, rcond=None)[0]), 2))

# ESTER: J(p) = 1 / ||W_up - W_down Phi||^2; keep the largest p with J(p) >= 0.1 max J
J = np.array([1 / esprit(p)[1] ** 2 for p in range(1, 8)])
print("ESTER J(p), p = 1..7:", np.round(J).astype(int), "-> order", 1 + np.max(np.where(J >= 0.1 * J.max())[0]))

# MUSIC pseudo-spectrum on a fine grid
Wn = U[:, 3:]                                                 # noise subspace
grid = np.linspace(0.15, 0.35, 4001)
steer = np.exp(2j * np.pi * np.outer(np.arange(L), grid))
P = 1 / np.linalg.norm(Wn.conj().T @ steer, axis=0) ** 2
pk = [i for i in range(1, len(grid) - 1) if P[i] > P[i - 1] and P[i] > P[i + 1]]
best = sorted(pk, key=lambda i: -P[i])[:3]
print("MUSIC peaks:", np.round(np.sort(grid[best]), 4))

At 30 dB, ESPRIT finds 0.1997, 0.2083 and 0.2999 cycles per sample (true 0.200, 0.208, 0.300), the damping −0.021 (true −0.02), and amplitudes 0.86, 0.84 and 0.51 (true 1.0, 0.8, 0.5); the two close amplitudes are the least accurate.

ESTER gives J = 15277 at p = 3, against 1243 and 93 at p = 1, 2 and at most 32 above, and selects K = 3. MUSIC peaks at 0.1998, 0.2083 and 0.3000.

11. Complements for self-study

This chapter collects proofs and methods that extend the in-class chapters. Each section names the chapter it complements; none is needed to follow the sessions, and all are useful for the projects.

EM for Gaussian mixtures (complement to ch. 2)

The GMM classifier of Definition 2.2 fits one mixture per class to the training frames of that class.

In practice the covariances are often diagonal, initialised by k-means, and floored to avoid degenerate components collapsing on single frames.

Feedback delay networks and convolution reverberation (complement to ch. 4)

Unlike Schroeder’s parallel combs, the network mixes every line into every other one, so the echo density grows quickly and the modes are denser. Frequency-dependent decay is obtained by replacing each gain \gamma^{D_j} by a low-pass filter with that gain at each frequency.

Convolution reverberation convolves the signal with a measured impulse response by FFT blocks (Proposition 0.2), partitioned to keep the latency low (Gardner, 1995).

Decay times are measured on Schroeder’s backward integral of a response h, E(t) = \int_t^\infty h^2(s)\, ds, which averages out the fluctuations of h^2. Its level in dB is close to a straight line, whose slope between −5 and −35 dB, extrapolated to 60 dB, gives T_{60}; the code of ch. 4 does exactly this.

Linear prediction (complement to ch. 5)

The Toeplitz system is solved in O(P^2) operations by the Levinson–Durbin recursion, and the autocorrelation method yields a stable all-pole model. An order of about f_s / 1000 + 2 suits speech.

Inverting a modified STFT (complement to ch. 1 and 5)

Masking (ch. 8) or re-phasing (ch. 5) produces an array Y(m, k) that is usually not the STFT of any signal. Overlap-add still returns the best possible signal in the following sense.

When only the magnitude |Y| is known, the Griffin–Lim algorithm alternates two steps: give |Y| the phases of the current estimate’s STFT, then invert by Proposition 11.4. Griffin and Lim (1984) show that the distance between the magnitude constraint and the set of consistent STFTs does not increase along the iterations.

Amplitude panning and Ambisonics (complement to ch. 6)

Between the perceptual rendering of binaural signals and the physical rendering of WFS, two loudspeaker techniques are widely used.

First-order Ambisonics (Gerzon, 1973) encodes a plane wave s from azimuth \theta and elevation \phi into four channels that do not depend on the loudspeaker layout: W = \frac{s}{\sqrt{2}}, \qquad X = s \cos\theta \cos\phi, \qquad Y = s \sin\theta \cos\phi, \qquad Z = s \sin\phi . These are the first spherical harmonics of the sound field at the listening point. A decoder computes loudspeaker feeds whose re-encoding matches these channels, for example by a pseudo-inverse. Higher orders sharpen the spatial resolution (Zotter and Frank, 2019).

Maximum-likelihood ICA (complement to ch. 7)

Maximum-likelihood ICA, equivalent to Infomax (Bell and Sejnowski, 1995), updates the whole demixing matrix by the natural gradient, where \phi is the score function of the assumed source density: \mathbf{W} \leftarrow \mathbf{W} + \mu \big(\mathbf{I} - \mathbb{E}[\phi(\mathbf{y})\, \mathbf{y}^\top]\big) \mathbf{W} .

It follows from the log-likelihood of one sample, \log |\det \mathbf{W}| + \sum_i \log p_i(y_i), whose gradient with respect to \mathbf{W}, multiplied on the right by \mathbf{W}^\top \mathbf{W} (the relative or natural gradient), gives the update. With \phi = \tanh it suits super-Gaussian sources such as speech.

Nonnegative matrix factorisation (complement to ch. 8)

The Wiener masks of Proposition 8.6 need the source variances v_j(m, k). Nonnegative matrix factorisation (NMF) models a power spectrogram as a sum of a few spectral templates with time-varying gains.

For separation, templates are learned per source (supervised NMF) or grouped after learning. Then v_j = \mathbf{W}_j \mathbf{H}_j feeds the Wiener masks. The Itakura–Saito divergence, which is scale-invariant and corresponds to the Gaussian model of Proposition 8.6, suits audio better (Févotte, Bertin and Durrieu, 2009), and multichannel NMF extends the model to convolutive mixtures (Ozerov and Févotte, 2010).

The code factorises a synthetic spectrogram made of two harmonic “notes” with sparse activations.

import numpy as np

rng = np.random.default_rng(0)
F, T, K = 64, 100, 2
W0 = np.zeros((F, K)); W0[[5, 10, 15, 20], 0] = 1; W0[[7, 14, 21, 28], 1] = 1   # two 'notes'
H0 = (rng.random((K, T)) < 0.3) * rng.random((K, T))                         # sparse activations
V = W0 @ H0 + 1e-3                                                            # power spectrogram
W, H = rng.random((F, K)) + 0.1, rng.random((K, T)) + 0.1
kl = lambda V, Vh: np.sum(V * np.log(V / Vh) - V + Vh)
hist = []
for it in range(300):                                  # Lee-Seung multiplicative updates (KL)
    H *= (W.T @ (V / (W @ H))) / W.sum(0)[:, None]
    W *= ((V / (W @ H)) @ H.T) / H.sum(1)[None, :]
    hist.append(kl(V, W @ H))
print("KL divergence at iterations 1, 10, 100, 300:", np.round([hist[0], hist[9], hist[99], hist[-1]], 3))
print("monotone decrease:", bool(np.all(np.diff(hist) <= 1e-9)))
print("partials of each learned template (4 largest bins):",
      [sorted(int(i) for i in np.argsort(W[:, k])[-4:]) for k in range(K)])

The divergence falls from 68.9 to 9.32 in 10 iterations and then stabilises at 9.319, decreasing at every step. The two learned templates recover the partials of the two notes (bins 5, 10, 15, 20 and 7, 14, 21, 28). The remaining divergence comes from the constant floor added to the spectrogram.

Harmonic–percussive separation (complement to ch. 8)

In a spectrogram, sustained harmonic sounds form horizontal lines and percussive sounds vertical lines. FitzGerald (2010) median-filters the spectrogram along time, which keeps the horizontal structures, and along frequency, which keeps the vertical ones. The two filtered spectrograms define soft or binary masks, a fast and useful front end for the rhythm analysis of ch. 3 and the pitch analysis of ch. 9.

The Cramér–Rao bound for a tone (complement to ch. 9)

For one tone, the maximum of the periodogram is the maximum-likelihood estimator of the frequency, and it reaches the bound above a threshold SNR; the high-resolution methods matter when several components are close. The code compares ESPRIT and the zero-padded periodogram peak with the bound, over 300 noise draws.

import numpy as np

rng = np.random.default_rng(0)
N, f0, trials = 64, 0.1234, 300
n = np.arange(N)
for snr_db in (0, 10, 20):
    s2 = 10 ** (-snr_db / 10)                                    # noise variance, |alpha| = 1
    crb = 6 * s2 / (N * (N ** 2 - 1)) / (2 * np.pi) ** 2         # in (cycles/sample)^2
    e_esp, e_per = [], []
    for _ in range(trials):
        x = np.exp(2j * np.pi * f0 * n + 1j * rng.uniform(0, 2 * np.pi)) \
            + np.sqrt(s2 / 2) * (rng.standard_normal(N) + 1j * rng.standard_normal(N))
        L = N // 3
        Hm = np.array([x[i:i + N - L + 1] for i in range(L)])
        u = np.linalg.svd(Hm)[0][:, :1]                          # ESPRIT with K = 1
        zhat = (np.linalg.pinv(u[:-1]) @ u[1:])[0, 0]
        e_esp.append(np.angle(zhat) / (2 * np.pi) - f0)
        Xf = np.abs(np.fft.fft(x, 64 * N))                      # grid step 1/(64 N)
        e_per.append(np.argmax(Xf) / (64 * N) - f0)
    print(f"SNR {snr_db:2d} dB: RMSE/sqrt(CRB) ESPRIT {np.sqrt(np.mean(np.square(e_esp)) / crb):.2f}, "
          f"periodogram peak {np.sqrt(np.mean(np.square(e_per)) / crb):.2f}")

The ratio of the root mean squared error to the square root of the bound is 1.16, 1.01 and 1.06 for ESPRIT at 0, 10 and 20 dB. For the periodogram peak it is 1.04 and 0.98 at 0 and 10 dB (the bound reached within Monte Carlo error), and 1.60 at 20 dB, where the frequency grid of step 1/(64N) becomes the limit.

Pisarenko’s method (1973), the ancestor of the subspace methods, takes L = K + 1: the noise subspace of Proposition 9.7 is then a single vector, and the roots of the polynomial whose coefficients are its entries are the poles. MUSIC and ESPRIT use larger L, which averages the noise over more dimensions.

Neural networks for indexing and separation (complement to ch. 2, 3 and 8)

Since the mid-2010s, convolutional networks applied to log-mel spectrograms have replaced hand-made features in large-scale audio classification (Hershey et al., 2017; Kong et al., 2020). Recurrent and convolutional networks estimate beat and downbeat activations, decoded by the dynamic programmes and hidden Markov models of ch. 3 (Böck, Krebs and Widmer, 2016).

In separation, networks estimate the masks of ch. 8 from the mixture spectrogram, or work directly on the waveform with a learned encoder (Luo and Mesgarani, 2019; Défossez et al., 2019). The evaluation pitfalls of Proposition 2.7 apply to them unchanged.

12. Research directions and open problems

This chapter is for self-study. It maps active questions that start from chapters 2–10, with pointers to papers listed in the references. Each item says what the question is and where to start reading.

Learned representations for machine listening

• Large labelled corpora such as AudioSet (Gemmeke et al., 2017) made deep audio classifiers possible (Hershey et al., 2017), and pretrained networks now serve as general feature extractors (Kong et al., 2020).
• Self-supervised learning removes the need for labels: wav2vec 2.0 (Baevski, Zhou, Mohamed and Auli, 2020) learns speech representations by a contrastive task on masked latent features. Contrastive language–audio pretraining (Elizalde, Deshmukh, Al Ismail and Wang, 2023) aligns sounds with text descriptions and enables zero-shot classification.
• Open questions: which representations transfer to music, whose structure (harmony, meter) differs from speech; and how to evaluate them without the leakage of Proposition 2.7.

Separation without isolated sources

• Supervised separation needs isolated sources to build training mixtures. Universal sound separation (Kavalerov et al., 2019) targets arbitrary sound classes, and mixture invariant training (Wisdom et al., 2020) learns from mixtures alone by remixing them.
• Evaluation is itself a research question: Le Roux, Wisdom, Erdogan and Hershey (2019) show that the SDR of Definition 7.8 can be gamed by filtering, and propose a scale-invariant version.
• Open questions: reverberant and moving sources, where the narrowband approximation of Proposition 7.7 fails; the gap between SDR and perceived quality; and how much of the spatial models of ch. 8 neural separators still need.

Signal processing inside learning

• Differentiable digital signal processing (Engel, Hantrakul, Gu and Roberts, 2020) writes oscillators, filters and reverberators as differentiable layers, so that the models of ch. 4–5 are trained by gradient descent from audio.
• Open questions: how far such structured models can replace black-box generators, and how to keep their parameters interpretable as pitch, envelope and room.

Super-resolution beyond subspace methods

• Off-the-grid sparse methods recover spectral lines by convex optimisation over measures: Candès and Fernandez-Granda (2014) prove exact recovery when lines are separated by a few times 1/N, and Tang, Bhaskar, Shah and Recht (2013) use the atomic norm with random samples.
• Subspace methods have been extended to exponentials modulated by polynomials, a model of amplitude modulation, and the performance of ESPRIT on them has been analysed (Badeau, Richard and David, 2008).
• Open questions: inharmonic and time-varying partials (piano, voice), the lowest separation at which close components can be resolved at a given SNR, and combinations of high-resolution front ends with learned models.

Rhythm, meter and structure

• Joint beat and downbeat tracking with neural activations and probabilistic decoding (Böck, Krebs and Widmer, 2016; Durand et al., 2017) is now standard.
• Open questions: expressive tempo and rubato, meters outside the Western 3/4 and 4/4, and the genuine ambiguity of the metrical level, which an F-measure against one annotation cannot capture.

Spatial audio

• Individualised HRTFs, predicted from morphology or images, and the capture and rendering of sound fields with higher-order Ambisonics (Zotter and Frank, 2019) and analytic sound field synthesis (Ahrens, 2012) are active topics.
• Open questions: perceptual evaluation protocols, and rendering for listeners who move.

Questions to explore

1. Can DUET be extended to reverberant rooms by replacing the single delay per source with a distribution of interchannel phase differences learned from data?
2. At which frequency separation, as a function of N and SNR, does ESPRIT stop resolving two equal-amplitude tones, and how does this compare with the Cramér–Rao bound for two tones?
3. Does a classifier trained on music features keep its accuracy when evaluated with strict artist filters on a new data set?
4. Can a phase vocoder with learned phase propagation remove phasiness without introducing transient smearing?

Exercises

Levels: ★ apply the course, ★★ combine results or prove a short statement, ★★★ open-ended or research-flavoured. Exercises 16, 17 and 18 are computational; exercise 19 is open-ended.

1. ★ Resolution and leakage. Two tones at 440 and 446 Hz are sampled at 8 kHz. (a) How many samples must a rectangular frame have for the gap to reach the Rayleigh limit? (b) And a Hann frame, judged by its main-lobe half-width? (c) Why does zero-padding the frame not help?

   Solution: (a) the limit is f_s / N \le 6 Hz, so N \ge 1334 samples (0.17 s). (b) The Hann main lobe is twice as wide, so N \ge 2667 samples (0.33 s).

   3. Zero-padding samples the same DTFT more densely; the width of the main lobe is set by the number of signal samples (Proposition 0.4).

2. ★ STFT parameters and overlap-add. At 44.1 kHz, take a periodic Hann window of N = 2048 samples and a hop R = 512. (a) Give the frame duration, the bin spacing and the number of frames per second. (b) Show that the shifted windows sum to 2 and their squares to 3/2. (c) Which synthesis window reconstructs exactly?

   Solution: (a) 46.4 ms, 21.5 Hz and 86.1 frames per second.

   2. w(n) = \tfrac12 - \tfrac12 \cos(2\pi n / N). Shifting by N/4 rotates the phase of the cosine by \pi/2, so the four shifted cosines cancel and the sum is 4 \times \tfrac12 = 2. Similarly w^2 = \tfrac38 - \tfrac12 \cos(2\pi n/N) + \tfrac18 \cos(4\pi n/N), and both cosines cancel over four shifts: the sum is 4 \times \tfrac38 = \tfrac32.

   3. By Proposition 1.2, w_s = \tfrac23 w.

3. ★ A mel filterbank by hand. (a) Compute the mel values of 1000, 4000 and 8000 Hz. (b) For 40 filters between 0 and 8000 Hz, find the centres of the first and last filters, and the widths between consecutive centres there. (c) Why is c_0 often discarded from the MFCCs?

   Solution: (a) 1000, 2146 and 2840 mel. (b) The 42 edges are 69.27 mel apart. The first centre is at 69.27 mel, that is 44.4 Hz, and the last at 2770.8 mel, that is 7481 Hz. Consecutive centres are 47 Hz apart at the bottom and 488 Hz apart at the top.

   3. c_0 is the sum of the log energies, which depends on the recording level rather than on timbre.

4. ★ Tempo from an autocorrelation lag. An onset function sampled at 100 frames per second has its weighted autocorrelation maximum at lag 47. (a) Give the tempo and the tempi of lags 46 and 48. (b) The unweighted autocorrelation also peaks at lag 94: which tempo does it suggest, and how does the prior of Definition 3.2 decide?

   Solution: (a) 127.7 BPM; lags 46 and 48 give 130.4 and 125.0 BPM, so the lag grid has a resolution of about 2.7 BPM here. Parabolic interpolation of the peak refines it.

   2. Lag 94 gives 63.8 BPM, the half tempo (Proposition 3.3). The log-Gaussian prior centred on 120 BPM weights lag 47 more, which selects 127.7 BPM.

5. ★ Flanger notches. A flanger at 44.1 kHz uses g = 0.7 and a delay swept between 1 and 5 ms. (a) Between which frequencies does the first notch move? (b) Give the depth of the notches and the height of the peaks in dB. (c) How many notches lie below 5 kHz when D = 5 ms?

   Solution: (a) the first notch is at f_s / 2D = 1 / (2 \times \text{delay}): from 500 Hz (1 ms) down to 100 Hz (5 ms).

   2. 20 \log_{10}(1 - g) = -10.5 dB and 20 \log_{10}(1 + g) = +4.6 dB (Proposition 4.2).

   3. Notches at 100, 300, 500, … Hz, odd multiples of 100 Hz, so 25 below 5 kHz.

6. ★ Tuning a Schroeder reverberator. At 44.1 kHz, compute the comb gains for T_{60} = 2 s and the delays 1557, 1617, 1491 and 1422 samples. Which mean absorption coefficient gives the same T_{60} in the room of ch. 4 (320 m³, 304 m² of walls)?

   Solution: by Proposition 4.9, g = 10^{-3D / (2 \times 44100)}, that is 0.885, 0.881, 0.890 and 0.895: longer delays need smaller gains.

   Sabine’s formula gives A = 0.161 \times 320 / 2 = 25.8 m², hence a mean absorption coefficient of \bar\alpha = 25.8 / 304 = 0.085 .

7. ★ Interaural time differences. With a = 8.75 cm and c = 343 m/s, compute the ITD at 45° in milliseconds and in samples at 48 kHz. Above which frequency does the interaural phase become ambiguous for a source at 90°?

   Solution: Definition 6.1 gives \mathrm{ITD}(45^\circ) = \frac{0.0875}{343} \Big(\frac{\pi}{4} + \frac{\sqrt2}{2}\Big) = 0.381 \text{ ms}, that is 18.3 samples. At 90° the ITD is 0.656 ms; the phase difference exceeds \pi when half a period is shorter than the ITD, above 1 / (2 \times 0.656 \text{ ms}) = 762 Hz. This is one reason why the ITD is used at low frequencies.

8. ★ Spacing a WFS array. (a) Which loudspeaker spacing avoids spatial aliasing up to 1500 Hz for every direction? (b) For frontal plane waves only? (c) How many loudspeakers does a 4 m array need in case (a)?

   Solution: by Proposition 6.5, (a) \Delta x \le c / (2 \times 1500) = 11.4 cm; (b) \Delta x \le c / 1500 = 22.9 cm. (c) 4 / 0.114 = 35 intervals, so 36 loudspeakers.

9. ★★ The Bayes rule over a texture window. (a) Show that summing frame log-likelihoods (Definition 2.2) is the Bayes rule when frames are independent given the class. (b) Two equiprobable classes, 100 frames, mean log-likelihood ratio 0.05 per frame: compute the posterior. (c) If the frames are in fact 25 independent groups of 4 identical frames, what is the correct posterior? Comment.

   Solution: (a) Under conditional independence the likelihood of the window is the product of the frame likelihoods p(\mathbf{x}_m \mid c); take the logarithm and add \log p(c).

   2. The log-odds are 100 \times 0.05 = 5, so the posterior is e^5 / (1 + e^5) = 0.993. (c) Only 25 independent observations: log-odds 1.25 and posterior 0.78.

   The decision is the same, but treating overlapping frames as independent makes the classifier overconfident, which matters when posteriors are combined or thresholded.

10. ★★ Distortion, harmonics and aliasing. (a) Show that f(x) = x + x^2 / 2 adds a DC offset and a second harmonic to A\cos\theta. (b) A 9 kHz tone passes through x \mapsto x^3 at 44.1 kHz: where does its third harmonic appear? (c) Which smallest integer oversampling factor guarantees that, for every input tone below 20 kHz, the alias of the third harmonic lands above 20 kHz, where the decimation filter removes it?

    Solution: (a) The square term gives \tfrac12 A^2 \cos^2\theta = \tfrac14 A^2 + \tfrac14 A^2 \cos 2\theta, a DC offset and a second harmonic: the even part of a nonlinearity produces even harmonics, consistent with Proposition 4.7.

    2. At 27 kHz, which aliases to 44.1 - 27 = 17.1 kHz, an inharmonic component (1.9 times the fundamental).

    3. At rate M \times 44.1 kHz, a third harmonic of at most 60 kHz aliases to 44.1 M - 60 kHz. This is at least 20 kHz when M \ge 1.81, so M = 2 suffices.

11. ★★ The phase vocoder hop. With a Hann window of N = 2048 samples, (a) give the largest analysis hop allowed by Proposition 5.4. (b) With R_a = 1024, a sinusoid lies 1.5 bins above the centre of bin k. What frequency does the phase vocoder estimate for bin k?

    Solution: (a) R_a \le N / 4 = 512.

    2. The phase offset is 1.5 \times (2\pi / 2048) \times 1024 = 1.5\pi, which \operatorname{princarg} maps to -0.5\pi. The estimate is therefore \omega_k - 0.5\pi / 1024, that is 0.5 bin below the centre instead of 1.5 bins above: an error of 2 bins, 43 Hz at 44.1 kHz.

12. ★★ Whitening and kurtosis. (a) Prove that \mathbf{Q} = \mathbf{B}\mathbf{A} of Proposition 7.3 is orthogonal using \mathbf{R}_x = \mathbf{A}\mathbf{A}^\top. (b) Compute the kurtosis of a uniform variable on [-\sqrt3, \sqrt3] and of a Laplace variable of unit variance. (c) Which sign of kurtosis do speech and a sinusoid with random phase have?

    Solution: (a) Whiteness of \mathbf{z} means \mathbf{B}\mathbf{R}_x\mathbf{B}^\top = \mathbf{I}, that is \mathbf{Q}\mathbf{Q}^\top = \mathbf{I}.

    2. Uniform: \mathbb{E}[y^4] = 9/5, so \kappa = -1.2. Laplace with unit variance: \mathbb{E}[y^4] = 6, so \kappa = 3.

    3. Speech is super-Gaussian (\kappa > 0). A sinusoid \sqrt2 \cos\Theta has \mathbb{E}[y^4] = 3/2, so \kappa = -1.5: sub-Gaussian. FastICA handles both signs (Proposition 7.6).

13. ★★ DUET and microphone spacing. Two microphones 5 cm apart record at 16 kHz. (a) What is the largest delay in samples? (b) Above which frequency can delays wrap? (c) Compute the symmetric attenuation for a = 2 and for a = 1/2, and explain why DUET uses it.

    Solution: (a) d f_s / c = 0.05 \times 16000 / 343 = 2.33 samples, more than the bound of Proposition 8.5.

    2. The phase \omega\delta exceeds \pi above c / (2d) = 3430 Hz. Estimate delays below that frequency only, or use microphones less than 2.1 cm apart.

    3. \alpha = 2 - 1/2 = 1.5 and \alpha = 1/2 - 2 = -1.5. Exchanging the microphones changes a into 1/a and only flips the sign of \alpha, so the histogram treats both channels symmetrically.

14. ★★ Crosstalk cancellation in free field. Model the contralateral path as G_c = \rho\, G_i\, e^{-i\omega\tau} with 0 < \rho < 1. (a) Compute \det \mathbf{G} and the frequencies where its modulus is minimal. (b) For \rho = 0.95, by how much must the inverse filter amplify the signal there?

    Solution: (a) \det \mathbf{G} = G_i^2 (1 - \rho^2 e^{-2i\omega\tau}). Its modulus is minimal, equal to |G_i|^2 (1 - \rho^2), when 2\omega\tau is a multiple of 2\pi, that is at f = m / (2\tau), including f = 0.

    2. At those frequencies the matrix has the eigenvalues G_i (1 \pm \rho): \mathbf{G} = G_i \begin{pmatrix} 1 & \pm\rho \\ \pm\rho & 1 \end{pmatrix} .

    The inverse therefore amplifies some binaural signals by 1 / ((1 - \rho)|G_i|). That is 1 / (1 - \rho) = 20 times, or 26 dB, more than the gain 1 / |G_i| needed without crosstalk. The determinant alone contributes the factor 1 / (1 - \rho^2) = 10.3, about 20 dB.

    The code of ch. 6 (\rho \approx 0.96) finds only 20–23 dB of extra gain at 100 and 1900 Hz, because its Tikhonov regularisation caps the gain.

15. ★★ ESPRIT by hand. Let x(n) = \cos(\pi n / 2) for n = 0, \dots, 4. Build the 3 \times 3 Hankel matrix, an orthonormal basis of its column space, the matrix \boldsymbol{\Phi} of Theorem 10.4 and its eigenvalues; then recover the amplitudes.

    Solution: x = (1, 0, -1, 0, 1), and the Hankel matrix has rows (1, 0, -1), (0, -1, 0) and (-1, 0, 1), of rank 2. Its column space has the orthonormal basis (1, 0, -1)/\sqrt2 and (0, 1, 0).

    Removing the last row, then the first, gives \mathbf{W}_\downarrow = \begin{pmatrix} 1/\sqrt2 & 0 \\ 0 & 1 \end{pmatrix}, \quad \mathbf{W}_\uparrow = \begin{pmatrix} 0 & 1 \\ -1/\sqrt2 & 0 \end{pmatrix}, \quad \boldsymbol{\Phi} = \mathbf{W}_\downarrow^{-1} \mathbf{W}_\uparrow = \begin{pmatrix} 0 & \sqrt2 \\ -1/\sqrt2 & 0 \end{pmatrix} . Its characteristic polynomial is \lambda^2 + 1, so the poles are \pm i = e^{\pm i\pi/2}. From x(0) = \alpha_1 + \alpha_2 = 1 and x(1) = i(\alpha_1 - \alpha_2) = 0, both amplitudes are 1/2.

16. ★★★ Beat tracking through a tempo change (Python). Modify the code of ch. 3 so that the tempo glides linearly from 110 to 140 BPM over 48 beats. Run the global-tempo tracker, then estimate the beat period on 4 s windows every second and let \tau in Proposition 3.5 depend on the frame.

    Hints and expected result: the global autocorrelation is smeared by the glide and the prior then picks the eighth-note level, 240 BPM, at the edge of the search range. The tracker returns about twice too many beats (96–98 for 48, F-measure about 0.65).

    With the local beat period (lags falling from about 53 to 44 frames), it finds 48–49 beats with an F-measure of about 0.98. The recursion of Proposition 3.5 remains exact with a frame-dependent \tau.

17. ★★★ The resolution threshold of ESPRIT (Python). For two equal-amplitude complex tones with random phases, N = 64, L = 21, and gaps of 0.25, 0.5, 0.75, 1 and 1.5 times 1/N, estimate over 200 draws how often ESPRIT and the two highest peaks of the zero-padded periodogram place both frequencies within a quarter of the gap. Use SNRs of 20 and 40 dB.

    Hints and expected result: at 20 dB, ESPRIT succeeds in about 95 % of the draws at a gap of 0.5/N, and the periodogram needs about 1.5/N (57 % at 1/N).

    At 40 dB, ESPRIT resolves 0.25/N in about 98 % of the draws, while the periodogram does not improve. Its limit is resolution, not noise (Proposition 9.2); the limit of ESPRIT moves with the SNR.

18. ★★★ NMF and Wiener masks (Python). Using Definition 11.6 and Proposition 8.6, simulate two sources with complex Gaussian STFT coefficients whose variances are one harmonic template times sparse random activations. Fix the templates, estimate the activations on the mixture by the multiplicative update for \mathbf{H}, apply Wiener masks, and compute the SDR of source 1.

    Take the sizes of the NMF code of ch. 11: 64 bins, 100 frames, templates floored at 10^{-3}, activations non-zero in 30 % of the frames. Compare templates with distinct partials (bins 5, 10, 15, 20 and 7, 14, 21, 28) with an octave (bins 10, 20, 30, 40 for source 2).

    Hints and expected result: with distinct partials, the SDR is about 25 dB; with the octave, about 12 dB, close to the 13 dB obtained with the true activations. Denser activations lower both values (about 22 and 7 dB at 80 %). Shared bins receive a mask between 0 and 1 even when the template gains are perfect, because the Wiener filter can only split their energy in proportion to the variances.

19. ★★★ DUET in a reverberant room (open). Replace the anechoic mixture of the DUET code by convolutions with two different synthetic room responses per source (for instance Schroeder reverberators from ch. 4 with different delays), for T_{60} from 0 to 0.6 s. Measure how the (\alpha, \delta) clusters and the SDRs degrade, and propose a remedy.

    Hints: as T_{60} grows, the ratios X_2 / X_1 at single-source points no longer have one modulus and one phase slope, so the clusters spread and the masks mix sources. Possible remedies are longer windows (Proposition 7.7), probabilistic models of interchannel phase and level differences fitted by EM (Mandel, Weiss and Ellis, 2010), or the multichannel local Gaussian model of Duong, Vincent and Gribonval (2010). This is a natural extension of the first solo project below.

Projects and assessment

The official assessment, per master-mva.com, is the reading and analysis of research papers, with a written report and an oral presentation; the Université Paris-Saclay course detail classifies it as “Evaluation Continue non Intégrale” (partial continuous assessment). The solo and group subjects below are paper studies sized for that format. The written exam is not part of the official assessment: it is a self-test proposal, and none of these proposals is official material.

Solo projects

1. [Solo · ≈ 20 h] DUET and W-disjoint orthogonality: a paper study. Reproduce the main experiments of Yilmaz and Rickard (2004).
   • Implement: DUET from Algorithm 8.4, with histogram peak picking, and the approximate W-disjoint orthogonality measures of the paper.
   • Data: anechoic mixtures of 2 to 5 speakers simulated from public speech recordings such as LibriSpeech (Panayotov et al., 2015), with chosen attenuations and sub-sample delays.
   • Success criterion: SDR, SIR and SAR (Definition 7.8) as a function of the number of sources and of the window length, reproducing the paper’s trends.
   • Deliverables: code, a 4–6 page report, a 10-minute talk. Chapters 1, 7, 8.
2. [Solo · ≈ 20 h] ESPRIT and ESTER on musical notes: a paper study. Study Roy and Kailath (1989) and the ESTER criterion of Badeau, David and Richard (2006).
   • Implement: Algorithm 10.5, ESTER and an MDL baseline; a Monte Carlo study against the Cramér–Rao bound (Theorem 9.3).
   • Data: synthetic damped sinusoids, then isolated acoustic keyboard notes from the public NSynth data set (Engel et al., 2017); for piano notes the partial frequencies follow f_k = k f_0 \sqrt{1 + B k^2}.
   • Success criterion: correct order in at least 90 % of synthetic draws at 20 dB, frequency RMSE within a factor 1.5 of the bound, and a fitted inharmonicity coefficient B for several piano notes.
   • Deliverables: code, report, talk. Chapters 9, 10, 11.
3. [Solo · ≈ 15 h] Beat tracking by dynamic programming: a paper study. Reimplement Ellis (2007) and evaluate it.
   • Data: the public Ballroom data set of the tempo comparison of Gouyon et al. (2006), with the beat annotations released by Krebs, Böck and Widmer (2013).
   • Work: the onset function, tempo estimation with and without the prior, and the effect of the tightness \alpha; F-measures computed with mir_eval (Raffel et al., 2014).
   • Success criterion: a table of F-measures per dance style, with and without the tempo prior; a higher mean F-measure with the prior; and an error analysis that separates metrical-level errors (half or double tempo) from phase errors (off-beat tracking).
   • Deliverables: code, report, talk. Chapter 3.

Group projects

1. [Group of 2–3 · ≈ 30 h each] Instrument recognition with an honest evaluation. Build and compare instrument-family classifiers on the public NSynth data set (Engel et al., 2017).
   • Work split: one member writes the features of ch. 1 (MFCC statistics, descriptors of Definition 2.6) and the GMM and SVM baselines; one trains a small convolutional network on log-mel spectrograms; one builds the evaluation, with splits that keep each instrument on one side (Proposition 2.7), confusion matrices and error bars.
   • Success criterion: the accuracy gap between random and instrument-disjoint splits, measured for every model, and a model that beats the baselines on the disjoint split.
   • Deliverables: code, a 6-page report, a 15-minute talk. Chapters 1, 2, 11.
2. [Group of 2–3 · ≈ 30 h each] A binaural and transaural renderer. Render moving virtual sources from the public CIPIC HRTF database (Algazi et al., 2001).
   • Work split: binaural synthesis with interpolation between measured directions; a crosstalk canceller with the regularisation of Kirkeby et al. (1998), studied against \beta and head displacement; a listening test with at least five listeners.
   • Success criterion: simulated channel separation of at least 20 dB between 200 Hz and 6 kHz at the sweet spot, and a localisation-error table from the listening test.
   • Deliverables: code, report, demo, talk. Chapters 6, 11.
3. [Group of 2–3 · ≈ 30 h each] An effects and reverberation toolkit: a paper study. Reproduce the Schroeder reverberator (Schroeder, 1962) and a feedback delay network (Jot and Chaigne, 1991), and add the effects of ch. 4–5.
   • Work split: flanger, phaser and an oversampled distortion with measured aliasing; an FDN reverberator with frequency-dependent decay; a phase vocoder with the phase locking of Laroche and Dolson (1999).
   • Success criterion: reverberation times within 5 % of the target in octave bands; aliasing products below −60 dB with oversampling; a blind listening comparison.
   • Deliverables: code, report, sound examples, talk. Chapters 4, 5, 11.

Written exam (proposal)

Format: 2 hours, one handwritten A4 sheet allowed, no electronic device; four independent exercises.

1. Course question (4 points). Define the STFT. State and prove the overlap-add reconstruction condition, and give a window and a hop that satisfy it.

   Expected: Definition 1.1; Proposition 1.2 with its proof (each frame contributes x(n)\, w\, w_s at time n); the periodic Hann window with R = N/4, whose shifted squares sum to 3/2, so w_s = \tfrac23 w.

2. Computation (6 points). (a) A feedforward comb at 48 kHz has g = 0.8 and a delay of 2.5 ms: give its notch frequencies, notch depth and peak height. (b) Which gain gives a feedback comb of 1500 samples a reverberation time of 1.2 s? (c) Which loudspeaker spacing makes WFS alias-free up to 1 kHz in every direction?

   Expected: (a) notches at 200, 600, 1000, … Hz, depth −14.0 dB, peaks +5.1 dB (Proposition 4.2). (b) g = 10^{-3 \times 1500 / (1.2 \times 48000)} = 0.835 (Proposition 4.9). (c) c / (2 \times 1000) = 17.2 cm (Proposition 6.5).

3. Short proof (5 points). Let \mathbf{W} be a basis of the signal subspace of a noiseless sum of K distinct complex exponentials, with L \ge K + 1 rows. Prove that the eigenvalues of the matrix \boldsymbol{\Phi} solving the shift equation are the poles.

   Expected: Theorem 10.4: the shift invariance of the Vandermonde matrix, \mathbf{W} = \mathbf{V}_L \mathbf{G}, the full column rank of the truncated Vandermonde matrix, and the similarity of \boldsymbol{\Phi} with the diagonal matrix of the poles.

4. Problem (5 points). In an anechoic two-microphone mixture, at a point where \omega_k = \pi/4 rad per sample, the STFT ratio is X_2 / X_1 = 0.6\, e^{-i\pi/8}. (a) Assuming one active source, give its attenuation, delay and symmetric attenuation. (b) Can this delay wrap below the Nyquist frequency? (c) If two Gaussian sources of variances 3 and 1 are active instead, give the MMSE estimate of the first.

   Expected: (a) a = 0.6, \delta = (\pi/8) / (\pi/4) = 0.5 sample, \alpha = 0.6 - 1/0.6 = -1.07 (Proposition 8.3). (b) No: |\omega\delta| \le \pi/2 up to Nyquist (Proposition 8.5). (c) The Wiener estimate \tfrac34 X_1 (Proposition 8.6).

Grading rubric (proposal)

For the paper studies, which are the official assessment format:

References

Official references (master-mva.com)

• Y. Grenier, R. Badeau and G. Richard. Polycopiés de cours sur le traitement du signal audio (course notes on audio signal processing, in French). Télécom ParisTech (no year given on the course page).
• M. Müller, D. P. W. Ellis, A. Klapuri and G. Richard. Signal processing for music analysis. IEEE Journal of Selected Topics in Signal Processing, 5(6):1088–1110, October 2011.
• Course page linked from the official description (“Plus d’information”): https://ecampus.paris-saclay.fr/course/view.php?id=21210.

Textbooks

• J. Ahrens. Analytic Methods of Sound Field Synthesis. Springer, 2012.
• J. Blauert. Spatial Hearing: The Psychophysics of Human Sound Localization. Revised edition, MIT Press, 1997.
• P. Comon and C. Jutten (eds.). Handbook of Blind Source Separation: Independent Component Analysis and Applications. Academic Press, 2010.
• A. Hyvärinen, J. Karhunen and E. Oja. Independent Component Analysis. Wiley, 2001.
• S. M. Kay. Modern Spectral Estimation: Theory and Application. Prentice Hall, 1988.
• A. Klapuri and M. Davy (eds.). Signal Processing Methods for Music Transcription. Springer, 2006.
• M. Müller. Fundamentals of Music Processing. Springer, 2015.
• A. V. Oppenheim and R. W. Schafer. Discrete-Time Signal Processing. 3rd edition, Pearson, 2010.
• P. Stoica and R. Moses. Spectral Analysis of Signals. Prentice Hall, 2005.
• E. Vincent, T. Virtanen and S. Gannot (eds.). Audio Source Separation and Speech Enhancement. Wiley, 2018.
• F. Zotter and M. Frank. Ambisonics: A Practical 3D Audio Theory for Recording, Studio Production, Sound Reinforcement, and Virtual Reality. Springer, 2019.
• U. Zölzer (ed.). DAFX: Digital Audio Effects. 2nd edition, Wiley, 2011.

Features and indexing (ch. 1–2)

• S. B. Davis and P. Mermelstein. Comparison of parametric representations for monosyllabic word recognition in continuously spoken sentences. IEEE Transactions on Acoustics, Speech, and Signal Processing, 28(4):357–366, 1980.
• S. Essid, G. Richard and B. David. Instrument recognition in polyphonic music based on automatic taxonomies. IEEE Transactions on Audio, Speech, and Language Processing, 14(1):68–80, 2006.
• D. W. Griffin and J. S. Lim. Signal estimation from modified short-time Fourier transform. IEEE Transactions on Acoustics, Speech, and Signal Processing, 32(2):236–243, 1984.
• G. Richard, S. Sundaram and S. Narayanan. An overview on perceptually motivated audio indexing and classification. Proceedings of the IEEE, 101(9):1939–1954, 2013.
• B. L. Sturm. The state of the art ten years after a state of the art: future research in music information retrieval. Journal of New Music Research, 43(2):147–172, 2014.
• G. Tzanetakis and P. Cook. Musical genre classification of audio signals. IEEE Transactions on Speech and Audio Processing, 10(5):293–302, 2002.

Rhythm (ch. 3)

• J. P. Bello, L. Daudet, S. Abdallah, C. Duxbury, M. Davies and M. B. Sandler. A tutorial on onset detection in music signals. IEEE Transactions on Speech and Audio Processing, 13(5):1035–1047, 2005.
• S. Böck, F. Krebs and G. Widmer. Joint beat and downbeat tracking with recurrent neural networks. Proceedings of the International Society for Music Information Retrieval Conference (ISMIR), 2016.
• S. Durand, J. P. Bello, B. David and G. Richard. Robust downbeat tracking using an ensemble of convolutional networks. IEEE/ACM Transactions on Audio, Speech, and Language Processing, 25(1):76–89, 2017.
• D. P. W. Ellis. Beat tracking by dynamic programming. Journal of New Music Research, 36(1):51–60, 2007.
• F. Gouyon, A. Klapuri, S. Dixon, M. Alonso, G. Tzanetakis, C. Uhle and P. Cano. An experimental comparison of audio tempo induction algorithms. IEEE Transactions on Audio, Speech, and Language Processing, 14(5):1832–1844, 2006.
• P. Grosche and M. Müller. Extracting predominant local pulse information from music recordings. IEEE Transactions on Audio, Speech, and Language Processing, 19(6):1688–1701, 2011.
• A. P. Klapuri, A. J. Eronen and J. T. Astola. Analysis of the meter of acoustic musical signals. IEEE Transactions on Audio, Speech, and Language Processing, 14(1):342–355, 2006.
• F. Krebs, S. Böck and G. Widmer. Rhythmic pattern modeling for beat and downbeat tracking in musical audio. Proceedings of ISMIR, 2013.
• C. Raffel, B. McFee, E. J. Humphrey, J. Salamon, O. Nieto, D. Liang and D. P. W. Ellis. mir_eval: a transparent implementation of common MIR metrics. Proceedings of ISMIR, 2014.
• E. D. Scheirer. Tempo and beat analysis of acoustic musical signals. Journal of the Acoustical Society of America, 103(1):588–601, 1998.

Effects and modifications (ch. 4–5)

• J. L. Flanagan and R. M. Golden. Phase vocoder. Bell System Technical Journal, 45(9):1493–1509, 1966.
• W. G. Gardner. Efficient convolution without input–output delay. Journal of the Audio Engineering Society, 43(3):127–136, 1995.
• J.-M. Jot and A. Chaigne. Digital delay networks for designing artificial reverberators. 90th Convention of the Audio Engineering Society, Paris, 1991.
• J. Laroche and M. Dolson. Improved phase vocoder time-scale modification of audio. IEEE Transactions on Speech and Audio Processing, 7(3):323–332, 1999.
• J. Makhoul. Linear prediction: a tutorial review. Proceedings of the IEEE, 63(4):561–580, 1975.
• E. Moulines and F. Charpentier. Pitch-synchronous waveform processing techniques for text-to-speech synthesis using diphones. Speech Communication, 9(5–6):453–467, 1990.
• M. R. Schroeder. Natural sounding artificial reverberation. Journal of the Audio Engineering Society, 10(3):219–223, 1962.
• V. Välimäki, J. D. Parker, L. Savioja, J. O. Smith and J. S. Abel. Fifty years of artificial reverberation. IEEE Transactions on Audio, Speech, and Language Processing, 20(5):1421–1448, 2012.
• W. Verhelst and M. Roelands. An overlap-add technique based on waveform similarity (WSOLA) for high quality time-scale modification of speech. Proceedings of ICASSP, 1993.

Spatial audio (ch. 6)

• V. R. Algazi, R. O. Duda, D. M. Thompson and C. Avendano. The CIPIC HRTF database. Proceedings of the IEEE Workshop on Applications of Signal Processing to Audio and Acoustics (WASPAA), 2001.
• A. J. Berkhout, D. de Vries and P. Vogel. Acoustic control by wave field synthesis. Journal of the Acoustical Society of America, 93(5):2764–2778, 1993.
• W. G. Gardner and K. D. Martin. HRTF measurements of a KEMAR. Journal of the Acoustical Society of America, 97(6):3907–3908, 1995.
• M. A. Gerzon. Periphony: with-height sound reproduction. Journal of the Audio Engineering Society, 21(1):2–10, 1973.
• O. Kirkeby, P. A. Nelson, H. Hamada and F. Orduña-Bustamante. Fast deconvolution of multichannel systems using regularization. IEEE Transactions on Speech and Audio Processing, 6(2):189–194, 1998.
• V. Pulkki. Virtual sound source positioning using vector base amplitude panning. Journal of the Audio Engineering Society, 45(6):456–466, 1997.
• S. Spors, R. Rabenstein and J. Ahrens. The theory of wave field synthesis revisited. 124th Convention of the Audio Engineering Society, Amsterdam, 2008.
• R. S. Woodworth. Experimental Psychology. Holt, 1938.

Source separation (ch. 7–8)

• A. J. Bell and T. J. Sejnowski. An information-maximization approach to blind separation and blind deconvolution. Neural Computation, 7(6):1129–1159, 1995.
• P. Bofill and M. Zibulevsky. Underdetermined blind source separation using sparse representations. Signal Processing, 81(11):2353–2362, 2001.
• P. Comon. Independent component analysis, a new concept? Signal Processing, 36(3):287–314, 1994.
• N. Q. K. Duong, E. Vincent and R. Gribonval. Under-determined reverberant audio source separation using a full-rank spatial covariance model. IEEE Transactions on Audio, Speech, and Language Processing, 18(7):1830–1840, 2010.
• A. Hyvärinen and E. Oja. Independent component analysis: algorithms and applications. Neural Networks, 13(4–5):411–430, 2000.
• A. Jourjine, S. Rickard and Ö. Yilmaz. Blind separation of disjoint orthogonal signals: demixing N sources from 2 mixtures. Proceedings of ICASSP, 2000.
• M. I. Mandel, R. J. Weiss and D. P. W. Ellis. Model-based expectation-maximization source separation and localization. IEEE Transactions on Audio, Speech, and Language Processing, 18(2):382–394, 2010.
• P. Smaragdis. Blind separation of convolved mixtures in the frequency domain. Neurocomputing, 22(1–3):21–34, 1998.
• E. Vincent, R. Gribonval and C. Févotte. Performance measurement in blind audio source separation. IEEE Transactions on Audio, Speech, and Language Processing, 14(4):1462–1469, 2006.
• Ö. Yilmaz and S. Rickard. Blind separation of speech mixtures via time-frequency masking. IEEE Transactions on Signal Processing, 52(7):1830–1847, 2004.

High-resolution spectral analysis (ch. 9–10)

• R. Badeau, B. David and G. Richard. Fast approximated power iteration subspace tracking. IEEE Transactions on Signal Processing, 53(8):2931–2941, 2005.
• R. Badeau, B. David and G. Richard. A new perturbation analysis for signal enumeration in rotational invariance techniques. IEEE Transactions on Signal Processing, 54(2):450–458, 2006.
• R. Badeau, G. Richard and B. David. Performance of ESPRIT for estimating mixtures of complex exponentials modulated by polynomials. IEEE Transactions on Signal Processing, 56(2):492–504, 2008.
• R. J. McAulay and T. F. Quatieri. Speech analysis/synthesis based on a sinusoidal representation. IEEE Transactions on Acoustics, Speech, and Signal Processing, 34(4):744–754, 1986.
• V. F. Pisarenko. The retrieval of harmonics from a covariance function. Geophysical Journal of the Royal Astronomical Society, 33(3):347–366, 1973.
• D. C. Rife and R. R. Boorstyn. Single-tone parameter estimation from discrete-time observations. IEEE Transactions on Information Theory, 20(5):591–598, 1974.
• R. Roy and T. Kailath. ESPRIT — estimation of signal parameters via rotational invariance techniques. IEEE Transactions on Acoustics, Speech, and Signal Processing, 37(7):984–995, 1989.
• R. O. Schmidt. Multiple emitter location and signal parameter estimation. IEEE Transactions on Antennas and Propagation, 34(3):276–280, 1986.
• X. Serra and J. O. Smith. Spectral modeling synthesis: a sound analysis/synthesis system based on a deterministic plus stochastic decomposition. Computer Music Journal, 14(4):12–24, 1990.
• M. Wax and T. Kailath. Detection of signals by information theoretic criteria. IEEE Transactions on Acoustics, Speech, and Signal Processing, 33(2):387–392, 1985.

Complements (ch. 11)

• A. Défossez, N. Usunier, L. Bottou and F. Bach. Music source separation in the waveform domain. arXiv:1911.13254, 2019.
• C. Févotte, N. Bertin and J.-L. Durrieu. Nonnegative matrix factorization with the Itakura–Saito divergence: with application to music analysis. Neural Computation, 21(3):793–830, 2009.
• D. FitzGerald. Harmonic/percussive separation using median filtering. Proceedings of the International Conference on Digital Audio Effects (DAFx), 2010.
• S. Hershey, S. Chaudhuri, D. P. W. Ellis, J. F. Gemmeke, A. Jansen, R. C. Moore, M. Plakal, D. Platt, R. A. Saurous, B. Seybold, M. Slaney, R. J. Weiss and K. Wilson. CNN architectures for large-scale audio classification. Proceedings of ICASSP, 2017.
• Q. Kong, Y. Cao, T. Iqbal, Y. Wang, W. Wang and M. D. Plumbley. PANNs: large-scale pretrained audio neural networks for audio pattern recognition. IEEE/ACM Transactions on Audio, Speech, and Language Processing, 28:2880–2894, 2020.
• D. D. Lee and H. S. Seung. Algorithms for non-negative matrix factorization. Advances in Neural Information Processing Systems 13, 2001.
• Y. Luo and N. Mesgarani. Conv-TasNet: surpassing ideal time–frequency magnitude masking for speech separation. IEEE/ACM Transactions on Audio, Speech, and Language Processing, 27(8):1256–1266, 2019.
• A. Ozerov and C. Févotte. Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation. IEEE Transactions on Audio, Speech, and Language Processing, 18(3):550–563, 2010.

Research directions (ch. 12) and data sets

• A. Baevski, Y. Zhou, A. Mohamed and M. Auli. wav2vec 2.0: a framework for self-supervised learning of speech representations. Advances in Neural Information Processing Systems 33, 2020.
• E. J. Candès and C. Fernandez-Granda. Towards a mathematical theory of super-resolution. Communications on Pure and Applied Mathematics, 67(6):906–956, 2014.
• B. Elizalde, S. Deshmukh, M. Al Ismail and H. Wang. CLAP: learning audio concepts from natural language supervision. Proceedings of ICASSP, 2023.
• J. Engel, L. Hantrakul, C. Gu and A. Roberts. DDSP: differentiable digital signal processing. International Conference on Learning Representations (ICLR), 2020.
• J. Engel, C. Resnick, A. Roberts, S. Dieleman, D. Eck, K. Simonyan and M. Norouzi. Neural audio synthesis of musical notes with WaveNet autoencoders. Proceedings of the International Conference on Machine Learning (ICML), 2017 (the NSynth data set).
• J. F. Gemmeke, D. P. W. Ellis, D. Freedman, A. Jansen, W. Lawrence, R. C. Moore, M. Plakal and M. Ritter. Audio Set: an ontology and human-labeled dataset for audio events. Proceedings of ICASSP, 2017.
• I. Kavalerov, S. Wisdom, H. Erdogan, B. Patton, K. Wilson, J. Le Roux and J. R. Hershey. Universal sound separation. Proceedings of WASPAA, 2019.
• J. Le Roux, S. Wisdom, H. Erdogan and J. R. Hershey. SDR — half-baked or well done? Proceedings of ICASSP, 2019.
• V. Panayotov, G. Chen, D. Povey and S. Khudanpur. LibriSpeech: an ASR corpus based on public domain audio books. Proceedings of ICASSP, 2015.
• G. Tang, B. N. Bhaskar, P. Shah and B. Recht. Compressed sensing off the grid. IEEE Transactions on Information Theory, 59(11):7465–7490, 2013.
• S. Wisdom, E. Tzinis, H. Erdogan, R. J. Weiss, K. Wilson and J. R. Hershey. Unsupervised sound separation using mixture invariant training. Advances in Neural Information Processing Systems 33, 2020.