linear_predictive_coding/
lib.rs

1//! # lpc-rs
2//!
3//! `lpc-rs` is a library for calculating Linear Predictive Coding (LPC) coefficients.
4//! It provides three methods to calculate LPC coefficients.
5//! - Low speed method (Temporarily commented out due to suspension of updates to dependent libraries)
6//! - High speed method
7//! - Burg method
8
9use ndarray::prelude::*;
10// use ndarray_inverse::Inverse;
11
12/// find the correlation of the input array
13/// # Arguments
14/// [;N]
15///
16/// # Returns
17/// [;N]
18pub fn correlate(a: ArrayView1<f64>) -> Array1<f64> {
19    a.iter()
20        .enumerate()
21        .map(|(n, _)| {
22            a.slice(s![n..])
23                .iter()
24                .zip(a.iter())
25                .map(|(x, y)| x * y)
26                .sum::<f64>()
27        })
28        .collect()
29}
30
31// pub fn calc_lpc_by_low_speed(a: ArrayView1<f64>, depth: usize) -> Array1<f64> {
32//     let r = correlate(a);
33
34//     let mut large_r: Array2<f64> = Array2::zeros((depth, depth));
35//     for i in 0..depth {
36//         for j in 0..depth {
37//             large_r[[i, j]] = r[(i as isize - j as isize).abs() as usize];
38//         }
39//     }
40
41//     println!("{:?}", large_r);
42
43//     let r = r.slice(s![1..=depth]);
44
45//     println!("{:?}", r);
46
47//     let inverse_large_r = large_r.inv().unwrap();
48
49//     let a = inverse_large_r.dot(&r);
50
51//     let minus_a = a.mapv(|x| -x);
52
53//     minus_a
54// }
55
56/// https://qiita.com/hirokisince1998/items/fd50c0515c7788458fce
57/// Levinson-Durbin recursion
58pub fn calc_lpc_by_levinson_durbin(a: ArrayView1<f64>, depth: usize) -> Option<(Array1<f64>, f64)> {
59    if a.len() < depth {
60        return None;
61    }
62
63    let r = correlate(a);
64    // println!("{:?}", r);
65    let r = r.slice(s![..=depth]);
66    // println!("{:?}", r);
67
68    fn calc_lpc_by_high_speed_inner(
69        a: ArrayView1<f64>,
70        depth: usize,
71        r: ArrayView1<f64>,
72    ) -> (Array1<f64>, f64) {
73        if depth == 1 {
74            let a = Array1::from_iter(vec![1.0, -r[1] as f64 / r[0]]);
75            let e = a.dot(&r.slice(s![..2]));
76
77            // println!("{:?}", a);
78
79            (a, e)
80        } else {
81            let (aa, ee) = calc_lpc_by_high_speed_inner(a, depth - 1, r);
82
83            let kk = -aa.dot(&r.slice(s![1..=depth; -1])) / ee;
84
85            let large_u = ndarray::concatenate![Axis(0), aa.view(), Array1::from_elem(1, 0.0)];
86
87            let large_v = large_u.slice(s![..; -1]);
88
89            let a = large_u.clone() + large_v.mapv(|x| x * kk);
90
91            let e = ee * (1.0 - kk * kk);
92
93            (a, e)
94        }
95    }
96
97    let (a, e) = calc_lpc_by_high_speed_inner(a, depth, r.view());
98
99    Some((a.slice_move(s![1..]), e))
100}
101
102/// Burg method
103pub fn calc_lpc_by_burg(x: ArrayView1<f64>, depth: usize) -> Option<Array1<f64>> {
104    if x.len() < depth {
105        return None;
106    }
107
108    let mut a = Array1::<f64>::zeros(depth + 1);
109
110    let mut k = Array1::<f64>::zeros(depth);
111
112    a[0] = 1.0;
113
114    let mut f = x.to_owned();
115
116    let mut b = x.to_owned();
117
118    let n = x.len();
119
120    for p in 0..depth {
121        let kf = f.slice(s![p + 1..]);
122        let kb = b.slice(s![..n - p - 1]);
123        // element-wise sum of squares
124        let d = kf.iter().map(|x| x * x).sum::<f64>() + kb.iter().map(|x| x * x).sum::<f64>();
125        k[p] = -2.0 * kf.iter().zip(kb.iter()).map(|(x, y)| x * y).sum::<f64>() / d;
126        let u = a.slice(s![..=p + 1]);
127        let v = u.slice(s![..; -1]);
128
129        // println!("u: {:?}", u);
130        // println!("v: {:?}", v);
131
132        let added = &u + &v.mapv(|x| x * k[p]);
133        a.slice_mut(s![..=p + 1]).assign(&added.view());
134        let fu = b.slice(s![..n - p - 1]).mapv(|x| x * k[p]);
135        let bu = f.slice(s![p + 1..]).mapv(|x| x * k[p]);
136        f.slice_mut(s![p + 1..])
137            .iter_mut()
138            .zip(fu.iter())
139            .for_each(|(x, fu)| *x += *fu);
140        b.slice_mut(s![..n - p - 1])
141            .iter_mut()
142            .zip(bu.iter())
143            .for_each(|(x, bu)| *x += bu);
144    }
145
146    Some(a.slice_move(s![1..]))
147}
148
149#[cfg(test)]
150mod tests {
151    use super::*;
152
153    #[test]
154    fn test_correlate() {
155        let a = Array1::from(vec![2., 3., -1., -2., 1., 4., 1.]);
156        let expected = Array1::from(vec![36.0, 11.0, -16.0, -7.0, 13.0, 11.0, 2.0]);
157        assert_eq!(correlate(a.view()), expected);
158    }
159
160    // #[test]
161    // fn test_calc_lpc_by_low_speed() {
162    //     let a = Array1::from(vec![2., 3., -1., -2., 1., 4., 1.]);
163    //     let depth = 3;
164    //     let expected = Array1::from(vec![
165    //         -0.6919053749597682,
166    //         0.7615062761506275,
167    //         -0.34575152880592214,
168    //     ]);
169    //     assert_eq!(calc_lpc_by_low_speed(a.view(), depth), expected);
170    // }
171
172    #[test]
173    fn test_calc_lpc_by_high_speed() {
174        let a = Array1::from(vec![2., 3., -1., -2., 1., 4., 1.]);
175        let depth = 3;
176        let expected = Array1::from(vec![
177            -0.6919053749597684,
178            0.7615062761506278,
179            -0.3457515288059223,
180        ]);
181        assert_eq!(
182            calc_lpc_by_levinson_durbin(a.view(), depth).unwrap().0,
183            expected
184        );
185    }
186
187    #[test]
188    fn test_calc_lpc_by_burg() {
189        let a = Array1::from(vec![2., 3., -1., -2., 1., 4., 1.]);
190        let depth = 3;
191        let expected = Array1::from(vec![
192            -1.0650404360323664,
193            1.157238171254371,
194            -0.5771692748969812,
195        ]);
196        assert_eq!(calc_lpc_by_burg(a.view(), depth), Some(expected));
197    }
198}
linear_predictive_coding/lib.rs

linear_predictive_coding/
lib.rs
