Classical Audio Classification

Chapter 1: Framing — Chopping Sound Into Stillness

Chapter 0 left us needing a fixed-length, phase-insensitive description of a sound. Before we can compute any such description, we have to answer a question that sounds trivial and is not: over what stretch of time?

Compute the average brightness of a whole four-second clip and you get one number describing a sentence, a siren sweep, and the bus that drove past — all averaged together into mush. Compute it over one sample and it is meaningless (a single number has no spectrum). The answer is in between, and the reason a specific answer exists is physical.

Sound holds still, briefly

A vocal tract is a tube of flesh. Its shape changes when you move your tongue and jaw, and flesh has mass, so it cannot change quickly — the articulators move at roughly ten significant changes per second. A guitar string’s resonance, a car engine’s cycle, a room’s reverberation: all change on a scale of tens of milliseconds or slower.

This gives us quasi-stationarity: over a window of about 20 to 40 milliseconds, the statistical character of a sound — its spectrum — is approximately constant. Not the samples (those oscillate thousands of times per second), the spectrum. That is the assumption the whole classical pipeline stands on, and it is a good one for speech, most music, and most machine sounds.

So we chop. We cut the signal into overlapping frames, treat each frame as one still photograph of the sound, and compute features from it. A four-second clip becomes a few hundred photographs. This step is called framing or blocking, and it converts a one-dimensional stream into a two-dimensional array of shape (frames, samples-per-frame).

The three numbers

Framing has exactly three parameters, and every audio toolkit exposes them under slightly different names:

The 25 ms / 10 ms pair is so standard that you will see it in Whisper’s front end, in Kaldi, in every ASR paper from 1985 to now. Environmental sound systems often use longer frames (40–50 ms) because the events are less rapid than phonemes; music systems sometimes use 100 ms or more.

Worked by hand: framing an 8-sample signal

Take an absurdly small signal so every index is visible. Eight samples:

Choose a frame length of 4 samples and a hop of 2 samples (that is 50% overlap). Frame i starts at sample index i × hop and covers hop-length onward:

Enumerate them:

Three frames. Notice that samples x₂ and x₃ appear in frames 0 and 1: with 50% overlap, every interior sample is used twice. That is deliberate and we will justify it in a moment.

The frame-count formula, derived

How many frames does a signal of N samples give? The last valid frame must start at an index s with s + L ≤ N, that is s ≤ N − L. Starts are multiples of h, so the largest usable one is h × ⌊(N − L)/h⌋, and counting from i = 0 gives

Check against our toy: 1 + ⌊(8 − 4)/2⌋ = 1 + ⌊2⌋ = 3. Correct.

Now the real case. A 4-second clip at 16 kHz is N = 64,000 samples; L = 400, h = 160:

So the array shape after framing is (398, 400) — 159,200 numbers, more than the 64,000 we started with, because of the overlap. Framing expands the data; the features that come next are what shrink it, from 400 numbers per frame down to about 39.

The frame rate is sr / h = 16000 / 160 = 100 frames per second. Every classical speech system runs at exactly this rate, which is why “frames” and “centiseconds” are used interchangeably in that literature.

Why overlap, really

The obvious objection: overlapping frames duplicate samples and cost 2.5× the computation. Why not hop by a full frame?

Two reasons, and the second is the real one. First, boundary events. A plosive consonant lasts about 5 ms. With non-overlapping 25 ms frames, a plosive that straddles a boundary is split into two halves, and neither half looks like a plosive. With 50–60% overlap, some frame always contains it whole.

Second, windowing throws away the edges. The window function we are about to introduce multiplies the first and last samples of every frame by zero. Without overlap, those samples contribute nothing to any frame — you have literally deleted parts of the signal. With 50% overlap, a sample that sits at the tapered edge of one frame sits near the centre of the next, so nothing is lost.

The window function, and the discontinuity it fixes

Cutting a frame out of a signal with a hard edge is itself a violent operation. The Fourier transform assumes the frame repeats forever; if the last sample is +0.9 and the first is −0.5, the implied repetition has a jump discontinuity every 400 samples. A jump is a click, and a click is broadband — it splatters energy across every frequency bin. This artefact is called spectral leakage, and it makes a clean tone look like a smeared hill.

The fix is to taper the frame to zero at both ends before transforming, by multiplying element-wise with a window function. The most common is the Hann window:

Worked by hand: an 8-point Hann window

With L = 8, the argument is 2πn/7. Compute each cosine and each weight:

Symmetric, zero at both ends, one at the centre. Now apply it to our 8-sample frame, multiplying element-wise:

The energy of the windowed frame is 0.0000 + 0.0128 + 0.0598 + 0.5782 + 0.0361 + 0.3026 + 0.0004 + 0.0000 = 0.990, versus 2.280 for the raw frame. Windowing removed 57% of the energy. That is expected and harmless as long as you are consistent — but it is exactly why you must never compare a windowed energy feature against an unwindowed one, and why libraries expose a normalization option.

Frequency resolution falls out of the frame length

Once you fix L you have fixed how finely you can resolve frequency. An L-sample frame at sample rate sr spans L/sr seconds, and the Fourier transform of that frame has bins spaced

For our speech setting: 16000 / 400 = 40 Hz per bin. Two tones 20 Hz apart cannot be separated in a 25 ms frame no matter what you do afterwards — not by zero-padding, not by a bigger FFT. Zero padding to 512 points makes the bin spacing 31.25 Hz, which interpolates the same underlying 40 Hz resolution more smoothly; it does not create information.

Code: three forms

python
# FORM 1 - the loop you would write on a whiteboard. Explicit, obviously correct.
import numpy as np

def frame_signal(x, L, h):
    n_frames = 1 + (len(x) - L) // h        # the formula we derived
    out = np.zeros((n_frames, L))
    for i in range(n_frames):
        out[i] = x[i*h : i*h + L]              # slice, copy
    return out

# FORM 2 - the same thing with zero copying, using stride tricks.
# Each row is a VIEW into x, offset by h samples. Memory: O(len(x)), not O(frames*L).
def frame_strided(x, L, h):
    n_frames = 1 + (len(x) - L) // h
    s = x.strides[0]
    return np.lib.stride_tricks.as_strided(x, shape=(n_frames, L), strides=(h*s, s))

# Apply the window to every frame at once - broadcasting, one line.
w = np.hanning(400)                              # (400,)
frames = frame_signal(x, 400, 160) * w             # (398, 400) * (400,) -> (398, 400)

# FORM 3 - the library one-liner. Same array, decades of edge cases handled.
import librosa
frames = librosa.util.frame(x, frame_length=400, hop_length=160).T * w

Run form 1 and form 3 on the same signal and subtract: the maximum absolute difference is zero to floating-point precision. The library is not doing anything mysterious. It is doing the loop, faster, with the boundary conventions written down.

Chapter 2: Time-Domain Features — Energy and Zero Crossings

We have a frame: 400 windowed samples, one still photograph of the sound. Now we describe it with a few numbers. The first two cost no Fourier transform at all — a sum and a comparison — which is exactly why they were the first features anybody used, on hardware that could not have afforded an FFT per frame.

Do not dismiss them as historical. Short-time energy and zero-crossing rate still sit in production voice-activity detectors today, because they are two orders of magnitude cheaper than anything spectral and they answer the two coarsest questions about a frame: is there anything here? and is it buzzy or smooth?

Feature 1: short-time energy

The short-time energy of a frame is just the sum of squared samples:

Squaring does two jobs. It makes every contribution positive (a pressure dip is as much sound as a pressure peak, and we must not let them cancel), and it weights loud excursions more than quiet ones, which matches the physics: acoustic power really does go with amplitude squared.

Energy grows with frame length, so we usually report the root-mean-square instead, which does not:

Worked by hand: energy of our 8-sample frame

Frame x = [ 0.1, 0.6, −0.4, −0.8, 0.2, 0.9, −0.1, −0.5 ]. Square each sample:

Add them left to right so you can check every partial sum: 0.01, then 0.37, then 0.53, then 1.17, then 1.21, then 2.02, then 2.03, then 2.28. So E = 2.28.

Mean square = 2.28 / 8 = 0.285. RMS = √0.285. Since 0.53² = 0.2809 and 0.54² = 0.2916, the answer is between them; refine: 0.534² = 0.28516, slightly high, so RMS = 0.5339.

Why we take the log

Audio dynamic range is enormous: a whisper and a jackhammer differ by a factor of about a million in amplitude, a trillion in energy. Feed raw energy into a distance-based classifier and the loud classes occupy a region a million times wider than the quiet ones. Every classical pipeline therefore uses log energy, usually in decibels:

For our frame: E/L = 0.285. Now, log₁₀(0.285) = log₁₀(2.85) − 1. Since log₁₀(2.85) ≈ 0.4548, we get −0.5452. Multiply by 10: E_dB = −5.45 dB relative to a full-scale sine.

The log also converts a multiplicative nuisance into an additive one. Halve the microphone gain and every sample halves, so E drops by a factor of 4 — but E_dB drops by exactly 6.02 dB regardless of what the sound was. A constant offset is something a classifier can be made blind to (subtract the per-clip mean); a multiplicative factor is not. This trick, turn multiplication into addition so it becomes removable, is going to reappear as the entire justification for the log in MFCCs.

Feature 2: zero-crossing rate

The zero-crossing rate (ZCR) counts how often the waveform changes sign inside the frame. A signal that oscillates rapidly crosses zero often; a slow, smooth signal crosses rarely. It is a frequency measurement made without any frequency mathematics.

You will also meet the equivalent textbook form using the absolute difference of signs, which counts each crossing as 2 and divides by 2L:

Worked by hand: ZCR of the same frame

Write out the signs: +   +   −   −   +   +   −   −

Now walk the seven adjacent pairs and mark the changes:

Three crossings out of seven pairs, so ZCR = 3/7 = 0.4286 crossings per sample pair. In the alternative convention: each of the 3 crossings contributes |+1 − (−1)| = 2, so the sum is 6, and 6 / (2×8) = 0.375. Same information, different bookkeeping — always check which convention a library uses before comparing numbers across papers.

