Linear predictive coding (LPC)

Deep dive

How it works

In: a frame of samples and an order p · Out: p predictor (or reflection) coefficients and the residual signal · State: none for the autocorrelation method; the decoder keeps the last p output samples · Cost: autocorrelation N(p+1) MACs, Levinson-Durbin O(p^2), filtering p MACs/sample

1. Window a frame (Window functions); 10-30 ms for speech, 4096 samples typical for FLAC.
2. Autocorrelation r[0..=p], r[l] = sum x[n] x[n+l].
3. Levinson-Durbin solves the Toeplitz system R a = -r order by order: k_i = -(r[i] + sum a[j] r[i-j]) / E_{i-1}, a_j += k_i a_{i-j}, a_i = k_i, E_i = E_{i-1} (1 - k_i^2). Lower orders come free.
4. Analysis filter A(z) = 1 + sum a[k] z^-k gives e[n] = x[n] + sum a[k] x[n-k]. Synthesis 1/A(z) inverts it.
5. Quantize and transmit: lossless coders use integer coefficients plus a shift and entropy-code e (Rice codes); lossy speech coders convert to reflection coefficients or line spectral frequencies (LSF), vector-quantize them, and code the excitation with a codebook.

x[n] --> [ A(z): FIR, p taps ] --> e[n] --> (entropy coder / quantizer) --> bits
                                                                             |
x^[n] <-- [ 1/A(z): all-pole ] <-- e^[n] <-- (decoder) <---------------------

Code

use std::f64::consts::PI;

/// Autocorrelation r[0..=p] of a frame (apply a window first).
pub fn autocorr(x: &[f64], p: usize) -> Vec<f64> {
    (0..=p).map(|lag| x[lag..].iter().zip(x).map(|(a, b)| a * b).sum()).collect()
}

pub fn hamming(x: &mut [f64]) {
    let n = x.len() as f64 - 1.0;
    for (i, v) in x.iter_mut().enumerate() {
        *v *= 0.54 - 0.46 * (2.0 * PI * i as f64 / n).cos();
    }
}

/// Levinson-Durbin on r[0..=p]. Returns (a, k, err): A(z) = 1 + a[1] z^-1 + ... (a[0] = 1),
/// reflection coefficients k[1..=p] and the final prediction-error energy.
pub fn levinson(r: &[f64]) -> (Vec<f64>, Vec<f64>, f64) {
    let p = r.len() - 1;
    let (mut a, mut k) = (vec![0.0; p + 1], vec![0.0; p + 1]);
    a[0] = 1.0;
    let mut err = r[0];
    for i in 1..=p {
        if err <= 0.0 { break; } // digital silence: keep the lower order
        let acc: f64 = r[i] + (1..i).map(|j| a[j] * r[i - j]).sum::<f64>();
        let ki = -acc / err;
        let prev = a.clone();
        for j in 1..i { a[j] = prev[j] + ki * prev[i - j]; }
        a[i] = ki;
        k[i] = ki;
        err *= 1.0 - ki * ki; // |ki| < 1 => err stays > 0 and 1/A(z) is stable
    }
    (a, k, err)
}

/// Bandwidth expansion: a[i] *= gamma^i pulls the poles toward the origin.
pub fn bw_expand(a: &mut [f64], gamma: f64) {
    let mut g = 1.0;
    for c in a.iter_mut() { *c *= g; g *= gamma; }
}

/// Whitening filter A(z): e[n] = x[n] + sum a[k] x[n-k] (history before n = 0 is zero).
pub fn residual(x: &[f64], a: &[f64]) -> Vec<f64> {
    (0..x.len()).map(|n| (0..a.len()).filter(|&k| k <= n).map(|k| a[k] * x[n - k]).sum()).collect()
}

/// Synthesis filter 1/A(z): x[n] = e[n] - sum a[k] x[n-k].
pub fn synthesize(e: &[f64], a: &[f64]) -> Vec<f64> {
    let mut x = vec![0.0; e.len()];
    for n in 0..e.len() {
        let fb: f64 = (1..a.len()).filter(|&k| k <= n).map(|k| a[k] * x[n - k]).sum();
        x[n] = e[n] - fb;
    }
    x
}

