Audio signal processing: spectral analysis, filters, and synthesis

Why this subject matters now

For a decade it looked as if deep learning would swallow signal processing whole: feed a network raw samples, let it learn its own filters, and retire the textbook. That is not what happened. The most effective audio models are the ones that keep the DSP scaffolding and make it differentiable. Neural vocoders predict a mel spectrogram, an explicitly windowed short-time Fourier magnitude, and invert it. Differentiable DSP (DDSP, from Google's Magenta group) puts a harmonic oscillator bank and a time-varying filtered-noise source, the century-old source-filter model, inside the autograd graph and learns their controls. Neural audio codecs quantize a learned latent but are trained against multi-resolution STFT losses. Room-acoustics simulators seed reverberators with the Schroeder topology. The practitioner who understands why a Hann window has a wider main lobe than a rectangular one, why an IIR filter can go unstable when a pole leaves the unit circle, and why the phase vocoder needs to unwrap phase, reads all of these systems as small variations on a fixed set of ideas rather than as a pile of unrelated tricks.

The fundamentals also decide correctness in ways a network cannot paper over. A one-line off-by-one in a hop size shifts every frame; a window applied without overlap-add compensation puts a comb filter in the output; a filter designed by naive impulse-response truncation ripples in the passband because nobody thought about the window. These are the bugs that ship, and they are invisible unless the person debugging knows the transform pair that predicts them. What a practitioner is expected to know today that they might not have five years ago is precisely how to differentiate through this machinery: how to write an STFT that backpropagates, why the phase is the hard part, and which loss functions on spectra actually correlate with what a listener hears.

  acoustic          A/D                  analysis                       synthesis
  pressure  --->  sample + quantize --> DFT / STFT / LPC / cepstrum --> filters, oscillators,
  wave x(t)       x[n]  (this page)     spectra, pitch, formants        delay lines, reverb
                    |                        |                              |
              Nyquist, SNR            leakage, resolution,           phase vocoder, Karplus-Strong,
              (sec. 1-2)              poles/zeros (sec. 3-6)         Schroeder, coding (sec. 7-10)
          

Sampling, aliasing, and quantization

The sampling theorem and where aliasing comes from

A microphone produces a continuous signal \( x(t) \). To store it we read its value every \( T \) seconds, producing the sequence \( x[n] = x(nT) \), where \( f_s = 1/T \) is the sample rate. The question the Nyquist-Shannon theorem answers is which continuous signals survive this: for which \( x(t) \) can the samples \( x[n] \) be inverted back to \( x(t) \) with no loss. The answer, proved by Shannon in 1949 building on Nyquist's 1928 analysis, is that a signal band-limited to frequencies below \( f_s/2 \) is reconstructed exactly by sinc interpolation of its samples, and no signal with content above \( f_s/2 \) is. The frequency \( f_s/2 \) is the Nyquist frequency, and \( f_s/2 \) is also the folding frequency, for the reason we now derive.

Consider two pure tones, \( x_1(t) = \cos(2\pi f_1 t) \) and \( x_2(t) = \cos(2\pi f_2 t) \), and suppose \( f_2 = f_1 + m f_s \) for some integer \( m \). Sample both at \( t = nT = n/f_s \):

because \( 2\pi m n \) is an integer multiple of \( 2\pi \) and cosine has period \( 2\pi \). The two distinct continuous tones produce byte-for-byte identical sample sequences. The same cancellation with \( \cos(-\theta)=\cos\theta \) shows that \( f_2 = -f_1 + m f_s \) also aliases onto \( f_1 \). So every frequency \( f \) is indistinguishable from the set \( \{\pm f + m f_s\} \), and the apparent frequency after sampling is whichever member of that set lands in the baseband \( [0, f_s/2] \). Geometrically the spectrum folds like an accordion about the multiples of \( f_s/2 \): a tone climbing past \( f_s/2 \) appears to descend again, which is the folding frequency.

Making this concrete, take \( f_s = 8000 \) Hz. Then \( 1000 \), \( 7000 \) (which is \( -1000 + 8000 \)), and \( 9000 \) (which is \( 1000 + 8000 \)) Hz all alias to \( 1000 \) Hz. Sampling \( \cos(2\pi\cdot 1000\, n/8000) \) and \( \cos(2\pi\cdot 9000\, n/8000) \) for \( n = 0,\dots,7 \) gives \( (1, 0.7071, 0, -0.7071, -1, -0.7071, 0, 0.7071) \) in both cases; the maximum difference between the two sequences is \( 3.2\times 10^{-15} \), pure floating-point noise (measured in NumPy 2.2.6 on this machine). This is why every A/D converter puts an analog anti-alias lowpass filter before the sampler: energy above \( f_s/2 \) must be removed in the continuous domain, because once it aliases there is no transform that can separate the \( 9000 \) Hz tone from the \( 1000 \) Hz tone. They are the same numbers.

Quantization noise and the six-decibels-per-bit rule

Sampling discretizes time; quantization discretizes amplitude. A \( b \)-bit converter maps the continuous amplitude onto \( 2^b \) levels spaced by a step \( \Delta \). Rounding to the nearest level introduces an error \( e[n] = \hat{x}[n] - x[n] \) bounded by \( |e| \le \Delta/2 \). The standard model, accurate when the signal is busy relative to \( \Delta \), treats \( e \) as a random variable uniform on \( [-\Delta/2, \Delta/2] \), independent of the signal. Its mean is zero and its power is the variance

Now put a full-scale sinusoid through the converter: amplitude \( A \), so it spans the full input range of width \( 2A \), and the step is \( \Delta = 2A / 2^b \). The signal power of a sinusoid is \( A^2/2 \). The signal-to-quantization-noise ratio is

In decibels, \( 10\log_{10}(\tfrac{3}{2} 2^{2b}) = 10\log_{10} 1.5 + 2b\cdot 10\log_{10} 2 \). Since \( 10\log_{10} 1.5 = 1.76 \) and \( 20\log_{10} 2 = 6.0206 \),

Each added bit halves the quantization step and buys \( 6.02 \) dB. Sixteen-bit CD audio gives a theoretical \( 98.08 \) dB; a full-scale \( 44.1 \) kHz sine quantized to 16 bits in NumPy measured \( 98.53 \) dB, the small surplus coming from the sine spending time near its peaks where the uniform-error assumption slightly under-counts. Eight bits give the model's \( 49.92 \) dB and measured \( 50.00 \) dB; twenty-four bits give \( 146.24 \) dB, comfortably below the thermal noise floor of any analog front end, which is why 24-bit converters are specified by their analog noise, not their bit depth. The \( 1.76 \) dB constant is specific to a full-scale sinusoid; a signal with a lower crest factor or one that does not reach full scale gives a smaller constant, which is why real converters quote SNR for a \( -1 \) dBFS or \( -3 \) dBFS tone.

The discrete Fourier transform and the FFT

From the continuous transform to the DFT

The continuous Fourier transform of \( x(t) \) is \( X(f) = \int_{-\infty}^{\infty} x(t) e^{-j2\pi f t}\, dt \). Three approximations turn this into something a computer can evaluate on a finite array. First, replace \( x(t) \) by its samples \( x[n] = x(nT) \) and the integral by a Riemann sum over samples; this gives the discrete-time Fourier transform (DTFT), \( X(e^{j\omega}) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n} \), a continuous, \( 2\pi \)-periodic function of the normalized angular frequency \( \omega = 2\pi f / f_s \). Second, keep only \( N \) samples, which as we will see multiplies the true DTFT by a window. Third, evaluate the DTFT only at the \( N \) equally spaced frequencies \( \omega_k = 2\pi k/N \), \( k = 0,\dots,N-1 \). What survives all three steps is the discrete Fourier transform,