Turning ZCR into hertz

ZCR is only interesting once you convert it to a physical rate. Seven sample intervals at 16 kHz span 7 / 16000 = 0.4375 ms. Three crossings in that time is

A pure sine crosses zero twice per cycle, so the implied dominant frequency is 6857 / 2 ≈ 3429 Hz. That is a bright, hissy frame — consistent with a fricative or noise, not with a vowel. Our eight-sample toy is barely enough to estimate anything, but the arithmetic is exactly what runs on a 400-sample frame.

What the two features mean together

Individually each is weak. Together they carve up the acoustic world surprisingly well, which is why the pair has survived since the 1970s:

The classic endpoint detector of Rabiner and Sambur (1975) is nothing more than thresholds on these two numbers: estimate energy and ZCR statistics from the first 100 ms (assumed to be background), then mark speech as starting where energy exceeds a threshold, and back the boundary up to include preceding high-ZCR low-energy frames — because a word starting with “s” has no energy onset to find. That backing-up rule is the entire reason ZCR was invented, and it is a beautiful piece of engineering: a feature designed to catch the specific failure of another feature.

Two more cheap time-domain numbers

Crest factor = peak / RMS tells you how spiky a frame is. For our frame: peak amplitude 0.9, RMS 0.5339, so crest = 0.9 / 0.5339 = 1.686. A pure sine has crest √2 = 1.414; Gaussian noise sits around 3–4; a click can exceed 10. Impact sounds and percussive events are exactly the high-crest ones.

Autocorrelation peak gives real pitch, unlike ZCR. Correlate the frame with a delayed copy of itself at lag τ; the lag with the largest peak (excluding τ = 0) is the period, and the pitch is sr / τ. Its height relative to the τ = 0 value measures periodicity, which is the honest version of the question ZCR was pretending to answer.

Push the DC-offset slider up with the mix at pure tone: the waveform lifts off the axis and the crossing ticks vanish one by one until the ZCR reads zero for a signal that has not changed audibly at all. That is the failure mode from the callout above, live.

From one frame to a contour

A single frame’s energy is a number; the sequence of frame energies is a contour, and the contour is where events live. Suppose five consecutive frames have RMS values

The first difference — frame-to-frame change, the discrete derivative — is [ +0.01, +0.38, −0.03, −0.26 ]. The huge positive jump between frames 1 and 2 is an onset: something started. The slow decay afterwards is a release. A classifier fed only the five raw numbers in some order could not tell an onset from an offset; fed the differences, it can.

This is the first appearance of delta features, and they will come back formally in Chapter 4. The idea is always the same: append the rate of change of each feature to the feature vector, because how a sound evolves is often more diagnostic than what it is at any instant. A door slam and a bass note can have identical average spectra and completely different envelopes.

Code: three forms

python
# FORM 1 - the arithmetic exactly as we did it by hand.
def energy_zcr_manual(frame):
    e = 0.0
    for s in frame:
        e += s * s                                  # sum of squares
    rms = (e / len(frame)) ** 0.5
    crossings = 0
    for n in range(1, len(frame)):
        if (frame[n] >= 0) != (frame[n-1] >= 0):
            crossings += 1
    return rms, crossings / (len(frame) - 1)

# FORM 2 - vectorised over ALL frames at once. Shape (n_frames, L) -> (n_frames,)
import numpy as np

def energy_zcr(frames):
    frames = frames - frames.mean(axis=1, keepdims=True)   # kill DC first!
    rms = np.sqrt((frames ** 2).mean(axis=1))
    sgn = np.signbit(frames)                                # True where negative
    zcr = (sgn[:, 1:] != sgn[:, :-1]).mean(axis=1)          # fraction of pairs
    return 20 * np.log10(rms + 1e-10), zcr             # +eps: log(0) is -inf

# FORM 3 - the library. Note librosa returns shape (1, n_frames).
import librosa
rms = librosa.feature.rms(y=x, frame_length=400, hop_length=160)
zcr = librosa.feature.zero_crossing_rate(x, frame_length=400, hop_length=160)

Chapter 3: Spectral Features — Centroid, Rolloff, Flux

Energy told us how much. ZCR told us how fast, badly. Now we take the Fourier transform of the windowed frame and describe the resulting spectrum with a handful of numbers — and this is where the features start to sound like the words a person would use: bright, thin, noisy, changing.

Concretely: our 400-sample windowed frame goes into an FFT and comes out as 201 complex numbers. Throw away the phase (Chapter 0 taught us why) and keep the magnitudes m_k, one per frequency bin k = 0…200, with bin k sitting at frequency f_k = k · sr / L = k × 40 Hz. Every feature in this chapter is a one-line summary of that list of 201 numbers.

The toy spectrum we will compute by hand

201 numbers is too many to do on paper, so shrink it to five bins spaced 2 kHz apart. This is a real spectrum shape — energy concentrated low-mid, tailing off — just coarsely sampled:

The total magnitude is 1.0 + 3.0 + 4.0 + 1.5 + 0.5 = 10.0. Keep that number; almost every feature below divides by it.

Feature 3: spectral centroid

The spectral centroid is the centre of mass of the spectrum: treat each magnitude as a weight sitting at its frequency and find the balance point.

Work the numerator term by term, so nothing is hidden:

Numerator = 0 + 6000 + 16000 + 9000 + 4000 = 35000. Divide by the total magnitude 10.0:

Sanity check the number against intuition: most of the weight sits at 2 and 4 kHz, and the small tail at 6–8 kHz drags the balance point up a little. 3500 Hz — just above the weighted middle — is exactly right.

Perceptually, centroid tracks brightness. A muted trumpet has a lower centroid than an open one; a dull thud lower than a hi-hat; the vowel “oo” lower than “ee”. Of all the hand-designed features, this is the one that correlates best with a word non-engineers actually use.

Feature 4: spectral spread (bandwidth)

Centroid is a mean, so the natural companion is a standard deviation: how far the energy is spread around the centre of mass.

With centroid = 3500, the deviations are −3500, −1500, +500, +2500, +4500. Square them:

Sum = 39,500,000. Divide by 10.0 → 3,950,000. Square root: 1987² = 3,948,169 and 1988² = 3,952,144, so

Interpretation: this frame’s energy is centred at 3.5 kHz and typically lands within about 2 kHz of it. A pure tone would have spread near zero; white noise across 0–8 kHz would have spread around 2.3 kHz. Ours is nearly as spread out as noise, which tells you the toy spectrum is fairly broadband.

Feature 5: spectral rolloff

Rolloff answers: below which frequency does a given fraction (usually 85%) of the total magnitude live? It is a percentile of the spectrum, and it captures skew — how much high-frequency tail there is — more robustly than the maximum frequency, which is dominated by the noise floor.

Target = 0.85 × 10.0 = 8.5. Walk the cumulative sum:

Notice how sharply this can move: at bin 2 the cumulative was 8.0, only 0.5 short. Nudge m₂ from 4.0 to 4.6 and the rolloff jumps from 6000 Hz down to 4000 Hz — a 33% change in the feature from a 6% change in one bin. Rolloff is a quantized feature; treat single-frame values with suspicion and rely on its average over many frames.

Feature 6: spectral flux

Every feature so far describes one frame in isolation. Spectral flux describes the change between consecutive frames — the speed at which the sound is transforming:

Let the next frame’s spectrum be m′ = [1.0, 2.0, 3.0, 3.0, 2.0] — energy has shifted upward in frequency. Difference bin by bin:

Sum of squares = 0 + 1.00 + 1.00 + 2.25 + 2.25 = 6.50, so flux = √6.50 = 2.550. The half-wave-rectified version keeps only the increases: 0 + 0 + 0 + 1.5 + 1.5 = 3.0.

Feature 7: spectral flatness

Flatness (also called Wiener entropy) is the ratio of the geometric mean to the arithmetic mean of the magnitudes. It answers “tonal or noisy?” on a scale from 0 to 1.

Product of our five magnitudes: 1.0 × 3.0 = 3.0; × 4.0 = 12.0; × 1.5 = 18.0; × 0.5 = 9.0. Fifth root: 9^1/5 = e^{(ln 9)/5} = e^2.1972/5 = e^0.4394 = 1.552. Arithmetic mean = 10.0 / 5 = 2.0. So

Check the extremes to calibrate. Pure tone, spectrum [0, 0, 10, 0, 0]: the product contains a zero, so the geometric mean is 0 and flatness = 0 — maximally tonal. White noise, spectrum [2, 2, 2, 2, 2]: geometric mean = 2, arithmetic mean = 2, flatness = 1 — maximally flat. Our 0.776 says “quite noisy, with some structure.”

What each feature buys you

Code: three forms

python
# FORM 1 - one frame, exactly the arithmetic we just did on paper.
import numpy as np

def spectral_feats_manual(m, freqs, m_prev=None, roll=0.85):
    total = m.sum() + 1e-10
    centroid = (freqs * m).sum() / total
    spread   = np.sqrt((m * (freqs - centroid) ** 2).sum() / total)
    cum      = np.cumsum(m)
    rolloff  = freqs[np.searchsorted(cum, roll * total)]      # first bin past 85%
    gm       = np.exp(np.log(m + 1e-10).mean())               # geometric mean, safely
    flatness = gm / (m.mean() + 1e-10)
    flux     = 0.0 if m_prev is None else np.sqrt(((m - m_prev) ** 2).sum())
    return centroid, spread, rolloff, flatness, flux

# Verify against the hand computation:
m = np.array([1.0, 3.0, 4.0, 1.5, 0.5])
f = np.array([0, 2000, 4000, 6000, 8000], dtype=float)
print(spectral_feats_manual(m, f))
# -> (3500.0, 1987.46, 6000.0, 0.7759, 0.0)   exactly our paper numbers

# FORM 2 - the whole clip at once. frames: (T, L) already windowed.
S = np.abs(np.fft.rfft(frames, axis=1))        # (T, L//2+1) magnitudes
freqs = np.fft.rfftfreq(L, 1/sr)                  # (L//2+1,) in Hz
tot = S.sum(1, keepdims=True) + 1e-10
centroid = (S * freqs).sum(1) / tot[:, 0]        # (T,)
flux = np.sqrt((np.diff(S, axis=0) ** 2).sum(1))    # (T-1,)

