add dithering to lp coeff quantization

diff --git a/src/libFLAC/lpc.c b/src/libFLAC/lpc.c
index ffcc36d..c21b6dd 100644 (file)
--- a/src/libFLAC/lpc.c
+++ b/src/libFLAC/lpc.c
@@ -152,10 +152,8 @@ void FLAC__lpc_compute_lp_coefficients(const FLAC__real autoc[], unsigned *max_o
 int FLAC__lpc_quantize_coefficients(const FLAC__real lp_coeff[], unsigned order, unsigned precision, FLAC__int32 qlp_coeff[], int *shift)
 {
        unsigned i;
-       FLAC__double d, cmax = -1e32;
+       FLAC__double cmax;
        FLAC__int32 qmax, qmin;
-       const int max_shiftlimit = (1 << (FLAC__SUBFRAME_LPC_QLP_SHIFT_LEN-1)) - 1;
-       const int min_shiftlimit = -max_shiftlimit - 1;
 
        FLAC__ASSERT(precision > 0);
        FLAC__ASSERT(precision >= FLAC__MIN_QLP_COEFF_PRECISION);
@@ -166,19 +164,21 @@ int FLAC__lpc_quantize_coefficients(const FLAC__real lp_coeff[], unsigned order,
        qmin = -qmax;
        qmax--;
 
+       /* calc cmax = max( |lp_coeff[i]| ) */
+       cmax = 0.0;
        for(i = 0; i < order; i++) {
-               if(lp_coeff[i] == 0.0)
-                       continue;
-               d = fabs(lp_coeff[i]);
+               const FLAC__double d = fabs(lp_coeff[i]);
                if(d > cmax)
                        cmax = d;
        }
-redo_it:
+
        if(cmax <= 0.0) {
                /* => coefficients are all 0, which means our constant-detect didn't work */
                return 2;
        }
        else {
+               const int max_shiftlimit = (1 << (FLAC__SUBFRAME_LPC_QLP_SHIFT_LEN-1)) - 1;
+               const int min_shiftlimit = -max_shiftlimit - 1;
                int log2cmax;
 
                (void)frexp(cmax, &log2cmax);
@@ -186,6 +186,9 @@ redo_it:
                *shift = (int)precision - log2cmax - 1;
 
                if(*shift < min_shiftlimit || *shift > max_shiftlimit) {
+#ifdef FLAC__OVERFLOW_DETECT
+                       fprintf(stderr,"FLAC__lpc_quantize_coefficients: shift out of limit: shift=%d cmax=%f precision=%u\n",q,qmax,*shift,cmax,precision+1);
+#endif
 #if 0
                        /*@@@ this does not seem to help at all, but was not extensively tested either: */
                        if(*shift > max_shiftlimit)
@@ -197,36 +200,53 @@ redo_it:
        }
 
        if(*shift >= 0) {
+               FLAC__double error = 0.0;
+               FLAC__int32 q;
                for(i = 0; i < order; i++) {
-                       qlp_coeff[i] = (FLAC__int32)floor((FLAC__double)lp_coeff[i] * (FLAC__double)(1 << *shift));
-
-                       /* double-check the result */
-                       if(qlp_coeff[i] > qmax || qlp_coeff[i] < qmin) {
+                       error += lp_coeff[i] * (1 << *shift);
+                       q = lround(error); /* round() is also suitable */
 #ifdef FLAC__OVERFLOW_DETECT
-                               fprintf(stderr,"FLAC__lpc_quantize_coefficients: compensating for overflow, qlp_coeff[%u]=%d, lp_coeff[%u]=%f, cmax=%f, precision=%u, shift=%d, q=%f, f(q)=%f\n", i, qlp_coeff[i], i, lp_coeff[i], cmax, precision, *shift, (FLAC__double)lp_coeff[i] * (FLAC__double)(1 << *shift), floor((FLAC__double)lp_coeff[i] * (FLAC__double)(1 << *shift)));
+                       if(q > qmax)
+                               fprintf(stderr,"FLAC__lpc_quantize_coefficients: quantizer overflow: q>qmax %d>%d shift=%d cmax=%f precision=%u lpc[%u]=%f\n",q,qmax,*shift,cmax,precision+1,i,lp_coeff[i]);
+                       else if(q < qmin)
+                               fprintf(stderr,"FLAC__lpc_quantize_coefficients: quantizer overflow: q<qmin %d<%d shift=%d cmax=%f precision=%u lpc[%u]=%f\n",q,qmin,*shift,cmax,precision+1,i,lp_coeff[i]);
 #endif
-                               cmax *= 2.0;
-                               goto redo_it;
-                       }
+                       if(q > qmax)
+                               q = qmax;
+                       else if(q < qmin)
+                               q = qmin;
+                       error -= q;
+                       qlp_coeff[i] = q;
                }
        }
-       else { /* (*shift < 0) */
+       /* negative shift is very rare but due to design flaw, negative shift is
+        * a NOP in the decoder, so it must be handled specially by scaling down
+        * coeffs
+        */
+       else {
                const int nshift = -(*shift);
+               FLAC__double error = 0.0;
+               FLAC__int32 q;
 #ifdef DEBUG
                fprintf(stderr,"FLAC__lpc_quantize_coefficients: negative shift = %d\n", *shift);
 #endif
                for(i = 0; i < order; i++) {
-                       qlp_coeff[i] = (FLAC__int32)floor((FLAC__double)lp_coeff[i] / (FLAC__double)(1 << nshift));
-
-                       /* double-check the result */
-                       if(qlp_coeff[i] > qmax || qlp_coeff[i] < qmin) {
+                       error += lp_coeff[i] / (1 << nshift);
+                       q = lround(error); /* round() is also suitable */
 #ifdef FLAC__OVERFLOW_DETECT
-                               fprintf(stderr,"FLAC__lpc_quantize_coefficients: compensating for overflow, qlp_coeff[%u]=%d, lp_coeff[%u]=%f, cmax=%f, precision=%u, shift=%d, q=%f, f(q)=%f\n", i, qlp_coeff[i], i, lp_coeff[i], cmax, precision, *shift, (FLAC__double)lp_coeff[i] / (FLAC__double)(1 << nshift), floor((FLAC__double)lp_coeff[i] / (FLAC__double)(1 << nshift)));
+                       if(q > qmax)
+                               fprintf(stderr,"FLAC__lpc_quantize_coefficients: quantizer overflow: q>qmax %d>%d shift=%d cmax=%f precision=%u lpc[%u]=%f\n",q,qmax,*shift,cmax,precision+1,i,lp_coeff[i]);
+                       else if(q < qmin)
+                               fprintf(stderr,"FLAC__lpc_quantize_coefficients: quantizer overflow: q<qmin %d<%d shift=%d cmax=%f precision=%u lpc[%u]=%f\n",q,qmin,*shift,cmax,precision+1,i,lp_coeff[i]);
 #endif
-                               cmax *= 2.0;
-                               goto redo_it;
-                       }
+                       if(q > qmax)
+                               q = qmax;
+                       else if(q < qmin)
+                               q = qmin;
+                       error -= q;
+                       qlp_coeff[i] = q;
                }
+               *shift = 0;
        }
 
        return 0;