Various optimizations for LC3

src/lc3.c

@@ -483,8 +483,10 @@ static void store_s16(
     float *xs = decoder->x + decoder->xs_off;
     int ns = lc3_ns(dt, sr);
 
+    /* The `copysignf` form is branchless, and lets the compiler vectorize
+     * the unit stride case, while rounding identically to `x +/- 0.5`. */
     for ( ; ns > 0; ns--, xs++, pcm += stride) {
-        int32_t s = *xs >= 0 ? (int)(*xs + 0.5f) : (int)(*xs - 0.5f);
+        const int32_t s = (int32_t)(*xs + copysignf(0.5f, *xs));
         *pcm = LC3_SAT16(s);
     }
 }
@@ -506,8 +508,7 @@ static void store_s24(
     int ns = lc3_ns(dt, sr);
 
     for ( ; ns > 0; ns--, xs++, pcm += stride) {
-        int32_t s = *xs >= 0 ? (int32_t)(lc3_ldexpf(*xs, 8) + 0.5f)
-                             : (int32_t)(lc3_ldexpf(*xs, 8) - 0.5f);
+        const int32_t s = (int32_t)(lc3_ldexpf(*xs, 8) + copysignf(0.5f, *xs));
         *pcm = LC3_SAT24(s);
     }
 }