/// FLAC-style integer predictor: c[k-1] = round(-a[k] * 2^shift), |c| < 2^(prec-1).
pub fn quantize(a: &[f64], prec: u32) -> (Vec<i32>, u32) {
    let cmax = a[1..].iter().fold(1e-9f64, |m, v| m.max(v.abs()));
    let shift = ((prec - 1) as i32 - cmax.log2().floor() as i32 - 1).clamp(0, 31) as u32;
    let lim = (1i64 << (prec - 1)) - 1;
    let c = a[1..].iter().map(|v| ((-v * (1u64 << shift) as f64).round() as i64).clamp(-lim, lim) as i32).collect();
    (c, shift)
}

/// Bit-exact integer residual for n >= order (the first `order` samples are sent verbatim).
/// i64 accumulator: 24-bit samples x 15-bit coefficients x 32 taps overflows i32.
pub fn int_residual(x: &[i32], c: &[i32], shift: u32) -> Vec<i32> {
    (c.len()..x.len())
        .map(|n| x[n] - (c.iter().enumerate().map(|(k, &ck)| ck as i64 * x[n - 1 - k] as i64).sum::<i64>() >> shift) as i32)
        .collect()
}

pub fn int_restore(warmup: &[i32], res: &[i32], c: &[i32], shift: u32) -> Vec<i32> {
    let mut x = warmup.to_vec();
    for &e in res {
        let n = x.len();
        let p = c.iter().enumerate().map(|(k, &ck)| ck as i64 * x[n - 1 - k] as i64).sum::<i64>() >> shift;
        x.push(e + p as i32);
    }
    x
}

#[cfg(test)]
mod tests {
    use super::*;
    #[test]
    fn ar2_roundtrip() {
        let mut s = 12345u64; // LCG noise in [-0.5, 0.5)
        let e: Vec<f64> = (0..20000).map(|_| {
            s = s.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
            (s >> 33) as f64 / (1u64 << 31) as f64 - 0.5
        }).collect();
        let x = synthesize(&e, &[1.0, -1.6, 0.8]); // true A(z) = 1 - 1.6 z^-1 + 0.8 z^-2
        let (a, k, err) = levinson(&autocorr(&x, 4)); // fit order 4 > true order 2
        assert!((a[1] + 1.6).abs() < 0.05 && (a[2] - 0.8).abs() < 0.05);
        assert!(a[3].abs() < 0.05 && a[4].abs() < 0.05 && k[1..].iter().all(|v| v.abs() < 1.0));
        let er: f64 = residual(&x, &a).iter().map(|v| v * v).sum();
        assert!((err - er).abs() / er < 0.01); // Levinson error = residual energy
        let back = synthesize(&residual(&x, &a), &a);
        assert!(back.iter().zip(&x).all(|(p, q)| (p - q).abs() < 1e-6));
        let mut fr = x[..240].to_vec();
        hamming(&mut fr);
        let (mut a10, k10, _) = levinson(&autocorr(&fr, 10));
        assert!(k10[1..].iter().all(|v| v.abs() < 1.0));
        bw_expand(&mut a10, 0.95);
        // lossless integer path: 16-bit-ish signal, order 4, 12-bit coefficients
        let xi: Vec<i32> = x.iter().map(|v| (v * 2000.0).round() as i32).collect();
        let (c, sh) = quantize(&a, 12);
        let ri = int_residual(&xi, &c, sh);
        assert_eq!(int_restore(&xi[..4], &ri, &c, sh), xi);
        let mag = |v: &[i32]| v.iter().map(|&q| q.unsigned_abs() as f64).sum::<f64>() / v.len() as f64;
        assert!(mag(&ri) * 3.0 < mag(&xi)); // residual much smaller than signal
        println!("shift {sh}, mean |x| {:.0}, mean |res| {:.0}", mag(&xi), mag(&ri));
    }
}

quantize follows the FLAC convention (prediction = sum c[k] x[n-k] >> shift): integer-only, bit-exact decoding. The test covers AR(2) recovery, filter round trip, stability and a lossless round trip.

Complexity and constants

• Autocorrelation: N(p+1) multiply-adds, the dominant analysis cost. Measured here (M4 Pro, scalar f64, order of): 4096 samples, p=12, about 30 us. Levinson-Durbin p=12: about 0.3 us, so the solve is free.
• Per-sample filtering: p multiply-adds, i.e. 12 MAC/sample at p=12 (a 4096-sample residual took about 19 us in the naive loop above).
• Side info: FLAC sends p coefficients of 1-15 bits each per subframe (RFC 9639), e.g. 12 x 12 = 144 bits per 4096-sample block.
• Returns diminish quickly with order; FLAC's presets stop at 8-12.