# FORM 3 - library one-liners. Each returns (1, T).
import librosa
centroid = librosa.feature.spectral_centroid(y=x, sr=sr, n_fft=400, hop_length=160)
rolloff  = librosa.feature.spectral_rolloff(y=x, sr=sr, roll_percent=0.85)
flatness = librosa.feature.spectral_flatness(y=x)
onset    = librosa.onset.onset_strength(y=x, sr=sr)   # rectified flux, mel-scaled

Chapter 4: MFCCs — The Feature Vector That Ran the World

Chapter 3 gave us seven summary statistics of the spectrum. They cannot distinguish “ah” from “eh”, because both are voiced, similarly bright, similarly noisy — they differ in the shape of the spectral envelope, in where the resonant bumps sit. To capture shape you need a handful of numbers that describe a curve, not a curve’s average.

The answer, from Davis and Mermelstein in 1980, is the Mel-Frequency Cepstral Coefficient vector. Thirteen numbers per frame. It became so dominant that for thirty years, “the features” in a speech paper meant MFCCs without further comment. If you learn one thing from this lesson, learn why each of its stages exists — because every stage is a decision, and every decision has a cost that Chapter 11 will collect.

Why the spectrum has two things mixed in it

The source-filter model: a voiced sound is made by the vocal folds producing a buzzy pulse train (the source, which sets the pitch), and that buzz then passes through the throat, mouth and nose, which act as a resonant tube (the filter, which sets the vowel). In the frequency domain, filtering is multiplication:

E(f) is a comb of sharp harmonics spaced at the pitch (say every 120 Hz for a male voice). H(f) is a smooth curve with a few broad bumps — the formants — and it is H(f) that determines which vowel. The observed spectrum is the comb with the smooth curve as its envelope.

Almost everything you want for classification is in H. The pitch in E is a nuisance: the same word spoken by a child and an adult has utterly different combs and nearly the same envelope. So the task is to separate a product into its factors, given only the product. That sounds impossible until you take a logarithm.

The transform of a log spectrum is called the cepstrum (“spectrum” with the first syllable reversed — Bogert, Healy and Tukey, 1963, who also gave us quefrency for its axis and liftering for filtering along it). Low quefrency = slowly varying = the envelope we want. High quefrency = the pitch comb we want to drop. Keeping the first ~13 coefficients is exactly “keep the envelope, discard the pitch.”

The five stages

Stage 2 recap: the mel scale

We built this in Audio Representations (and it is derived in detail in EE269-13: Cepstrum & MFCC), so here is just the arithmetic you need. The mel scale warps hertz to match how finely we can hear differences:

Check a few values with a calculator: m(1000) = 2595 × log₁₀(1 + 1.4286) = 2595 × log₁₀(2.4286) = 2595 × 0.3854 = 1000.0 mel — the constant 2595 is chosen so that 1000 Hz maps to 1000 mel. Then m(2000) = 1521, m(4000) = 2146, m(8000) = 2840.

Now compare two stretches of the frequency axis that are worlds apart in width. The first 100 Hz — from 0 to 100 Hz — spans m(100) = 150 mel. The stretch from 4000 to 8000 Hz is 4000 Hz wide, forty times bigger, and spans only 2840 − 2146 = 694 mel, under five times more perceptual distance. Per octave the mel cost is nearly constant (1→2 kHz gains 521 mel, 2→4 kHz gains 625, 4→8 kHz gains 694) even though the hertz cost doubles every time.

Spacing the filterbank evenly in mel therefore gives narrow, closely packed triangles down low and wide, sparse triangles up high — resolution spent where the ear has it, and 201 FFT bins collapsed into 40 band energies along the way.

Worked by hand: log-mel to MFCC on four bands

Forty bands is too many for paper; use four. Suppose the mel filterbank outputs the band energies

Stage 3, take natural logs: ln 2 = 0.693, ln 8 = 2.079, ln 4 = 1.386, ln 1 = 0.000. So the log-mel vector is

Stage 4, the DCT-II. Its definition, for N values:

With N = 4, the cosine table is small enough to write out completely. The angles are πk(n+0.5)/4 in radians, shown here in degrees:

Now multiply the row by s and add. c₀ is the easy one — all cosines are 1, so it is just the sum:

c₁, term by term: 0.693 × 0.9239 = 0.640; 2.079 × 0.3827 = 0.796; 1.386 × (−0.3827) = −0.530; 0.000 × (−0.9239) = 0.000.

c₂: every cosine is ±0.7071, so factor it out: 0.7071 × (0.693 − 2.079 − 1.386 + 0.000) = 0.7071 × (−2.772) =

c₃: 0.693 × 0.3827 = 0.265; 2.079 × (−0.9239) = −1.921; 1.386 × 0.9239 = 1.281; 0.000 × (−0.3827) = 0.000.

So the four-band MFCC vector is [4.158, 0.906, −1.960, −0.375]. (If you check this in SciPy you will get [8.318, 1.811, −3.921, −0.750] — exactly twice ours, because scipy.fftpack.dct(type=2, norm=None) includes a factor of 2. Conventions differ by constant factors everywhere in this corner of DSP; what matters is that you use one consistently.)

Reading the coefficients

Each c_k is the amount of a particular cosine ripple present in the log-spectrum shape:

Confirm against the data: our band energies were [2, 8, 4, 1] — a peak at band 1, decreasing to the right. That is exactly “tilted low” (c₁ > 0) with “a middle bump” (c₂ < 0). The coefficients are not a black box; you can read a rough picture of the spectrum straight off them.

Why the DCT, really: a parameter-count argument

The perceptual story explains the mel scale and the log. The DCT’s justification is different and purely statistical, and it is the most under-explained step in the whole pipeline.

Neighbouring mel bands are strongly correlated: a broad formant lights up four or five adjacent triangles together. The classifier we are heading for (Chapter 6) is a Gaussian mixture, and a Gaussian over correlated dimensions needs a full covariance matrix. For a 39-dimensional feature that is 39 × 40 / 2 = 780 parameters per component, and you need to estimate every one of them from your data. With a diagonal covariance you need 39. That is a twentyfold reduction — the difference between a model you can fit with an hour of speech and one you cannot.

The DCT is (approximately) the eigenbasis of the covariance of log-mel vectors, so after applying it the dimensions are nearly uncorrelated and the diagonal assumption becomes defensible. In other words: MFCCs exist so that a diagonal-covariance GMM is legitimate. When the classifier changed — when neural nets arrived and correlated inputs stopped being a problem — the DCT immediately became unnecessary, which is why every modern audio model eats log-mel, not MFCC. The feature was co-designed with its classifier, and it retired with it.

Deltas: 13 becomes 39

A static MFCC vector says nothing about motion, and Chapter 2 already showed motion is diagnostic. The standard fix is to append the first and second time derivatives, estimated by a regression over a window of ±M frames:

With M = 2 the denominator is 2(1 + 4) = 10, so Δc_t = [ (c_t+1 − c_t−1) + 2(c_t+2 − c_t−2) ] / 10. Apply it to the c₁ sequence [0.4, 0.6, 0.9, 1.3, 1.6] at the centre frame: (1.3 − 0.6) + 2(1.6 − 0.4) = 0.7 + 2.4 = 3.1, divided by 10 gives 0.31 per frame. Applying the same formula to the delta stream gives the delta-delta. Result: 13 static + 13 delta + 13 delta-delta = 39 dimensions per frame, the number stamped on twenty years of speech systems.

Cepstral mean normalization, and why it works

Record the same voice on a cheap phone and on a studio microphone. The channel is a fixed filter G(f), so the observed spectrum is S(f)·G(f) — another multiplication. In the log domain it is an additive constant, identical in every frame. In the cepstral domain, a DCT of a constant is still a constant. Therefore: subtract the per-clip mean of each coefficient and the channel is gone.

This one line of code — cepstral mean normalization, CMN — was worth several percent absolute in every ASR system of the era, and it is worth understanding why: it works only because the log turned a convolution into an addition. Every stage of MFCC is paying off some other stage. That is what a well-designed hand pipeline looks like.

Code: three forms

python
# FORM 1 - the DCT exactly as computed on paper. O(N^2), obviously correct.
import numpy as np

def dct2_manual(s):
    N = len(s)
    c = np.zeros(N)
    for k in range(N):
        for n in range(N):
            c[k] += s[n] * np.cos(np.pi * k * (n + 0.5) / N)
    return c

print(dct2_manual(np.log([2., 8., 4., 1.])))
# -> [ 4.159  0.906 -1.961 -0.375 ]   the paper answer

# FORM 2 - full MFCC front end from scratch, shapes annotated.
def mfcc_from_scratch(frames, sr, n_mels=40, n_mfcc=13):
    S    = np.abs(np.fft.rfft(frames, axis=1)) ** 2     # (T, 201) power
    fb   = mel_filterbank(sr, frames.shape[1], n_mels)  # (40, 201) triangles
    mel  = S @ fb.T                                    # (T, 40)  band energies
    logm = np.log(mel + 1e-10)                          # (T, 40)  product -> sum
    c    = dct2_manual_batch(logm)[:, :n_mfcc]         # (T, 13)  keep envelope
    c   -= c.mean(axis=0, keepdims=True)                # CMN: channel removed
    d1   = delta(c); d2 = delta(d1)
    return np.hstack([c, d1, d2])                      # (T, 39)

# FORM 3 - the one-liner. Identical maths, three decades of tuning.
import librosa
m = librosa.feature.mfcc(y=x, sr=16000, n_mfcc=13, n_fft=400, hop_length=160, n_mels=40)
feats = np.vstack([m, librosa.feature.delta(m), librosa.feature.delta(m, order=2)])  # (39, T)

Chapter 5: k-Nearest Neighbours — The Laziest Classifier That Works

Every clip in the training set is now a point: a fixed-length vector of pooled features with a label attached. A new clip arrives and becomes another point. The simplest possible thing you could do is look at which labelled points it landed near.

