|
| 1 | +use super::coefficients::Coefficients; |
| 2 | + |
| 3 | +/// N == pd + 1 |
| 4 | +#[derive(Debug, Clone)] |
| 5 | +pub struct MelLogSpectrumApproximation<const N: usize> { |
| 6 | + d11: [f64; N], |
| 7 | + d12: [f64; N], |
| 8 | + d21: [Vec<f64>; N], |
| 9 | + d22: [f64; N], |
| 10 | +} |
| 11 | + |
| 12 | +pub trait Pade<const N: usize> { |
| 13 | + const PPADE: [f64; N]; |
| 14 | +} |
| 15 | +macro_rules! impl_pade { |
| 16 | + ($($i:literal: $ppade:expr),* $(,)?) => { |
| 17 | + $( |
| 18 | + impl Pade<$i> for MelLogSpectrumApproximation<$i> { |
| 19 | + const PPADE: [f64; $i] = $ppade; |
| 20 | + } |
| 21 | + )* |
| 22 | + }; |
| 23 | +} |
| 24 | + |
| 25 | +impl_pade!( |
| 26 | + 1: [1.00000000000f64], |
| 27 | + 2: [1.00000000000f64, 0.00000000000f64], |
| 28 | + 3: [1.00000000000f64, 0.00000000000f64, 0.00000000000f64], |
| 29 | + 4: [1.00000000000f64, 0.00000000000f64, 0.00000000000f64, 0.00000000000f64], |
| 30 | + 5: [1.00000000000f64, 0.49992730000f64, 0.10670050000f64, 0.01170221000f64, 0.00056562790f64], |
| 31 | + 6: [1.00000000000f64, 0.49993910000f64, 0.11070980000f64, 0.01369984000f64, 0.00095648530f64, 0.00003041721f64], |
| 32 | +); |
| 33 | + |
| 34 | +impl<const N: usize> MelLogSpectrumApproximation<N> |
| 35 | +where |
| 36 | + Self: Pade<N>, |
| 37 | +{ |
| 38 | + pub fn new(nmcp: usize) -> Self { |
| 39 | + Self { |
| 40 | + d11: [0.0; N], |
| 41 | + d12: [0.0; N], |
| 42 | + d21: std::array::from_fn(|_| vec![0.0; nmcp]), |
| 43 | + d22: [0.0; N], |
| 44 | + } |
| 45 | + } |
| 46 | + |
| 47 | + #[inline(always)] |
| 48 | + pub fn df(&mut self, x: &mut f64, alpha: f64, coefficients: &'_ Coefficients) { |
| 49 | + self.df1(x, alpha, coefficients); |
| 50 | + self.df2(x, alpha, coefficients); |
| 51 | + } |
| 52 | + |
| 53 | + #[inline(always)] |
| 54 | + fn df1(&mut self, x: &mut f64, alpha: f64, coefficients: &'_ Coefficients) { |
| 55 | + let aa = 1.0 - alpha * alpha; |
| 56 | + let mut out = 0.0; |
| 57 | + for i in (1..N).rev() { |
| 58 | + self.d11[i] = aa * self.d12[i - 1] + alpha * self.d11[i]; |
| 59 | + self.d12[i] = self.d11[i] * coefficients[1]; |
| 60 | + let v = self.d12[i] * Self::PPADE[i]; |
| 61 | + *x += if i & 1 != 0 { v } else { -v }; |
| 62 | + out += v; |
| 63 | + } |
| 64 | + self.d12[0] = *x; |
| 65 | + *x += out; |
| 66 | + } |
| 67 | + |
| 68 | + #[inline(always)] |
| 69 | + fn df2(&mut self, x: &mut f64, alpha: f64, coefficients: &'_ Coefficients) { |
| 70 | + let mut out = 0.0; |
| 71 | + for i in (1..N).rev() { |
| 72 | + self.d22[i] = fir(&mut self.d21[i - 1], self.d22[i - 1], alpha, coefficients); |
| 73 | + let v = self.d22[i] * Self::PPADE[i]; |
| 74 | + *x += if i & 1 != 0 { v } else { -v }; |
| 75 | + out += v; |
| 76 | + } |
| 77 | + self.d22[0] = *x; |
| 78 | + *x += out; |
| 79 | + } |
| 80 | +} |
| 81 | + |
| 82 | +#[inline(always)] |
| 83 | +fn fir(d: &mut [f64], x: f64, alpha: f64, coefficients: &[f64]) -> f64 { |
| 84 | + d[0] = x; |
| 85 | + |
| 86 | + let a = alpha; |
| 87 | + let iaa = 1.0 - a * a; |
| 88 | + let mut rem = 0.0; |
| 89 | + |
| 90 | + for di in &mut d[..] { |
| 91 | + (*di, rem) = (a * *di + rem, iaa * *di - a * rem); |
| 92 | + } |
| 93 | + |
| 94 | + let mut y = 0.0; |
| 95 | + let c = &coefficients[..d.len()]; |
| 96 | + for i in 2..d.len() { |
| 97 | + y += d[i] * c[i]; |
| 98 | + } |
| 99 | + y |
| 100 | +} |
0 commit comments