Efficiency

• Cost: measured here (M4 Pro, scalar f64, order of): autocorrelation of 4096 samples at p=12 about 30 us; Levinson-Durbin p=12 about 0.3 us; filtering p MACs per sample (a 4096-sample residual about 19 us in the naive loop). FLAC sends 12 x 12 = 144 coefficient bits per 4096-sample block. Bound: the encoder is compute-bound in autocorrelation (N*p MACs); FLAC encodes and decodes without floating point (RFC 9639).
• Levers:
  • order search from one recursion — Levinson gives E_i for every i <= p; libFLAC's FLAC__lpc_compute_lp_coefficients returns all orders at once — free
  • fewer apodization windows — each extra window costs another autocorrelation, solve and residual pass; flac presets -6 to -8 use subdivide_tukey with 2-3 windows (flac man page) — slower presets only
  • SIMD over lags — libFLAC ships intrinsic variants; f32 only for short, well-conditioned frames — long frames
  • FLAC fixed predictors (orders 0-4) — no analysis at all — fastest presets
• Naive → standard → state of the art:
  • fixed predictors, orders 0-4 · no analysis, weaker prediction
  • windowed autocorrelation plus Levinson-Durbin, order 8-12 · about 30 us per 4096 samples here
  • several windows with order search, or neural vocoder on LPC features (LPCNet) · more analysis for fewer bits or better speech

Implementation notes

• Speech: update every 10-20 ms, interpolate in the LSF or reflection domain, keep filter memory across frames.
• Opus SILK uses a modified Burg analysis (no window needed, short subframes) and a two-stage vector quantizer for normalized LSFs; NLSF-to-LPC conversion applies up to 10 rounds of bandwidth expansion to fit Q12 coefficients (RFC 6716).
• Speech order rule of thumb: sample rate in kHz plus 2-4.
• White-noise correction (r[0] *= 1 + 1e-9) or a lag window cures ill-conditioning on tonal or band-limited input.
• Lossless: i64 accumulator for 24-bit audio; the first p samples go out verbatim (FLAC warm-up); FLAC lists coefficients newest-sample first.

Pitfalls

• No window: edge leakage distorts the estimate.
• r[0] = 0 (digital silence) divides by zero; guard it as in the code.
• Quantizing or interpolating a[k] directly can push poles outside the unit circle.
• Sign conventions differ (1 + sum here, 1 - sum in many textbooks, prediction sum c x in FLAC); mixing them gives a plausible but wrong filter.
• Float decoders can differ between platforms: lossless needs the integer path.
• High order on short frames overfits; coefficient bits eat the gain.

Variants and alternatives

In production

• Lossless audio: FLAC (order up to 32, Subset limits it to 12 up to 48 kHz) (RFC 9639); FFmpeg's FLAC and ALAC encoders (libavcodec/lpc.c).
• Speech coding and telephony: the Opus SILK layer, a linear-prediction layer for speech (RFC 6716, Opus); GSM and AMR speech codecs; Codec2.
• Speech synthesis: LPCNet feeds linear prediction into a neural vocoder (Sources).
• Statistics and time series: statsmodels offers the Levinson-Durbin recursion for autoregressive processes (docs).

Implementations

Sources

• John Makhoul, Linear prediction: A tutorial review, Proceedings of the IEEE 63(4), 1975. paper
• B. S. Atal, S. L. Hanauer, Speech Analysis and Synthesis by Linear Prediction of the Speech Wave, J. Acoust. Soc. Am. 50(2B), 1971. paper
• Norman Levinson, The Wiener RMS error criterion in filter design and prediction, J. Math. Phys. 25, 1947; J. Durbin, The fitting of time-series models, Rev. Int. Stat. Inst. 28, 1960.
• M. van Beurden, A. Weaver, Free Lossless Audio Codec (FLAC), RFC 9639, 2024. RFC · code
• J.-M. Valin et al., Definition of the Opus Audio Codec, RFC 6716, 2012. RFC · code
• J.-M. Valin, J. Skoglund, LPCNet: Improving Neural Speech Synthesis Through Linear Prediction, ICASSP 2019. paper