That is k-nearest-neighbour classification in full. Store everything. To classify, compute the distance to every stored point, take the k closest, and let them vote. There is no training step, no parameters to fit, no optimizer. It is the right first classifier for any new problem because it tells you, immediately, whether your features are any good — if k-NN cannot separate the classes, the geometry is not there and no fancier model will conjure it.

The setup, with real audio features

Use two features so we can draw and compute everything: zero-crossing rate and spectral centroid in kHz, each averaged over the clip. Four labelled training clips:

Worked by hand: raw distances

Euclidean distance is √((ΔZCR)² + (Δcentroid)²). Compute each, showing both squared terms so you can see which one matters:

With k = 1 the answer is C, so we predict music. With k = 3 the three nearest are C (0.204), D (0.800) and B (0.802): two music, one speech, so again music. Fine — but look at the columns.

The (ΔZCR)² column runs from 0.0001 to 0.0025. The (Δcent)² column runs from 0.04 to 1.69. The centroid contributes between 16 and 17,000 times more to every distance. The ZCR feature is, arithmetically, not being used at all. We computed a two-dimensional distance and got a one-dimensional answer.

Worked by hand: standardization (z-scoring)

Give every feature the same voice by converting it to “standard deviations away from the training mean”:

ZCR statistics. Mean: (0.10 + 0.14 + 0.05 + 0.07)/4 = 0.36/4 = 0.090. Deviations: +0.01, +0.05, −0.04, −0.02. Squares: 0.0001, 0.0025, 0.0016, 0.0004; sum 0.0046; divide by 4 to get variance 0.00115; square root σ = 0.0339.

Centroid statistics. Mean: (1.5 + 2.0 + 3.0 + 3.6)/4 = 10.1/4 = 2.525. Deviations: −1.025, −0.525, +0.475, +1.075. Squares: 1.0506, 0.2756, 0.2256, 1.1556; sum 2.7075; variance 0.6769; σ = 0.8227.

Now transform every point. For A: z₁ = (0.10 − 0.090)/0.0339 = +0.295; z₂ = (1.5 − 2.525)/0.8227 = −1.246. The rest, the same way:

Recompute the four distances in this new space. To A: Δ = (0.295 − 0,  −1.246 − 0.334) = (0.295, −1.580); squares 0.0870 and 2.4964; sum 2.5834; √ = 1.607. The full set:

The nearest neighbour changed, from C to D. The predicted label happens to survive here (both are music), but the ranking is different, and on a harder query it would flip the answer. Same data, same algorithm, different result — purely because we chose units. Euclidean distance is only meaningful when the axes are commensurate, and making them commensurate is your job, not the classifier’s.

Choosing k

k = 1 gives a decision boundary that wraps tightly around every training point, including mislabelled and freak ones: zero training error, high variance. Large k averages over a wide neighbourhood: smooth, stable, but it will smear away genuinely small classes. In the limit k = N, every query gets the majority label of the whole dataset.

Pick k by cross-validation, and prefer odd k with two classes so votes cannot tie. A refinement worth knowing: distance weighting, where neighbour i votes with weight 1/d_i, so a neighbour at distance 0.2 outweighs one at 0.8 by four to one. That makes k less critical and is usually a free improvement.

The curse of dimensionality, measured

k-NN degrades as dimensionality grows, and the reason is worth seeing numerically rather than as a slogan. Scatter 2000 random points uniformly in a d-dimensional unit cube, pick a query, and look at how far away the nearest and farthest points are:

In two dimensions the nearest neighbour is meaningfully nearer than everything else: the notion of “close” has teeth. In 78 dimensions — the size of our pooled feature vector — the farthest point is only about 1.6× the distance of the nearest. Everything is roughly the same distance from everything, and the word “neighbour” loses most of its meaning. That is why classical systems worked so hard to keep feature vectors small (13 MFCCs, not 40), and why k-NN is usually the worst of the three classifiers in this lesson on real audio benchmarks.

Other distances you will meet

Note that standardizing and then using Euclidean distance is exactly Mahalanobis with a diagonal Σ. The two ideas are the same idea; full Mahalanobis additionally undoes correlations between features, which is what the DCT in Chapter 4 was already doing for us.

Worked by hand: distance-weighted voting

Return to the raw-distance ranking: C at 0.204 (music), D at 0.800 (music), B at 0.802 (speech). A plain k = 3 vote gives music 2, speech 1 — a 67% “confidence” that ignores how much nearer C is than the other two. Weight each vote by 1/d instead:

Music scores 4.902 + 1.250 = 6.152; speech scores 1.247. Normalizing, music takes 6.152 / 7.399 = 83% of the vote rather than 67%. The weighting also makes the result far less sensitive to k: adding a fourth neighbour at distance 1.5 would contribute only 0.67, barely moving the total. This is why weights='distance' is usually a free improvement, and why it is the first knob to turn when k-NN results look unstable.

The cost, and why k-NN stayed at the clip level

Brute-force k-NN costs O(N · D) per query. With N = 2000 training clips and D = 86 features, that is 172,000 multiply-adds — microseconds. Perfectly fine. Now consider using it per frame: a hundred hours of audio is 36 million frames, so each query costs 36 million × 39 operations and you have 100 queries per second of test audio. Hopeless.

Two classical mitigations are worth knowing. KD-trees and ball trees partition the space so a query examines only a fraction of the points; they help enormously in low dimensions and, thanks to the concentration effect above, help almost not at all beyond about 20 dimensions. Condensed nearest neighbour instead prunes the training set itself, keeping only the points near class boundaries — typically 10–20% of the data, with nearly identical accuracy. That is a striking foreshadowing of the next chapter but one: an SVM does the same thing, but chooses which points to keep by optimization rather than by heuristic.

Code: three forms

python
# FORM 1 - the loop, exactly the table we computed by hand.
import numpy as np
from collections import Counter

def knn_manual(X_train, y_train, q, k=3):
    dists = []
    for i, p in enumerate(X_train):
        d = 0.0
        for j in range(len(q)):
            d += (p[j] - q[j]) ** 2          # sum of squared differences
        dists.append((d ** 0.5, y_train[i]))
    dists.sort()                              # nearest first
    return Counter(l for _, l in dists[:k]).most_common(1)[0][0]

# FORM 2 - vectorised, all queries at once. No Python loop at all.
def knn_batch(X_train, y_train, Q, k=3):
    # ||a-b||^2 = ||a||^2 + ||b||^2 - 2 a.b   -- the expansion that makes this fast
    d2 = (Q**2).sum(1)[:, None] + (X_train**2).sum(1)[None, :] - 2 * Q @ X_train.T
    idx = np.argsort(d2, axis=1)[:, :k]        # (n_query, k)
    votes = y_train[idx]                        # (n_query, k)
    return np.array([np.bincount(v).argmax() for v in votes])

# FORM 3 - the library, with the scaler fitted on TRAIN ONLY.
from sklearn.preprocessing import StandardScaler
from sklearn.neighbors import KNeighborsClassifier
from sklearn.pipeline import make_pipeline

clf = make_pipeline(StandardScaler(), KNeighborsClassifier(n_neighbors=5, weights='distance'))
clf.fit(X_train, y_train)        # the pipeline fits the scaler on train only - no leakage
print(clf.score(X_test, y_test))

Chapter 6: Gaussian Mixtures & One EM Iteration By Hand

k-NN memorised the training set. What we want instead is a compact description of each class: a probability density over feature space that says “dog barks live around here, with this much spread.” Then classifying is scoring — ask each class density how likely it thinks the new clip is, and take the winner.

This was the dominant approach of the classical era, and not by a small margin. Every speech recognizer from 1990 to 2012 was a bank of Gaussian mixtures. Speaker verification, language identification, the DCASE 2013 acoustic-scene baseline: Gaussian mixtures. It is worth doing the arithmetic yourself once, which is what this chapter is.

Start with one Gaussian

In one dimension, a Gaussian says “values cluster around μ with a typical wander of σ”:

Two parameters, both trivially estimated: μ is the sample mean, σ² the sample variance. For a D-dimensional feature vector with a diagonal covariance you simply do this per dimension and multiply the results — which is precisely the independence assumption the DCT in Chapter 4 bought us.

But one Gaussian is a single blob, and a sound class is rarely one blob. “Dog” contains yaps and deep woofs. The phoneme /t/ has a silent closure followed by a burst. A single Gaussian fitted to a two-lump cloud puts its mean in the empty valley between the lumps and assigns high probability to a region where no data ever occurs. That is not a small inaccuracy; it is a density that is confidently wrong exactly where it matters.

A mixture of Gaussians

So use several, and weight them:

The generative story is what makes the maths tractable: to produce a data point, first roll a die weighted by π to pick a component k, then draw from that component’s Gaussian. If we could see the die roll for each training point, fitting would be trivial — sort the points by component and compute means and variances within each group.

We cannot see it. The component identity is a latent variable. And that is the entire difficulty, dissolved by one idea: if we cannot observe the die roll, we will estimate the probability of each possible roll, and use those probabilities as soft group memberships.

The EM loop

The E-step is just Bayes’ rule. The probability that component k produced point x is the probability of picking k times the probability k would emit x, normalized over components:

The M-step is the ordinary mean/variance formulas with each point counted r_nk times instead of once:

Worked by hand: the complete first iteration

Four data points — think of them as a one-dimensional feature (say c₁) from four frames:

Initialize deliberately badly, so there is something to fix: μ₁ = 4, μ₂ = 7, σ₁² = σ₂² = 4, π₁ = π₂ = 0.5.

E-step

Because the weights and variances are equal, the normalizing constants and the πs cancel in the ratio, and the responsibility depends only on the exponents. With 2σ² = 8:

Point x = 2. Squared distances: (2−4)² = 4 and (2−7)² = 25. Exponents: −4/8 = −0.500 and −25/8 = −3.125. Exponentials: e^−0.500 = 0.60653, e^−3.125 = 0.04394. Their sum is 0.65047. So

Point x = 3. (3−4)² = 1, (3−7)² = 16. Exponents −0.125 and −2.000. Exponentials 0.88250 and 0.13534, sum 1.01784.

