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