source/modules/juce_audio_formats/codecs/flac/libFLAC/fixed.c

   1 /* libFLAC - Free Lossless Audio Codec library
   2  * Copyright (C) 2000-2009  Josh Coalson
   3  * Copyright (C) 2011-2014  Xiph.Org Foundation
   4  *
   5  * Redistribution and use in source and binary forms, with or without
   6  * modification, are permitted provided that the following conditions
   7  * are met:
   8  *
   9  * - Redistributions of source code must retain the above copyright
  10  * notice, this list of conditions and the following disclaimer.
  11  *
  12  * - Redistributions in binary form must reproduce the above copyright
  13  * notice, this list of conditions and the following disclaimer in the
  14  * documentation and/or other materials provided with the distribution.
  15  *
  16  * - Neither the name of the Xiph.org Foundation nor the names of its
  17  * contributors may be used to endorse or promote products derived from
  18  * this software without specific prior written permission.
  19  *
  20  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  21  * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  22  * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  23  * A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
  24  * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
  25  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
  26  * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
  27  * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
  28  * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
  29  * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  30  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  31  */
  32
  33 #ifdef HAVE_CONFIG_H
  34 #  include <config.h>
  35 #endif
  36
  37 #include <math.h>
  38 #include <string.h>
  39 #include "../compat.h"
  40 #include "include/private/bitmath.h"
  41 #include "include/private/fixed.h"
  42 #include "../assert.h"
  43
  44 #ifdef local_abs
  45 #undef local_abs
  46 #endif
  47 #define local_abs(x) ((unsigned)((x)<0? -(x) : (x)))
  48
  49 #ifdef FLAC__INTEGER_ONLY_LIBRARY
  50 /* rbps stands for residual bits per sample
  51  *
  52  *             (ln(2) * err)
  53  * rbps = log  (-----------)
  54  *           2 (     n     )
  55  */
  56 static FLAC__fixedpoint local__compute_rbps_integerized(FLAC__uint32 err, FLAC__uint32 n)
  57 {
  58         FLAC__uint32 rbps;
  59         unsigned bits; /* the number of bits required to represent a number */
  60         int fracbits; /* the number of bits of rbps that comprise the fractional part */
  61
  62         FLAC__ASSERT(sizeof(rbps) == sizeof(FLAC__fixedpoint));
  63         FLAC__ASSERT(err > 0);
  64         FLAC__ASSERT(n > 0);
  65
  66         FLAC__ASSERT(n <= FLAC__MAX_BLOCK_SIZE);
  67         if(err <= n)
  68                 return 0;
  69         /*
  70          * The above two things tell us 1) n fits in 16 bits; 2) err/n > 1.
  71          * These allow us later to know we won't lose too much precision in the
  72          * fixed-point division (err<<fracbits)/n.
  73          */
  74
  75         fracbits = (8*sizeof(err)) - (FLAC__bitmath_ilog2(err)+1);
  76
  77         err <<= fracbits;
  78         err /= n;
  79         /* err now holds err/n with fracbits fractional bits */
  80
  81         /*
  82          * Whittle err down to 16 bits max.  16 significant bits is enough for
  83          * our purposes.
  84          */
  85         FLAC__ASSERT(err > 0);
  86         bits = FLAC__bitmath_ilog2(err)+1;
  87         if(bits > 16) {
  88                 err >>= (bits-16);
  89                 fracbits -= (bits-16);
  90         }
  91         rbps = (FLAC__uint32)err;
  92
  93         /* Multiply by fixed-point version of ln(2), with 16 fractional bits */
  94         rbps *= FLAC__FP_LN2;
  95         fracbits += 16;
  96         FLAC__ASSERT(fracbits >= 0);
  97
  98         /* FLAC__fixedpoint_log2 requires fracbits%4 to be 0 */
  99         {
 100                 const int f = fracbits & 3;
 101                 if(f) {
 102                         rbps >>= f;
 103                         fracbits -= f;
 104                 }
 105         }
 106
 107         rbps = FLAC__fixedpoint_log2(rbps, fracbits, (unsigned)(-1));
 108
 109         if(rbps == 0)
 110                 return 0;
 111
 112         /*
 113          * The return value must have 16 fractional bits.  Since the whole part
 114          * of the base-2 log of a 32 bit number must fit in 5 bits, and fracbits
 115          * must be >= -3, these assertion allows us to be able to shift rbps
 116          * left if necessary to get 16 fracbits without losing any bits of the
 117          * whole part of rbps.
 118          *
 119          * There is a slight chance due to accumulated error that the whole part
 120          * will require 6 bits, so we use 6 in the assertion.  Really though as
 121          * long as it fits in 13 bits (32 - (16 - (-3))) we are fine.
 122          */
 123         FLAC__ASSERT((int)FLAC__bitmath_ilog2(rbps)+1 <= fracbits + 6);
 124         FLAC__ASSERT(fracbits >= -3);
 125
 126         /* now shift the decimal point into place */
 127         if(fracbits < 16)
 128                 return rbps << (16-fracbits);
 129         else if(fracbits > 16)
 130                 return rbps >> (fracbits-16);
 131         else
 132                 return rbps;
 133 }
 134
 135 static FLAC__fixedpoint local__compute_rbps_wide_integerized(FLAC__uint64 err, FLAC__uint32 n)
 136 {
 137         FLAC__uint32 rbps;
 138         unsigned bits; /* the number of bits required to represent a number */
 139         int fracbits; /* the number of bits of rbps that comprise the fractional part */
 140
 141         FLAC__ASSERT(sizeof(rbps) == sizeof(FLAC__fixedpoint));
 142         FLAC__ASSERT(err > 0);
 143         FLAC__ASSERT(n > 0);
 144
 145         FLAC__ASSERT(n <= FLAC__MAX_BLOCK_SIZE);
 146         if(err <= n)
 147                 return 0;
 148         /*
 149          * The above two things tell us 1) n fits in 16 bits; 2) err/n > 1.
 150          * These allow us later to know we won't lose too much precision in the
 151          * fixed-point division (err<<fracbits)/n.
 152          */
 153
 154         fracbits = (8*sizeof(err)) - (FLAC__bitmath_ilog2_wide(err)+1);
 155
 156         err <<= fracbits;
 157         err /= n;
 158         /* err now holds err/n with fracbits fractional bits */
 159
 160         /*
 161          * Whittle err down to 16 bits max.  16 significant bits is enough for
 162          * our purposes.
 163          */
 164         FLAC__ASSERT(err > 0);
 165         bits = FLAC__bitmath_ilog2_wide(err)+1;
 166         if(bits > 16) {
 167                 err >>= (bits-16);
 168                 fracbits -= (bits-16);
 169         }
 170         rbps = (FLAC__uint32)err;
 171
 172         /* Multiply by fixed-point version of ln(2), with 16 fractional bits */
 173         rbps *= FLAC__FP_LN2;
 174         fracbits += 16;
 175         FLAC__ASSERT(fracbits >= 0);
 176
 177         /* FLAC__fixedpoint_log2 requires fracbits%4 to be 0 */
 178         {
 179                 const int f = fracbits & 3;
 180                 if(f) {
 181                         rbps >>= f;
 182                         fracbits -= f;
 183                 }
 184         }
 185
 186         rbps = FLAC__fixedpoint_log2(rbps, fracbits, (unsigned)(-1));
 187
 188         if(rbps == 0)
 189                 return 0;
 190
 191         /*
 192          * The return value must have 16 fractional bits.  Since the whole part
 193          * of the base-2 log of a 32 bit number must fit in 5 bits, and fracbits
 194          * must be >= -3, these assertion allows us to be able to shift rbps
 195          * left if necessary to get 16 fracbits without losing any bits of the
 196          * whole part of rbps.
 197          *
 198          * There is a slight chance due to accumulated error that the whole part
 199          * will require 6 bits, so we use 6 in the assertion.  Really though as
 200          * long as it fits in 13 bits (32 - (16 - (-3))) we are fine.
 201          */
 202         FLAC__ASSERT((int)FLAC__bitmath_ilog2(rbps)+1 <= fracbits + 6);
 203         FLAC__ASSERT(fracbits >= -3);
 204
 205         /* now shift the decimal point into place */
 206         if(fracbits < 16)
 207                 return rbps << (16-fracbits);
 208         else if(fracbits > 16)
 209                 return rbps >> (fracbits-16);
 210         else
 211                 return rbps;
 212 }
 213 #endif
 214
 215 #ifndef FLAC__INTEGER_ONLY_LIBRARY
 216 unsigned FLAC__fixed_compute_best_predictor(const FLAC__int32 data[], unsigned data_len, FLAC__float residual_bits_per_sample[FLAC__MAX_FIXED_ORDER+1])
 217 #else
 218 unsigned FLAC__fixed_compute_best_predictor(const FLAC__int32 data[], unsigned data_len, FLAC__fixedpoint residual_bits_per_sample[FLAC__MAX_FIXED_ORDER+1])
 219 #endif
 220 {
 221         FLAC__int32 last_error_0 = data[-1];
 222         FLAC__int32 last_error_1 = data[-1] - data[-2];
 223         FLAC__int32 last_error_2 = last_error_1 - (data[-2] - data[-3]);
 224         FLAC__int32 last_error_3 = last_error_2 - (data[-2] - 2*data[-3] + data[-4]);
 225         FLAC__int32 error, save;
 226         FLAC__uint32 total_error_0 = 0, total_error_1 = 0, total_error_2 = 0, total_error_3 = 0, total_error_4 = 0;
 227         unsigned i, order;
 228
 229         for(i = 0; i < data_len; i++) {
 230                 error  = data[i]     ; total_error_0 += local_abs(error);                      save = error;
 231                 error -= last_error_0; total_error_1 += local_abs(error); last_error_0 = save; save = error;
 232                 error -= last_error_1; total_error_2 += local_abs(error); last_error_1 = save; save = error;
 233                 error -= last_error_2; total_error_3 += local_abs(error); last_error_2 = save; save = error;
 234                 error -= last_error_3; total_error_4 += local_abs(error); last_error_3 = save;
 235         }
 236
 237         if(total_error_0 < flac_min(flac_min(flac_min(total_error_1, total_error_2), total_error_3), total_error_4))
 238                 order = 0;
 239         else if(total_error_1 < flac_min(flac_min(total_error_2, total_error_3), total_error_4))
 240                 order = 1;
 241         else if(total_error_2 < flac_min(total_error_3, total_error_4))
 242                 order = 2;
 243         else if(total_error_3 < total_error_4)
 244                 order = 3;
 245         else
 246                 order = 4;
 247
 248         /* Estimate the expected number of bits per residual signal sample. */
 249         /* 'total_error*' is linearly related to the variance of the residual */
 250         /* signal, so we use it directly to compute E(|x|) */
 251         FLAC__ASSERT(data_len > 0 || total_error_0 == 0);
 252         FLAC__ASSERT(data_len > 0 || total_error_1 == 0);
 253         FLAC__ASSERT(data_len > 0 || total_error_2 == 0);
 254         FLAC__ASSERT(data_len > 0 || total_error_3 == 0);
 255         FLAC__ASSERT(data_len > 0 || total_error_4 == 0);
 256 #ifndef FLAC__INTEGER_ONLY_LIBRARY
 257         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);
 258         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);
 259         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);
 260         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);
 261         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);
 262 #else
 263         residual_bits_per_sample[0] = (total_error_0 > 0) ? local__compute_rbps_integerized(total_error_0, data_len) : 0;
 264         residual_bits_per_sample[1] = (total_error_1 > 0) ? local__compute_rbps_integerized(total_error_1, data_len) : 0;
 265         residual_bits_per_sample[2] = (total_error_2 > 0) ? local__compute_rbps_integerized(total_error_2, data_len) : 0;
 266         residual_bits_per_sample[3] = (total_error_3 > 0) ? local__compute_rbps_integerized(total_error_3, data_len) : 0;
 267         residual_bits_per_sample[4] = (total_error_4 > 0) ? local__compute_rbps_integerized(total_error_4, data_len) : 0;
 268 #endif
 269
 270         return order;
 271 }
 272
 273 #ifndef FLAC__INTEGER_ONLY_LIBRARY
 274 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])
 275 #else
 276 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])
 277 #endif
 278 {
 279         FLAC__int32 last_error_0 = data[-1];
 280         FLAC__int32 last_error_1 = data[-1] - data[-2];
 281         FLAC__int32 last_error_2 = last_error_1 - (data[-2] - data[-3]);
 282         FLAC__int32 last_error_3 = last_error_2 - (data[-2] - 2*data[-3] + data[-4]);
 283         FLAC__int32 error, save;
 284         /* total_error_* are 64-bits to avoid overflow when encoding
 285          * erratic signals when the bits-per-sample and blocksize are
 286          * large.
 287          */
 288         FLAC__uint64 total_error_0 = 0, total_error_1 = 0, total_error_2 = 0, total_error_3 = 0, total_error_4 = 0;
 289         unsigned i, order;
 290
 291         for(i = 0; i < data_len; i++) {
 292                 error  = data[i]     ; total_error_0 += local_abs(error);                      save = error;
 293                 error -= last_error_0; total_error_1 += local_abs(error); last_error_0 = save; save = error;
 294                 error -= last_error_1; total_error_2 += local_abs(error); last_error_1 = save; save = error;
 295                 error -= last_error_2; total_error_3 += local_abs(error); last_error_2 = save; save = error;
 296                 error -= last_error_3; total_error_4 += local_abs(error); last_error_3 = save;
 297         }
 298
 299         if(total_error_0 < flac_min(flac_min(flac_min(total_error_1, total_error_2), total_error_3), total_error_4))
 300                 order = 0;
 301         else if(total_error_1 < flac_min(flac_min(total_error_2, total_error_3), total_error_4))
 302                 order = 1;
 303         else if(total_error_2 < flac_min(total_error_3, total_error_4))
 304                 order = 2;
 305         else if(total_error_3 < total_error_4)
 306                 order = 3;
 307         else
 308                 order = 4;
 309
 310         /* Estimate the expected number of bits per residual signal sample. */
 311         /* 'total_error*' is linearly related to the variance of the residual */
 312         /* signal, so we use it directly to compute E(|x|) */
 313         FLAC__ASSERT(data_len > 0 || total_error_0 == 0);
 314         FLAC__ASSERT(data_len > 0 || total_error_1 == 0);
 315         FLAC__ASSERT(data_len > 0 || total_error_2 == 0);
 316         FLAC__ASSERT(data_len > 0 || total_error_3 == 0);
 317         FLAC__ASSERT(data_len > 0 || total_error_4 == 0);
 318 #ifndef FLAC__INTEGER_ONLY_LIBRARY
 319         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);
 320         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);
 321         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);
 322         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);
 323         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);
 324 #else
 325         residual_bits_per_sample[0] = (total_error_0 > 0) ? local__compute_rbps_wide_integerized(total_error_0, data_len) : 0;
 326         residual_bits_per_sample[1] = (total_error_1 > 0) ? local__compute_rbps_wide_integerized(total_error_1, data_len) : 0;
 327         residual_bits_per_sample[2] = (total_error_2 > 0) ? local__compute_rbps_wide_integerized(total_error_2, data_len) : 0;
 328         residual_bits_per_sample[3] = (total_error_3 > 0) ? local__compute_rbps_wide_integerized(total_error_3, data_len) : 0;
 329         residual_bits_per_sample[4] = (total_error_4 > 0) ? local__compute_rbps_wide_integerized(total_error_4, data_len) : 0;
 330 #endif
 331
 332         return order;
 333 }
 334
 335 void FLAC__fixed_compute_residual(const FLAC__int32 data[], unsigned data_len, unsigned order, FLAC__int32 residual[])
 336 {
 337         const int idata_len = (int)data_len;
 338         int i;
 339
 340         switch(order) {
 341                 case 0:
 342                         FLAC__ASSERT(sizeof(residual[0]) == sizeof(data[0]));
 343                         memcpy(residual, data, sizeof(residual[0])*data_len);
 344                         break;
 345                 case 1:
 346                         for(i = 0; i < idata_len; i++)
 347                                 residual[i] = data[i] - data[i-1];
 348                         break;
 349                 case 2:
 350                         for(i = 0; i < idata_len; i++)
 351 #if 1 /* OPT: may be faster with some compilers on some systems */
 352                                 residual[i] = data[i] - (data[i-1] << 1) + data[i-2];
 353 #else
 354                                 residual[i] = data[i] - 2*data[i-1] + data[i-2];
 355 #endif
 356                         break;
 357                 case 3:
 358                         for(i = 0; i < idata_len; i++)
 359 #if 1 /* OPT: may be faster with some compilers on some systems */
 360                                 residual[i] = data[i] - (((data[i-1]-data[i-2])<<1) + (data[i-1]-data[i-2])) - data[i-3];
 361 #else
 362                                 residual[i] = data[i] - 3*data[i-1] + 3*data[i-2] - data[i-3];
 363 #endif
 364                         break;
 365                 case 4:
 366                         for(i = 0; i < idata_len; i++)
 367 #if 1 /* OPT: may be faster with some compilers on some systems */
 368                                 residual[i] = data[i] - ((data[i-1]+data[i-3])<<2) + ((data[i-2]<<2) + (data[i-2]<<1)) + data[i-4];
 369 #else
 370                                 residual[i] = data[i] - 4*data[i-1] + 6*data[i-2] - 4*data[i-3] + data[i-4];
 371 #endif
 372                         break;
 373                 default:
 374                         FLAC__ASSERT(0);
 375         }
 376 }
 377
 378 void FLAC__fixed_restore_signal(const FLAC__int32 residual[], unsigned data_len, unsigned order, FLAC__int32 data[])
 379 {
 380         int i, idata_len = (int)data_len;
 381
 382         switch(order) {
 383                 case 0:
 384                         FLAC__ASSERT(sizeof(residual[0]) == sizeof(data[0]));
 385                         memcpy(data, residual, sizeof(residual[0])*data_len);
 386                         break;
 387                 case 1:
 388                         for(i = 0; i < idata_len; i++)
 389                                 data[i] = residual[i] + data[i-1];
 390                         break;
 391                 case 2:
 392                         for(i = 0; i < idata_len; i++)
 393 #if 1 /* OPT: may be faster with some compilers on some systems */
 394                                 data[i] = residual[i] + (data[i-1]<<1) - data[i-2];
 395 #else
 396                                 data[i] = residual[i] + 2*data[i-1] - data[i-2];
 397 #endif
 398                         break;
 399                 case 3:
 400                         for(i = 0; i < idata_len; i++)
 401 #if 1 /* OPT: may be faster with some compilers on some systems */
 402                                 data[i] = residual[i] + (((data[i-1]-data[i-2])<<1) + (data[i-1]-data[i-2])) + data[i-3];
 403 #else
 404                                 data[i] = residual[i] + 3*data[i-1] - 3*data[i-2] + data[i-3];
 405 #endif
 406                         break;
 407                 case 4:
 408                         for(i = 0; i < idata_len; i++)
 409 #if 1 /* OPT: may be faster with some compilers on some systems */
 410                                 data[i] = residual[i] + ((data[i-1]+data[i-3])<<2) - ((data[i-2]<<2) + (data[i-2]<<1)) - data[i-4];
 411 #else
 412                                 data[i] = residual[i] + 4*data[i-1] - 6*data[i-2] + 4*data[i-3] - data[i-4];
 413 #endif
 414                         break;
 415                 default:
 416                         FLAC__ASSERT(0);
 417         }
 418 }