be more lenient with negative lpc_errors due to accuracy problems with float
authorJosh Coalson <jcoalson@users.sourceforce.net>
Thu, 24 May 2001 19:27:08 +0000 (19:27 +0000)
committerJosh Coalson <jcoalson@users.sourceforce.net>
Thu, 24 May 2001 19:27:08 +0000 (19:27 +0000)
src/libFLAC/lpc.c

index db5d2c4..d85902c 100644 (file)
@@ -85,7 +85,7 @@ void FLAC__lpc_compute_lp_coefficients(const real autoc[], unsigned max_order, r
 
        for(i = 0; i < max_order; i++) {
                /* Sum up this iteration's reflection coefficient. */
-               r =autoc[i+1];
+               r = -autoc[i+1];
                for(j = 0; j < i; j++)
                        r -= lpc[j] * autoc[i-j];
                ref[i] = (r/=err);
@@ -101,6 +101,7 @@ void FLAC__lpc_compute_lp_coefficients(const real autoc[], unsigned max_order, r
                        lpc[j] += lpc[j] * r;
 
                err *= (1.0 - r * r);
+if(err<0.0)fprintf(stderr,"@@@ err went negative, order=%u (%f,%f)\n",i,err,r);
 
                /* save this order */
                for(j = 0; j <= i; j++)
@@ -112,7 +113,7 @@ void FLAC__lpc_compute_lp_coefficients(const real autoc[], unsigned max_order, r
 int FLAC__lpc_quantize_coefficients(const real lp_coeff[], unsigned order, unsigned precision, unsigned bits_per_sample, int32 qlp_coeff[], int *shift)
 {
        unsigned i;
-       real d, cmax = -1e99;
+       real d, cmax = -1e10;
 
        assert(bits_per_sample > 0);
        assert(bits_per_sample <= sizeof(int32)*8);
@@ -248,20 +249,39 @@ void FLAC__lpc_restore_signal(const int32 residual[], unsigned data_len, const i
 
 real FLAC__lpc_compute_expected_bits_per_residual_sample(real lpc_error, unsigned total_samples)
 {
-       real escale;
+       real error_scale;
 
-       assert(lpc_error >= 0.0); /* the error can never be negative */
        assert(total_samples > 0);
 
-       escale = 0.5 * M_LN2 * M_LN2 / (real)total_samples;
+       error_scale = 0.5 * M_LN2 * M_LN2 / (real)total_samples;
 
        if(lpc_error > 0.0) {
-               real bps = 0.5 * log(escale * lpc_error) / M_LN2;
+               real bps = 0.5 * log(error_scale * lpc_error) / M_LN2;
                if(bps >= 0.0)
                        return bps;
                else
                        return 0.0;
        }
+       else if(lpc_error < 0.0) { /* error should not be negative but can happen due to inadequate float resolution */
+               return 1e10;
+       }
+       else {
+               return 0.0;
+       }
+}
+
+real FLAC__lpc_compute_expected_bits_per_residual_sample_with_error_scale(real lpc_error, real error_scale)
+{
+       if(lpc_error > 0.0) {
+               real bps = 0.5 * log(error_scale * lpc_error) / M_LN2;
+               if(bps >= 0.0)
+                       return bps;
+               else
+                       return 0.0;
+       }
+       else if(lpc_error < 0.0) { /* error should not be negative but can happen due to inadequate float resolution */
+               return 1e10;
+       }
        else {
                return 0.0;
        }
@@ -271,14 +291,18 @@ unsigned FLAC__lpc_compute_best_order(const real lpc_error[], unsigned max_order
 {
        unsigned order, best_order;
        real best_bits, tmp_bits;
+       const error_scale;
 
        assert(max_order > 0);
+       assert(total_samples > 0);
+
+       error_scale = 0.5 * M_LN2 * M_LN2 / (real)total_samples;
 
        best_order = 0;
-       best_bits = FLAC__lpc_compute_expected_bits_per_residual_sample(lpc_error[0], total_samples) * (real)total_samples;
+       best_bits = FLAC__lpc_compute_expected_bits_per_residual_sample_with_error_scale(lpc_error[0], error_scale) * (real)total_samples;
 
        for(order = 1; order < max_order; order++) {
-               tmp_bits = FLAC__lpc_compute_expected_bits_per_residual_sample(lpc_error[order], total_samples) * (real)(total_samples - order) + (real)(order * bits_per_signal_sample);
+               tmp_bits = FLAC__lpc_compute_expected_bits_per_residual_sample_with_error_scale(lpc_error[order], error_scale) * (real)(total_samples - order) + (real)(order * bits_per_signal_sample);
                if(tmp_bits < best_bits) {
                        best_order = order;
                        best_bits = tmp_bits;