nihav-codec-support/src/dsp/lpc.rs

   1 //! LPC filter calculation.
   2
   3 /// Calculated LPC filter for the provided input data, filter order is defined by `filter` length.
   4 pub fn calc_lpc_filter<T: Copy>(data: &[T], filter: &mut [f64])
   5 where f64: From<T>
   6 {
   7     let order = filter.len();
   8
   9     let mut acorr = vec![0.0f64; order + 1];
  10
  11     for (i, sum) in acorr.iter_mut().enumerate() {
  12         let src1 = &data[i..];
  13         *sum = data.iter().zip(src1.iter()).fold(0.0, |acc, (&a, &b)| acc + f64::from(a) * f64::from(b));
  14     }
  15
  16     for el in filter.iter_mut() {
  17         *el = 0.0;
  18     }
  19
  20     if acorr[0].abs() < 1.0e-6 {
  21         return;
  22     }
  23
  24     let mut bwd0 = vec![0.0; order];
  25     let mut bwd1 = vec![0.0; order];
  26
  27     bwd0[0] = 1.0 / acorr[0];
  28     filter[order - 1] = acorr[1] / acorr[0];
  29
  30     for i in 2..=order {
  31         let eb = acorr[1..i].iter().zip(bwd0.iter()).fold(0.0, |acc, (&a, &b)| acc + a * b);
  32
  33         let scale = 1.0 / (1.0 - eb * eb);
  34         for j in 0..i {
  35             bwd1[j] = if j > 0 { bwd0[j - 1] * scale } else { 0.0 };
  36             if j < i - 1 {
  37                 bwd1[j] -= bwd0[i - 2 - j] * eb * scale;
  38             }
  39         }
  40
  41         let ex = acorr[1..i].iter().rev().zip(filter.iter().rev()).fold(0.0, |acc, (&a, &b)| acc + a * b);
  42         let scale = acorr[i] - ex;
  43         for (filt, &add) in filter.iter_mut().rev().zip(bwd1.iter()).take(i) {
  44             *filt += add * scale;
  45         }
  46
  47         std::mem::swap(&mut bwd0, &mut bwd1);
  48     }
  49 }