| 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 | } |