X-Git-Url: http://review.tizen.org/git/?a=blobdiff_plain;f=src%2FlibFLAC%2Flpc.c;h=4c10f289421d49e0e3cce93baa38caeb4b37f133;hb=afd8107872c6a877b66957ee192b43530782c6ec;hp=27cfc4b8c203b390497dfbffcf8d42d6db3d24c5;hpb=c625f9033ce47e737e33178b7a27113c904a2455;p=platform%2Fupstream%2Fflac.git diff --git a/src/libFLAC/lpc.c b/src/libFLAC/lpc.c index 27cfc4b..4c10f28 100644 --- a/src/libFLAC/lpc.c +++ b/src/libFLAC/lpc.c @@ -1,62 +1,104 @@ /* libFLAC - Free Lossless Audio Codec library - * Copyright (C) 2000,2001 Josh Coalson + * Copyright (C) 2000,2001,2002,2003 Josh Coalson * - * This library is free software; you can redistribute it and/or - * modify it under the terms of the GNU Library General Public - * License as published by the Free Software Foundation; either - * version 2 of the License, or (at your option) any later version. + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: * - * This library is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Library General Public License for more details. + * - Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. * - * You should have received a copy of the GNU Library General Public - * License along with this library; if not, write to the - * Free Software Foundation, Inc., 59 Temple Place - Suite 330, - * Boston, MA 02111-1307, USA. + * - Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * - Neither the name of the Xiph.org Foundation nor the names of its + * contributors may be used to endorse or promote products derived from + * this software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR + * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ -#include #include -#include +#include "FLAC/assert.h" #include "FLAC/format.h" +#include "private/bitmath.h" #include "private/lpc.h" +#if defined DEBUG || defined FLAC__OVERFLOW_DETECT || defined FLAC__OVERFLOW_DETECT_VERBOSE +#include +#endif #ifndef M_LN2 /* math.h in VC++ doesn't seem to have this (how Microsoft is that?) */ #define M_LN2 0.69314718055994530942 #endif -void FLAC__lpc_compute_autocorrelation(const real data[], unsigned data_len, unsigned lag, real autoc[]) +void FLAC__lpc_compute_autocorrelation(const FLAC__real data[], unsigned data_len, unsigned lag, FLAC__real autoc[]) { - real d; + /* a readable, but slower, version */ +#if 0 + FLAC__real d; unsigned i; - assert(lag > 0); - assert(lag <= data_len); + FLAC__ASSERT(lag > 0); + FLAC__ASSERT(lag <= data_len); while(lag--) { for(i = lag, d = 0.0; i < data_len; i++) d += data[i] * data[i - lag]; autoc[lag] = d; } +#endif + + /* + * this version tends to run faster because of better data locality + * ('data_len' is usually much larger than 'lag') + */ + FLAC__real d; + unsigned sample, coeff; + const unsigned limit = data_len - lag; + + FLAC__ASSERT(lag > 0); + FLAC__ASSERT(lag <= data_len); + + for(coeff = 0; coeff < lag; coeff++) + autoc[coeff] = 0.0; + for(sample = 0; sample <= limit; sample++) { + d = data[sample]; + for(coeff = 0; coeff < lag; coeff++) + autoc[coeff] += d * data[sample+coeff]; + } + for(; sample < data_len; sample++) { + d = data[sample]; + for(coeff = 0; coeff < data_len - sample; coeff++) + autoc[coeff] += d * data[sample+coeff]; + } } -void FLAC__lpc_compute_lp_coefficients(const real autoc[], unsigned max_order, real lp_coeff[][FLAC__MAX_LPC_ORDER], real error[]) +void FLAC__lpc_compute_lp_coefficients(const FLAC__real autoc[], unsigned max_order, FLAC__real lp_coeff[][FLAC__MAX_LPC_ORDER], FLAC__real error[]) { unsigned i, j; - real r, err, ref[FLAC__MAX_LPC_ORDER], lpc[FLAC__MAX_LPC_ORDER]; + double r, err, ref[FLAC__MAX_LPC_ORDER], lpc[FLAC__MAX_LPC_ORDER]; - assert(0 < max_order); - assert(max_order <= FLAC__MAX_LPC_ORDER); - assert(autoc[0] != 0.0); + FLAC__ASSERT(0 < max_order); + FLAC__ASSERT(max_order <= FLAC__MAX_LPC_ORDER); + FLAC__ASSERT(autoc[0] != 0.0); err = autoc[0]; 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); @@ -64,7 +106,7 @@ void FLAC__lpc_compute_lp_coefficients(const real autoc[], unsigned max_order, r /* Update LPC coefficients and total error. */ lpc[i]=r; for(j = 0; j < (i>>1); j++) { - real tmp = lpc[j]; + double tmp = lpc[j]; lpc[j] += r * lpc[i-1-j]; lpc[i-1-j] += r * tmp; } @@ -75,71 +117,27 @@ void FLAC__lpc_compute_lp_coefficients(const real autoc[], unsigned max_order, r /* save this order */ for(j = 0; j <= i; j++) - lp_coeff[i][j] = -lpc[j]; /* N.B. why do we have to negate here? */ - error[i] = err; + lp_coeff[i][j] = (FLAC__real)(-lpc[j]); /* negate FIR filter coeff to get predictor coeff */ + error[i] = (FLAC__real)err; } } -#if 0 -int FLAC__lpc_quantize_coefficients(const real lp_coeff[], unsigned order, unsigned precision, unsigned bits_per_sample, int32 qlp_coeff[], int *shift) +int FLAC__lpc_quantize_coefficients(const FLAC__real lp_coeff[], unsigned order, unsigned precision, FLAC__int32 qlp_coeff[], int *shift) { unsigned i; - real d, rprecision = (real)precision, maxlog = -1e99, minlog = 1e99; - - assert(bits_per_sample > 0); - assert(bits_per_sample <= sizeof(int32)*8); - assert(precision >= FLAC__MIN_QLP_COEFF_PRECISION); - assert(precision + bits_per_sample < sizeof(int32)*8); -#ifdef NDEBUG - (void)bits_per_sample; /* silence compiler warning about unused parameter */ -#endif - - for(i = 0; i < order; i++) { - if(lp_coeff[i] == 0.0) - continue; - d = log(fabs(lp_coeff[i])) / M_LN2; - if(d > maxlog) - maxlog = d; - if(d < minlog) - minlog = d; - } - if(maxlog < minlog) - return 2; - else if(maxlog - minlog >= (real)(precision+1)) - return 1; - else if((rprecision-1.0) - maxlog >= (real)(precision+1)) - rprecision = (real)precision + maxlog + 1.0; - - *shift = (int)floor((rprecision-1.0) - maxlog); /* '-1' because *shift can be negative and the sign bit costs 1 bit */ - if(*shift > (int)precision || *shift <= -(int)precision) { - fprintf(stderr, "@@@ FLAC__lpc_quantize_coefficients(): ERROR: *shift=%d, maxlog=%f, minlog=%f, precision=%u, rprecision=%f\n", *shift, maxlog, minlog, precision, rprecision); - return 1; - } - - if(*shift != 0) { /* just to avoid wasting time... */ - for(i = 0; i < order; i++) - qlp_coeff[i] = (int32)floor(lp_coeff[i] * (real)(1 << *shift)); - } - return 0; -} -#endif + double d, cmax = -1e32; + FLAC__int32 qmax, qmin; + const int max_shiftlimit = (1 << (FLAC__SUBFRAME_LPC_QLP_SHIFT_LEN-1)) - 1; + const int min_shiftlimit = -max_shiftlimit - 1; -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;//@@@, cmin = 1e99; - - assert(bits_per_sample > 0); - assert(bits_per_sample <= sizeof(int32)*8); - assert(precision > 0); - assert(precision >= FLAC__MIN_QLP_COEFF_PRECISION); - assert(precision + bits_per_sample < sizeof(int32)*8); -#ifdef NDEBUG - (void)bits_per_sample; /* silence compiler warning about unused parameter */ -#endif + FLAC__ASSERT(precision > 0); + FLAC__ASSERT(precision >= FLAC__MIN_QLP_COEFF_PRECISION); /* drop one bit for the sign; from here on out we consider only |lp_coeff[i]| */ precision--; + qmax = 1 << precision; + qmin = -qmax; + qmax--; for(i = 0; i < order; i++) { if(lp_coeff[i] == 0.0) @@ -147,65 +145,99 @@ int FLAC__lpc_quantize_coefficients(const real lp_coeff[], unsigned order, unsig d = fabs(lp_coeff[i]); if(d > cmax) cmax = d; -//@@@ if(d < cmin) -//@@@ cmin = d; } -//@@@ if(cmax < cmin) - if(cmax < 0) { - /* => coeffients are all 0, which means our constant-detect didn't work */ -fprintf(stderr,"@@@ LPCQ ERROR, all lpc_coeffs are 0\n"); +redo_it: + if(cmax <= 0.0) { + /* => coefficients are all 0, which means our constant-detect didn't work */ return 2; } else { -//@@@ const int minshift = (int)precision - floor(log(cmin) / M_LN2) - 1; - const int maxshift = (int)precision - floor(log(cmax) / M_LN2) - 1; -//@@@ assert(maxshift >= minshift); - const int max_shiftlimit = (1 << (FLAC__SUBFRAME_LPC_QLP_SHIFT_LEN-1)) - 1; - const int min_shiftlimit = -max_shiftlimit - 1; + int log2cmax; - *shift = maxshift; + (void)frexp(cmax, &log2cmax); + log2cmax--; + *shift = (int)precision - log2cmax - 1; if(*shift < min_shiftlimit || *shift > max_shiftlimit) { -fprintf(stderr,"@@@ LPCQ ERROR, shift is outside shiftlimit\n"); - return 1; +#if 0 + /*@@@ this does not seem to help at all, but was not extensively tested either: */ + if(*shift > max_shiftlimit) + *shift = max_shiftlimit; + else +#endif + return 1; } } - if(*shift != 0) { /* just to avoid wasting time... */ - for(i = 0; i < order; i++) - qlp_coeff[i] = (int32)floor(lp_coeff[i] * (real)(1 << *shift)); + if(*shift >= 0) { + for(i = 0; i < order; i++) { + qlp_coeff[i] = (FLAC__int32)floor((double)lp_coeff[i] * (double)(1 << *shift)); + + /* double-check the result */ + if(qlp_coeff[i] > qmax || qlp_coeff[i] < qmin) { +#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, (double)lp_coeff[i] * (double)(1 << *shift), floor((double)lp_coeff[i] * (double)(1 << *shift))); +#endif + cmax *= 2.0; + goto redo_it; + } + } + } + else { /* (*shift < 0) */ + const int nshift = -(*shift); +#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((double)lp_coeff[i] / (double)(1 << nshift)); + + /* double-check the result */ + if(qlp_coeff[i] > qmax || qlp_coeff[i] < qmin) { +#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, (double)lp_coeff[i] / (double)(1 << nshift), floor((double)lp_coeff[i] / (double)(1 << nshift))); +#endif + cmax *= 2.0; + goto redo_it; + } + } } + return 0; } -void FLAC__lpc_compute_residual_from_qlp_coefficients(const int32 data[], unsigned data_len, const int32 qlp_coeff[], unsigned order, int lp_quantization, int32 residual[]) +void FLAC__lpc_compute_residual_from_qlp_coefficients(const FLAC__int32 data[], unsigned data_len, const FLAC__int32 qlp_coeff[], unsigned order, int lp_quantization, FLAC__int32 residual[]) { -#ifdef FLAC_OVERFLOW_DETECT - int64 sumo; +#ifdef FLAC__OVERFLOW_DETECT + FLAC__int64 sumo; #endif unsigned i, j; - int32 sum; - const int32 *history; + FLAC__int32 sum; + const FLAC__int32 *history; -#ifdef FLAC_OVERFLOW_DETECT_VERBOSE +#ifdef FLAC__OVERFLOW_DETECT_VERBOSE fprintf(stderr,"FLAC__lpc_compute_residual_from_qlp_coefficients: data_len=%d, order=%u, lpq=%d",data_len,order,lp_quantization); for(i=0;i 0); + FLAC__ASSERT(order > 0); for(i = 0; i < data_len; i++) { -#ifdef FLAC_OVERFLOW_DETECT +#ifdef FLAC__OVERFLOW_DETECT sumo = 0; #endif sum = 0; history = data; for(j = 0; j < order; j++) { sum += qlp_coeff[j] * (*(--history)); -#ifdef FLAC_OVERFLOW_DETECT - sumo += (int64)qlp_coeff[j] * (int64)(*history); - if(sumo > 2147483647ll || sumo < -2147483648ll) { +#ifdef FLAC__OVERFLOW_DETECT + sumo += (FLAC__int64)qlp_coeff[j] * (FLAC__int64)(*history); +#if defined _MSC_VER + if(sumo > 2147483647I64 || sumo < -2147483648I64) +#else + if(sumo > 2147483647ll || sumo < -2147483648ll) +#endif + { fprintf(stderr,"FLAC__lpc_compute_residual_from_qlp_coefficients: OVERFLOW, i=%u, j=%u, c=%d, d=%d, sumo=%lld\n",i,j,qlp_coeff[j],*history,sumo); } #endif @@ -223,34 +255,72 @@ void FLAC__lpc_compute_residual_from_qlp_coefficients(const int32 data[], unsign */ } -void FLAC__lpc_restore_signal(const int32 residual[], unsigned data_len, const int32 qlp_coeff[], unsigned order, int lp_quantization, int32 data[]) +void FLAC__lpc_compute_residual_from_qlp_coefficients_wide(const FLAC__int32 data[], unsigned data_len, const FLAC__int32 qlp_coeff[], unsigned order, int lp_quantization, FLAC__int32 residual[]) { -#ifdef FLAC_OVERFLOW_DETECT - int64 sumo; + unsigned i, j; + FLAC__int64 sum; + const FLAC__int32 *history; + +#ifdef FLAC__OVERFLOW_DETECT_VERBOSE + fprintf(stderr,"FLAC__lpc_compute_residual_from_qlp_coefficients_wide: data_len=%d, order=%u, lpq=%d",data_len,order,lp_quantization); + for(i=0;i 0); + + for(i = 0; i < data_len; i++) { + sum = 0; + history = data; + for(j = 0; j < order; j++) + sum += (FLAC__int64)qlp_coeff[j] * (FLAC__int64)(*(--history)); +#ifdef FLAC__OVERFLOW_DETECT + if(FLAC__bitmath_silog2_wide(sum >> lp_quantization) > 32) { + fprintf(stderr,"FLAC__lpc_compute_residual_from_qlp_coefficients_wide: OVERFLOW, i=%u, sum=%lld\n", i, sum >> lp_quantization); + break; + } + if(FLAC__bitmath_silog2_wide((FLAC__int64)(*data) - (sum >> lp_quantization)) > 32) { + fprintf(stderr,"FLAC__lpc_compute_residual_from_qlp_coefficients_wide: OVERFLOW, i=%u, data=%d, sum=%lld, residual=%lld\n", i, *data, sum >> lp_quantization, (FLAC__int64)(*data) - (sum >> lp_quantization)); + break; + } +#endif + *(residual++) = *(data++) - (FLAC__int32)(sum >> lp_quantization); + } +} + +void FLAC__lpc_restore_signal(const FLAC__int32 residual[], unsigned data_len, const FLAC__int32 qlp_coeff[], unsigned order, int lp_quantization, FLAC__int32 data[]) +{ +#ifdef FLAC__OVERFLOW_DETECT + FLAC__int64 sumo; #endif unsigned i, j; - int32 sum; - const int32 *history; + FLAC__int32 sum; + const FLAC__int32 *history; -#ifdef FLAC_OVERFLOW_DETECT_VERBOSE +#ifdef FLAC__OVERFLOW_DETECT_VERBOSE fprintf(stderr,"FLAC__lpc_restore_signal: data_len=%d, order=%u, lpq=%d",data_len,order,lp_quantization); for(i=0;i 0); + FLAC__ASSERT(order > 0); for(i = 0; i < data_len; i++) { -#ifdef FLAC_OVERFLOW_DETECT +#ifdef FLAC__OVERFLOW_DETECT sumo = 0; #endif sum = 0; history = data; for(j = 0; j < order; j++) { sum += qlp_coeff[j] * (*(--history)); -#ifdef FLAC_OVERFLOW_DETECT - sumo += (int64)qlp_coeff[j] * (int64)(*history); - if(sumo > 2147483647ll || sumo < -2147483648ll) { +#ifdef FLAC__OVERFLOW_DETECT + sumo += (FLAC__int64)qlp_coeff[j] * (FLAC__int64)(*history); +#if defined _MSC_VER + if(sumo > 2147483647I64 || sumo < -2147483648I64) +#else + if(sumo > 2147483647ll || sumo < -2147483648ll) +#endif + { fprintf(stderr,"FLAC__lpc_restore_signal: OVERFLOW, i=%u, j=%u, c=%d, d=%d, sumo=%lld\n",i,j,qlp_coeff[j],*history,sumo); } #endif @@ -268,39 +338,83 @@ 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) +void FLAC__lpc_restore_signal_wide(const FLAC__int32 residual[], unsigned data_len, const FLAC__int32 qlp_coeff[], unsigned order, int lp_quantization, FLAC__int32 data[]) { - real escale; + unsigned i, j; + FLAC__int64 sum; + const FLAC__int32 *history; - assert(lpc_error >= 0.0); /* the error can never be negative */ - assert(total_samples > 0); +#ifdef FLAC__OVERFLOW_DETECT_VERBOSE + fprintf(stderr,"FLAC__lpc_restore_signal_wide: data_len=%d, order=%u, lpq=%d",data_len,order,lp_quantization); + for(i=0;i 0); + + for(i = 0; i < data_len; i++) { + sum = 0; + history = data; + for(j = 0; j < order; j++) + sum += (FLAC__int64)qlp_coeff[j] * (FLAC__int64)(*(--history)); +#ifdef FLAC__OVERFLOW_DETECT + if(FLAC__bitmath_silog2_wide(sum >> lp_quantization) > 32) { + fprintf(stderr,"FLAC__lpc_restore_signal_wide: OVERFLOW, i=%u, sum=%lld\n", i, sum >> lp_quantization); + break; + } + if(FLAC__bitmath_silog2_wide((FLAC__int64)(*residual) + (sum >> lp_quantization)) > 32) { + fprintf(stderr,"FLAC__lpc_restore_signal_wide: OVERFLOW, i=%u, residual=%d, sum=%lld, data=%lld\n", i, *residual, sum >> lp_quantization, (FLAC__int64)(*residual) + (sum >> lp_quantization)); + break; + } +#endif + *(data++) = *(residual++) + (FLAC__int32)(sum >> lp_quantization); + } +} + +FLAC__real FLAC__lpc_compute_expected_bits_per_residual_sample(FLAC__real lpc_error, unsigned total_samples) +{ + double error_scale; + + FLAC__ASSERT(total_samples > 0); - escale = 0.5 * M_LN2 * M_LN2 / (real)total_samples; + error_scale = 0.5 * M_LN2 * M_LN2 / (FLAC__real)total_samples; + return FLAC__lpc_compute_expected_bits_per_residual_sample_with_error_scale(lpc_error, error_scale); +} + +FLAC__real FLAC__lpc_compute_expected_bits_per_residual_sample_with_error_scale(FLAC__real lpc_error, double error_scale) +{ if(lpc_error > 0.0) { - real bps = 0.5 * log(escale * lpc_error) / M_LN2; + FLAC__real bps = (FLAC__real)((double)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 (FLAC__real)1e32; + } else { return 0.0; } } -unsigned FLAC__lpc_compute_best_order(const real lpc_error[], unsigned max_order, unsigned total_samples, unsigned bits_per_signal_sample) +unsigned FLAC__lpc_compute_best_order(const FLAC__real lpc_error[], unsigned max_order, unsigned total_samples, unsigned bits_per_signal_sample) { unsigned order, best_order; - real best_bits, tmp_bits; + FLAC__real best_bits, tmp_bits; + double error_scale; + + FLAC__ASSERT(max_order > 0); + FLAC__ASSERT(total_samples > 0); - assert(max_order > 0); + error_scale = 0.5 * M_LN2 * M_LN2 / (FLAC__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) * (FLAC__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) * (FLAC__real)(total_samples - order) + (FLAC__real)(order * bits_per_signal_sample); if(tmp_bits < best_bits) { best_order = order; best_bits = tmp_bits;