Writing \( W_N = e^{-j2\pi/N} \) for the primitive \( N \)-th root of unity, \( X[k] = \sum_n x[n] W_N^{kn} \). Bin \( k \) corresponds to the physical frequency \( f_k = k f_s / N \), so the spacing between adjacent bins, the frequency resolution, is exactly

At \( f_s = 16 \) kHz an \( N = 512 \) transform resolves \( 31.25 \) Hz per bin, \( N = 1024 \) gives \( 15.625 \) Hz, and \( N = 2048 \) gives \( 7.8125 \) Hz (all exact). More samples in the analysis window means finer frequency resolution; that this must cost time resolution is the subject of the STFT section below.

Three properties that do all the work

Linearity is immediate from the sum: \( a x[n] + b y[n] \) transforms to \( a X[k] + b Y[k] \). The shift theorem says a delay is a linear phase ramp. Substituting \( m = n - n_0 \) into the transform of \( x[n - n_0] \) (indices modulo \( N \)),

A time shift leaves every magnitude \( |X[k]| \) untouched and rotates each phase by an amount proportional to both the bin and the delay. This is why the magnitude spectrogram is shift-tolerant and why phase carries the timing information, a fact the phase vocoder later exploits and neural vocoders struggle with. The convolution theorem is the reason the FFT matters for audio at all. Circular convolution in time is pointwise multiplication in frequency:

To convolve a length-\( L_x \) signal with a length-\( L_h \) filter linearly, zero-pad both to at least \( N \ge L_x + L_h - 1 \) so the circular wrap does not corrupt the result, transform, multiply, and invert. Convolving \( x = (1,2,3,4) \) with \( h = (1,-1,2) \) by direct summation gives \( (1,1,3,5,2,8) \); the length-6 FFT product inverts to the same vector to machine precision (verified). Direct convolution costs \( O(L_x L_h) \); the FFT route costs \( O(N \log N) \), which for long filters and reverb impulse responses tens of thousands of taps long is the difference between real time and not.

The radix-2 FFT and why the audio use is the point

