feat!: Remove unstable SIMD optimization (#99)

Author: cm-ayf

Cargo.toml

@@ -16,7 +16,6 @@ documentation = "https://docs.rs/jbonsai"
 default = ["htsvoice"]
 binary = ["htsvoice", "dep:hound"]
 htsvoice = ["dep:nom"]
-simd = []
 
 [package.metadata.docs.rs]
 all-features = true

README-ja.md

@@ -29,21 +29,6 @@ jbonsai = "0.3.0"
 
 <!-- x-release-please-end -->
 
-### SIMD (experimental)
-
-jbonsaiは、feature [portable_simd](https://github.com/rust-lang/portable-simd)による高速化をサポートしています。SIMD高速化を有効にするには、
-
-- nightly ツールチェーンを使用する必要があります。
-- 次のように`features = ["simd"]`を指定する必要があります。
-  <!-- x-release-please-start-version -->
-  ```toml
-  [dependencies]
-  jbonsai = { version = "0.3.0", features = ["simd"] }
-  ```
-  <!-- x-release-please-end -->
-
-SIMDサポートは非常に実験的であり、いつでも変更される可能性があります。
-
 ## 使用例
 
 以下の例は，「盆栽」と読み上げる音声を生成し，`speech`変数にモノラル, 48000 HzのPCMとして格納します．

README.md

@@ -29,21 +29,6 @@ jbonsai = "0.3.0"
 
 <!-- x-release-please-end -->
 
-### SIMD (experimental)
-
-jbonsai supports acceleration provided by feature [portable_simd](https://github.com/rust-lang/portable-simd). In order to enable SIMD acceleration,
-
-- you must use nightly toolchain.
-- you have to specify `features = ["simd"]` as follows:
-  <!-- x-release-please-start-version -->
-  ```toml
-  [dependencies]
-  jbonsai = { version = "0.3.0", features = ["simd"] }
-  ```
-  <!-- x-release-please-end -->
-
-The SIMD support is highly experimental and may change at any time.
-
 ## Example
 
 This example produces a mono, 48,000 Hz (typically) PCM data saying 「盆栽」(ぼんさい; bonsai) in `speech` variable.

src/lib.rs

@@ -1,6 +1,5 @@
 #![deny(missing_docs)]
 #![cfg_attr(docsrs, feature(doc_cfg, doc_auto_cfg))]
-#![cfg_attr(feature = "simd", feature(portable_simd))]
 #![doc = include_str!("../README.md")]
 
 mod constants;

src/vocoder/mlsa.rs

@@ -0,0 +1,119 @@
+use super::coefficients::Coefficients;
+
+/// N == pd + 1
+#[derive(Debug, Clone)]
+pub struct MelLogSpectrumApproximation<const N: usize> {
+    d11: [f64; N],
+    d12: [f64; N],
+    d21: [Vec<f64>; N],
+    d22: [f64; N],
+}
+
+pub trait Pade<const N: usize> {
+    const PPADE: [f64; N];
+}
+macro_rules! impl_pade {
+    ($($i:literal: $ppade:expr),* $(,)?) => {
+        $(
+            impl Pade<$i> for MelLogSpectrumApproximation<$i> {
+                const PPADE: [f64; $i] = $ppade;
+            }
+        )*
+    };
+}
+
+impl_pade!(
+    1: [1.00000000000f64],
+    2: [1.00000000000f64, 0.00000000000f64],
+    3: [1.00000000000f64, 0.00000000000f64, 0.00000000000f64],
+    4: [1.00000000000f64, 0.00000000000f64, 0.00000000000f64, 0.00000000000f64],
+    5: [1.00000000000f64, 0.49992730000f64, 0.10670050000f64, 0.01170221000f64, 0.00056562790f64],
+    6: [1.00000000000f64, 0.49993910000f64, 0.11070980000f64, 0.01369984000f64, 0.00095648530f64, 0.00003041721f64],
+);
+
+impl<const N: usize> MelLogSpectrumApproximation<N>
+where
+    Self: Pade<N>,
+{
+    pub fn new(nmcp: usize) -> Self {
+        Self {
+            d11: [0.0; N],
+            d12: [0.0; N],
+            d21: std::array::from_fn(|_| vec![0.0; nmcp]),
+            d22: [0.0; N],
+        }
+    }
+
+    #[inline(always)]
+    pub fn df(&mut self, x: &mut f64, alpha: f64, coefficients: &'_ Coefficients) {
+        self.df1(x, alpha, coefficients);
+        self.df2(x, alpha, coefficients);
+    }
+
+    #[inline(always)]
+    fn df1(&mut self, x: &mut f64, alpha: f64, coefficients: &'_ Coefficients) {
+        let aa = 1.0 - alpha * alpha;
+        let mut out = 0.0;
+        for i in (1..N).rev() {
+            self.d11[i] = aa * self.d12[i - 1] + alpha * self.d11[i];
+            self.d12[i] = self.d11[i] * coefficients[1];
+            let v = self.d12[i] * Self::PPADE[i];
+            *x += if i & 1 != 0 { v } else { -v };
+            out += v;
+        }
+        self.d12[0] = *x;
+        *x += out;
+    }
+
+    #[inline(always)]
+    fn df2(&mut self, x: &mut f64, alpha: f64, coefficients: &'_ Coefficients) {
+        let mut out = 0.0;
+        for i in (1..N).rev() {
+            self.d22[i] = fir(&mut self.d21[i - 1], self.d22[i - 1], alpha, coefficients);
+            let v = self.d22[i] * Self::PPADE[i];
+            *x += if i & 1 != 0 { v } else { -v };
+            out += v;
+        }
+        self.d22[0] = *x;
+        *x += out;
+    }
+}
+
+fn fir(d: &mut [f64], x: f64, alpha: f64, coefficients: &[f64]) -> f64 {
+    d[0] = x;
+
+    let a = alpha;
+    let aa = a * a;
+    let aaaa = aa * aa;
+    let iaa = 1.0 - aa;
+    let mut rem = 0.0;
+
+    let mut chunks = d.chunks_exact_mut(4);
+    for chunk in chunks.by_ref() {
+        (chunk[0], chunk[1], chunk[2], chunk[3], rem) = (
+            iaa * rem + a * chunk[0],
+            iaa * (-a * rem + chunk[0]) + a * chunk[1],
+            iaa * (aa * rem + (-a * chunk[0] + chunk[1])) + a * chunk[2],
+            iaa * (aa * (-a * rem + chunk[0]) + (-a * chunk[1] + chunk[2])) + a * chunk[3],
+            aaaa * rem + (aa * (-a * chunk[0] + chunk[1]) + (-a * chunk[2] + chunk[3])),
+        );
+    }
+    for di in chunks.into_remainder() {
+        (*di, rem) = (iaa * rem + a * *di, -a * rem + *di);
+    }
+
+    let mut y = [0.0; 8];
+    let mut c = coefficients[2..d.len()].chunks_exact(8);
+    let mut d = d[2..d.len()].chunks_exact(8);
+
+    use std::iter::zip;
+    for (c, d) in zip(&mut c, &mut d) {
+        for (y, (c, d)) in zip(&mut y, zip(c, d)) {
+            *y += c * d;
+        }
+    }
+    for (y, (c, d)) in zip(&mut y, zip(c.remainder(), d.remainder())) {
+        *y += c * d;
+    }
+    ((y[0] + y[1]) + (y[2] + y[3])) + ((y[4] + y[5]) + (y[6] + y[7]))
+}