Point x = 8. (8−4)² = 16, (8−7)² = 1 — the mirror image of x = 3, so r₁(8) = 0.1330, r₂(8) = 0.8670.

Point x = 9. Mirror of x = 2: r₁(9) = 0.0675, r₂(9) = 0.9325.

Every point belongs partly to both components. Nothing has been assigned; everything has been weighted. That softness is the whole reason EM converges smoothly where hard assignment (k-means) jumps.

M-step

Effective counts. N₁ = 0.9325 + 0.8670 + 0.1330 + 0.0675 = 2.000. By symmetry N₂ = 2.000. Total 4 — as it must be, since each point’s responsibilities sum to 1.

New means. Weight each point by its responsibility:

Term by term: 0.9325 × 2 = 1.8649; 0.8670 × 3 = 2.6011; 0.1330 × 8 = 1.0637; 0.0675 × 9 = 0.6079. Sum = 6.1376. Divide by 2.000:

The means moved from (4, 7) to (3.07, 7.93) — outward, towards the real clusters at (2.5, 8.5). Check the symmetry: 3.069 + 7.931 = 11.0, and the data is symmetric about 5.5, so the means must stay symmetric about 5.5 too. They do.

New variances. Using the updated μ₁ = 3.069:

Sum = 1.0651 + 0.0041 + 3.2333 + 2.3764 = 6.6789. Divide by N₁ = 2.000:

New weights. π₁ = 2.000 / 4 = 0.5, π₂ = 0.5 — unchanged, because the data is perfectly balanced.

Did it actually improve?

EM promises the log-likelihood never decreases. Verify it. The mixture density at a point is π₁N(x|μ₁,σ₁²) + π₂N(x|μ₂,σ₂²).

Before. The normalizing constant is 1/√(2π×4) = 1/5.0133 = 0.19947. At x = 2: 0.5 × 0.19947 × (0.60653 + 0.04394) = 0.099735 × 0.65047 = 0.06487, whose natural log is −2.7355. At x = 3: 0.099735 × 1.01784 = 0.10151, log −2.2876. Points 8 and 9 mirror them, so

After. Now σ² = 3.3394, so the constant is 1/√(2π×3.3394) = 1/4.5806 = 0.21831. At x = 2: distances 1.1423 and 35.179; exponents −1.1423/6.6788 = −0.17103 and −5.2672; exponentials 0.84280 and 0.00517; sum 0.84797; density 0.5 × 0.21831 × 0.84797 = 0.09256, log −2.3796. At x = 3: exponentials 0.99929 and 0.02626, sum 1.02555, density 0.11194, log −2.1899. So

It rose by 0.907 nats in one iteration. Run it to convergence and the means walk out to roughly 2.5 and 8.5 with small variances, exactly where a human would put them.

From one dimension to a real audio system

Everything above holds in D dimensions with the scalar variance replaced by a diagonal covariance vector. The bookkeeping for a speech system:

Classification is then Bayes: pick the class c maximizing log p(X | class c) + log P(class c). Because the log-likelihood is a sum over frames, a longer clip produces a more negative score — so if clip lengths differ, divide by the frame count before comparing, or your classifier will develop opinions about duration.

Code: two forms

python
# FORM 1 - EM from scratch, 1-D, exactly the arithmetic above.
import numpy as np

def em_1d(x, mu, var, pi, n_iter=20, var_floor=1e-3):
    for it in range(n_iter):
        # E-step: (N, K) responsibilities
        comp = pi * np.exp(-(x[:, None] - mu)**2 / (2*var)) / np.sqrt(2*np.pi*var)
        ll   = np.log(comp.sum(1) + 1e-300).sum()      # total log-likelihood
        r    = comp / comp.sum(1, keepdims=True)
        # M-step
        Nk  = r.sum(0)                                    # (K,) effective counts
        mu  = (r * x[:, None]).sum(0) / Nk
        var = (r * (x[:, None] - mu)**2).sum(0) / Nk
        var = np.maximum(var, var_floor)               # the floor that prevents collapse
        pi  = Nk / len(x)
        print(it, ll.round(4), mu.round(4), var.round(4))
    return mu, var, pi

em_1d(np.array([2., 3., 8., 9.]), np.array([4., 7.]), np.array([4., 4.]), np.array([.5, .5]))
# iteration 0 prints ll = -10.0457, then mu = [3.0688 7.9312], var = [3.3394 3.3394]

# FORM 2 - one GMM per class on real pooled features, then Bayes classification.
from sklearn.mixture import GaussianMixture

models = {}
for c in classes:
    models[c] = GaussianMixture(n_components=32, covariance_type='diag',
                                 reg_covar=1e-4, n_init=3).fit(X_train[y_train == c])
# reg_covar IS the variance floor; n_init=3 restarts to escape bad local optima

scores = np.stack([models[c].score_samples(X_test) for c in classes])   # (C, n_test)
pred   = np.array(classes)[scores.argmax(0)]

Chapter 7: Support Vector Machines & the Kernel Trick

A GMM describes what each class looks like: it models the density, then compares. That is a generative approach, and it spends its capacity describing regions of feature space you may never need to distinguish. If all you want is the answer, model the boundary.

That is the discriminative view, and its champion in this era is the support vector machine. On the pooled clip-level vectors we built in Chapter 5, an SVM with an RBF kernel was — and honestly still is — a strong baseline you should beat before claiming anything.

Many boundaries separate; one is widest

Take two classes that are linearly separable in feature space. There are infinitely many straight lines that split them, and the training error of every one of them is zero. Training error cannot choose. So choose by a different principle: pick the line with the widest empty corridor around it.

Why is width the right criterion? Because your features are noisy. Centroid estimates wobble by tens of hertz between recordings of the same event; pooled statistics shift with clip length. A boundary that skims a training point will misclassify that same sound on a slightly different day. A boundary with a 2-unit corridor tolerates 1 unit of feature noise on either side. Margin is a robustness budget, purchased in advance.

Worked by hand: the two-point SVM

Reduce to the smallest problem with an answer. One positive point at (3, 3) and one negative at (1, 1). The classifier is

The SVM convention scales w and b so the closest points sit exactly at f = ±1 (this is a choice of units, and it costs nothing since scaling w and b together does not move the boundary). So we need

By symmetry w must point along (1,1) — the direction from the negative point to the positive one — so write w = (a, a). Substitute:

Subtract the second equation from the first: (6a + b) − (2a + b) = 1 − (−1), so 4a = 2 and a = 0.5. Substitute back into 6a + b = 1: 3 + b = 1, so b = −2. Therefore

Check both constraints: f(3,3) = 1.5 + 1.5 − 2 = +1 ✓. f(1,1) = 0.5 + 0.5 − 2 = −1 ✓. The boundary passes through (2,2), the midpoint. Exactly what a human would draw.

The margin, and why the objective is what it is

The distance from a point to the hyperplane is |f(x)| / ‖w‖. Our support points have |f| = 1, so each sits 1/‖w‖ from the boundary and the full corridor is

Here ‖w‖ = √(0.25 + 0.25) = √0.5 = 0.7071, so the margin is 2 / 0.7071 = 2.828. Sanity check: the two points are √((3−1)² + (3−1)²) = √8 = 2.828 apart, and with only two points the widest corridor is obviously the whole gap between them. The formula agrees with the picture.

Maximizing 2/‖w‖ is the same as minimizing ‖w‖, which is the same as minimizing ½‖w‖² (nicer derivative). Hence the canonical statement:

Read it as English: make the weight vector as small as possible while still pushing every training point at least one unit onto its own side. Small w means a gentle slope means a wide corridor.

Soft margins: real data overlaps

Real audio classes overlap — a shout and a scream share feature space. With the hard constraint above there is no solution. So allow violations and pay for them, using the hinge loss

Compute a few hinge values to feel it. A point with y = +1 and f = 2.5 is correct and beyond the margin: loss max(0, 1 − 2.5) = 0 — it contributes nothing and could be deleted. A point with f = 0.4 is on the right side but inside the corridor: loss max(0, 1 − 0.4) = 0.6. A point with f = −0.3 is misclassified: loss max(0, 1 + 0.3) = 1.3.

C is the exchange rate between corridor width and violations. Small C: a wide corridor, many points allowed inside, a smooth boundary that ignores individuals — high bias. Large C: violations are expensive, so the boundary contorts to get every training point right — high variance, and with noisy audio features that means memorizing the noise. C is chosen by cross-validation, always, typically over a log grid from 0.01 to 1000.

Non-linear boundaries without non-linear maths

Audio classes are not linearly separable in a pooled feature space. The classical fix is to map the features into a higher-dimensional space where they are separable, and draw a straight line there — which is a curve back in the original space.

Doing that literally is expensive: a degree-2 map of a 78-dimensional vector has 3081 dimensions; degree 3 has over 90,000. The kernel trick avoids building the space at all. The key observation is that the SVM solution can be written so the data appears only inside dot products:

So if some cheap function K(u, v) happens to equal φ(u) · φ(v) for a rich map φ, we can substitute K for the dot product and get the rich space for free.

Worked by hand: the kernel trick is not magic

Take the degree-2 polynomial kernel K(u, v) = (u · v)² and the two vectors u = (1, 2), v = (3, 1).

The cheap way. u · v = 1×3 + 2×1 = 3 + 2 = 5. Square it: K(u, v) = 5² = 25. Two multiplications, one addition, one squaring.

The expensive way. Expand the map explicitly. For two dimensions,

Compute both images: φ(u) = (1², √2×1×2, 2²) = (1, 2.8284, 4) and φ(v) = (3², √2×3×1, 1²) = (9, 4.2426, 1). Their dot product, term by term:

Identical. We obtained the dot product in a three-dimensional expanded space without ever building the three-dimensional vectors. In 78 dimensions the saving is between 3081 numbers and one. The √2 in the map is exactly what makes the cross-term appear twice in the expansion of (u·v)² — the kernel dictates the map, not the other way round.

The RBF kernel, and its one knob

For our vectors with γ = 0.5: ‖u − v‖² = (1−3)² + (2−1)² = 4 + 1 = 5, so K = e^−2.5 = 0.0821. Two identical points give K = 1; distant points give K → 0. The RBF kernel is a similarity that decays with distance, and its implicit feature space is infinite-dimensional.

