refactor: make silk_LPC_fit safe

Author: DCNick3

src/silk/LPC_fit.rs

@@ -1,84 +1,59 @@
 use crate::silk::bwexpander_32::silk_bwexpander_32;
+use crate::silk::SigProc_FIX::{silk_RSHIFT_ROUND, silk_SAT16, SILK_FIX_CONST};
+use ndarray::azip;
 
-pub mod typedef_h {
-    pub const silk_int16_MIN: i32 = i16::MIN as i32;
-    pub const silk_int16_MAX: i32 = i16::MAX as i32;
-}
-
-pub use self::typedef_h::{silk_int16_MAX, silk_int16_MIN};
+/// Convert int32 coefficients to int16 coefs and make sure there's no wrap-around
+///
+/// ```text
+/// a_QOUT   O     Output signal
+/// a_QIN    I/O   Input signal
+/// QOUT     I     Input Q domain
+/// QIN      I     Input Q domain
+/// d        I     Filter order
+/// ```
+pub fn silk_LPC_fit(a_QOUT: &mut [i16], a_QIN: &mut [i32], QOUT: i32, QIN: i32) {
+    let d = a_QOUT.len();
+    assert_eq!(a_QIN.len(), d);
 
-pub unsafe fn silk_LPC_fit(a_QOUT: *mut i16, a_QIN: *mut i32, QOUT: i32, QIN: i32, d: i32) {
-    let mut i: i32 = 0;
-    let mut k: i32 = 0;
-    let mut idx: i32 = 0;
-    let mut maxabs: i32 = 0;
-    let mut absval: i32 = 0;
-    let mut chirp_Q16: i32 = 0;
-    i = 0;
+    /* Limit the maximum absolute value of the prediction coefficients, so that they'll fit in int16 */
+    let mut i = 0;
     while i < 10 {
-        maxabs = 0;
-        k = 0;
+        /* Find maximum absolute value and its index */
+        let mut maxabs = 0;
+        let mut idx = 0;
+        let mut k = 0;
         while k < d {
-            absval = if *a_QIN.offset(k as isize) > 0 {
-                *a_QIN.offset(k as isize)
-            } else {
-                -*a_QIN.offset(k as isize)
-            };
+            let absval = a_QIN[k].abs();
             if absval > maxabs {
                 maxabs = absval;
                 idx = k;
             }
             k += 1;
         }
-        maxabs = if QIN - QOUT == 1 {
-            (maxabs >> 1) + (maxabs & 1)
+        maxabs = silk_RSHIFT_ROUND(maxabs, QIN - QOUT);
+
+        if maxabs > i16::MAX as i32 {
+            /* Reduce magnitude of prediction coefficients */
+            maxabs = std::cmp::min(maxabs, 163838); /* ( silk_int32_MAX >> 14 ) + silk_int16_MAX = 163838 */
+            let chirp_Q16 = SILK_FIX_CONST(0.999f64, 16)
+                - ((maxabs - i16::MAX as i32) << 14) / ((maxabs * (idx as i32 + 1)) >> 2);
+            silk_bwexpander_32(a_QIN, chirp_Q16);
         } else {
-            (maxabs >> QIN - QOUT - 1) + 1 >> 1
-        };
-        if !(maxabs > silk_int16_MAX) {
             break;
         }
-        maxabs = if maxabs < 163838 { maxabs } else { 163838 };
-        chirp_Q16 = (0.999f64 * ((1) << 16) as f64 + 0.5f64) as i32
-            - (((maxabs - 0x7fff) as u32) << 14) as i32 / (maxabs * (idx + 1) >> 2);
-        silk_bwexpander_32(std::slice::from_raw_parts_mut(a_QIN, d as usize), chirp_Q16);
+
         i += 1;
     }
+
     if i == 10 {
-        k = 0;
-        while k < d {
-            *a_QOUT.offset(k as isize) = (if (if QIN - QOUT == 1 {
-                (*a_QIN.offset(k as isize) >> 1) + (*a_QIN.offset(k as isize) & 1)
-            } else {
-                (*a_QIN.offset(k as isize) >> QIN - QOUT - 1) + 1 >> 1
-            }) > silk_int16_MAX
-            {
-                silk_int16_MAX
-            } else if (if QIN - QOUT == 1 {
-                (*a_QIN.offset(k as isize) >> 1) + (*a_QIN.offset(k as isize) & 1)
-            } else {
-                (*a_QIN.offset(k as isize) >> QIN - QOUT - 1) + 1 >> 1
-            }) < silk_int16_MIN
-            {
-                silk_int16_MIN
-            } else if QIN - QOUT == 1 {
-                (*a_QIN.offset(k as isize) >> 1) + (*a_QIN.offset(k as isize) & 1)
-            } else {
-                (*a_QIN.offset(k as isize) >> QIN - QOUT - 1) + 1 >> 1
-            }) as i16;
-            *a_QIN.offset(k as isize) =
-                ((*a_QOUT.offset(k as isize) as i32 as u32) << QIN - QOUT) as i32;
-            k += 1;
-        }
+        /* Reached the last iteration, clip the coefficients */
+        azip!((out in a_QOUT, input in a_QIN) {
+            *out = silk_SAT16(silk_RSHIFT_ROUND(*input, QIN - QOUT)) as i16;
+            *input = (*out as i32) << (QIN - QOUT);
+        });
     } else {
-        k = 0;
-        while k < d {
-            *a_QOUT.offset(k as isize) = (if QIN - QOUT == 1 {
-                (*a_QIN.offset(k as isize) >> 1) + (*a_QIN.offset(k as isize) & 1)
-            } else {
-                (*a_QIN.offset(k as isize) >> QIN - QOUT - 1) + 1 >> 1
-            }) as i16;
-            k += 1;
-        }
+        azip!((out in a_QOUT, &mut input in a_QIN) {
+            *out = silk_RSHIFT_ROUND(input, QIN - QOUT) as i16;
+        });
     };
 }

src/silk/NLSF2A.rs

@@ -84,7 +84,12 @@ pub unsafe fn silk_NLSF2A(a_Q12: *mut i16, NLSF: *const i16, d: i32, _arch: i32)
         a32_QA1[(d - k - 1) as usize] = Qtmp - Ptmp;
         k += 1;
     }
-    silk_LPC_fit(a_Q12, a32_QA1.as_mut_ptr(), 12, QA + 1, d);
+    silk_LPC_fit(
+        std::slice::from_raw_parts_mut(a_Q12, d as usize),
+        &mut a32_QA1[..d as usize],
+        12,
+        QA + 1,
+    );
     i = 0;
     while silk_LPC_inverse_pred_gain_c(a_Q12, d) == 0 && i < MAX_LPC_STABILIZE_ITERATIONS {
         silk_bwexpander_32(&mut a32_QA1[..d as usize], 65536 - ((2) << i) as i32);

src/silk/SigProc_FIX.rs

@@ -1,6 +1,17 @@
 /// max order of the LPC analysis in schur() and k2a()
 pub const SILK_MAX_ORDER_LPC: usize = 24;
 
+#[inline]
+pub fn silk_SAT16(a: i32) -> i32 {
+    if a > i16::MAX as i32 {
+        i16::MAX as i32
+    } else if a < i16::MIN as i32 {
+        i16::MIN as i32
+    } else {
+        a
+    }
+}
+
 #[inline]
 pub fn silk_ROR32(a32: i32, rot: i32) -> i32 {
     let x: u32 = a32 as u32;
@@ -16,17 +27,27 @@ pub fn silk_ROR32(a32: i32, rot: i32) -> i32 {
 }
 
 #[inline]
-pub fn silk_max_16(a: i16, b: i16) -> i16 {
-    a.max(b)
+pub fn silk_RSHIFT_ROUND(a: i32, shift: i32) -> i32 {
+    if shift == 1 {
+        ((a) >> 1) + ((a) & 1)
+    } else {
+        (((a) >> ((shift) - 1)) + 1) >> 1
+    }
 }
 
 #[inline]
-pub fn silk_min_32(a: i32, b: i32) -> i32 {
-    a.min(b)
+pub fn silk_RSHIFT_ROUND64(a: i64, shift: i32) -> i64 {
+    if shift == 1 {
+        ((a) >> 1) + ((a) & 1)
+    } else {
+        (((a) >> ((shift) - 1)) + 1) >> 1
+    }
 }
+
+/// Convert floating-point constants to fixed-point
 #[inline]
-pub fn silk_max_32(a: i32, b: i32) -> i32 {
-    a.max(b)
+pub fn SILK_FIX_CONST(C: f64, Q: i32) -> i32 {
+    (C * (1i64 << Q) as f64 + 0.5f64) as i32
 }
 
 #[inline]
@@ -39,32 +60,17 @@ pub fn silk_max_int(a: i32, b: i32) -> i32 {
 }
 
 #[inline]
-pub fn silk_RSHIFT_ROUND(a: i32, shift: i32) -> i32 {
-    if shift == 1 {
-        ((a) >> 1) + ((a) & 1)
-    } else {
-        (((a) >> ((shift) - 1)) + 1) >> 1
-    }
+pub fn silk_max_16(a: i16, b: i16) -> i16 {
+    a.max(b)
 }
 
 #[inline]
-pub fn silk_RSHIFT_ROUND64(a: i64, shift: i32) -> i64 {
-    if shift == 1 {
-        ((a) >> 1) + ((a) & 1)
-    } else {
-        (((a) >> ((shift) - 1)) + 1) >> 1
-    }
+pub fn silk_min_32(a: i32, b: i32) -> i32 {
+    a.min(b)
 }
-
 #[inline]
-pub fn silk_SAT16(a: i32) -> i32 {
-    if a > i16::MAX as i32 {
-        i16::MAX as i32
-    } else if a < i16::MIN as i32 {
-        i16::MIN as i32
-    } else {
-        a
-    }
+pub fn silk_max_32(a: i32, b: i32) -> i32 {
+    a.max(b)
 }
 
 #[inline]