| 1 | /* libFLAC - Free Lossless Audio Codec library |
| 2 | * Copyright (C) 2000-2009 Josh Coalson |
| 3 | * Copyright (C) 2011-2016 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, 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] = (float)((total_error_0 > 0) ? log(M_LN2 * (double)total_error_0 / (double)data_len) / M_LN2 : 0.0); |
| 259 | residual_bits_per_sample[1] = (float)((total_error_1 > 0) ? log(M_LN2 * (double)total_error_1 / (double)data_len) / M_LN2 : 0.0); |
| 260 | residual_bits_per_sample[2] = (float)((total_error_2 > 0) ? log(M_LN2 * (double)total_error_2 / (double)data_len) / M_LN2 : 0.0); |
| 261 | residual_bits_per_sample[3] = (float)((total_error_3 > 0) ? log(M_LN2 * (double)total_error_3 / (double)data_len) / M_LN2 : 0.0); |
| 262 | residual_bits_per_sample[4] = (float)((total_error_4 > 0) ? log(M_LN2 * (double)total_error_4 / (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, 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] = (float)((total_error_0 > 0) ? log(M_LN2 * (double)total_error_0 / (double)data_len) / M_LN2 : 0.0); |
| 321 | residual_bits_per_sample[1] = (float)((total_error_1 > 0) ? log(M_LN2 * (double)total_error_1 / (double)data_len) / M_LN2 : 0.0); |
| 322 | residual_bits_per_sample[2] = (float)((total_error_2 > 0) ? log(M_LN2 * (double)total_error_2 / (double)data_len) / M_LN2 : 0.0); |
| 323 | residual_bits_per_sample[3] = (float)((total_error_3 > 0) ? log(M_LN2 * (double)total_error_3 / (double)data_len) / M_LN2 : 0.0); |
| 324 | residual_bits_per_sample[4] = (float)((total_error_4 > 0) ? log(M_LN2 * (double)total_error_4 / (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 | residual[i] = data[i] - 2*data[i-1] + data[i-2]; |
| 353 | break; |
| 354 | case 3: |
| 355 | for(i = 0; i < idata_len; i++) |
| 356 | residual[i] = data[i] - 3*data[i-1] + 3*data[i-2] - data[i-3]; |
| 357 | break; |
| 358 | case 4: |
| 359 | for(i = 0; i < idata_len; i++) |
| 360 | residual[i] = data[i] - 4*data[i-1] + 6*data[i-2] - 4*data[i-3] + data[i-4]; |
| 361 | break; |
| 362 | default: |
| 363 | FLAC__ASSERT(0); |
| 364 | } |
| 365 | } |
| 366 | |
| 367 | void FLAC__fixed_restore_signal(const FLAC__int32 residual[], unsigned data_len, unsigned order, FLAC__int32 data[]) |
| 368 | { |
| 369 | int i, idata_len = (int)data_len; |
| 370 | |
| 371 | switch(order) { |
| 372 | case 0: |
| 373 | FLAC__ASSERT(sizeof(residual[0]) == sizeof(data[0])); |
| 374 | memcpy(data, residual, sizeof(residual[0])*data_len); |
| 375 | break; |
| 376 | case 1: |
| 377 | for(i = 0; i < idata_len; i++) |
| 378 | data[i] = residual[i] + data[i-1]; |
| 379 | break; |
| 380 | case 2: |
| 381 | for(i = 0; i < idata_len; i++) |
| 382 | data[i] = residual[i] + 2*data[i-1] - data[i-2]; |
| 383 | break; |
| 384 | case 3: |
| 385 | for(i = 0; i < idata_len; i++) |
| 386 | data[i] = residual[i] + 3*data[i-1] - 3*data[i-2] + data[i-3]; |
| 387 | break; |
| 388 | case 4: |
| 389 | for(i = 0; i < idata_len; i++) |
| 390 | data[i] = residual[i] + 4*data[i-1] - 6*data[i-2] + 4*data[i-3] - data[i-4]; |
| 391 | break; |
| 392 | default: |
| 393 | FLAC__ASSERT(0); |
| 394 | } |
| 395 | } |