huff2/quantize: optimize scale-factor and codebook processing

libfaac/huff2.c

@@ -50,21 +50,19 @@ static hcode16_t * const hmap[12] = {
     book06, book07, book08, book09, book10, book11
 };
 
-/* Bitwise branchless non-zero check: returns 1 if x != 0, else 0. */
-static inline int is_nonzero(int x)
-{
-    return (int)(((unsigned int)x | (unsigned int)-x) >> 31);
-}
 
 /* Both books of a pair share the index expression and the sign-bit count; only
  * the table differs. One walk, two lookups.
  *
+ * ISO 14496-3 multidimensional Huffman section tuple indexing. Constant dimensions
+ * (DIM_S4, DIM_M4, DIM_S2, DIM_M2_7, DIM_M2_12) allow constant folding into shift-adds.
+ *
  * bnum is always a pair base -- huffbook takes HCB_ESC without sizing it -- so
  * there is deliberately no escape case. */
-static void huffcode_size_pair(int *qs, int len, int bnum, int *bits_a, int *bits_b)
+static void huffcode_size_pair(const int * __restrict qs, int len, int bnum, int *bits_a, int *bits_b)
 {
-    hcode16_t *booka = hmap[bnum];
-    hcode16_t *bookb = hmap[bnum + 1];
+    const hcode16_t *booka = hmap[bnum];
+    const hcode16_t *bookb = hmap[bnum + 1];
     int a = 0, b = 0;
     int i;
 
@@ -78,9 +76,9 @@ static void huffcode_size_pair(int *qs, int len, int bnum, int *bits_a, int *bit
         break;
     case HCB_3:
         for (i = 0; i < len; i += 4) {
-            int idx = DIM_M4*DIM_M4*DIM_M4 * abs(qs[i]) + DIM_M4*DIM_M4 * abs(qs[i+1]) + DIM_M4 * abs(qs[i+2]) + abs(qs[i+3]);
-            /* Branchless sign-bit counting using bitwise logic */
-            int sign = is_nonzero(qs[i]) + is_nonzero(qs[i+1]) + is_nonzero(qs[i+2]) + is_nonzero(qs[i+3]);
+            int a0 = abs(qs[i]), a1 = abs(qs[i+1]), a2 = abs(qs[i+2]), a3 = abs(qs[i+3]);
+            int idx = DIM_M4*DIM_M4*DIM_M4 * a0 + DIM_M4*DIM_M4 * a1 + DIM_M4 * a2 + a3;
+            int sign = (a0 != 0) + (a1 != 0) + (a2 != 0) + (a3 != 0);
             a += booka[idx].len + sign;
             b += bookb[idx].len + sign;
         }
@@ -94,16 +92,18 @@ static void huffcode_size_pair(int *qs, int len, int bnum, int *bits_a, int *bit
         break;
     case HCB_7:
         for (i = 0; i < len; i += 2) {
-            int idx = DIM_M2_7 * abs(qs[i]) + abs(qs[i+1]);
-            int sign = is_nonzero(qs[i]) + is_nonzero(qs[i+1]);
+            int a0 = abs(qs[i]), a1 = abs(qs[i+1]);
+            int idx = DIM_M2_7 * a0 + a1;
+            int sign = (a0 != 0) + (a1 != 0);
             a += booka[idx].len + sign;
             b += bookb[idx].len + sign;
         }
         break;
     case HCB_9:
         for (i = 0; i < len; i += 2) {
-            int idx = DIM_M2_12 * abs(qs[i]) + abs(qs[i+1]);
-            int sign = is_nonzero(qs[i]) + is_nonzero(qs[i+1]);
+            int a0 = abs(qs[i]), a1 = abs(qs[i+1]);
+            int idx = DIM_M2_12 * a0 + a1;
+            int sign = (a0 != 0) + (a1 != 0);
             a += booka[idx].len + sign;
             b += bookb[idx].len + sign;
         }
@@ -117,10 +117,10 @@ static void huffcode_size_pair(int *qs, int len, int bnum, int *bits_a, int *bit
 }
 
 /* Bitstream mutation function, called once per finalized frame. */
-static void huffcode_write(int *qs, int len, int bnum, CoderInfo *coder)
+static void huffcode_write(const int * __restrict qs, int len, int bnum, CoderInfo *coder)
 {
-    hcode16_t *book = hmap[bnum];
-    int i, j;
+    const hcode16_t *book = hmap[bnum];
+    int i;
     int datacnt = coder->datacnt;
 
     switch (bnum) {
@@ -135,15 +135,15 @@ static void huffcode_write(int *qs, int len, int bnum, CoderInfo *coder)
     case HCB_3:
     case HCB_4:
         for (i = 0; i < len; i += 4) {
-            int idx = DIM_M4*DIM_M4*DIM_M4 * abs(qs[i]) + DIM_M4*DIM_M4 * abs(qs[i+1]) + DIM_M4 * abs(qs[i+2]) + abs(qs[i+3]);
+            int q0 = qs[i], q1 = qs[i+1], q2 = qs[i+2], q3 = qs[i+3];
+            int a0 = abs(q0), a1 = abs(q1), a2 = abs(q2), a3 = abs(q3);
+            int idx = DIM_M4*DIM_M4*DIM_M4 * a0 + DIM_M4*DIM_M4 * a1 + DIM_M4 * a2 + a3;
             int blen = book[idx].len;
             int data = book[idx].data;
-            for (j = 0; j < 4; j++) {
-                if (qs[i+j]) {
-                    blen++;
-                    data = (data << 1) | (qs[i+j] < 0);
-                }
-            }
+            if (q0) { blen++; data = (data << 1) | (q0 < 0); }
+            if (q1) { blen++; data = (data << 1) | (q1 < 0); }
+            if (q2) { blen++; data = (data << 1) | (q2 < 0); }
+            if (q3) { blen++; data = (data << 1) | (q3 < 0); }
             coder->s[datacnt].data = data;
             coder->s[datacnt++].len = blen;
         }
@@ -159,31 +159,27 @@ static void huffcode_write(int *qs, int len, int bnum, CoderInfo *coder)
     case HCB_7:
     case HCB_8:
         for (i = 0; i < len; i += 2) {
-            int idx = DIM_M2_7 * abs(qs[i]) + abs(qs[i+1]);
+            int q0 = qs[i], q1 = qs[i+1];
+            int a0 = abs(q0), a1 = abs(q1);
+            int idx = DIM_M2_7 * a0 + a1;
             int blen = book[idx].len;
             int data = book[idx].data;
-            for (j = 0; j < 2; j++) {
-                if (qs[i+j]) {
-                    blen++;
-                    data = (data << 1) | (qs[i+j] < 0);
-                }
-            }
+            if (q0) { blen++; data = (data << 1) | (q0 < 0); }
+            if (q1) { blen++; data = (data << 1) | (q1 < 0); }
             coder->s[datacnt].data = data;
             coder->s[datacnt++].len = blen;
         }
         break;
     case HCB_9:
     case HCB_10:
         for (i = 0; i < len; i += 2) {
-            int idx = DIM_M2_12 * abs(qs[i]) + abs(qs[i+1]);
+            int q0 = qs[i], q1 = qs[i+1];
+            int a0 = abs(q0), a1 = abs(q1);
+            int idx = DIM_M2_12 * a0 + a1;
             int blen = book[idx].len;
             int data = book[idx].data;
-            for (j = 0; j < 2; j++) {
-                if (qs[i+j]) {
-                    blen++;
-                    data = (data << 1) | (qs[i+j] < 0);
-                }
-            }
+            if (q0) { blen++; data = (data << 1) | (q0 < 0); }
+            if (q1) { blen++; data = (data << 1) | (q1 < 0); }
             coder->s[datacnt].data = data;
             coder->s[datacnt++].len = blen;
         }
@@ -228,16 +224,10 @@ static void huffcode_write(int *qs, int len, int bnum, CoderInfo *coder)
 }
 
 /* Pick the codebook that minimizes the bit cost for a given band. */
-int huffbook(CoderInfo *coder, int *qs, int len)
+int huffbook(CoderInfo *coder, const int *qs, int len, int maxq)
 {
-    int i, maxq = 0;
     int bookmin = HCB_ZERO;
 
-    for (i = 0; i < len; i++) {
-        int q = abs(qs[i]);
-        if (maxq < q) maxq = q;
-    }
-
     if (maxq > 0) {
         /* Each spectral book covers values up to its LAV; select the range-pair
          * whose lower book just fits maxq, then pick the partner if it costs fewer

libfaac/huff2.h

@@ -94,7 +94,7 @@ static inline int clamp_sf_diff(int diff)
 /* Forward declaration for CoderInfo */
 struct CoderInfo;
 
-int huffbook(struct CoderInfo *coder, int *qs, int len);
+int huffbook(struct CoderInfo *coder, const int *qs, int len, int maxq);
 int writebooks(struct CoderInfo *coder, BitStream *stream, int writeFlag);
 int writesf(struct CoderInfo *coder, BitStream *bitStream, int writeFlag);
 

libfaac/quantize.c

@@ -22,34 +22,44 @@
 #include "huff2.h"
 #include "cpu_compute.h"
 
-typedef void (*QuantizeFunc)(const float * __restrict xr, int * __restrict xi, int n, float sfacfix);
+typedef int (*QuantizeFunc)(const float * __restrict xr, int * __restrict xi, int n, float sfacfix);
 
 #if defined(HAVE_SSE2)
-extern void quantize_sse2(const float * __restrict xr, int * __restrict xi, int n, float sfacfix);
+extern int quantize_sse2(const float * __restrict xr, int * __restrict xi, int n, float sfacfix);
 #endif
 
-static void quantize_scalar(const float * __restrict xr, int * __restrict xi, int n, float sfacfix)
+static int quantize_scalar(const float * __restrict xr, int * __restrict xi, int n, float sfacfix)
 {
     const float magic = MAGIC_NUMBER;
-    int i;
+    int i, maxq = 0;
     for (i = 0; i < n; i++)
     {
         float val = xr[i];
         float tmp = fabsf(val) * sfacfix;
         tmp = sqrtf(tmp * sqrtf(tmp));
         int q = (int)(tmp + magic);
+        if (q > maxq) maxq = q;
         xi[i] = (val < 0) ? -q : q;
     }
+    return maxq;
 }
 
 static QuantizeFunc qfunc = quantize_scalar;
 static float sfstep;
 static float max_quant_limit;
 
+#define GAIN_LUT_SIZE 512
+#define GAIN_LUT_BIAS 256
+/* ISO 14496-3 §8.3.4 scalefactors use quarter-dB steps (1 unit = 2^(1/4) in amplitude).
+ * Precomputed 2^(sfac/4) LUT eliminates repeated transcendental powf calls during gain coupling. */
+static float gain_lut[GAIN_LUT_SIZE];
+static float log10_width_sf_lut[128];
+
 #define SF_CHAIN_UNSET INT_MIN
 
 void QuantizeInit(void)
 {
+    int i;
 #if defined(HAVE_SSE2)
     CPUCaps caps = get_cpu_caps();
     if (caps & CPU_CAP_SSE2)
@@ -59,20 +69,36 @@ void QuantizeInit(void)
         qfunc = quantize_scalar;
 
     sfstep = SF_STEP_AMPL;
+    /* Pre-calculate lookup table for 10^(sfac / sfstep) = 2^(sfac / 4) */
+    for (i = -GAIN_LUT_BIAS; i < GAIN_LUT_SIZE - GAIN_LUT_BIAS; i++)
+        gain_lut[i + GAIN_LUT_BIAS] = powf(10.0f, (float)i / sfstep);
+
+    /* Pre-multiply width logarithm by SF_STEP_ENRG (= sfstep / 2) */
+    for (i = 1; i < 128; i++)
+        log10_width_sf_lut[i] = log10f((float)i) * SF_STEP_ENRG;
+
     /* One-time constant: computed in double so the stored float is
      * correctly rounded, at zero runtime cost. */
     max_quant_limit = (float)pow((double)MAX_HUFF_ESC_VAL + 1.0 - (double)MAGIC_NUMBER, 4.0/3.0);
 }
 
+static inline float sfac_to_gain(int sfac)
+{
+    unsigned int idx = (unsigned int)(sfac + GAIN_LUT_BIAS);
+    if (idx < GAIN_LUT_SIZE)
+        return gain_lut[idx];
+    return powf(10.0f, (float)sfac / sfstep);
+}
+
 /* sfac and gain are coupled; clamping one forces a recompute of the other. */
 static float gain_with_overflow_clamp(int *sfac, float band_peak)
 {
-    float gain = powf(10, *sfac / sfstep);
+    float gain = sfac_to_gain(*sfac);
     if (band_peak > 0.0f && gain * band_peak > max_quant_limit)
     {
         gain = max_quant_limit / band_peak;
         *sfac = (int)floorf(log10f(gain) * sfstep);
-        gain = powf(10, *sfac / sfstep);
+        gain = sfac_to_gain(*sfac);
     }
     return gain;
 }
@@ -101,14 +127,15 @@ static void measure_band_energy(const CoderInfo * __restrict ci, const float * _
     for (sfb = 0; sfb < ci->sfbn; sfb++)
     {
         int lo = ci->sfb_offset[sfb], hi = ci->sfb_offset[sfb + 1];
+        int len = hi - lo;
         float sum = 0.0f, peak = 0.0f;
         int w;
 
         for (w = 0; w < gsize; w++)
         {
-            const float *line = xr0 + w * BLOCK_LEN_SHORT + lo;
+            const float * __restrict line = xr0 + w * BLOCK_LEN_SHORT + lo;
             int k;
-            for (k = 0; k < hi - lo; k++)
+            for (k = 0; k < len; k++)
             {
                 float e = line[k] * line[k];
                 sum += e;
@@ -253,17 +280,22 @@ static void assign_band_codebooks(CoderInfo * __restrict ci, const float * __res
             continue;
         }
 
+        /* Log-domain identity: log10(target/rms) = log10(target) - 0.5*log10(avg) + 0.5*log10(width).
+         * Reuses sf_enrg_avg (log10(avg) * SF_STEP_ENRG) shared with PNS to avoid division and sqrtf. */
+        float sf_enrg_avg = log10f(avg_per_window) * SF_STEP_ENRG;
+
         /* PNS is fine inside TNS-covered bands -- the decoder's inverse
          * TNS filter shapes the substituted noise too. */
         if (target[sb] < pns_threshold)
         {
             ci->book[band] = HCB_PNS;
-            ci->sf[band] += lrintf(log10f(avg_per_window) * SF_STEP_ENRG);
+            ci->sf[band] += lrintf(sf_enrg_avg);
             ci->bandcnt++;
             continue;
         }
 
-        int sfac = lrintf(log10f(target[sb] / rms) * sfstep);
+        float log10_w_sf = (width < 128) ? log10_width_sf_lut[width] : log10f((float)width) * SF_STEP_ENRG;
+        int sfac = lrintf(log10f(target[sb]) * sfstep - sf_enrg_avg + log10_w_sf);
         int sf_rel = SF_OFFSET - sfac;
         int sf_bias = ci->sf[band];
 
@@ -276,11 +308,14 @@ static void assign_band_codebooks(CoderInfo * __restrict ci, const float * __res
             int sf_abs;
             float gain = resolve_band_gain(sfac, sf_bias, bandpeak[sb], *p_last_abs, &sf_rel, &sf_abs);
             int xi[FRAME_LEN];
-            int win;
+            int win, maxq = 0;
 
             for (win = 0; win < gsize; win++)
-                qfunc(xr0 + win * BLOCK_LEN_SHORT + lo, xi + win * width, width, gain);
-            huffbook(ci, xi, gsize * width);
+            {
+                int qm = qfunc(xr0 + win * BLOCK_LEN_SHORT + lo, xi + win * width, width, gain);
+                if (qm > maxq) maxq = qm;
+            }
+            huffbook(ci, xi, gsize * width, maxq);
             *p_last_abs = sf_abs;
         }
 

libfaac/quantize_sse.c

@@ -21,13 +21,14 @@
 #include <math.h>
 #include "quantize.h"
 
-void quantize_sse2(const float * __restrict xr, int * __restrict xi, int n, float sfacfix)
+int quantize_sse2(const float * __restrict xr, int * __restrict xi, int n, float sfacfix)
 {
     const __m128 zero = _mm_setzero_ps();
     const __m128 sfac = _mm_set1_ps(sfacfix);
     const __m128 magic = _mm_set1_ps(MAGIC_NUMBER);
     // Mask to strip the sign bit (0x7FFFFFFF)
     const __m128 abs_mask = _mm_castsi128_ps(_mm_set1_epi32(0x7FFFFFFF));
+    __m128i max_vec = _mm_setzero_si128();
     int cnt = 0;
 
     // Process 4 elements per iteration
@@ -47,6 +48,8 @@ void quantize_sse2(const float * __restrict xr, int * __restrict xi, int n, floa
 
         // Convert to integer
         __m128i xi_vec = _mm_cvttps_epi32(x);
+        __m128i mask = _mm_cmpgt_epi32(xi_vec, max_vec);
+        max_vec = _mm_or_si128(_mm_and_si128(mask, xi_vec), _mm_andnot_si128(mask, max_vec));
 
         // Bitwise Sign Fix: (val ^ mask) - mask
         __m128i m_int = _mm_castps_si128(sign_mask);
@@ -55,14 +58,24 @@ void quantize_sse2(const float * __restrict xr, int * __restrict xi, int n, floa
         _mm_storeu_si128((__m128i*)&xi[cnt], xi_vec);
     }
 
+    int maxq_arr[4];
+    _mm_storeu_si128((__m128i*)maxq_arr, max_vec);
+    int maxq = maxq_arr[0];
+    if (maxq_arr[1] > maxq) maxq = maxq_arr[1];
+    if (maxq_arr[2] > maxq) maxq = maxq_arr[2];
+    if (maxq_arr[3] > maxq) maxq = maxq_arr[3];
+
     // Safe scalar remainder loop for widths not multiple of 4
     for (; cnt < n; cnt++)
     {
-	float val = xr[cnt];
+        float val = xr[cnt];
         float tmp = fabsf(val);
         tmp *= sfacfix;
         tmp = sqrtf(tmp * sqrtf(tmp));
         int q = (int)(tmp + (float)MAGIC_NUMBER);
+        if (q > maxq) maxq = q;
         xi[cnt] = (val < 0) ? -q : q;
     }
+
+    return maxq;
 }