Read the classifier f(x) = ∑α_iy_iK(x_i, x) + b and you will see something familiar: it is a weighted vote of the support vectors, with weights that fall off with distance. An RBF SVM is a k-NN that learned which neighbours matter and how much.

Fifty classes from a two-class machine

An SVM is binary. ESC-50 has fifty classes. Two standard reductions:

LIBSVM, the library behind almost every classical result you will read, uses one-vs-one. It sounds absurd — 1225 classifiers — but each sees only 2/50 of the data, and training cost grows faster than linearly with dataset size, so the total is usually cheaper than 50 full-data problems.

Why SVMs sat at the clip level and GMMs at the frame level

This is the structural fact that explains the whole era’s architecture. Kernel SVM training solves a quadratic program over an n × n kernel matrix: memory grows as n², time as roughly n² to n³. With n = 2000 pooled clip vectors, the kernel matrix is 4 million entries — trivial. With n = 10 million frames (a hundred hours of speech), it is 10¹⁴ entries — impossible.

GMM training, by contrast, is linear in the number of points and streams beautifully. So the division of labour was forced by computation: GMMs for frame-level acoustic modelling, where the data is enormous, and SVMs for clip-level classification, where the data is small and the boundary is what you want. Not fashion — arithmetic.

Code: two forms

python
# FORM 1 - the kernel trick, verified against the hand computation.
import numpy as np

u, v = np.array([1., 2.]), np.array([3., 1.])
K_cheap = (u @ v) ** 2                                  # 25.0
phi = lambda x: np.array([x[0]**2, np.sqrt(2)*x[0]*x[1], x[1]**2])
K_explicit = phi(u) @ phi(v)                            # 25.000000000000004
print(K_cheap, K_explicit)                            # the trick, demonstrated

# The full RBF kernel matrix, vectorised (no loops):
def rbf(X, Y, gamma):
    d2 = (X**2).sum(1)[:, None] + (Y**2).sum(1)[None, :] - 2 * X @ Y.T
    return np.exp(-gamma * d2)                          # (n, m)

# FORM 2 - the practical recipe: standardize, grid-search, cross-validate.
from sklearn.svm import SVC
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import GridSearchCV, StratifiedKFold

pipe = make_pipeline(StandardScaler(), SVC(kernel='rbf'))
grid = {'svc__C': [0.1, 1, 10, 100], 'svc__gamma': [1e-3, 1e-2, 1e-1, 1]}
cv   = GridSearchCV(pipe, grid, cv=StratifiedKFold(5, shuffle=True, random_state=0))
cv.fit(X_train, y_train)
print(cv.best_params_, cv.score(X_test, y_test))
print(len(cv.best_estimator_[-1].support_), 'support vectors of', len(X_train))

Chapter 8: Hidden Markov Models & the Forward Trellis

Classify each frame independently with your best GMM and watch the output: speech, speech, music, speech, speech, music, music… It flickers. No physical sound alternates between categories at 100 Hz. The classifier is treating consecutive frames as unrelated events when they are anything but.

Two things are missing. First, a prior on continuity: whatever was happening 10 ms ago is overwhelmingly likely to still be happening. Second, the ability to score and decode a whole sequence rather than a bag of frames. The hidden Markov model supplies both, and it is the single most important model of the classical era.

The model, in three tables

An HMM says: there is a hidden state that evolves over time, and each state emits an observation. We see the observations, never the states — hence “hidden.” It is specified by three things:

In a real system b_i is a Gaussian mixture over the 39-dimensional MFCC vector — the GMM from Chapter 6, slotted in as the emission model. That composite is the GMM-HMM, and it was the entire field of speech recognition for two decades. Here we use two discrete symbols so the arithmetic fits on a page.

Our worked example

Two hidden states, S (speech) and M (music). Two possible observations per frame: H (high ZCR) and L (low ZCR). The parameters:

Read the transitions: speech stays speech 70% of the time, so the model expects runs, not flicker. Read the emissions: speech usually shows high ZCR (consonant-rich), music usually low. Observed sequence:

The question — how likely is this sequence under this model? — is what lets us compare a “speech/music” HMM against a “traffic” HMM and pick a winner.

The naive answer, and why it is impossible

The sequence could have been produced by any state path. There are 2³ = 8 paths here, and we could compute the probability of each and add. Let us actually do two of them so the structure is concrete.

Path S→S→S: start in S (0.6) and emit H (0.8); stay in S (0.7) and emit L (0.2); stay in S (0.7) and emit L (0.2).

Path S→M→M: start in S (0.6), emit H (0.8); move to M (0.3), emit L (0.7); stay in M (0.6), emit L (0.7).

All eight, for completeness:

Their sum is 0.10248. That is the answer — and the method is useless, because the number of paths is N^T. For a modest 3-state model over a one-second clip (T = 100 frames) that is 3¹⁰⁰ ≈ 5 × 10⁴⁷ paths. You cannot enumerate them, ever.

The forward algorithm: share the work

Look at the eight paths again. Paths S→S→S and S→S→M share their first two steps and recompute them independently. That waste is the whole problem, and the fix is to compute each shared prefix once.

Define the forward variable

It bundles every path that reaches state i at time t into a single number. The recursion follows from the Markov property — given the current state, the past is irrelevant to the future:

In words: to be in state j now, you must have been in some state before (sum over i), moved here (a_ij), and then emitted what we saw (b_j). The final answer is P(O) = ∑_i α_T(i).

Worked by hand: every alpha

t = 1, observation H.

t = 2, observation L. Incoming to S: from S with 0.7, from M with 0.4.

Incoming to M: from S with 0.3, from M with 0.6.

Already interesting: at t = 1 speech was four times more likely than music; after one low-ZCR frame, music has overtaken it. The evidence flipped the belief, and the transition prior slowed the flip down — that is the smoothing we wanted.

t = 3, observation L.

Exactly the brute-force sum — but computed with 8 multiply-adds instead of 8 full path expansions, and the gap explodes with T. The cost is O(T · N²) instead of O(N^T): for the 3-state, 100-frame case, 900 operations instead of 5 × 10⁴⁷.

Decoding: which path, not just how likely

Replace the sum with a max and remember where the max came from, and the same trellis gives you the single most likely state sequence — the Viterbi algorithm:

For our example the best path is S→M→M with probability 0.042336, which you can read straight off the eight-path table. Note that its probability, 0.042, is far below the total 0.102: the single best path accounts for only 41% of the sequence probability. That difference is exactly why the forward algorithm (sum) and Viterbi (max) are different algorithms answering different questions — “how likely is this sound at all” versus “what was happening when.”

How a real speech recognizer is wired

The 1990s architecture, in numbers, so the scale is concrete:

Recognition is then a single enormous Viterbi search through this composite graph, with beam pruning to keep it tractable. Training uses the Baum-Welch algorithm, which is exactly EM from Chapter 6 with the responsibilities computed by a forward pass and a matching backward pass. Everything in this lesson composes.

Underflow: the practical detail that bites everyone

α is a product of probabilities. After 1000 frames of a real system, values around 0.1 per frame give α ≈ 10⁻¹⁰⁰⁰ — which is exactly 0.0 in double precision, whose smallest normal value is about 10⁻³⁰⁸. Every implementation must either work in logs (replacing multiplication with addition and the sum with a log-sum-exp) or rescale α to sum to 1 at each step and accumulate the log of the scale factors. This is not an optimization; without it the algorithm returns zero and your classifier compares zero to zero.

Code: three forms

python
# FORM 1 - the recursion exactly as computed by hand.
import numpy as np

pi = np.array([0.6, 0.4])
A  = np.array([[0.7, 0.3], [0.4, 0.6]])
B  = {'H': np.array([0.8, 0.3]), 'L': np.array([0.2, 0.7])}

def forward(obs):
    a = pi * B[obs[0]]                     # alpha_1 -> [0.48, 0.12]
    for o in obs[1:]:
        a = (a @ A) * B[o]                 # sum over i, then emit
    return a.sum()

print(forward(['H', 'L', 'L']))            # 0.10248 - the paper answer

# FORM 2 - log domain, the version you must ship. No underflow, ever.
from scipy.special import logsumexp

def forward_log(obs, logpi, logA, logB):
    a = logpi + logB[obs[0]]
    for o in obs[1:]:
        a = logsumexp(a[:, None] + logA, axis=0) + logB[o]
    return logsumexp(a)                       # log P(O)

# FORM 3 - Viterbi: same trellis, max instead of sum, plus backpointers.
def viterbi(obs):
    d = pi * B[obs[0]]; back = []
    for o in obs[1:]:
        scores = d[:, None] * A                # (from, to)
        back.append(scores.argmax(0))
        d = scores.max(0) * B[o]
    path = [int(d.argmax())]
    for bp in reversed(back):
        path.append(int(bp[path[-1]]))
    return path[::-1], d.max()

print(viterbi(['H', 'L', 'L']))            # ([0, 1, 1], 0.042336)  =  S M M

Chapter 9: The Pipeline & the Baseline Era

Everything is built. Now we bolt it together, run it on the datasets people actually competed on, and put honest numbers next to it — including the numbers that show where it stopped.

The missing link: statistics pooling

After Chapter 4 a clip is an array of shape (T, 39), and T depends on how long the recording was. The classifiers of Chapters 5–7 need a single fixed-length vector. The classical answer is statistics pooling: summarize each feature dimension over time.

The mean says where this clip typically sits; the standard deviation says how much it moves. A steady drone and a rapid alternation can have the same mean and very different spread — so the second statistic is doing real work, not decoration.

Worked by hand: pooling four frames

Take the ZCR values of four consecutive frames:

Mean. 0.42 + 0.30 = 0.72; + 0.55 = 1.27; + 0.33 = 1.60. Divide by 4:

Deviations. 0.42 − 0.40 = +0.02; 0.30 − 0.40 = −0.10; 0.55 − 0.40 = +0.15; 0.33 − 0.40 = −0.07. (Check: they must sum to zero. +0.02 − 0.10 + 0.15 − 0.07 = 0.00 ✓.)