Evaluated as written the DFT is a matrix-vector product with an \( N\times N \) dense matrix, \( O(N^2) \) operations. The Cooley-Tukey factorization (1965, and in Gauss's notebooks of 1805) splits the sum by the parity of the time index. With \( N \) even, separate the even-indexed and odd-indexed samples:

using \( W_N^{2rk} = e^{-j2\pi\cdot 2 rk/N} = e^{-j2\pi rk/(N/2)} = W_{N/2}^{rk} \). Here \( E \) and \( O \) are the \( (N/2) \)-point DFTs of the even and odd sub-sequences. The twiddle factor \( W_N^k \) has period \( N \), while \( E \) and \( O \) have period \( N/2 \), so the upper half of the spectrum reuses the same sub-transforms with a sign flip:

This is the decimation-in-time butterfly. Each of the \( N/2 \) values of \( k \) produces two outputs from one complex multiply and two adds. The work at size \( N \) obeys \( T(N) = 2 T(N/2) + \Theta(N) \), which unrolls to \( \Theta(N \log_2 N) \); the recursion tree has \( \log_2 N \) levels, each doing \( \Theta(N) \) butterfly work. A length-4096 FFT is roughly \( 4096\cdot 12 \approx 49{,}000 \) versus \( 4096^2 \approx 1.7\times 10^7 \) operations, a factor near 340. The derivation of the recurrence, the bit-reversal permutation, and the numerical-stability analysis belong to algorithm design and analysis, which owns the FFT as a divide-and-conquer algorithm. What matters for audio is the consequence: because the FFT made the DFT cheap, every real-time analysis, the spectrogram, the phase vocoder, the perceptual coder's filterbank, fast convolution reverb, is built on repeated transforms of short frames. A direct-DFT spectrogram would not run in real time; an FFT one runs on a phone. As a sanity check, a length-4096 FFT of white noise agreed with the explicit matrix DFT \( W x \) to a maximum error of \( 1.4\times 10^{-10} \) (measured), confirming the fast route computes the same transform, only faster.

Windows, spectral leakage, and the short-time transform

Why a finite window leaks

Taking \( N \) samples of a signal is multiplying it by a rectangular window \( w[n] \) that is one on \( [0, N-1] \) and zero elsewhere. Multiplication in time is convolution in frequency, so the observed spectrum is the true spectrum convolved with the window's transform. The rectangular window's DTFT is the Dirichlet kernel

a narrow main lobe of width \( 4\pi/N \) (first nulls at \( \omega = \pm 2\pi/N \), i.e. two DFT bins on each side of center) flanked by side lobes whose first peak sits only \( 13.3 \) dB below the main lobe and rolls off slowly at \( 6 \) dB per octave. A single sinusoid at a frequency that lands exactly on a bin is captured cleanly, because the Dirichlet nulls fall on all the other bins. A sinusoid that lands between bins is smeared across all of them by these side lobes: this is spectral leakage. Concretely, at \( f_s = 1000 \) Hz with \( N = 64 \) (bin spacing \( 15.625 \) Hz), a \( 125 \) Hz tone lands exactly on bin 8 and its energy piles into that one bin, the ratio of the peak bin to its neighbour is astronomical. Move the tone to \( 128 \) Hz, bin \( 8.19 \), and the energy spills so that the peak bin exceeds its neighbour by a factor of only \( 4.48 \) (both measured). A weak spectral line near a strong off-bin line can be completely buried under the strong line's side lobes.

The main-lobe / side-lobe tradeoff

Tapering the window to zero at its edges suppresses the side lobes at the cost of a wider main lobe. The three classic tapers, with \( n = 0,\dots,N-1 \), are

The Hamming coefficients are chosen to cancel the first side lobe of the Dirichlet kernel, which is why its nearest side lobe is so low; Blackman adds a third cosine to push the side lobes lower still. Harris's 1978 survey tabulates dozens of windows by exactly these two numbers. Measuring the DTFT of a length-1024 window (zero-padded 8x) gives:

The side-lobe numbers here are measured (\( -13.4, -31.5, -42.7, -58.2 \) dB) and match the textbook values of roughly \( -13, -31, -43, -58 \) dB. The pattern is the whole story: there is no window that is both narrow and low, and the choice is dictated by whether the task is separating two close tones of similar strength (want a narrow main lobe, Hamming) or seeing a quiet tone beside a loud one (want low side lobes, Blackman). The Hann window is the default for the STFT because its main lobe and its \( 18 \) dB/octave roll-off are a good compromise and because it overlap-adds to a constant at \( 50\% \) or \( 75\% \) hop, which matters for reconstruction.

The STFT and the time-frequency uncertainty tradeoff

Music and speech are not stationary; their spectra change moment to moment. The short-time Fourier transform slides a window of length \( L \) along the signal in hops of \( H \) samples and transforms each frame:

The magnitude \( |X_t[k]|^2 \), laid out with frame index \( t \) on one axis and bin \( k \) on the other, is the spectrogram. Two resolutions are in tension. Frequency resolution is set by the window length, \( \Delta f \approx f_s / L \): a longer window packs more cycles of a tone and separates nearby frequencies better. Time resolution is also set by the window length, \( \Delta t = L / f_s \): a longer window blurs a transient over its whole span. Their product is fixed at one,

the discrete image of the Gabor-Heisenberg uncertainty principle, whose continuous form \( \Delta t\,\Delta f \ge \tfrac{1}{4\pi} \) says no signal can be arbitrarily localized in both domains at once, with equality only for the Gaussian window. A \( 46 \) ms window at \( 44.1 \) kHz (\( L = 2048 \)) resolves \( 21.5 \) Hz but smears a drum hit across \( 46 \) ms; an \( 8 \) ms window catches the transient but cannot tell a bass note from its neighbour a semitone away. This is why constant-Q and wavelet transforms exist, and why analysis pipelines often run two window sizes and fuse them. The hop \( H \) controls redundancy and reconstruction, not the resolution tradeoff: with a constant-overlap-add window, summing the inverse transforms of the frames recovers the signal exactly, so the STFT is an invertible, if overcomplete, representation.

Digital filters: the z-transform, poles, and zeros

FIR versus IIR, and the transfer function

A linear time-invariant filter is defined by a difference equation relating output to input,

If all \( a_k = 0 \) for \( k \ge 1 \), the output is a finite weighted sum of past inputs: a finite impulse response (FIR) filter, whose impulse response is just \( (b_0, \dots, b_L) \) and stops. If some \( a_k \neq 0 \), past outputs feed back and the impulse response rings on forever: an infinite impulse response (IIR) filter. The z-transform turns the difference equation into algebra. Define \( X(z) = \sum_n x[n] z^{-n} \); because a unit delay \( x[n-1] \) transforms to \( z^{-1} X(z) \) (by the same index-shift argument as the DFT shift theorem), the difference equation becomes \( A(z) Y(z) = B(z) X(z) \), and the transfer function is the ratio of two polynomials in \( z^{-1} \):

The roots of \( B \) are the zeros of the filter, frequencies it nulls; the roots of \( A \) are the poles, resonances it boosts. Evaluating \( H(z) \) on the unit circle \( z = e^{j\omega} \) gives the frequency response \( H(e^{j\omega}) \), whose magnitude is the gain and whose angle is the phase at each frequency. A pole near the unit circle at angle \( \theta \) produces a sharp resonant peak at \( \omega = \theta \); a zero on the circle produces a notch. This pole-zero picture is the working language of filter design: place poles where you want resonance, zeros where you want rejection, and read off the response.

Stability and the unit circle

A causal LTI filter is bounded-input bounded-output stable if and only if every pole lies strictly inside the unit circle, \( |p_i| < 1 \). The reason is direct: partial-fraction expansion writes the impulse response as a sum of terms \( C_i\, p_i^{\,n} \) for \( n \ge 0 \), and \( p_i^{\,n} \) decays to zero exactly when \( |p_i| < 1 \), grows without bound when \( |p_i| > 1 \), and rings forever when \( |p_i| = 1 \). FIR filters have all their poles at the origin (\( A(z) = 1 \)) and are therefore unconditionally stable, one of their two great virtues. IIR filters are cheaper for a given selectivity but must be checked: a coefficient rounded in fixed-point can push a pole across the circle and turn a filter into an oscillator. Consider the one-pole feedback \( y[n] = x[n] + a\, y[n-1] \), with pole at \( z = a \). For \( a = 0.9 \) the pole is at \( 0.9 \), inside the circle, and the filter is a stable leaky integrator. For \( a = 1.2 \) the pole is at \( 1.2 \), outside, and the output diverges geometrically. A second-order resonator with poles at radius \( r = 0.95 \) and angle \( \theta = \pi/4 \) has poles \( 0.672 \pm 0.672 j \), magnitude \( 0.95 < 1 \), stable, and rings with a decay set by how close \( r \) is to one (all three verified by root-finding on the denominator).

FIR design by windowed sinc, and linear phase

The ideal lowpass filter with cutoff \( f_c \) has a brick-wall frequency response and, by inverse DTFT, the impulse response

which is infinitely long and non-causal, its energy trailing off in both directions. The windowed-sinc method makes it realizable in two moves. Truncate it to \( M+1 \) taps and multiply by a window \( w[n] \) to control the passband ripple and stopband attenuation, then shift by \( M/2 \) to make it causal. The truncation alone (a rectangular window) gives the Gibbs ripple, roughly \( 9\% \) overshoot at the band edge no matter how long the filter; a Hamming or Blackman window trades a slightly wider transition for a stopband tens of decibels deeper. A \( 65 \)-tap Hamming-windowed sinc with \( f_c = 2000 \) Hz at \( f_s = 16 \) kHz measured its \( -6 \) dB point at exactly \( 2000 \) Hz with a worst-case stopband rejection of \( -60 \) dB beyond \( 3 \) kHz, and it agrees tap-for-tap with \( \texttt{scipy.signal. firwin} \) (verified).

The decisive property of a symmetric FIR filter, \( h[n] = h[M - n] \), is exactly linear phase. Factoring the symmetry out of the frequency response,

the bracket is real, so the phase is exactly \( -\omega M/2 \), a straight line, and the group delay \( -d\phi/d\omega = M/2 \) is constant across all frequencies. Every frequency is delayed by the same \( M/2 \) samples, so the waveform shape is preserved and no frequency-dependent smearing occurs. This is why FIR filters are used wherever phase matters, crossovers, hearing aids, anything where transients must not disperse, and it is a property IIR filters cannot have exactly. The price is length: a sharp FIR filter needs many taps, and its constant delay of \( M/2 \) samples can be substantial.

IIR design by the bilinear transform, and warping

The efficient way to get a sharp filter with few coefficients is to design a classical analog prototype, a Butterworth, Chebyshev, or elliptic filter with its poles laid out for a maximally flat or equiripple response, and map it to the digital domain. The bilinear transform is the substitution

It has exactly the properties a mapping needs. It sends the entire imaginary axis of the analog \( s \)-plane onto the unit circle of the \( z \)-plane once, so the whole analog frequency axis fits in the digital band without aliasing, and it sends the stable left half-plane \( \mathrm{Re}(s) < 0 \) into the interior of the unit circle, so a stable analog filter always maps to a stable digital one. Setting \( s = j\Omega \) (analog frequency) and \( z = e^{j\omega} \) (digital frequency) and simplifying,

This relation is the frequency warping: the analog axis \( \Omega \in (0,\infty) \) is compressed nonlinearly onto the digital axis \( \omega \in (0,\pi) \) by the tangent. The map is nearly linear for small \( \omega \) but stretches severely as \( \omega \to \pi \), so a filter's critical frequencies land in the wrong place unless they are pre-warped. The fix is to design the analog prototype not at the desired digital frequency but at \( \Omega_c = \tfrac{2}{T}\tan(\omega_c/2) \), so that after the inverse warp it comes back to \( \omega_c \). For a digital cutoff of \( 2000 \) Hz at \( f_s = 16 \) kHz, \( \omega_c = 2\pi\cdot 2000/16000 = \pi/4 \), and the pre-warped analog cutoff is \( \tfrac{2}{T}\tan(\pi/8) \), which corresponds to \( 2109.6 \) Hz (measured): the analog design is aimed \( 5.5\% \) high so the digital filter lands on target. A fourth-order Butterworth designed this way has poles at radii \( 0.758 \) and \( 0.458 \), all inside the unit circle (verified), confirming the transform preserved stability. The tradeoff against FIR is exactly the linear-phase property: the bilinear IIR filter reaches a given selectivity in a quarter the coefficients but has a nonlinear, frequency-dependent phase.

The phase vocoder: independent time and pitch scaling

Slowing music down without dropping its pitch, or transposing a voice without slowing it, are the same operation viewed two ways, and the phase vocoder (Flanagan and Golden 1966; Portnoff 1976 for the FFT implementation) is the classic solution. Its difficulty is entirely in the phase. The STFT gives each frame a complex value \( X_t[k] = |X_t[k]| e^{j\phi_t[k]} \) per bin. Resynthesizing frames at a different hop rate changes the timing, but naively reusing the analysis phases makes successive frames disagree about where each sinusoid is in its cycle, and the result is the smeared, phasey artifact of a broken time-stretch. The fix is to estimate the true instantaneous frequency in each bin and lay down synthesis phases that are consistent with it.

The estimate comes from phase unwrapping. A sinusoid sitting in bin \( k \) advances in phase by \( \omega_k H = \tfrac{2\pi k}{N} H \) radians per analysis hop of \( H \) samples if it is exactly at the bin center. The measured phase increment between consecutive frames is \( \Delta\phi = \phi_t[k] - \phi_{t-1}[k] \). Subtract the expected bin-center advance and wrap the remainder into \( (-\pi, \pi] \) with the principal-argument operator, since phase is only known modulo \( 2\pi \):

The wrapped residual is the deviation of the true frequency from the bin center, accumulated over \( H \) samples, so the instantaneous frequency in that bin is

Testing this on a pure \( 1100 \) Hz tone at \( f_s = 16 \) kHz with \( N = 1024, H = 256 \): the nearest bin is number \( 70 \), whose center is \( 1093.75 \) Hz, but the phase-unwrap estimate recovers \( 1100.00 \) Hz to two decimals (measured). The bin index alone is off by \( 6 \) Hz; the phase correction nails it. To time-stretch by a factor \( \alpha \), resynthesize with an analysis hop \( H_a \) and a synthesis hop \( H_s = \alpha H_a \), advancing the synthesis phase by the instantaneous frequency over the synthesis hop, \( \phi^{\mathrm{syn}}_t[k] = \phi^{\mathrm{syn}}_{t-1}[k] + \hat{\omega}_k H_s \), and inverting each frame with overlap-add. The magnitudes are copied unchanged; only the phase progression is rewritten. Pitch-shifting is then time-stretching by \( \alpha \) followed by resampling by \( 1/\alpha \): the stretch changes duration at constant pitch, the resample changes both back to the original duration, leaving pitch scaled. Laroche and Dolson (1999) improved the basic method by locking the phases of bins around each spectral peak to the peak's phase, which removes the residual phasiness by keeping each sinusoid's bins coherent, and phase-locked variants are what production time-stretchers use.

Pitch, formants, and linear prediction

Autocorrelation and cepstral pitch detection

A periodic signal repeats every \( P \) samples, so its autocorrelation \( r[\tau] = \sum_n x[n] x[n+\tau] \) has a strong peak at \( \tau = P \) (and at multiples of \( P \)). Finding the pitch is finding that lag: search \( r[\tau] \) over the plausible range of periods and take the peak, then \( f_0 = f_s / P \). On a synthetic glottal-like signal, a \( 200 \) Hz fundamental plus its second and third harmonics at \( f_s = 8000 \) Hz, the autocorrelation peaks at lag \( 40 \) samples, giving \( 8000/40 = 200.00 \) Hz exactly (measured). The failure mode is octave errors: because \( r[2P] \) is also large, a detector can lock onto twice the period, and practical algorithms (YIN, pYIN, the cumulative-mean-normalized difference function) are engineered specifically to suppress that.

The cepstrum offers a complementary route by separating source from filter multiplicatively. Voiced speech is a periodic glottal source convolved with the vocal-tract filter, so its log-magnitude spectrum is the sum of a slowly varying filter envelope and a rapidly oscillating harmonic comb. Taking the inverse transform of the log-magnitude spectrum, the real cepstrum

maps the slow envelope to low quefrency (small \( n \)) and the fast harmonic ripple to a sharp peak at the quefrency equal to the pitch period (Noll 1967). Reading the pitch off the high-quefrency peak and the formants off the low-quefrency part (liftering) is the cepstral method, and the same low-quefrency envelope, warped onto the mel scale, becomes the MFCCs used as classical speech features, which are derived on the speech page.

Linear prediction and the source-filter model

The source-filter model of speech (Fant's theory, realized digitally by Atal and Hanauer 1971 and reviewed by Makhoul 1975) says the vocal tract is an all-pole resonator excited by a source: a periodic pulse train for voiced sounds, noise for unvoiced. Linear predictive coding estimates that resonator by predicting each sample from a linear combination of its predecessors,

and choosing the coefficients \( a_k \) to minimize the total squared prediction error \( E = \sum_n e[n]^2 \). Setting \( \partial E / \partial a_j = 0 \) for each \( j \):

Writing the autocorrelation \( r[m] = \sum_n x[n] x[n-m] = r[-m] \), this is the set of normal equations

a \( p\times p \) linear system whose matrix is symmetric Toeplitz (constant along diagonals), because entry \( (j,k) \) depends only on \( |j-k| \). That structure is what makes LPC cheap: the Levinson-Durbin recursion solves a Toeplitz system in \( O(p^2) \) instead of the \( O(p^3) \) of a general solve, and it does so while guaranteeing the resulting filter is stable. Fitting an order-2 model to a signal generated by the AR process \( x[n] = 1.3\,x[n-1] - 0.6\,x[n-2] + \varepsilon[n] \) recovered coefficients \( (1.297, -0.598) \) against the true \( (1.3, -0.6) \), with the normal equations satisfied to a residual of \( 3.6\times 10^{-12} \) (measured, and the Levinson-Durbin solve of the Toeplitz system agrees with a direct general solve). The vocal-tract resonances, the formants, are then the pole angles of the estimated all-pole filter \( 1/A(z) \). On a synthetic two-formant vowel with resonances placed at \( 700 \) and \( 1220 \) Hz, an order-4 LPC recovered formants at \( 700.3 \) and \( 1216.6 \) Hz (measured), an error under \( 0.3\% \). LPC is the backbone of low-bit-rate speech codecs (the vocal-tract filter plus a coded excitation, as in CELP) precisely because a dozen coefficients capture the spectral envelope that carries phonetic identity.

Physical modeling: the Karplus-Strong string

Instead of adding sinusoids (additive synthesis) or shaping a spectrum (subtractive synthesis), physical modeling simulates the object that makes the sound. The simplest and most famous example is the Karplus-Strong plucked string (1983), which sounds startlingly like a real string for how little it computes. Fill a delay line of length \( L \) with random noise, the pluck, then recirculate it through a one-tap averaging lowpass filter:

Two things happen. The delay line of length \( L \) is a comb filter whose resonances fall at the harmonics of \( f_0 \approx f_s / (L + \tfrac{1}{2}) \), the extra half-sample coming from the group delay of the averaging filter at low frequencies, so the noise burst is immediately shaped into a harmonic tone at that pitch. And the averaging filter, being a lowpass, attenuates high frequencies more each pass around the loop, so the bright noisy attack decays into a mellow sustain exactly as a plucked string does, high partials dying first. Julius Smith later showed that this is the crudest member of the digital waveguide family: a real string is a pair of delay lines carrying left- and right-going traveling waves, and Karplus-Strong is that model collapsed to one loop with a lumped loss filter (Smith 1992). The same delay-line-plus-filter template, with the filter chosen to match the physics, models wind instruments, bowed strings, and drum membranes.

The pitch quantization is the catch. \( L \) must be an integer, so the achievable fundamentals are \( f_s/(L + \tfrac{1}{2}) \) for integer \( L \), and the spacing is coarse at high pitch. At \( f_s = 44.1 \) kHz, \( A_4 = 440 \) Hz wants \( L = f_s/440 - 0.5 = 99.7 \), rounded to \( L = 100 \), which yields \( f_s/(100.5) = 438.8 \) Hz, noticeably flat. Synthesizing a \( 220 \) Hz string uses \( L = \operatorname{round}(44100/220 - 0.5) = 200 \), a nominal fundamental of \( f_s/(L+0.5) = 219.95 \) Hz; recovering the period by autocorrelation of the decaying tone returned a peak at lag \( 199 \), i.e. \( 221.6 \) Hz (measured), the small discrepancy being the imperfection of lag-based pitch detection on a comb that is losing its high partials. Getting the tuning right therefore requires a fractional delay: an allpass or Lagrange interpolation filter inserted in the loop supplies the sub-sample part of \( L \) without disturbing the loop gain, which is exactly the Jaffe-Smith (1983) extension that made Karplus-Strong usable for real music.

Reverberation and perceptual coding

Schroeder reverberators and convolution reverb

A room turns a single click into a dense, decaying cloud of reflections. Schroeder's insight (1961, 1962) was that this cloud can be built from two elementary recirculating filters. A comb filter, \( y[n] = x[n] + g\, y[n-L] \), produces an infinite train of echoes spaced \( L \) samples apart, each \( g \) times the last, so it decays exponentially; the trouble is that its magnitude response is a comb, peaks at the harmonics of \( f_s/L \), which colors the sound metallically. Running several combs of mutually prime lengths in parallel fills in the spectral gaps and thickens the echo density. The comb's reverberation time, the \( -60 \) dB decay \( T_{60} \), follows from the per-loop attenuation: \( |g|^{T_{60} f_s / L} = 10^{-3} \) gives \( T_{60} = \tfrac{-3 L}{f_s \log_{10}|g|} \). A comb with \( L = 1500 \) (\( 34 \) ms) and \( g = 0.7 \) gives \( T_{60} = 0.659 \) s (measured). The second element is the allpass filter,

whose magnitude is exactly one at every frequency, \( |H_{\mathrm{ap}}(e^{j\omega})| = 1 \) (measured: the response is flat to within \( 10^{-4} \)), so it disperses the signal in time, smearing echoes into a smooth tail, without coloring it. Schroeder's design cascades a few allpasses after the parallel combs: the combs supply the exponential decay and echo density, the allpasses make it colorless. This topology is still the skeleton of algorithmic reverb, refined into feedback delay networks (Jot). The alternative, convolution reverb, measures the actual impulse response of a real hall by playing a swept sine and deconvolving, then convolves any dry signal with that response. It is exact for the measured room, and it is only affordable because of the FFT: a \( 3 \)-second impulse response at \( 48 \) kHz is \( 144{,}000 \) taps, hopeless by direct convolution but routine by partitioned overlap-add FFT convolution, the audio payoff of the \( O(N\log N) \) transform derived above.

Perceptual coding: throwing away the inaudible

Lossless coding of CD audio saves maybe a factor of two. MP3, AAC, and their successors reach ten to one or more by discarding sound the listener cannot hear, and the science that says what is inaudible is psychoacoustics (Zwicker and Fastl). Two facts do the work. The cochlea analyzes sound in critical bands, roughly 24 frequency groups (the Bark scale) each about a third of an octave wide, and within a band the ear integrates energy rather than resolving fine structure. And a loud tone raises the hearing threshold for nearby quieter sounds, simultaneous masking in frequency and temporal masking just before and after in time, so a soft component sitting under a loud one's masking skirt is simply not perceived. A perceptual coder computes, for each short frame, a masking threshold as a function of frequency from a psychoacoustic model: it finds the tonal and noise maskers, spreads each across its critical band, and combines them with the absolute threshold of hearing into a per-band noise floor that will be inaudible. It then quantizes each band coarsely enough that the quantization noise, whose power is set by the number of bits as \( \sigma_e^2 = \Delta^2/12 \) from the SNR derivation above, stays just below that masking threshold. Bits go where the ear can hear them and are starved from bands that are masked. The transform, the quantizer, and the arithmetic are ordinary DSP; the intelligence is the masking model deciding how many bits each critical band deserves. The same principle, keep only what a perceptual loss says is audible, reappears in modern neural codecs trained against perceptual and multi-resolution STFT objectives.

Implementation

The first block is the analysis pipeline in NumPy and SciPy: it computes a windowed spectrum, quantifies leakage, and measures the frequency resolution, reproducing the numbers quoted above. Everything here runs on CPU and is deterministic.

import numpy as np
from scipy import signal

fs = 1000.0            # sample rate, Hz
N  = 64               # DFT length
print("frequency resolution df = fs/N =", fs / N, "Hz/bin")   # 15.625

# on-bin vs off-bin leakage -----------------------------------------------
n = np.arange(N)
for f in (125.0, 128.0):                 # 125 Hz lands on bin 8 exactly
    x = np.cos(2 * np.pi * f * n / fs)   # rectangular (implicit) window
    X = np.abs(np.fft.rfft(x))
    top = np.sort(X)[-2:]
    print(f"f={f:5.1f} Hz  bin={f/(fs/N):5.2f}  peak/next = {top[1]/top[0]:.2f}")
# f=125.0 -> peak/next huge (energy in one bin);  f=128.0 -> 4.48 (leakage)

# window main-lobe / side-lobe tradeoff -----------------------------------
for name in ("boxcar", "hann", "hamming", "blackman"):
    w  = signal.get_window(name, 1024, fftbins=True)     # analysis window
    W  = np.abs(np.fft.rfft(w, 8 * 1024)); W /= W.max()
    dB = 20 * np.log10(W + 1e-12)
    i  = 1
    while dB[i + 1] < dB[i]:                              # descend the main lobe
        i += 1
    print(f"{name:8s} peak side lobe = {dB[i:].max():6.1f} dB")
# boxcar -13.4 | hann -31.5 | hamming -42.7 | blackman -58.2  dB

# convolution theorem: linear convolution via the FFT --------------------
x = np.array([1., 2., 3., 4.]); h = np.array([1., -1., 2.])
Nc  = len(x) + len(h) - 1                     # zero-pad to avoid circular wrap
fast = np.real(np.fft.ifft(np.fft.fft(x, Nc) * np.fft.fft(h, Nc)))
print("direct :", np.convolve(x, h))          # [1 1 3 5 2 8]
print("via FFT:", np.round(fast, 6))          # [1 1 3 5 2 8]

The second block designs both filter families and reads off their poles and stability, then confirms the FIR filter is linear phase and the IIR filter's poles lie inside the unit circle.

import numpy as np
from scipy import signal

fs, fc, M = 16000.0, 2000.0, 64          # sample rate, cutoff, FIR order

# FIR: windowed-sinc lowpass, symmetric -> exactly linear phase -----------
n = np.arange(M + 1)
h = np.sinc(2 * fc / fs * (n - M / 2)) * (2 * fc / fs)
h *= signal.get_window("hamming", M + 1, fftbins=False)   # symmetric window
h /= h.sum()                                              # unity DC gain
print("FIR symmetric (linear phase)?",
      np.max(np.abs(h - h[::-1])) < 1e-12)                # True; group delay = M/2 = 32
w, H = signal.freqz(h, worN=8000, fs=fs)
mag  = 20 * np.log10(np.abs(H) + 1e-12)
print("FIR -6 dB point:", round(w[np.argmin(np.abs(mag + 6))], 1), "Hz")   # 2000.0
print("FIR stopband  :", round(mag[w > 3000].max(), 1), "dB")             # -60.0

# IIR: 4th-order Butterworth via bilinear transform ----------------------
b, a   = signal.butter(4, fc, fs=fs)
_, p, _ = signal.tf2zpk(b, a)
print("IIR pole radii:", np.round(np.abs(p), 4), "stable:", np.all(np.abs(p) < 1))
# [0.7577 0.7577 0.4579 0.4579]  stable: True

# bilinear frequency warping: digital fc -> prewarped analog frequency ----
wd = 2 * np.pi * fc / fs
print("prewarped analog cutoff:", round(2 * fs * np.tan(wd / 2) / (2 * np.pi), 2), "Hz")  # 2109.57

# pole/zero stability of hand-written recursions -------------------------
for coef in ([1, -0.9], [1, -1.2], [1, -2*0.95*np.cos(np.pi/4), 0.95**2]):
    poles = np.roots(coef)
    print("poles", np.round(poles, 3), "stable:", np.all(np.abs(poles) < 1))
# 0.9 stable | 1.2 UNSTABLE | 0.672+-0.672j (r=0.95) stable

The third block does the synthesis and analysis half: LPC by the autocorrelation method solved through the Levinson-Durbin recursion, formant extraction from the pole angles, and a Karplus-Strong string whose pitch is measured back out by autocorrelation.

import numpy as np
from scipy import signal
from scipy.linalg import solve_toeplitz

# LPC by the autocorrelation method (normal equations, Toeplitz) ---------
rng = np.random.default_rng(1)
a_true = [1, -1.3, 0.6]                                   # AR(2): x = 1.3 x[-1] - 0.6 x[-2] + e
x = signal.lfilter([1], a_true, rng.standard_normal(4000))
p = 2
r = np.array([np.sum(x[:len(x) - k] * x[k:]) for k in range(p + 1)])   # autocorrelation
a = solve_toeplitz(r[:p], r[1:p + 1])                     # Levinson-Durbin, O(p^2)
print("LPC coeffs:", np.round(a, 4), "  (true 1.3, -0.6)")   # [1.297 -0.5978]

# formants = pole angles of the all-pole vocal-tract filter 1/A(z) -------
def resonator(f, bw):                                     # one formant, given fs=8000
    rr, th = np.exp(-np.pi * bw / 8000), 2 * np.pi * f / 8000
    return [1, -2 * rr * np.cos(th), rr * rr]
A_syn = np.convolve(resonator(700, 80), resonator(1220, 90))
src   = (np.arange(2000) % 80 == 0).astype(float)         # 100 Hz glottal pulse train
v     = signal.lfilter([1], A_syn, src)
rr    = np.array([np.sum(v[:len(v) - k] * v[k:]) for k in range(5)])
acoef = solve_toeplitz(rr[:4], rr[1:5])
roots = np.roots(np.concatenate([[1], -acoef]))
fmts  = np.sort(np.angle(roots[np.imag(roots) >= 0]) * 8000 / (2 * np.pi))
print("estimated formants:", np.round(fmts[fmts > 50], 1), "Hz")   # [700.3 1216.6]

# Karplus-Strong plucked string, pitch measured by autocorrelation ------
def karplus_strong(f0, fs=44100, dur=0.5, seed=0):
    L   = int(round(fs / f0 - 0.5))                       # loop length
    buf = list(np.random.default_rng(seed).uniform(-1, 1, L))
    out = []
    for _ in range(int(fs * dur)):
        first = buf.pop(0)
        buf.append(0.5 * (first + buf[0]))                # averaging lowpass in the loop
        out.append(first)
    return np.asarray(out)

y   = karplus_strong(220.0)
ac  = np.correlate(y, y, "full")[len(y) - 1:]
lag = np.argmax(ac[50:400]) + 50
print("KS 220 Hz string -> measured", round(44100 / lag, 2), "Hz")   # 221.61 (lag 199)

The last block is a differentiable STFT and a multi-resolution spectral loss, the objective that trains neural vocoders and DDSP-style synthesizers. Because the STFT is a linear map followed by a magnitude, it backpropagates: a gradient on the spectrogram flows back to the waveform. PyTorch and JAX express the same computation; the batched forward is identical, only the autodiff plumbing differs.

import torch
import torch.nn.functional as F

def stft_mag(x, n_fft, hop):                        # x: (B, T) waveform
    win = torch.hann_window(n_fft, device=x.device, dtype=x.dtype)
    X = torch.stft(x, n_fft=n_fft, hop_length=hop, window=win,
                   center=True, return_complex=True)   # (B, F, frames) complex
    return X.abs()                                      # magnitude spectrogram

def multires_stft_loss(pred, target, ffts=(512, 1024, 2048)):
    # sum of L1 magnitude errors at several resolutions -> a perceptual proxy
    loss = pred.new_zeros(())
    for n_fft in ffts:
        hop = n_fft // 4
        P = stft_mag(pred,   n_fft, hop)
        T = stft_mag(target, n_fft, hop)
        loss = loss + F.l1_loss(torch.log(P + 1e-5), torch.log(T + 1e-5))
    return loss

# the spectrogram is differentiable: optimize a waveform to match a target
target = torch.sin(2 * torch.pi * 440 * torch.arange(16000) / 16000)[None]  # (1, 16000)
x = torch.randn(1, 16000, requires_grad=True)
opt = torch.optim.Adam([x], lr=1e-2)
for step in range(200):
    opt.zero_grad()
    loss = multires_stft_loss(x, target)
    loss.backward()                                    # gradient flows through the STFT
    opt.step()
print("final multi-resolution STFT loss:", float(loss))

import jax, jax.numpy as jnp
from jax import lax

def stft_mag(x, n_fft, hop):                            # x: (T,) waveform
    win    = jnp.hanning(n_fft)
    starts = jnp.arange(0, x.shape[0] - n_fft + 1, hop)
    frames = jax.vmap(lambda s: lax.dynamic_slice(x, (s,), (n_fft,)))(starts)
    return jnp.abs(jnp.fft.rfft(frames * win, axis=-1))     # (frames, F) magnitude

def multires_stft_loss(pred, target, ffts=(512, 1024, 2048)):
    loss = 0.0
    for n_fft in ffts:
        P = stft_mag(pred,   n_fft, n_fft // 4)
        T = stft_mag(target, n_fft, n_fft // 4)
        loss = loss + jnp.mean(jnp.abs(jnp.log(P + 1e-5) - jnp.log(T + 1e-5)))
    return loss

target = jnp.sin(2 * jnp.pi * 440 * jnp.arange(16000) / 16000)
loss_and_grad = jax.value_and_grad(lambda x: multires_stft_loss(x, target))

x = jax.random.normal(jax.random.PRNGKey(0), (16000,))
lr = 1e-2
for step in range(200):
    loss, g = loss_and_grad(x)                          # grad through the STFT
    x = x - lr * g
print("final multi-resolution STFT loss:", float(loss))

How it is done in practice

Production audio DSP diverges from the textbook in the direction of latency, throughput, and numerical care. Real-time systems cannot afford a full-signal FFT; they run block-based processing with overlap-add or overlap-save so that a long convolution (a reverb impulse response, an FIR EQ) is broken into partitions whose per-block cost stays bounded and whose first output appears within one block of latency. Uniformly partitioned convolution keeps latency low for the early part of an impulse response and switches to larger, cheaper partitions for the diffuse tail, the standard architecture of a convolution-reverb plugin. Sample rates in music production are \( 44.1 \) or \( 48 \) kHz, oversampled internally to \( 2\times \) or \( 4\times \) around nonlinear stages (saturation, clipping) because a nonlinearity generates harmonics that would otherwise alias; the oversampler is itself a steep FIR or polyphase filter.

Fixed-point and low-precision arithmetic reintroduce quantization inside the algorithm, not just at the converter. Recursive IIR filters accumulate rounding error in their feedback path, and a filter that is stable in exact arithmetic can limit-cycle or go unstable in 16-bit fixed point when a pole coefficient is rounded; the mitigation is to factor high-order filters into cascaded second-order sections (biquads), each with well-conditioned coefficients, which is why every parametric EQ is a stack of biquads rather than one big polynomial. On the analysis side, the spectrograms feeding neural models are computed once and cached: a \( 16 \) kHz utterance framed at \( 25 \) ms windows and \( 10 \) ms hops yields about \( 100 \) frames per second, each a short FFT, and libraries like torchaudio and librosa compute these on GPU in batches so the front end is never the bottleneck. The tolerance for phase error differs sharply by task: a magnitude-only mel spectrogram discards phase entirely and relies on a neural vocoder or Griffin-Lim to invent a plausible one, whereas a time-stretch or a lossless codec must preserve phase to the sample, which is why the phase vocoder's unwrap step is where its engineering effort concentrates.

The current research frontier

The dominant recent theme is making classical DSP differentiable so it can be trained end to end. Differentiable Digital Signal Processing (DDSP, Engel and colleagues at Google Magenta, ICLR 2020) puts a harmonic-plus-noise synthesizer, a bank of sinusoidal oscillators and a time-varying filtered noise source, directly in the autograd graph, so a network predicts interpretable controls (pitch, loudness, filter coefficients) and the audio is rendered by the source-filter model this page derived. The payoff is data efficiency and control: the model cannot waste capacity relearning what a sinusoid is. Around it has grown a family of differentiable filters, differentiable all-pole and IIR layers (work from Queen Mary University of London and others on backpropagating through recursive filters), and differentiable room acoustics.

Neural audio codecs are the other major front. SoundStream (Google) and EnCodec (Meta AI) replace the hand-designed psychoacoustic quantizer with a learned encoder, a residual vector quantizer, and a decoder, trained against adversarial and multi-resolution STFT losses; the descendant discrete tokens (as in the DAC codec from Descript, and the audio tokens underlying music and speech language models) have become the interface between audio and transformer language models. The loss functions that make this work, multi-resolution STFT and mel losses, are collected in openly available libraries and are direct descendants of the perceptual-coding idea that error should be measured where the ear hears it. On the source-separation and enhancement side, time-domain models (Conv-TasNet and successors from groups at Columbia, Mitsubishi Electric Research Labs, and elsewhere) argued that a learned analysis filterbank can beat the fixed STFT for separation, while a competing line insists the STFT's interpretability and its exact invertibility are worth keeping; both camps agree the window and hop are now hyperparameters to be tuned, not constants. Physical modeling has its own neural revival, with differentiable waveguide and modal models learning instrument parameters from recordings, closing the loop back to Karplus-Strong.

Open source to read

• librosa/librosa: the reference Python library for music and audio analysis. Open librosa/core/spectrum.py to see the STFT, ISTFT, and Griffin-Lim exactly as derived here, and librosa/core/pitch.py for autocorrelation and YIN pitch tracking.
• scipy/scipy: scipy/signal is the canonical implementation of filter design (firwin, butter, bilinear), windows, and freqz. Reading signal/_filter_design.py shows the bilinear transform and pre-warping in production code.
• pytorch/audio: torchaudio. Start with torchaudio/functional/functional.py for the differentiable STFT, spectrogram, and Griffin-Lim that feed neural front ends on GPU.
• magenta/ddsp: differentiable oscillators and filtered noise. Open ddsp/synths.py and ddsp/core.py to see the source-filter model built as differentiable layers.
• csteinmetz1/auraloss: audio-focused loss functions in PyTorch. auraloss/freq.py is the multi-resolution STFT and mel-spectrogram losses used to train vocoders and codecs.
• LCAV/pyroomacoustics: room simulation and beamforming from EPFL. pyroomacoustics/room.py implements the image-source method for generating room impulse responses to convolve with dry audio.
• spatialaudio/python-sounddevice: real-time audio I/O binding to PortAudio, the practical way to run any of these filters on a live stream and hear the result.

Common misconceptions

"The Nyquist frequency is the highest frequency you can record." It is the highest you can record unambiguously. Content above \( f_s/2 \) is still recorded, it just aliases down onto a lower frequency and is then indistinguishable from a genuine tone there. The anti-alias filter's job is to remove that content in the analog domain before it aliases, because after sampling nothing can undo it.

"Zero-padding the FFT increases frequency resolution." Zero-padding interpolates the spectrum onto a finer grid but does not narrow the main lobe, which is set by the true window length. Two tones closer than the main-lobe width stay merged no matter how much you zero-pad; only a longer analysis window (more real samples) resolves them, and that costs time resolution.

"More bits always means better sound." Each bit adds \( 6.02 \) dB of headroom against quantization noise, but past the point where quantization noise drops below the analog noise floor (around \( 20 \)-\( 21 \) effective bits for the best converters), extra bits capture only thermal noise. Twenty-four bit converters are limited by their analog electronics, not their bit depth.

"FIR filters are always better than IIR because they are stable." FIR filters are unconditionally stable and can be exactly linear phase, but an IIR filter reaches the same selectivity with a fraction of the coefficients and a fraction of the delay. The right choice depends on whether phase linearity or computational cost dominates; both are used constantly.

"The spectrogram contains everything about the signal." The complex STFT is invertible and does contain everything, but the magnitude spectrogram most systems use throws away the phase. Reconstructing audio from magnitude alone requires guessing a consistent phase (Griffin-Lim, or a neural vocoder), and the guess is imperfect, which is a large part of why vocoding is hard.

"MP3 works by finding redundancy and compressing it like a zip file." Lossless entropy coding is only the last, minor stage. The compression comes from a psychoacoustic model deciding which spectral components are masked and therefore inaudible, then quantizing them coarsely or dropping them. It is lossy by design: the decoded signal differs from the original everywhere, just not audibly.

"The phase vocoder shifts pitch by shifting the FFT bins." Naively translating bins changes pitch but also corrupts the phase relationships and produces artifacts. Correct pitch-shifting estimates each partial's instantaneous frequency by phase unwrapping, time-stretches with phase-coherent resynthesis, and then resamples; the phase bookkeeping is the entire difficulty.

Self-check

References

1. Oppenheim, A. V., & Schafer, R. W. (2010). Discrete-Time Signal Processing (3rd ed.). Pearson. The standard graduate text for the DFT, z-transform, filter design, and the DTFT window analysis used throughout.
2. Proakis, J. G., & Manolakis, D. G. (2007). Digital Signal Processing: Principles, Algorithms, and Applications (4th ed.). Pearson. Companion reference with detailed FIR/IIR design and multirate material.
3. Smith, J. O. (2011). Spectral Audio Signal Processing. W3K Publishing / CCRMA. ccrma.stanford.edu/~jos/sasp. Free online textbook on the STFT, windows, and spectral modeling.
4. Smith, J. O. (2010). Physical Audio Signal Processing. W3K Publishing / CCRMA. ccrma.stanford.edu/~jos/pasp. Free online textbook on digital waveguides, Karplus-Strong, and reverberation.
5. Nyquist, H. (1928). Certain topics in telegraph transmission theory. Transactions of the AIEE, 47(2), 617-644.
6. Shannon, C. E. (1949). Communication in the presence of noise. Proceedings of the IRE, 37(1), 10-21. The sampling theorem.
7. Cooley, J. W., & Tukey, J. W. (1965). An algorithm for the machine calculation of complex Fourier series. Mathematics of Computation, 19(90), 297-301. doi:10.1090/S0025-5718-1965-0178586-1.
8. Harris, F. J. (1978). On the use of windows for harmonic analysis with the discrete Fourier transform. Proceedings of the IEEE, 66(1), 51-83. doi:10.1109/PROC.1978.10837. The definitive window survey.
9. Karplus, K., & Strong, A. (1983). Digital synthesis of plucked-string and drum timbres. Computer Music Journal, 7(2), 43-55. doi:10.2307/3680062.
10. Jaffe, D. A., & Smith, J. O. (1983). Extensions of the Karplus-Strong plucked-string algorithm. Computer Music Journal, 7(2), 56-69. doi:10.2307/3680063.
11. Smith, J. O. (1992). Physical modeling using digital waveguides. Computer Music Journal, 16(4), 74-91. doi:10.2307/3680470.
12. Schroeder, M. R. (1962). Natural sounding artificial reverberation. Journal of the Audio Engineering Society, 10(3), 219-223.
13. Makhoul, J. (1975). Linear prediction: A tutorial review. Proceedings of the IEEE, 63(4), 561-580. doi:10.1109/PROC.1975.9792.
14. Atal, B. S., & Hanauer, S. L. (1971). Speech analysis and synthesis by linear prediction of the speech wave. Journal of the Acoustical Society of America, 50(2B), 637-655. doi:10.1121/1.1912679.
15. Noll, A. M. (1967). Cepstrum pitch determination. Journal of the Acoustical Society of America, 41(2), 293-309. doi:10.1121/1.1910339.
16. Flanagan, J. L., & Golden, R. M. (1966). Phase vocoder. Bell System Technical Journal, 45(9), 1493-1509. doi:10.1002/j.1538-7305.1966.tb01706.x.
17. Portnoff, M. R. (1976). Implementation of the digital phase vocoder using the fast Fourier transform. IEEE Transactions on Acoustics, Speech, and Signal Processing, 24(3), 243-248. doi:10.1109/TASSP.1976.1162810.
18. Laroche, J., & Dolson, M. (1999). Improved phase vocoder time-scale modification of audio. IEEE Transactions on Speech and Audio Processing, 7(3), 323-332. doi:10.1109/89.759041.
19. Flanagan, J. L. (1972). Speech Analysis, Synthesis and Perception (2nd ed.). Springer. The classic source-filter reference.
20. Zwicker, E., & Fastl, H. (2007). Psychoacoustics: Facts and Models (3rd ed.). Springer. Critical bands, masking, and the Bark scale behind perceptual coding.
21. Painter, T., & Spanias, A. (2000). Perceptual coding of digital audio. Proceedings of the IEEE, 88(4), 451-515. doi:10.1109/5.842996.
22. Välimäki, V., Parker, J. D., Savioja, L., Smith, J. O., & Abel, J. S. (2012). Fifty years of artificial reverberation. IEEE Transactions on Audio, Speech, and Language Processing, 20(5), 1421-1448. doi:10.1109/TASL.2012.2189567.
23. Engel, J., Hantrakul, L., Gu, C., & Roberts, A. (2020). DDSP: Differentiable Digital Signal Processing. ICLR. arXiv:2001.04643 (Google Magenta).
24. McFee, B., Raffel, C., Liang, D., Ellis, D. P. W., McVicar, M., Battenberg, E., & Nieto, O. (2015). librosa: Audio and music signal analysis in Python. Proceedings of the 14th Python in Science Conference, 18-25. doi:10.25080/Majora-7b98e3ed-003.
25. Scheibler, R., Bezzam, E., & Dokmanić, I. (2018). Pyroomacoustics: A Python package for audio room simulation and array processing algorithms. ICASSP, 351-355. arXiv:1710.04196 (EPFL LCAV).