@@ -529,10 +530,8 @@ static void store_s24_3le(
     int ns = lc3_ns(dt, sr);
 
     for ( ; ns > 0; ns--, xs++, pcm += 3*stride) {
-        int32_t s = *xs >= 0 ? (int32_t)(lc3_ldexpf(*xs, 8) + 0.5f)
-                             : (int32_t)(lc3_ldexpf(*xs, 8) - 0.5f);
+        const int32_t s = LC3_SAT24((int32_t)(lc3_ldexpf(*xs, 8) + copysignf(0.5f, *xs)));
 
-        s = LC3_SAT24(s);
         pcm[0] = (s >>  0) & 0xff;
         pcm[1] = (s >>  8) & 0xff;
         pcm[2] = (s >> 16) & 0xff;

src/ltpf.c

@@ -597,11 +597,11 @@ static bool detect_pitch(struct lc3_ltpf_analysis *ltpf,
     const int16_t *x1 = x - (r0 + t1);
     const int16_t *x2 = x - (r0 + t2);
 
-    float nc1 = rm1 <= 0 ? 0 :
-        rm1 / sqrtf(dot(x, x, n) * dot(x1, x1, n));
+    const float e = dot(x, x, n);
 
-    float nc2 = rm2 <= 0 ? 0 :
-        rm2 / sqrtf(dot(x, x, n) * dot(x2, x2, n));
+    const float nc1 = rm1 <= 0 ? 0 : rm1 / sqrtf(e * dot(x1, x1, n));
+
+    const float nc2 = rm2 <= 0 ? 0 : rm2 / sqrtf(e * dot(x2, x2, n));
 
     int t1sel = nc2 <= 0.85f * nc1;
     ltpf->tc = (t1sel ? t1 : t2);

src/mdct.c

@@ -93,7 +93,8 @@ LC3_HOT static inline void fft_5(
 #ifndef fft_bf3
 LC3_HOT static inline void fft_bf3(
     const struct lc3_fft_bf3_twiddles *twiddles,
-    const struct lc3_complex *x, struct lc3_complex *y, int n)
+    const struct lc3_complex *x, struct lc3_complex *y,
+    int n)
 {
     int n3 = twiddles->n3;
     const struct lc3_complex (*w0)[2] = twiddles->t;
@@ -135,7 +136,8 @@ LC3_HOT static inline void fft_bf3(
 #ifndef fft_bf2
 LC3_HOT static inline void fft_bf2(
     const struct lc3_fft_bf2_twiddles *twiddles,
-    const struct lc3_complex *x, struct lc3_complex *y, int n)
+    const struct lc3_complex *x, struct lc3_complex *y,
+    int n)
 {
     int n2 = twiddles->n2;
     const struct lc3_complex *w = twiddles->t;
@@ -147,11 +149,20 @@ LC3_HOT static inline void fft_bf2(
 
         for (int j = 0; j < n2; j++, x0++, x1++) {
 
-            y0[j].re = x0->re + x1->re * w[j].re - x1->im * w[j].im;
-            y0[j].im = x0->im + x1->im * w[j].re + x1->re * w[j].im;
+            /* Load before storing, and compute the rotated value `x1 * w` once for both outputs. */
+
+            const float x0r = x0->re, x0i = x0->im;
+            const float x1r = x1->re, x1i = x1->im;
+            const float wr = w[j].re, wi = w[j].im;
+
+            const float tr = x1r * wr - x1i * wi;
+            const float ti = x1i * wr + x1r * wi;
+
+            y0[j].re = x0r + tr;
+            y0[j].im = x0i + ti;
 
-            y1[j].re = x0->re - x1->re * w[j].re + x1->im * w[j].im;
-            y1[j].im = x0->im - x1->im * w[j].re - x1->re * w[j].im;
+            y1[j].re = x0r - tr;
+            y1[j].im = x0i - ti;
         }
     }
 }

src/spec.c

@@ -246,25 +246,33 @@ LC3_HOT static int adjust_gain(
  * dt, sr          Duration and samplerate of the frame
  * g_int           Quantization gain value
  * x               Spectral coefficients, scaled as output
+ * xq              Return the quantized coefficients
  * n               Return count of significants
  */
 LC3_HOT static void quantize(
-    enum lc3_dt dt, enum lc3_srate sr, int g_int, float *x, int *n)
+    enum lc3_dt dt, enum lc3_srate sr, int g_int, float *x, int32_t *xq, int *n)
 {
     float g_inv = unquantize_gain(-g_int);
     int ne = lc3_ne(dt, sr);
 
-    *n = ne;
+    /* The quantized magnitude is derived once here, instead of being
+     * recomputed from the floating point values on each coding pass.
+     * Note that `|x| >= xq_min` is equivalent to `xq != 0`. */
 
-    for (int i = 0; i < ne; i += 2) {
-        float xq_min = lc3_hr(sr) ? 0.5f : 10.f/16;
+    const float xq_off = lc3_hr(sr) ? 0.5f : 6.f/16;
 
-        x[i+0] *= g_inv;
-        x[i+1] *= g_inv;
+    for (int i = 0; i < ne; i++) {
+        const float v = x[i] * g_inv;
+        const int32_t a = (int32_t)(fabsf(v) + xq_off);
 
-        *n = fabsf(x[i+0]) >= xq_min ||
-             fabsf(x[i+1]) >= xq_min   ? ne : *n - 2;
+        x[i] = v;
+        xq[i] = v < 0 ? -a : a;
     }
+
+    int i;
+    for (i = ne; i >= 2 && !xq[i-1] && !xq[i-2]; i -= 2);
+
+    *n = i;
 }
 
 /**
@@ -314,15 +322,15 @@ static bool resolve_modes(enum lc3_srate sr, int nbytes, bool *p_lsb_mode)
 /**
  * Bit consumption
  * dt, sr, nbytes  Duration, samplerate and size of the frame
- * x               Spectral quantized coefficients
+ * xq              Spectral quantized coefficients
  * n               Count of significant coefficients, updated on truncation
  * nbits_budget    Truncate to stay in budget, when not zero
  * p_lsb_mode      Return True when LSB's are not AC coded, or NULL
  * return          The number of bits coding the spectrum
  */
 LC3_HOT static int compute_nbits(
     enum lc3_dt dt, enum lc3_srate sr, int nbytes,
-    const float *x, int *n, int nbits_budget, bool *p_lsb_mode)
+    const int32_t *xq, int *n, int nbits_budget, bool *p_lsb_mode)
 {
     bool lsb_mode, high_rate = resolve_modes(sr, nbytes, &lsb_mode);
     int ne = lc3_ne(dt, sr);
@@ -343,9 +351,8 @@ LC3_HOT static int compute_nbits(
         for ( ; i < LC3_MIN(*n, (ne + 2) >> (1 - h))
                 && nbits <= nbits_budget; i += 2) {
 
-            float xq_off = lc3_hr(sr) ? 0.5f : 6.f/16;
-            uint32_t a = fabsf(x[i+0]) + xq_off;
-            uint32_t b = fabsf(x[i+1]) + xq_off;
+            uint32_t a = (uint32_t)LC3_ABS(xq[i+0]);
+            uint32_t b = (uint32_t)LC3_ABS(xq[i+1]);
 
             const uint8_t *lut = lut_coeff[state];
 
@@ -825,6 +832,12 @@ void lc3_spec_analyze(
 {
     bool reset_off;
 
+    /* Quantized spectrum, derived once and shared by the two bit
+     * consumption passes below. Scoped to the encoder path, so that
+     * the decoder stack is not affected. */
+
+    int32_t xq[LC3_MAX_NE];
+
     /* --- Bit budget --- */
 
     const int nbits_gain = 8;
@@ -847,9 +860,9 @@ void lc3_spec_analyze(
 
     /* --- Quantization --- */
 
-    quantize(dt, sr, g_int, x, &side->nq);
+    quantize(dt, sr, g_int, x, xq, &side->nq);
 
-    int nbits = compute_nbits(dt, sr, nbytes, x, &side->nq, 0, NULL);
+    const int nbits = compute_nbits(dt, sr, nbytes, xq, &side->nq, 0, NULL);
 
     spec->nbits_off = reset_off ? 0 : nbits_off;
     spec->nbits_spare = reset_off ? 0 : nbits_budget - nbits;
@@ -860,11 +873,10 @@ void lc3_spec_analyze(
         g_off + g_int, nbits, nbits_budget, g_off + g_min);
 
     if (g_adj)
-        quantize(dt, sr, g_adj, x, &side->nq);
+        quantize(dt, sr, g_adj, x, xq, &side->nq);
 
     side->g_idx = g_int + g_adj + g_off;
-    nbits = compute_nbits(dt, sr, nbytes,
-        x, &side->nq, nbits_budget, &side->lsb_mode);
+    (void)compute_nbits(dt, sr, nbytes, xq, &side->nq, nbits_budget, &side->lsb_mode);
 }
 
 /**

test/spec_py.c

@@ -88,7 +88,9 @@ static PyObject *quantize_py(PyObject *m, PyObject *args)
 
     CTYPES_CHECK("x", x_obj = to_1d_ptr(x_obj, NPY_FLOAT, ne, &x));
 
-    quantize(dt, sr, g_int, x, &nq);
+    int32_t xq[LC3_MAX_NE];
+
+    quantize(dt, sr, g_int, x, xq, &nq);
 
     return Py_BuildValue("Oi", x_obj, nq);
 }
@@ -113,8 +115,19 @@ static PyObject *compute_nbits_py(PyObject *m, PyObject *args)
 
     bool lsb_mode;
 
-    int nbits = compute_nbits(
-        dt, sr, nbytes, x, &nq, nbits_budget, &lsb_mode);
+    /* `compute_nbits()` now takes the quantized spectrum, which the
+     * encoder derives in `quantize()`. Rebuild it here, to keep the
+     * Python side interface on floating point values. */
+
+    int32_t xq[LC3_MAX_NE];
+    const float xq_off = lc3_hr(sr) ? 0.5f : 6.f/16;
+
+    for (int i = 0; i < ne; i++) {
+        const int32_t a = (int32_t)(fabsf(x[i]) + xq_off);
+        xq[i] = x[i] < 0 ? -a : a;
+    }
+
+    const int nbits = compute_nbits(dt, sr, nbytes, xq, &nq, nbits_budget, &lsb_mode);
 
     return Py_BuildValue("iii", nbits, nq, lsb_mode);
 }