Squares. 0.0004, 0.0100, 0.0225, 0.0049. Sum = 0.0378.

Variance and sd. Population variance divides by n: 0.0378 / 4 = 0.00945, so sd = √0.00945 = 0.0972. Sample variance divides by n − 1: 0.0378 / 3 = 0.0126, so sd = 0.1122. Both conventions appear in the literature; with T in the hundreds the difference is under half a percent, but with short clips it is not, and mixing conventions between training and test is a real (if embarrassing) source of degradation.

The complete shape trace

One 4-second clip, all the way through, with the transformation and the loss at each step:

64,000 numbers to 78: a compression of 820×. Each throw-away was defensible when it was made. Read the third column downward and you are reading, in advance, the list of things a convolutional network would later be praised for keeping.

Evaluating honestly

The protocol matters as much as the model, and this era learned that the hard way.

Stratified k-fold cross-validation. Split the clips into k folds preserving class proportions, train on k−1, test on the held-out one, rotate, average. With only 40 clips per class you cannot afford a single train/test split — the variance of the estimate would swamp the differences you are trying to measure.

Group by source recording. ESC-50’s 2000 clips are cut from a smaller number of original Freesound recordings, and the dataset ships five predefined folds constructed so that clips from the same source file never straddle a fold boundary. Ignore them, shuffle randomly, and your model can recognise the background hum of a particular recording rather than the sound class. Reported accuracies jump by double digits, and the system fails completely in the field. This is the single most important methodological rule in audio classification, and it is violated constantly.

Report the confusion matrix, not just accuracy. With 50 classes, an accuracy number hides which pairs are being confused — and the confusions are diagnostic: helicopter/chainsaw confusions mean your features cannot separate broadband periodic textures, which is a feature problem, not a classifier problem.

Worked by hand: why accuracy lies on unbalanced data

Suppose a three-class detector is tested on 100 clips: 80 of class A (background), 15 of class B (alarm), 5 of class C (glass break). The confusion matrix comes back as

Plain accuracy counts the diagonal: (72 + 3 + 0) / 100 = 75%. That sounds like a working system. Now compute the per-class recalls: A = 72/80 = 0.900, B = 3/15 = 0.200, C = 0/5 = 0.000. The balanced accuracy is their mean:

The same system is 75% or 37% depending on which number you report, and the second one is the truth about what it does: it never detects glass breaking, which is the only class anybody cared about. Always report per-class recall, and prefer balanced accuracy (or macro-F1) whenever the classes are unbalanced. ESC-50 is perfectly balanced — 40 clips per class — so plain accuracy is honest there, which is part of why it became the standard benchmark.

Pooling variants worth knowing

Every one of these is still permutation-invariant. Note the fifth row carefully: pooling the deltas smuggles a little temporal information past the pooling stage, because the deltas were computed before order was destroyed. That trick — compute time-sensitive quantities early, pool them late — is the best the bag-of-frames framework can do, and it is worth several points in practice.

The numbers, as reported

ESC-50 (Piczak, 2015): 2000 five-second clips, 50 environmental classes, 40 clips each, five predefined folds. The paper’s own baselines use MFCC and zero-crossing-rate statistics pooled per clip — exactly the pipeline of this lesson:

Sit with the first four rows. The best hand-feature classifier reached about 44% where humans reach about 81% — a gap of nearly a factor of two, on a task humans find easy. That gap is the honest summary of the classical era on general environmental sound.

DCASE acoustic scene classification. The 2013 challenge shipped an MFCC-plus-GMM baseline that scored around 55% on ten scene classes, with the best submitted systems reaching roughly 76%. By the 2016 edition (15 classes) the MFCC-GMM baseline was around 72–77% and the winning entries approached 90%, using i-vectors and early neural fusion. The baseline barely moved in three years; the winners moved a lot, and what changed in the winners was the representation.

Speech. Here the classical pipeline was not a weak baseline — it was the state of the art for twenty years, and GMM-HMM systems delivered usable dictation. The 2012 turn came when deep networks replaced the GMM emission model while keeping the HMM: across several groups and benchmarks, the reported word-error-rate reductions were on the order of a third relative. It is worth noticing which block was replaced first. Not the features, not the sequence model — the density estimator.

Code: the entire system, end to end

python
# A complete, runnable classical audio classifier. Every function is a chapter of this lesson.
import numpy as np, librosa, glob, os
from sklearn.svm import SVC
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline
from sklearn.model_selection import GroupKFold, cross_val_score

SR, N_FFT, HOP = 16000, 400, 160          # Ch 1: 25 ms frames, 10 ms hop

def clip_features(path):
    x, _ = librosa.load(path, sr=SR, mono=True)               # (n_samples,)
    # --- Ch 2-4: per-frame features ------------------------------------
    mf  = librosa.feature.mfcc(y=x, sr=SR, n_mfcc=13, n_fft=N_FFT, hop_length=HOP)
    mf  = mf - mf.mean(axis=1, keepdims=True)                 # CMN: remove the channel
    d1  = librosa.feature.delta(mf)
    d2  = librosa.feature.delta(mf, order=2)
    zcr = librosa.feature.zero_crossing_rate(x, frame_length=N_FFT, hop_length=HOP)
    cen = librosa.feature.spectral_centroid(y=x, sr=SR, n_fft=N_FFT, hop_length=HOP)
    rol = librosa.feature.spectral_rolloff(y=x, sr=SR, n_fft=N_FFT, hop_length=HOP)
    flt = librosa.feature.spectral_flatness(y=x, n_fft=N_FFT, hop_length=HOP)
    F   = np.vstack([mf, d1, d2, zcr, cen / SR, rol / SR, flt])       # (43, T)
    # --- Ch 9: statistics pooling --------------------------------------
    return np.concatenate([F.mean(axis=1), F.std(axis=1)])            # (86,)

files  = sorted(glob.glob('ESC-50/audio/*.wav'))
X      = np.stack([clip_features(f) for f in files])                  # (2000, 86)
y      = np.array([int(os.path.basename(f).split('-')[3][:-4]) for f in files])
folds  = np.array([int(os.path.basename(f).split('-')[0]) for f in files])   # official folds!

# Ch 5-7: standardize INSIDE the pipeline, evaluate on the official folds only.
clf    = make_pipeline(StandardScaler(), SVC(kernel='rbf', C=10, gamma='scale'))
scores = cross_val_score(clf, X, y, groups=folds, cv=GroupKFold(5))
print(scores.mean().round(3))     # expect the high 0.3s / low 0.4s - the era's ceiling

That script is about forty lines and it reproduces, to within a couple of points, the numbers people published. There is no hidden trick and no missing ingredient. The ceiling is structural.

One detail in it repays attention: groups=folds. That single keyword is the difference between a reproducible 40% and a fraudulent 60%. Everything else in the script — the features, the kernel, the pooling — is a modelling choice you can argue about. The grouping is not a choice; it is the definition of the question being asked.

Chapter 10: The 2008 Machine (showcase)

Here is the whole lesson, running. A sound is generated, framed, described by features frame by frame, pooled, and classified — live, with every intermediate stage on screen. Nothing in this simulation is faked: the features are the arithmetic from Chapters 2 and 3 applied to synthetic frames, the decision regions are the actual classifier evaluated on a grid, and the verdict is what that classifier says.

What you are looking at

Why these two features and not the other five? Because zero-crossing rate and spectral centroid are almost independent: ZCR responds to noisiness, centroid to brightness, and a sound can be any combination of the two. Plotting two correlated features would give a diagonal smear in which every class sits on one line, which teaches nothing. Choosing axes that vary independently is the same instinct the DCT formalized in Chapter 4, applied here to a picture instead of a model.

Five experiments to run right now

A simulation you only watch teaches nothing. Do these in order; each one demonstrates a claim made earlier in the lesson.

Experiment 1 — the trajectory matters, the pooled point does not know it. Select siren and press Play. Watch the small frame dots sweep upward across the plot as the siren rises in pitch: that is genuine temporal structure, plainly visible. Now look at the big pooled dot: it sits placidly in the middle of the sweep and never moves along it. Everything the classifier receives is that one dot. This is the Chapter 9 blindness, in the one place where you can see both things at once.

Experiment 2 — noise is a translation, not just a blur. Set the sound to speech and raise the noise slider from 0.1 to 0.7. The frame dots do not merely scatter — they migrate right and up, because added broadband noise raises both the zero-crossing rate and the spectral centroid. The pooled point walks out of the speech region and into the rain region. Your classifier has not been confused by randomness; it has been given features that genuinely moved. That distinction matters when you debug: this failure is fixed by a better feature (or noise-robust normalization), not by a bigger classifier.

Experiment 3 — how little data a Gaussian needs. Set the classifier to Gaussian and pull the training-clips slider down to 3 per class. The decision regions stay broadly sensible: a Gaussian estimates two numbers per dimension and can do that from three points. Now switch to nearest neighbour with the same three points and press “Resample training set” a few times: the regions lurch wildly between samples. That is the bias-variance trade-off, made visual. A parametric model imposes a shape and is stable when data is scarce; a memory-based model imposes nothing and is at the mercy of which points it happened to see.

Experiment 4 — the boundary between the two hardest classes. Select rain and then scream, watching where each lands. Rain is high-ZCR and high-centroid; a scream is loud and high-centroid but far more periodic, so its ZCR is lower. In this two-feature space they are adjacent, and with noise at 0.5 they overlap. In a real system you would fix this by adding a feature that separates them structurally — spectral flatness, which is near 1 for rain and near 0 for a scream’s harmonic stack. Note what “fixing it” meant: a human thought about the physics and designed a number. That is the era’s method and its limit.

Experiment 5 — break it on purpose. Set noise to 1.0. Every class collapses towards the same corner of the plot, and the verdict becomes essentially arbitrary while the score bars stay confidently tall. A generative classifier reports the relative likelihood of classes; nothing in it says “this clip is unlike everything I was trained on.” To get that, you would threshold the absolute log-likelihood — the open-set rejection from Chapter 6 — and this experiment is why you would bother.

