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 "share/compat.h" 40 #include "private/bitmath.h" 41 #include "private/fixed.h" 42 #include "private/macros.h" 43 #include "FLAC/assert.h" 44 45 #ifdef local_abs 46 #undef local_abs 47 #endif 48 #define local_abs(x) ((unsigned)((x)<0? -(x) : (x))) 49 50 #ifdef FLAC__INTEGER_ONLY_LIBRARY 51 /* rbps stands for residual bits per sample 52 * 53 * (ln(2) * err) 54 * rbps = log (-----------) 55 * 2 ( n ) 56 */ 57 static FLAC__fixedpoint local__compute_rbps_integerized(FLAC__uint32 err, FLAC__uint32 n) 58 { 59 FLAC__uint32 rbps; 60 unsigned bits; /* the number of bits required to represent a number */ 61 int fracbits; /* the number of bits of rbps that comprise the fractional part */ 62 63 FLAC__ASSERT(sizeof(rbps) == sizeof(FLAC__fixedpoint)); 64 FLAC__ASSERT(err > 0); 65 FLAC__ASSERT(n > 0); 66 67 FLAC__ASSERT(n <= FLAC__MAX_BLOCK_SIZE); 68 if(err <= n) 69 return 0; 70 /* 71 * The above two things tell us 1) n fits in 16 bits; 2) err/n > 1. 72 * These allow us later to know we won't lose too much precision in the 73 * fixed-point division (err<<fracbits)/n. 74 */ 75 76 fracbits = (8*sizeof(err)) - (FLAC__bitmath_ilog2(err)+1); 77 78 err <<= fracbits; 79 err /= n; 80 /* err now holds err/n with fracbits fractional bits */ 81 82 /* 83 * Whittle err down to 16 bits max. 16 significant bits is enough for 84 * our purposes. 85 */ 86 FLAC__ASSERT(err > 0); 87 bits = FLAC__bitmath_ilog2(err)+1; 88 if(bits > 16) { 89 err >>= (bits-16); 90 fracbits -= (bits-16); 91 } 92 rbps = (FLAC__uint32)err; 93 94 /* Multiply by fixed-point version of ln(2), with 16 fractional bits */ 95 rbps *= FLAC__FP_LN2; 96 fracbits += 16; 97 FLAC__ASSERT(fracbits >= 0); 98 99 /* FLAC__fixedpoint_log2 requires fracbits%4 to be 0 */ 100 { 101 const int f = fracbits & 3; 102 if(f) { 103 rbps >>= f; 104 fracbits -= f; 105 } 106 } 107 108 rbps = FLAC__fixedpoint_log2(rbps, fracbits, (unsigned)(-1)); 109 110 if(rbps == 0) 111 return 0; 112 113 /* 114 * The return value must have 16 fractional bits. Since the whole part 115 * of the base-2 log of a 32 bit number must fit in 5 bits, and fracbits 116 * must be >= -3, these assertion allows us to be able to shift rbps 117 * left if necessary to get 16 fracbits without losing any bits of the 118 * whole part of rbps. 119 * 120 * There is a slight chance due to accumulated error that the whole part 121 * will require 6 bits, so we use 6 in the assertion. Really though as 122 * long as it fits in 13 bits (32 - (16 - (-3))) we are fine. 123 */ 124 FLAC__ASSERT((int)FLAC__bitmath_ilog2(rbps)+1 <= fracbits + 6); 125 FLAC__ASSERT(fracbits >= -3); 126 127 /* now shift the decimal point into place */ 128 if(fracbits < 16) 129 return rbps << (16-fracbits); 130 else if(fracbits > 16) 131 return rbps >> (fracbits-16); 132 else 133 return rbps; 134 } 135 136 static FLAC__fixedpoint local__compute_rbps_wide_integerized(FLAC__uint64 err, FLAC__uint32 n) 137 { 138 FLAC__uint32 rbps; 139 unsigned bits; /* the number of bits required to represent a number */ 140 int fracbits; /* the number of bits of rbps that comprise the fractional part */ 141 142 FLAC__ASSERT(sizeof(rbps) == sizeof(FLAC__fixedpoint)); 143 FLAC__ASSERT(err > 0); 144 FLAC__ASSERT(n > 0); 145 146 FLAC__ASSERT(n <= FLAC__MAX_BLOCK_SIZE); 147 if(err <= n) 148 return 0; 149 /* 150 * The above two things tell us 1) n fits in 16 bits; 2) err/n > 1. 151 * These allow us later to know we won't lose too much precision in the 152 * fixed-point division (err<<fracbits)/n. 153 */ 154 155 fracbits = (8*sizeof(err)) - (FLAC__bitmath_ilog2_wide(err)+1); 156 157 err <<= fracbits; 158 err /= n; 159 /* err now holds err/n with fracbits fractional bits */ 160 161 /* 162 * Whittle err down to 16 bits max. 16 significant bits is enough for 163 * our purposes. 164 */ 165 FLAC__ASSERT(err > 0); 166 bits = FLAC__bitmath_ilog2_wide(err)+1; 167 if(bits > 16) { 168 err >>= (bits-16); 169 fracbits -= (bits-16); 170 } 171 rbps = (FLAC__uint32)err; 172 173 /* Multiply by fixed-point version of ln(2), with 16 fractional bits */ 174 rbps *= FLAC__FP_LN2; 175 fracbits += 16; 176 FLAC__ASSERT(fracbits >= 0); 177 178 /* FLAC__fixedpoint_log2 requires fracbits%4 to be 0 */ 179 { 180 const int f = fracbits & 3; 181 if(f) { 182 rbps >>= f; 183 fracbits -= f; 184 } 185 } 186 187 rbps = FLAC__fixedpoint_log2(rbps, fracbits, (unsigned)(-1)); 188 189 if(rbps == 0) 190 return 0; 191 192 /* 193 * The return value must have 16 fractional bits. Since the whole part 194 * of the base-2 log of a 32 bit number must fit in 5 bits, and fracbits 195 * must be >= -3, these assertion allows us to be able to shift rbps 196 * left if necessary to get 16 fracbits without losing any bits of the 197 * whole part of rbps. 198 * 199 * There is a slight chance due to accumulated error that the whole part 200 * will require 6 bits, so we use 6 in the assertion. Really though as 201 * long as it fits in 13 bits (32 - (16 - (-3))) we are fine. 202 */ 203 FLAC__ASSERT((int)FLAC__bitmath_ilog2(rbps)+1 <= fracbits + 6); 204 FLAC__ASSERT(fracbits >= -3); 205 206 /* now shift the decimal point into place */ 207 if(fracbits < 16) 208 return rbps << (16-fracbits); 209 else if(fracbits > 16) 210 return rbps >> (fracbits-16); 211 else 212 return rbps; 213 } 214 #endif 215 216 #ifndef FLAC__INTEGER_ONLY_LIBRARY 217 unsigned FLAC__fixed_compute_best_predictor(const FLAC__int32 data[], unsigned data_len, FLAC__float residual_bits_per_sample[FLAC__MAX_FIXED_ORDER+1]) 218 #else 219 unsigned FLAC__fixed_compute_best_predictor(const FLAC__int32 data[], unsigned data_len, FLAC__fixedpoint residual_bits_per_sample[FLAC__MAX_FIXED_ORDER+1]) 220 #endif 221 { 222 FLAC__int32 last_error_0 = data[-1]; 223 FLAC__int32 last_error_1 = data[-1] - data[-2]; 224 FLAC__int32 last_error_2 = last_error_1 - (data[-2] - data[-3]); 225 FLAC__int32 last_error_3 = last_error_2 - (data[-2] - 2*data[-3] + data[-4]); 226 FLAC__int32 error, save; 227 FLAC__uint32 total_error_0 = 0, total_error_1 = 0, total_error_2 = 0, total_error_3 = 0, total_error_4 = 0; 228 unsigned i, order; 229 230 for(i = 0; i < data_len; i++) { 231 error = data[i] ; total_error_0 += local_abs(error); save = error; 232 error -= last_error_0; total_error_1 += local_abs(error); last_error_0 = save; save = error; 233 error -= last_error_1; total_error_2 += local_abs(error); last_error_1 = save; save = error; 234 error -= last_error_2; total_error_3 += local_abs(error); last_error_2 = save; save = error; 235 error -= last_error_3; total_error_4 += local_abs(error); last_error_3 = save; 236 } 237 238 if(total_error_0 < flac_min(flac_min(flac_min(total_error_1, total_error_2), total_error_3), total_error_4)) 239 order = 0; 240 else if(total_error_1 < flac_min(flac_min(total_error_2, total_error_3), total_error_4)) 241 order = 1; 242 else if(total_error_2 < flac_min(total_error_3, total_error_4)) 243 order = 2; 244 else if(total_error_3 < total_error_4) 245 order = 3; 246 else 247 order = 4; 248 249 /* Estimate the expected number of bits per residual signal sample. */ 250 /* 'total_error*' is linearly related to the variance of the residual */ 251 /* signal, so we use it directly to compute E(|x|) */ 252 FLAC__ASSERT(data_len > 0 || total_error_0 == 0); 253 FLAC__ASSERT(data_len > 0 || total_error_1 == 0); 254 FLAC__ASSERT(data_len > 0 || total_error_2 == 0); 255 FLAC__ASSERT(data_len > 0 || total_error_3 == 0); 256 FLAC__ASSERT(data_len > 0 || total_error_4 == 0); 257 #ifndef FLAC__INTEGER_ONLY_LIBRARY 258 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); 259 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); 260 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); 261 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); 262 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); 263 #else 264 residual_bits_per_sample[0] = (total_error_0 > 0) ? local__compute_rbps_integerized(total_error_0, data_len) : 0; 265 residual_bits_per_sample[1] = (total_error_1 > 0) ? local__compute_rbps_integerized(total_error_1, data_len) : 0; 266 residual_bits_per_sample[2] = (total_error_2 > 0) ? local__compute_rbps_integerized(total_error_2, data_len) : 0; 267 residual_bits_per_sample[3] = (total_error_3 > 0) ? local__compute_rbps_integerized(total_error_3, data_len) : 0; 268 residual_bits_per_sample[4] = (total_error_4 > 0) ? local__compute_rbps_integerized(total_error_4, data_len) : 0; 269 #endif 270 271 return order; 272 } 273 274 #ifndef FLAC__INTEGER_ONLY_LIBRARY 275 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]) 276 #else 277 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]) 278 #endif 279 { 280 FLAC__int32 last_error_0 = data[-1]; 281 FLAC__int32 last_error_1 = data[-1] - data[-2]; 282 FLAC__int32 last_error_2 = last_error_1 - (data[-2] - data[-3]); 283 FLAC__int32 last_error_3 = last_error_2 - (data[-2] - 2*data[-3] + data[-4]); 284 FLAC__int32 error, save; 285 /* total_error_* are 64-bits to avoid overflow when encoding 286 * erratic signals when the bits-per-sample and blocksize are 287 * large. 288 */ 289 FLAC__uint64 total_error_0 = 0, total_error_1 = 0, total_error_2 = 0, total_error_3 = 0, total_error_4 = 0; 290 unsigned i, order; 291 292 for(i = 0; i < data_len; i++) { 293 error = data[i] ; total_error_0 += local_abs(error); save = error; 294 error -= last_error_0; total_error_1 += local_abs(error); last_error_0 = save; save = error; 295 error -= last_error_1; total_error_2 += local_abs(error); last_error_1 = save; save = error; 296 error -= last_error_2; total_error_3 += local_abs(error); last_error_2 = save; save = error; 297 error -= last_error_3; total_error_4 += local_abs(error); last_error_3 = save; 298 } 299 300 if(total_error_0 < flac_min(flac_min(flac_min(total_error_1, total_error_2), total_error_3), total_error_4)) 301 order = 0; 302 else if(total_error_1 < flac_min(flac_min(total_error_2, total_error_3), total_error_4)) 303 order = 1; 304 else if(total_error_2 < flac_min(total_error_3, total_error_4)) 305 order = 2; 306 else if(total_error_3 < total_error_4) 307 order = 3; 308 else 309 order = 4; 310 311 /* Estimate the expected number of bits per residual signal sample. */ 312 /* 'total_error*' is linearly related to the variance of the residual */ 313 /* signal, so we use it directly to compute E(|x|) */ 314 FLAC__ASSERT(data_len > 0 || total_error_0 == 0); 315 FLAC__ASSERT(data_len > 0 || total_error_1 == 0); 316 FLAC__ASSERT(data_len > 0 || total_error_2 == 0); 317 FLAC__ASSERT(data_len > 0 || total_error_3 == 0); 318 FLAC__ASSERT(data_len > 0 || total_error_4 == 0); 319 #ifndef FLAC__INTEGER_ONLY_LIBRARY 320 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); 321 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); 322 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); 323 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); 324 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); 325 #else 326 residual_bits_per_sample[0] = (total_error_0 > 0) ? local__compute_rbps_wide_integerized(total_error_0, data_len) : 0; 327 residual_bits_per_sample[1] = (total_error_1 > 0) ? local__compute_rbps_wide_integerized(total_error_1, data_len) : 0; 328 residual_bits_per_sample[2] = (total_error_2 > 0) ? local__compute_rbps_wide_integerized(total_error_2, data_len) : 0; 329 residual_bits_per_sample[3] = (total_error_3 > 0) ? local__compute_rbps_wide_integerized(total_error_3, data_len) : 0; 330 residual_bits_per_sample[4] = (total_error_4 > 0) ? local__compute_rbps_wide_integerized(total_error_4, data_len) : 0; 331 #endif 332 333 return order; 334 } 335 336 void FLAC__fixed_compute_residual(const FLAC__int32 data[], unsigned data_len, unsigned order, FLAC__int32 residual[]) 337 { 338 const int idata_len = (int)data_len; 339 int i; 340 341 switch(order) { 342 case 0: 343 FLAC__ASSERT(sizeof(residual[0]) == sizeof(data[0])); 344 memcpy(residual, data, sizeof(residual[0])*data_len); 345 break; 346 case 1: 347 for(i = 0; i < idata_len; i++) 348 residual[i] = data[i] - data[i-1]; 349 break; 350 case 2: 351 for(i = 0; i < idata_len; i++) 352 #if 1 /* OPT: may be faster with some compilers on some systems */ 353 residual[i] = data[i] - (data[i-1] << 1) + data[i-2]; 354 #else 355 residual[i] = data[i] - 2*data[i-1] + data[i-2]; 356 #endif 357 break; 358 case 3: 359 for(i = 0; i < idata_len; i++) 360 #if 1 /* OPT: may be faster with some compilers on some systems */ 361 residual[i] = data[i] - (((data[i-1]-data[i-2])<<1) + (data[i-1]-data[i-2])) - data[i-3]; 362 #else 363 residual[i] = data[i] - 3*data[i-1] + 3*data[i-2] - data[i-3]; 364 #endif 365 break; 366 case 4: 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]+data[i-3])<<2) + ((data[i-2]<<2) + (data[i-2]<<1)) + data[i-4]; 370 #else 371 residual[i] = data[i] - 4*data[i-1] + 6*data[i-2] - 4*data[i-3] + data[i-4]; 372 #endif 373 break; 374 default: 375 FLAC__ASSERT(0); 376 } 377 } 378 379 void FLAC__fixed_restore_signal(const FLAC__int32 residual[], unsigned data_len, unsigned order, FLAC__int32 data[]) 380 { 381 int i, idata_len = (int)data_len; 382 383 switch(order) { 384 case 0: 385 FLAC__ASSERT(sizeof(residual[0]) == sizeof(data[0])); 386 memcpy(data, residual, sizeof(residual[0])*data_len); 387 break; 388 case 1: 389 for(i = 0; i < idata_len; i++) 390 data[i] = residual[i] + data[i-1]; 391 break; 392 case 2: 393 for(i = 0; i < idata_len; i++) 394 #if 1 /* OPT: may be faster with some compilers on some systems */ 395 data[i] = residual[i] + (data[i-1]<<1) - data[i-2]; 396 #else 397 data[i] = residual[i] + 2*data[i-1] - data[i-2]; 398 #endif 399 break; 400 case 3: 401 for(i = 0; i < idata_len; i++) 402 #if 1 /* OPT: may be faster with some compilers on some systems */ 403 data[i] = residual[i] + (((data[i-1]-data[i-2])<<1) + (data[i-1]-data[i-2])) + data[i-3]; 404 #else 405 data[i] = residual[i] + 3*data[i-1] - 3*data[i-2] + data[i-3]; 406 #endif 407 break; 408 case 4: 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]+data[i-3])<<2) - ((data[i-2]<<2) + (data[i-2]<<1)) - data[i-4]; 412 #else 413 data[i] = residual[i] + 4*data[i-1] - 6*data[i-2] + 4*data[i-3] - data[i-4]; 414 #endif 415 break; 416 default: 417 FLAC__ASSERT(0); 418 } 419 } 420
