fixed.c

00001 /* libFLAC - Free Lossless Audio Codec library
00002  * Copyright (C) 2000-2009  Josh Coalson
00003  * Copyright (C) 2011-2014  Xiph.Org Foundation
00004  *
00005  * Redistribution and use in source and binary forms, with or without
00006  * modification, are permitted provided that the following conditions
00007  * are met:
00008  *
00009  * - Redistributions of source code must retain the above copyright
00010  * notice, this list of conditions and the following disclaimer.
00011  *
00012  * - Redistributions in binary form must reproduce the above copyright
00013  * notice, this list of conditions and the following disclaimer in the
00014  * documentation and/or other materials provided with the distribution.
00015  *
00016  * - Neither the name of the Xiph.org Foundation nor the names of its
00017  * contributors may be used to endorse or promote products derived from
00018  * this software without specific prior written permission.
00019  *
00020  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
00021  * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
00022  * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
00023  * A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
00024  * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00025  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00026  * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00027  * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00028  * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00029  * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00030  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00031  */
00032 
00033 #ifdef HAVE_CONFIG_H
00034 #  include <config.h>
00035 #endif
00036 
00037 #include <math.h>
00038 #include <string.h>
00039 #include "share/compat.h"
00040 #include "private/bitmath.h"
00041 #include "private/fixed.h"
00042 #include "private/macros.h"
00043 #include "FLAC/assert.h"
00044 
00045 #ifdef local_abs
00046 #undef local_abs
00047 #endif
00048 #define local_abs(x) ((unsigned)((x)<0? -(x) : (x)))
00049 
00050 #ifdef FLAC__INTEGER_ONLY_LIBRARY
00051 /* rbps stands for residual bits per sample
00052  *
00053  *             (ln(2) * err)
00054  * rbps = log  (-----------)
00055  *           2 (     n     )
00056  */
00057 static FLAC__fixedpoint local__compute_rbps_integerized(FLAC__uint32 err, FLAC__uint32 n)
00058 {
00059     FLAC__uint32 rbps;
00060     unsigned bits; /* the number of bits required to represent a number */
00061     int fracbits; /* the number of bits of rbps that comprise the fractional part */
00062 
00063     FLAC__ASSERT(sizeof(rbps) == sizeof(FLAC__fixedpoint));
00064     FLAC__ASSERT(err > 0);
00065     FLAC__ASSERT(n > 0);
00066 
00067     FLAC__ASSERT(n <= FLAC__MAX_BLOCK_SIZE);
00068     if(err <= n)
00069         return 0;
00070     /*
00071      * The above two things tell us 1) n fits in 16 bits; 2) err/n > 1.
00072      * These allow us later to know we won't lose too much precision in the
00073      * fixed-point division (err<<fracbits)/n.
00074      */
00075 
00076     fracbits = (8*sizeof(err)) - (FLAC__bitmath_ilog2(err)+1);
00077 
00078     err <<= fracbits;
00079     err /= n;
00080     /* err now holds err/n with fracbits fractional bits */
00081 
00082     /*
00083      * Whittle err down to 16 bits max.  16 significant bits is enough for
00084      * our purposes.
00085      */
00086     FLAC__ASSERT(err > 0);
00087     bits = FLAC__bitmath_ilog2(err)+1;
00088     if(bits > 16) {
00089         err >>= (bits-16);
00090         fracbits -= (bits-16);
00091     }
00092     rbps = (FLAC__uint32)err;
00093 
00094     /* Multiply by fixed-point version of ln(2), with 16 fractional bits */
00095     rbps *= FLAC__FP_LN2;
00096     fracbits += 16;
00097     FLAC__ASSERT(fracbits >= 0);
00098 
00099     /* FLAC__fixedpoint_log2 requires fracbits%4 to be 0 */
00100     {
00101         const int f = fracbits & 3;
00102         if(f) {
00103             rbps >>= f;
00104             fracbits -= f;
00105         }
00106     }
00107 
00108     rbps = FLAC__fixedpoint_log2(rbps, fracbits, (unsigned)(-1));
00109 
00110     if(rbps == 0)
00111         return 0;
00112 
00113     /*
00114      * The return value must have 16 fractional bits.  Since the whole part
00115      * of the base-2 log of a 32 bit number must fit in 5 bits, and fracbits
00116      * must be >= -3, these assertion allows us to be able to shift rbps
00117      * left if necessary to get 16 fracbits without losing any bits of the
00118      * whole part of rbps.
00119      *
00120      * There is a slight chance due to accumulated error that the whole part
00121      * will require 6 bits, so we use 6 in the assertion.  Really though as
00122      * long as it fits in 13 bits (32 - (16 - (-3))) we are fine.
00123      */
00124     FLAC__ASSERT((int)FLAC__bitmath_ilog2(rbps)+1 <= fracbits + 6);
00125     FLAC__ASSERT(fracbits >= -3);
00126 
00127     /* now shift the decimal point into place */
00128     if(fracbits < 16)
00129         return rbps << (16-fracbits);
00130     else if(fracbits > 16)
00131         return rbps >> (fracbits-16);
00132     else
00133         return rbps;
00134 }
00135 
00136 static FLAC__fixedpoint local__compute_rbps_wide_integerized(FLAC__uint64 err, FLAC__uint32 n)
00137 {
00138     FLAC__uint32 rbps;
00139     unsigned bits; /* the number of bits required to represent a number */
00140     int fracbits; /* the number of bits of rbps that comprise the fractional part */
00141 
00142     FLAC__ASSERT(sizeof(rbps) == sizeof(FLAC__fixedpoint));
00143     FLAC__ASSERT(err > 0);
00144     FLAC__ASSERT(n > 0);
00145 
00146     FLAC__ASSERT(n <= FLAC__MAX_BLOCK_SIZE);
00147     if(err <= n)
00148         return 0;
00149     /*
00150      * The above two things tell us 1) n fits in 16 bits; 2) err/n > 1.
00151      * These allow us later to know we won't lose too much precision in the
00152      * fixed-point division (err<<fracbits)/n.
00153      */
00154 
00155     fracbits = (8*sizeof(err)) - (FLAC__bitmath_ilog2_wide(err)+1);
00156 
00157     err <<= fracbits;
00158     err /= n;
00159     /* err now holds err/n with fracbits fractional bits */
00160 
00161     /*
00162      * Whittle err down to 16 bits max.  16 significant bits is enough for
00163      * our purposes.
00164      */
00165     FLAC__ASSERT(err > 0);
00166     bits = FLAC__bitmath_ilog2_wide(err)+1;
00167     if(bits > 16) {
00168         err >>= (bits-16);
00169         fracbits -= (bits-16);
00170     }
00171     rbps = (FLAC__uint32)err;
00172 
00173     /* Multiply by fixed-point version of ln(2), with 16 fractional bits */
00174     rbps *= FLAC__FP_LN2;
00175     fracbits += 16;
00176     FLAC__ASSERT(fracbits >= 0);
00177 
00178     /* FLAC__fixedpoint_log2 requires fracbits%4 to be 0 */
00179     {
00180         const int f = fracbits & 3;
00181         if(f) {
00182             rbps >>= f;
00183             fracbits -= f;
00184         }
00185     }
00186 
00187     rbps = FLAC__fixedpoint_log2(rbps, fracbits, (unsigned)(-1));
00188 
00189     if(rbps == 0)
00190         return 0;
00191 
00192     /*
00193      * The return value must have 16 fractional bits.  Since the whole part
00194      * of the base-2 log of a 32 bit number must fit in 5 bits, and fracbits
00195      * must be >= -3, these assertion allows us to be able to shift rbps
00196      * left if necessary to get 16 fracbits without losing any bits of the
00197      * whole part of rbps.
00198      *
00199      * There is a slight chance due to accumulated error that the whole part
00200      * will require 6 bits, so we use 6 in the assertion.  Really though as
00201      * long as it fits in 13 bits (32 - (16 - (-3))) we are fine.
00202      */
00203     FLAC__ASSERT((int)FLAC__bitmath_ilog2(rbps)+1 <= fracbits + 6);
00204     FLAC__ASSERT(fracbits >= -3);
00205 
00206     /* now shift the decimal point into place */
00207     if(fracbits < 16)
00208         return rbps << (16-fracbits);
00209     else if(fracbits > 16)
00210         return rbps >> (fracbits-16);
00211     else
00212         return rbps;
00213 }
00214 #endif
00215 
00216 #ifndef FLAC__INTEGER_ONLY_LIBRARY
00217 unsigned FLAC__fixed_compute_best_predictor(const FLAC__int32 data[], unsigned data_len, FLAC__float residual_bits_per_sample[FLAC__MAX_FIXED_ORDER+1])
00218 #else
00219 unsigned FLAC__fixed_compute_best_predictor(const FLAC__int32 data[], unsigned data_len, FLAC__fixedpoint residual_bits_per_sample[FLAC__MAX_FIXED_ORDER+1])
00220 #endif
00221 {
00222     FLAC__int32 last_error_0 = data[-1];
00223     FLAC__int32 last_error_1 = data[-1] - data[-2];
00224     FLAC__int32 last_error_2 = last_error_1 - (data[-2] - data[-3]);
00225     FLAC__int32 last_error_3 = last_error_2 - (data[-2] - 2*data[-3] + data[-4]);
00226     FLAC__int32 error, save;
00227     FLAC__uint32 total_error_0 = 0, total_error_1 = 0, total_error_2 = 0, total_error_3 = 0, total_error_4 = 0;
00228     unsigned i, order;
00229 
00230     for(i = 0; i < data_len; i++) {
00231         error  = data[i]     ; total_error_0 += local_abs(error);                      save = error;
00232         error -= last_error_0; total_error_1 += local_abs(error); last_error_0 = save; save = error;
00233         error -= last_error_1; total_error_2 += local_abs(error); last_error_1 = save; save = error;
00234         error -= last_error_2; total_error_3 += local_abs(error); last_error_2 = save; save = error;
00235         error -= last_error_3; total_error_4 += local_abs(error); last_error_3 = save;
00236     }
00237 
00238     if(total_error_0 < flac_min(flac_min(flac_min(total_error_1, total_error_2), total_error_3), total_error_4))
00239         order = 0;
00240     else if(total_error_1 < flac_min(flac_min(total_error_2, total_error_3), total_error_4))
00241         order = 1;
00242     else if(total_error_2 < flac_min(total_error_3, total_error_4))
00243         order = 2;
00244     else if(total_error_3 < total_error_4)
00245         order = 3;
00246     else
00247         order = 4;
00248 
00249     /* Estimate the expected number of bits per residual signal sample. */
00250     /* 'total_error*' is linearly related to the variance of the residual */
00251     /* signal, so we use it directly to compute E(|x|) */
00252     FLAC__ASSERT(data_len > 0 || total_error_0 == 0);
00253     FLAC__ASSERT(data_len > 0 || total_error_1 == 0);
00254     FLAC__ASSERT(data_len > 0 || total_error_2 == 0);
00255     FLAC__ASSERT(data_len > 0 || total_error_3 == 0);
00256     FLAC__ASSERT(data_len > 0 || total_error_4 == 0);
00257 #ifndef FLAC__INTEGER_ONLY_LIBRARY
00258     residual_bits_per_sample[0] = (FLAC__float)((total_error_0 > 0) ? log(M_LN2 * (FLAC__double)total_error_0 / (FLAC__double)data_len) / M_LN2 : 0.0);
00259     residual_bits_per_sample[1] = (FLAC__float)((total_error_1 > 0) ? log(M_LN2 * (FLAC__double)total_error_1 / (FLAC__double)data_len) / M_LN2 : 0.0);
00260     residual_bits_per_sample[2] = (FLAC__float)((total_error_2 > 0) ? log(M_LN2 * (FLAC__double)total_error_2 / (FLAC__double)data_len) / M_LN2 : 0.0);
00261     residual_bits_per_sample[3] = (FLAC__float)((total_error_3 > 0) ? log(M_LN2 * (FLAC__double)total_error_3 / (FLAC__double)data_len) / M_LN2 : 0.0);
00262     residual_bits_per_sample[4] = (FLAC__float)((total_error_4 > 0) ? log(M_LN2 * (FLAC__double)total_error_4 / (FLAC__double)data_len) / M_LN2 : 0.0);
00263 #else
00264     residual_bits_per_sample[0] = (total_error_0 > 0) ? local__compute_rbps_integerized(total_error_0, data_len) : 0;
00265     residual_bits_per_sample[1] = (total_error_1 > 0) ? local__compute_rbps_integerized(total_error_1, data_len) : 0;
00266     residual_bits_per_sample[2] = (total_error_2 > 0) ? local__compute_rbps_integerized(total_error_2, data_len) : 0;
00267     residual_bits_per_sample[3] = (total_error_3 > 0) ? local__compute_rbps_integerized(total_error_3, data_len) : 0;
00268     residual_bits_per_sample[4] = (total_error_4 > 0) ? local__compute_rbps_integerized(total_error_4, data_len) : 0;
00269 #endif
00270 
00271     return order;
00272 }
00273 
00274 #ifndef FLAC__INTEGER_ONLY_LIBRARY
00275 unsigned FLAC__fixed_compute_best_predictor_wide(const FLAC__int32 data[], unsigned data_len, FLAC__float residual_bits_per_sample[FLAC__MAX_FIXED_ORDER+1])
00276 #else
00277 unsigned FLAC__fixed_compute_best_predictor_wide(const FLAC__int32 data[], unsigned data_len, FLAC__fixedpoint residual_bits_per_sample[FLAC__MAX_FIXED_ORDER+1])
00278 #endif
00279 {
00280     FLAC__int32 last_error_0 = data[-1];
00281     FLAC__int32 last_error_1 = data[-1] - data[-2];
00282     FLAC__int32 last_error_2 = last_error_1 - (data[-2] - data[-3]);
00283     FLAC__int32 last_error_3 = last_error_2 - (data[-2] - 2*data[-3] + data[-4]);
00284     FLAC__int32 error, save;
00285     /* total_error_* are 64-bits to avoid overflow when encoding
00286      * erratic signals when the bits-per-sample and blocksize are
00287      * large.
00288      */
00289     FLAC__uint64 total_error_0 = 0, total_error_1 = 0, total_error_2 = 0, total_error_3 = 0, total_error_4 = 0;
00290     unsigned i, order;
00291 
00292     for(i = 0; i < data_len; i++) {
00293         error  = data[i]     ; total_error_0 += local_abs(error);                      save = error;
00294         error -= last_error_0; total_error_1 += local_abs(error); last_error_0 = save; save = error;
00295         error -= last_error_1; total_error_2 += local_abs(error); last_error_1 = save; save = error;
00296         error -= last_error_2; total_error_3 += local_abs(error); last_error_2 = save; save = error;
00297         error -= last_error_3; total_error_4 += local_abs(error); last_error_3 = save;
00298     }
00299 
00300     if(total_error_0 < flac_min(flac_min(flac_min(total_error_1, total_error_2), total_error_3), total_error_4))
00301         order = 0;
00302     else if(total_error_1 < flac_min(flac_min(total_error_2, total_error_3), total_error_4))
00303         order = 1;
00304     else if(total_error_2 < flac_min(total_error_3, total_error_4))
00305         order = 2;
00306     else if(total_error_3 < total_error_4)
00307         order = 3;
00308     else
00309         order = 4;
00310 
00311     /* Estimate the expected number of bits per residual signal sample. */
00312     /* 'total_error*' is linearly related to the variance of the residual */
00313     /* signal, so we use it directly to compute E(|x|) */
00314     FLAC__ASSERT(data_len > 0 || total_error_0 == 0);
00315     FLAC__ASSERT(data_len > 0 || total_error_1 == 0);
00316     FLAC__ASSERT(data_len > 0 || total_error_2 == 0);
00317     FLAC__ASSERT(data_len > 0 || total_error_3 == 0);
00318     FLAC__ASSERT(data_len > 0 || total_error_4 == 0);
00319 #ifndef FLAC__INTEGER_ONLY_LIBRARY
00320     residual_bits_per_sample[0] = (FLAC__float)((total_error_0 > 0) ? log(M_LN2 * (FLAC__double)total_error_0 / (FLAC__double)data_len) / M_LN2 : 0.0);
00321     residual_bits_per_sample[1] = (FLAC__float)((total_error_1 > 0) ? log(M_LN2 * (FLAC__double)total_error_1 / (FLAC__double)data_len) / M_LN2 : 0.0);
00322     residual_bits_per_sample[2] = (FLAC__float)((total_error_2 > 0) ? log(M_LN2 * (FLAC__double)total_error_2 / (FLAC__double)data_len) / M_LN2 : 0.0);
00323     residual_bits_per_sample[3] = (FLAC__float)((total_error_3 > 0) ? log(M_LN2 * (FLAC__double)total_error_3 / (FLAC__double)data_len) / M_LN2 : 0.0);
00324     residual_bits_per_sample[4] = (FLAC__float)((total_error_4 > 0) ? log(M_LN2 * (FLAC__double)total_error_4 / (FLAC__double)data_len) / M_LN2 : 0.0);
00325 #else
00326     residual_bits_per_sample[0] = (total_error_0 > 0) ? local__compute_rbps_wide_integerized(total_error_0, data_len) : 0;
00327     residual_bits_per_sample[1] = (total_error_1 > 0) ? local__compute_rbps_wide_integerized(total_error_1, data_len) : 0;
00328     residual_bits_per_sample[2] = (total_error_2 > 0) ? local__compute_rbps_wide_integerized(total_error_2, data_len) : 0;
00329     residual_bits_per_sample[3] = (total_error_3 > 0) ? local__compute_rbps_wide_integerized(total_error_3, data_len) : 0;
00330     residual_bits_per_sample[4] = (total_error_4 > 0) ? local__compute_rbps_wide_integerized(total_error_4, data_len) : 0;
00331 #endif
00332 
00333     return order;
00334 }
00335 
00336 void FLAC__fixed_compute_residual(const FLAC__int32 data[], unsigned data_len, unsigned order, FLAC__int32 residual[])
00337 {
00338     const int idata_len = (int)data_len;
00339     int i;
00340 
00341     switch(order) {
00342         case 0:
00343             FLAC__ASSERT(sizeof(residual[0]) == sizeof(data[0]));
00344             memcpy(residual, data, sizeof(residual[0])*data_len);
00345             break;
00346         case 1:
00347             for(i = 0; i < idata_len; i++)
00348                 residual[i] = data[i] - data[i-1];
00349             break;
00350         case 2:
00351             for(i = 0; i < idata_len; i++)
00352 #if 1 /* OPT: may be faster with some compilers on some systems */
00353                 residual[i] = data[i] - (data[i-1] << 1) + data[i-2];
00354 #else
00355                 residual[i] = data[i] - 2*data[i-1] + data[i-2];
00356 #endif
00357             break;
00358         case 3:
00359             for(i = 0; i < idata_len; i++)
00360 #if 1 /* OPT: may be faster with some compilers on some systems */
00361                 residual[i] = data[i] - (((data[i-1]-data[i-2])<<1) + (data[i-1]-data[i-2])) - data[i-3];
00362 #else
00363                 residual[i] = data[i] - 3*data[i-1] + 3*data[i-2] - data[i-3];
00364 #endif
00365             break;
00366         case 4:
00367             for(i = 0; i < idata_len; i++)
00368 #if 1 /* OPT: may be faster with some compilers on some systems */
00369                 residual[i] = data[i] - ((data[i-1]+data[i-3])<<2) + ((data[i-2]<<2) + (data[i-2]<<1)) + data[i-4];
00370 #else
00371                 residual[i] = data[i] - 4*data[i-1] + 6*data[i-2] - 4*data[i-3] + data[i-4];
00372 #endif
00373             break;
00374         default:
00375             FLAC__ASSERT(0);
00376     }
00377 }
00378 
00379 void FLAC__fixed_restore_signal(const FLAC__int32 residual[], unsigned data_len, unsigned order, FLAC__int32 data[])
00380 {
00381     int i, idata_len = (int)data_len;
00382 
00383     switch(order) {
00384         case 0:
00385             FLAC__ASSERT(sizeof(residual[0]) == sizeof(data[0]));
00386             memcpy(data, residual, sizeof(residual[0])*data_len);
00387             break;
00388         case 1:
00389             for(i = 0; i < idata_len; i++)
00390                 data[i] = residual[i] + data[i-1];
00391             break;
00392         case 2:
00393             for(i = 0; i < idata_len; i++)
00394 #if 1 /* OPT: may be faster with some compilers on some systems */
00395                 data[i] = residual[i] + (data[i-1]<<1) - data[i-2];
00396 #else
00397                 data[i] = residual[i] + 2*data[i-1] - data[i-2];
00398 #endif
00399             break;
00400         case 3:
00401             for(i = 0; i < idata_len; i++)
00402 #if 1 /* OPT: may be faster with some compilers on some systems */
00403                 data[i] = residual[i] + (((data[i-1]-data[i-2])<<1) + (data[i-1]-data[i-2])) + data[i-3];
00404 #else
00405                 data[i] = residual[i] + 3*data[i-1] - 3*data[i-2] + data[i-3];
00406 #endif
00407             break;
00408         case 4:
00409             for(i = 0; i < idata_len; i++)
00410 #if 1 /* OPT: may be faster with some compilers on some systems */
00411                 data[i] = residual[i] + ((data[i-1]+data[i-3])<<2) - ((data[i-2]<<2) + (data[i-2]<<1)) - data[i-4];
00412 #else
00413                 data[i] = residual[i] + 4*data[i-1] - 6*data[i-2] + 4*data[i-3] - data[i-4];
00414 #endif
00415             break;
00416         default:
00417             FLAC__ASSERT(0);
00418     }
00419 }