Worked by hand: what the score bars are computing

The bars at the bottom of the simulation are diagonal-Gaussian log-likelihoods. Let us compute two of them by hand so nothing is mysterious. Suppose the pooled point for the current clip is

and two of the class models are

With a diagonal covariance the log-likelihood is a sum over dimensions:

Speech, dimension 1. The standardized distance is (0.30 − 0.16)/0.05 = 2.800, so the quadratic term is 2.800²/2 = 3.920. The constant term is −½ log(2π × 0.0025) = −½ log(0.015708) = −½(−4.1533) = +2.077. Contribution: 2.077 − 3.920 = −1.843.

Speech, dimension 2. (3.4 − 2.0)/0.6 = 2.333; quadratic term 2.722; constant −½ log(2π × 0.36) = −½(0.8163) = −0.408. Contribution: −0.408 − 2.722 = −3.130.

Rain, dimension 1. (0.30 − 0.35)/0.06 = −0.833; quadratic 0.347; constant −½ log(2π × 0.0036) = +1.895. Contribution: +1.547. Dimension 2. (3.4 − 4.2)/0.9 = −0.889; quadratic 0.395; constant −0.814. Contribution: −1.209.

Rain wins by 0.339 − (−4.974) = 5.312 nats. With equal class priors, that is a posterior odds ratio of e^5.312 ≈ 203 to 1. Note that a log-likelihood can be positive: this is a density, not a probability, and a tight Gaussian (small σ) can exceed 1 at its peak. If you ever see a “probability” above 1 in an audio classifier, this is why, and it is not a bug.

How the playground is implemented

javascript
// The classifier behind the shaded regions - the same maths as the hand computation.
function logGauss(x, y, m) {                 // m = {mx, my, sx, sy}
  const a = -0.5 * Math.log(2 * Math.PI * m.sx * m.sx) - (x - m.mx) ** 2 / (2 * m.sx * m.sx);
  const b = -0.5 * Math.log(2 * Math.PI * m.sy * m.sy) - (y - m.my) ** 2 / (2 * m.sy * m.sy);
  return a + b;                              // diagonal covariance = sum over dims
}

// Fit: the M-step of a one-component GMM is just the sample mean and variance.
function fitClass(points) {
  const n = points.length;
  const mx = points.reduce((s, p) => s + p[0], 0) / n;
  const my = points.reduce((s, p) => s + p[1], 0) / n;
  const vx = points.reduce((s, p) => s + (p[0] - mx) ** 2, 0) / n;
  const vy = points.reduce((s, p) => s + (p[1] - my) ** 2, 0) / n;
  // variance floor - Chapter 6's singularity guard, needed here too with 3 points/class
  return { mx, my, sx: Math.sqrt(Math.max(vx, 1e-4)), sy: Math.sqrt(Math.max(vy, 1e-3)) };
}

// The decision region is just argmax evaluated on a coarse grid, one colour per winner.
for (let gx = 0; gx < GRID; gx++)
  for (let gy = 0; gy < GRID; gy++)
    region[gx][gy] = argmax(models.map(m => logGauss(xOf(gx), yOf(gy), m)));

Twenty lines for a working classifier, because a diagonal Gaussian has a closed-form fit — no gradient descent, no iterations, no learning rate. That closed form is a large part of why this family dominated when compute was scarce: you could fit a model faster than you could load the data.

The same budget argument governs inference. Per frame, the playground computes a few hundred arithmetic operations for the features and five Gaussian evaluations for the scores — call it a thousand operations at 100 frames per second, so 100,000 operations per second of audio. A 2008 laptop ran roughly a billion operations per second, so the pipeline used about one ten-thousandth of the machine and could comfortably run on dozens of live microphone streams at once. That is not a historical curiosity: it is why energy-and-ZCR voice detection still runs on the always-listening core of your phone while the neural model stays asleep.

Finally, notice what the playground does not need: a training loop, a validation curve, an early stopping criterion, a learning-rate schedule, or a GPU. Press “Resample training set” and the model is refitted instantly, from scratch, in front of you. Whatever else it lacked, this era had a development loop measured in seconds.

What a real submission added on top

The playground uses two features and one classifier for legibility. A competitive 2010 DCASE entry differed in scale, not in kind:

None of those additions changes the shape of the pipeline. They are all “more of the same, tuned harder,” which is exactly the profile of a mature technology approaching its ceiling. The next chapter is about what that ceiling was made of.

Ablations: what each stage is actually worth

If you build the real system, these are the moves that change the number, roughly in order of how much they change it. The magnitudes below are typical of what you will see on a small environmental-sound task; treat them as expectations to test, not constants.

A debugging checklist

When a classical audio classifier underperforms, the cause is almost always in this list, and almost never in the classifier:

Chapter 11: The Plateau & the Handoff

This pipeline was refined for a quarter of a century by thousands of very good engineers. Then, between roughly 2012 and 2016, it was overtaken across every audio task at once. It is worth being precise about why, because “deep learning is better” is not an explanation — it is a restatement of the observation.

There are four structural reasons, and each of them is a specific decision we made earlier in this lesson, coming due.

Reason 1: the representation cannot learn

In the classical stack, exactly one block adapts to your task: the classifier at the end. The framing, the filterbank, the log, the DCT, the truncation, the pooling — all fixed before any data is seen. If your task needs a feature nobody thought of, there is no mechanism by which the system can discover it. The only path is a human having an idea.

In a convolutional network on log-mel input, the first layer is a learned filterbank, the middle layers are learned pattern detectors, the pooling is learned, and the classifier is learned — and the error signal from a misclassified clip flows back through all of them. The difference is not that nets are “more powerful” in the abstract; it is that gradient reaches the representation.

Reason 2: MFCCs were designed for speech, then used everywhere

Every discard in the MFCC pipeline is correct for speech recognition:

A feature set is a statement about what does not matter. MFCCs are an excellent statement about speech and a poor one about sirens, and the field used them for both for fifteen years because they were what existed.

Reason 3: no local time-frequency patterns

This is the deepest one, and the simulation below is built to show it. Look at a spectrogram and the diagnostic evidence is often a two-dimensional shape: a rising diagonal streak (a chirp), a vertical line followed by horizontal bands (a struck note), a checkerboard of harmonics moving together (a voice). These are patterns in time and frequency jointly.

Per-frame features collapse each column to a few numbers before any of that shape can be measured, and statistics pooling then discards the column order. There is no stage in the classical pipeline at which a 2-D pattern could be detected even in principle. A convolutional kernel sliding over the spectrogram does nothing else.

Reason 4: it could not absorb more data

Give a GMM-plus-MFCC system ten times more training data and it improves a little: the parameter count is fixed by your choice of K, and beyond a few hours the estimates are already converged. Give a neural network ten times more data and you can also make it ten times larger, and it improves a lot. Once datasets like AudioSet (about two million clips, 2017) existed, the classical stack had no way to use them.

That is the honest ordering of events. The deep models did not merely have better inductive biases; they arrived at the same time as datasets and hardware that rewarded scale, and the classical pipeline had no scale dimension to grow along.

What the handoff kept, and what it dropped

Read the “kept” rows: everything up to and including the log survived unchanged. The front end of a 2024 audio model is the front end of a 1990 speech recognizer. What was replaced is precisely everything that was lossy and hand-decided — and everything that was lossless and perceptually motivated stayed. That is a remarkably clean verdict on twenty-five years of design work: the physics was right, the summarization was wrong.

The handoff, year by year

Notice the order in which the blocks fell: emission model, then sequence model, then the features. The front end — framing, FFT, mel, log — never fell at all. Technologies are not replaced whole; they are dismantled from whichever end is most obviously arbitrary.

The honest counterargument

It would be tidy to say the classical pipeline was simply worse. It was not, and pretending otherwise makes you a worse engineer.

For most of the 2010s, on tasks with a few hundred labelled clips, a well-tuned MFCC-plus-SVM system beat a from-scratch convolutional network, because the net had nowhere near enough data and the hand features encoded decades of correct prior knowledge for free. What changed the verdict was not architecture alone; it was transfer — the ability to pretrain on AudioSet or on unlabelled audio and arrive at your small task with a representation already formed. Before pretraining was routine, hand-crafted priors were competitive with learned ones, exactly as theory predicts when data is scarce.

The lasting lesson is therefore not “never design features.” It is: a hand-designed feature is a prior, priors dominate when data is scarce, and data stopped being scarce.

What is still classical in production today

That last row is the practical reason this lesson is not history. Before fine-tuning a pretrained audio transformer on your 400 labelled clips, spend an afternoon on the classical pipeline. If it gets 70%, your problem is easy and you have a deployable model already. If it gets chance, your labels or your recordings are broken, and you have learned that for the price of an afternoon instead of a week of GPU time.

Choosing among the three classifiers, for good

If you remember one row, remember the likelihood row. A system that can say “I have never heard anything like this” behaves completely differently in deployment from one that must pick a label, and that capability came from modelling the density rather than the boundary. Most modern classifiers gave it up, which is why open-set audio recognition is a live research problem again.

Classical or learned? A decision procedure

You will actually face this choice, so here it is as a procedure rather than an opinion. Work down the list and stop at the first row that matches your situation.

The last row is the one people skip. If you can describe in words what makes your two classes different, there is often a five-line feature that captures it, and that feature will be more robust, more debuggable and a thousand times cheaper than anything you could train. The classical era’s real bequest is not MFCCs; it is the habit of listening to the data first and asking what actually differs.

Cheat sheet

One last calibration before you leave. Everything in this lesson — twelve chapters, seven features, three classifiers, a sequence model — adds up to a system that reaches roughly 44% where humans reach 81% on general environmental sound. That is not a failure. It is a measurement of exactly how much of hearing is signal processing and how much is learned experience, taken by a generation of engineers who had no way to supply the second half. When you read a modern paper reporting 95%, you are reading what happened when the second half arrived.

Where to go next

“Far better an approximate answer to the right question, which is often vague, than an exact answer to the wrong question, which can always be made precise.”
— John Tukey, who also gave us the FFT, the cepstrum’s name, and most of this lesson’s vocabulary