]>
Commit | Line | Data |
---|---|---|
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 | } |