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 "FLAC/assert.h" 38#include "FLAC/format.h" 39#include "private/bitmath.h" 40#include "private/lpc.h" 41#if defined DEBUG || defined FLAC__OVERFLOW_DETECT || defined FLAC__OVERFLOW_DETECT_VERBOSE 42#include <stdio.h> 43#endif 44 45#ifndef FLAC__INTEGER_ONLY_LIBRARY 46 47#ifndef M_LN2 48/* math.h in VC++ doesn't seem to have this (how Microsoft is that?) */ 49#define M_LN2 0.69314718055994530942 50#endif 51 52/* OPT: #undef'ing this may improve the speed on some architectures */ 53#define FLAC__LPC_UNROLLED_FILTER_LOOPS 54 55 56void FLAC__lpc_window_data(const FLAC__int32 in[], const FLAC__real window[], FLAC__real out[], unsigned data_len) 57{ 58 unsigned i; 59 for(i = 0; i < data_len; i++) 60 out[i] = in[i] * window[i]; 61} 62 63void FLAC__lpc_compute_autocorrelation(const FLAC__real data[], unsigned data_len, unsigned lag, FLAC__real autoc[]) 64{ 65 /* a readable, but slower, version */ 66#if 0 67 FLAC__real d; 68 unsigned i; 69 70 FLAC__ASSERT(lag > 0); 71 FLAC__ASSERT(lag <= data_len); 72 73 /* 74 * Technically we should subtract the mean first like so: 75 * for(i = 0; i < data_len; i++) 76 * data[i] -= mean; 77 * but it appears not to make enough of a difference to matter, and 78 * most signals are already closely centered around zero 79 */ 80 while(lag--) { 81 for(i = lag, d = 0.0; i < data_len; i++) 82 d += data[i] * data[i - lag]; 83 autoc[lag] = d; 84 } 85#endif 86 87 /* 88 * this version tends to run faster because of better data locality 89 * ('data_len' is usually much larger than 'lag') 90 */ 91 FLAC__real d; 92 unsigned sample, coeff; 93 const unsigned limit = data_len - lag; 94 95 FLAC__ASSERT(lag > 0); 96 FLAC__ASSERT(lag <= data_len); 97 98 for(coeff = 0; coeff < lag; coeff++) 99 autoc[coeff] = 0.0; 100 for(sample = 0; sample <= limit; sample++) { 101 d = data[sample]; 102 for(coeff = 0; coeff < lag; coeff++) 103 autoc[coeff] += d * data[sample+coeff]; 104 } 105 for(; sample < data_len; sample++) { 106 d = data[sample]; 107 for(coeff = 0; coeff < data_len - sample; coeff++) 108 autoc[coeff] += d * data[sample+coeff]; 109 } 110} 111 112void FLAC__lpc_compute_lp_coefficients(const FLAC__real autoc[], unsigned *max_order, FLAC__real lp_coeff[][FLAC__MAX_LPC_ORDER], FLAC__double error[]) 113{ 114 unsigned i, j; 115 FLAC__double r, err, ref[FLAC__MAX_LPC_ORDER], lpc[FLAC__MAX_LPC_ORDER]; 116 117 FLAC__ASSERT(0 != max_order); 118 FLAC__ASSERT(0 < *max_order); 119 FLAC__ASSERT(*max_order <= FLAC__MAX_LPC_ORDER); 120 FLAC__ASSERT(autoc[0] != 0.0); 121 122 err = autoc[0]; 123 124 for(i = 0; i < *max_order; i++) { 125 /* Sum up this iteration's reflection coefficient. */ 126 r = -autoc[i+1]; 127 for(j = 0; j < i; j++) 128 r -= lpc[j] * autoc[i-j]; 129 ref[i] = (r/=err); 130 131 /* Update LPC coefficients and total error. */ 132 lpc[i]=r; 133 for(j = 0; j < (i>>1); j++) { 134 FLAC__double tmp = lpc[j]; 135 lpc[j] += r * lpc[i-1-j]; 136 lpc[i-1-j] += r * tmp; 137 } 138 if(i & 1) 139 lpc[j] += lpc[j] * r; 140 141 err *= (1.0 - r * r); 142 143 /* save this order */ 144 for(j = 0; j <= i; j++) 145 lp_coeff[i][j] = (FLAC__real)(-lpc[j]); /* negate FIR filter coeff to get predictor coeff */ 146 error[i] = err; 147 148 /* see SF bug #1601812 http://sourceforge.net/tracker/index.php?func=detail&aid=1601812&group_id=13478&atid=113478 */ 149 if(err == 0.0) { 150 *max_order = i+1; 151 return; 152 } 153 } 154} 155 156int FLAC__lpc_quantize_coefficients(const FLAC__real lp_coeff[], unsigned order, unsigned precision, FLAC__int32 qlp_coeff[], int *shift) 157{ 158 unsigned i; 159 FLAC__double cmax; 160 FLAC__int32 qmax, qmin; 161 162 FLAC__ASSERT(precision > 0); 163 FLAC__ASSERT(precision >= FLAC__MIN_QLP_COEFF_PRECISION); 164 165 /* drop one bit for the sign; from here on out we consider only |lp_coeff[i]| */ 166 precision--; 167 qmax = 1 << precision; 168 qmin = -qmax; 169 qmax--; 170 171 /* calc cmax = max( |lp_coeff[i]| ) */ 172 cmax = 0.0; 173 for(i = 0; i < order; i++) { 174 const FLAC__double d = fabs(lp_coeff[i]); 175 if(d > cmax) 176 cmax = d; 177 } 178 179 if(cmax <= 0.0) { 180 /* => coefficients are all 0, which means our constant-detect didn't work */ 181 return 2; 182 } 183 else { 184 const int max_shiftlimit = (1 << (FLAC__SUBFRAME_LPC_QLP_SHIFT_LEN-1)) - 1; 185 const int min_shiftlimit = -max_shiftlimit - 1; 186 int log2cmax; 187 188 (void)frexp(cmax, &log2cmax); 189 log2cmax--; 190 *shift = (int)precision - log2cmax - 1; 191 192 if(*shift > max_shiftlimit) 193 *shift = max_shiftlimit; 194 else if(*shift < min_shiftlimit) 195 return 1; 196 } 197 198 if(*shift >= 0) { 199 FLAC__double error = 0.0; 200 FLAC__int32 q; 201 for(i = 0; i < order; i++) { 202 error += lp_coeff[i] * (1 << *shift); 203#if 1 /* unfortunately lround() is C99 */ 204 if(error >= 0.0) 205 q = (FLAC__int32)(error + 0.5); 206 else 207 q = (FLAC__int32)(error - 0.5); 208#else 209 q = lround(error); 210#endif 211#ifdef FLAC__OVERFLOW_DETECT 212 if(q > qmax+1) /* we expect q==qmax+1 occasionally due to rounding */ 213 fprintf(stderr,"FLAC__lpc_quantize_coefficients: quantizer overflow: q>qmax %d>%d shift=%d cmax=%f precision=%u lpc[%u]=%f\n",q,qmax,*shift,cmax,precision+1,i,lp_coeff[i]); 214 else if(q < qmin) 215 fprintf(stderr,"FLAC__lpc_quantize_coefficients: quantizer overflow: q<qmin %d<%d shift=%d cmax=%f precision=%u lpc[%u]=%f\n",q,qmin,*shift,cmax,precision+1,i,lp_coeff[i]); 216#endif 217 if(q > qmax) 218 q = qmax; 219 else if(q < qmin) 220 q = qmin; 221 error -= q; 222 qlp_coeff[i] = q; 223 } 224 } 225 /* negative shift is very rare but due to design flaw, negative shift is 226 * a NOP in the decoder, so it must be handled specially by scaling down 227 * coeffs 228 */ 229 else { 230 const int nshift = -(*shift); 231 FLAC__double error = 0.0; 232 FLAC__int32 q; 233#ifdef DEBUG 234 fprintf(stderr,"FLAC__lpc_quantize_coefficients: negative shift=%d order=%u cmax=%f\n", *shift, order, cmax); 235#endif 236 for(i = 0; i < order; i++) { 237 error += lp_coeff[i] / (1 << nshift); 238#if 1 /* unfortunately lround() is C99 */ 239 if(error >= 0.0) 240 q = (FLAC__int32)(error + 0.5); 241 else 242 q = (FLAC__int32)(error - 0.5); 243#else 244 q = lround(error); 245#endif 246#ifdef FLAC__OVERFLOW_DETECT 247 if(q > qmax+1) /* we expect q==qmax+1 occasionally due to rounding */ 248 fprintf(stderr,"FLAC__lpc_quantize_coefficients: quantizer overflow: q>qmax %d>%d shift=%d cmax=%f precision=%u lpc[%u]=%f\n",q,qmax,*shift,cmax,precision+1,i,lp_coeff[i]); 249 else if(q < qmin) 250 fprintf(stderr,"FLAC__lpc_quantize_coefficients: quantizer overflow: q<qmin %d<%d shift=%d cmax=%f precision=%u lpc[%u]=%f\n",q,qmin,*shift,cmax,precision+1,i,lp_coeff[i]); 251#endif 252 if(q > qmax) 253 q = qmax; 254 else if(q < qmin) 255 q = qmin; 256 error -= q; 257 qlp_coeff[i] = q; 258 } 259 *shift = 0; 260 } 261 262 return 0; 263} 264 265void FLAC__lpc_compute_residual_from_qlp_coefficients(const FLAC__int32 *data, unsigned data_len, const FLAC__int32 qlp_coeff[], unsigned order, int lp_quantization, FLAC__int32 residual[]) 266#if defined(FLAC__OVERFLOW_DETECT) || !defined(FLAC__LPC_UNROLLED_FILTER_LOOPS) 267{ 268 FLAC__int64 sumo; 269 unsigned i, j; 270 FLAC__int32 sum; 271 const FLAC__int32 *history; 272 273#ifdef FLAC__OVERFLOW_DETECT_VERBOSE 274 fprintf(stderr,"FLAC__lpc_compute_residual_from_qlp_coefficients: data_len=%d, order=%u, lpq=%d",data_len,order,lp_quantization); 275 for(i=0;i<order;i++) 276 fprintf(stderr,", q[%u]=%d",i,qlp_coeff[i]); 277 fprintf(stderr,"\n"); 278#endif 279 FLAC__ASSERT(order > 0); 280 281 for(i = 0; i < data_len; i++) { 282 sumo = 0; 283 sum = 0; 284 history = data; 285 for(j = 0; j < order; j++) { 286 sum += qlp_coeff[j] * (*(--history)); 287 sumo += (FLAC__int64)qlp_coeff[j] * (FLAC__int64)(*history); 288#if defined _MSC_VER 289 if(sumo > 2147483647I64 || sumo < -2147483648I64) 290 fprintf(stderr,"FLAC__lpc_compute_residual_from_qlp_coefficients: OVERFLOW, i=%u, j=%u, c=%d, d=%d, sumo=%I64d\n",i,j,qlp_coeff[j],*history,sumo); 291#else 292 if(sumo > 2147483647ll || sumo < -2147483648ll) 293 fprintf(stderr,"FLAC__lpc_compute_residual_from_qlp_coefficients: OVERFLOW, i=%u, j=%u, c=%d, d=%d, sumo=%lld\n",i,j,qlp_coeff[j],*history,(long long)sumo); 294#endif 295 } 296 *(residual++) = *(data++) - (sum >> lp_quantization); 297 } 298 299 /* Here's a slower but clearer version: 300 for(i = 0; i < data_len; i++) { 301 sum = 0; 302 for(j = 0; j < order; j++) 303 sum += qlp_coeff[j] * data[i-j-1]; 304 residual[i] = data[i] - (sum >> lp_quantization); 305 } 306 */ 307} 308#else /* fully unrolled version for normal use */ 309{ 310 int i; 311 FLAC__int32 sum; 312 313 FLAC__ASSERT(order > 0); 314 FLAC__ASSERT(order <= 32); 315 316 /* 317 * We do unique versions up to 12th order since that's the subset limit. 318 * Also they are roughly ordered to match frequency of occurrence to 319 * minimize branching. 320 */ 321 if(order <= 12) { 322 if(order > 8) { 323 if(order > 10) { 324 if(order == 12) { 325 for(i = 0; i < (int)data_len; i++) { 326 sum = 0; 327 sum += qlp_coeff[11] * data[i-12]; 328 sum += qlp_coeff[10] * data[i-11]; 329 sum += qlp_coeff[9] * data[i-10]; 330 sum += qlp_coeff[8] * data[i-9]; 331 sum += qlp_coeff[7] * data[i-8]; 332 sum += qlp_coeff[6] * data[i-7]; 333 sum += qlp_coeff[5] * data[i-6]; 334 sum += qlp_coeff[4] * data[i-5]; 335 sum += qlp_coeff[3] * data[i-4]; 336 sum += qlp_coeff[2] * data[i-3]; 337 sum += qlp_coeff[1] * data[i-2]; 338 sum += qlp_coeff[0] * data[i-1]; 339 residual[i] = data[i] - (sum >> lp_quantization); 340 } 341 } 342 else { /* order == 11 */ 343 for(i = 0; i < (int)data_len; i++) { 344 sum = 0; 345 sum += qlp_coeff[10] * data[i-11]; 346 sum += qlp_coeff[9] * data[i-10]; 347 sum += qlp_coeff[8] * data[i-9]; 348 sum += qlp_coeff[7] * data[i-8]; 349 sum += qlp_coeff[6] * data[i-7]; 350 sum += qlp_coeff[5] * data[i-6]; 351 sum += qlp_coeff[4] * data[i-5]; 352 sum += qlp_coeff[3] * data[i-4]; 353 sum += qlp_coeff[2] * data[i-3]; 354 sum += qlp_coeff[1] * data[i-2]; 355 sum += qlp_coeff[0] * data[i-1]; 356 residual[i] = data[i] - (sum >> lp_quantization); 357 } 358 } 359 } 360 else { 361 if(order == 10) { 362 for(i = 0; i < (int)data_len; i++) { 363 sum = 0; 364 sum += qlp_coeff[9] * data[i-10]; 365 sum += qlp_coeff[8] * data[i-9]; 366 sum += qlp_coeff[7] * data[i-8]; 367 sum += qlp_coeff[6] * data[i-7]; 368 sum += qlp_coeff[5] * data[i-6]; 369 sum += qlp_coeff[4] * data[i-5]; 370 sum += qlp_coeff[3] * data[i-4]; 371 sum += qlp_coeff[2] * data[i-3]; 372 sum += qlp_coeff[1] * data[i-2]; 373 sum += qlp_coeff[0] * data[i-1]; 374 residual[i] = data[i] - (sum >> lp_quantization); 375 } 376 } 377 else { /* order == 9 */ 378 for(i = 0; i < (int)data_len; i++) { 379 sum = 0; 380 sum += qlp_coeff[8] * data[i-9]; 381 sum += qlp_coeff[7] * data[i-8]; 382 sum += qlp_coeff[6] * data[i-7]; 383 sum += qlp_coeff[5] * data[i-6]; 384 sum += qlp_coeff[4] * data[i-5]; 385 sum += qlp_coeff[3] * data[i-4]; 386 sum += qlp_coeff[2] * data[i-3]; 387 sum += qlp_coeff[1] * data[i-2]; 388 sum += qlp_coeff[0] * data[i-1]; 389 residual[i] = data[i] - (sum >> lp_quantization); 390 } 391 } 392 } 393 } 394 else if(order > 4) { 395 if(order > 6) { 396 if(order == 8) { 397 for(i = 0; i < (int)data_len; i++) { 398 sum = 0; 399 sum += qlp_coeff[7] * data[i-8]; 400 sum += qlp_coeff[6] * data[i-7]; 401 sum += qlp_coeff[5] * data[i-6]; 402 sum += qlp_coeff[4] * data[i-5]; 403 sum += qlp_coeff[3] * data[i-4]; 404 sum += qlp_coeff[2] * data[i-3]; 405 sum += qlp_coeff[1] * data[i-2]; 406 sum += qlp_coeff[0] * data[i-1]; 407 residual[i] = data[i] - (sum >> lp_quantization); 408 } 409 } 410 else { /* order == 7 */ 411 for(i = 0; i < (int)data_len; i++) { 412 sum = 0; 413 sum += qlp_coeff[6] * data[i-7]; 414 sum += qlp_coeff[5] * data[i-6]; 415 sum += qlp_coeff[4] * data[i-5]; 416 sum += qlp_coeff[3] * data[i-4]; 417 sum += qlp_coeff[2] * data[i-3]; 418 sum += qlp_coeff[1] * data[i-2]; 419 sum += qlp_coeff[0] * data[i-1]; 420 residual[i] = data[i] - (sum >> lp_quantization); 421 } 422 } 423 } 424 else { 425 if(order == 6) { 426 for(i = 0; i < (int)data_len; i++) { 427 sum = 0; 428 sum += qlp_coeff[5] * data[i-6]; 429 sum += qlp_coeff[4] * data[i-5]; 430 sum += qlp_coeff[3] * data[i-4]; 431 sum += qlp_coeff[2] * data[i-3]; 432 sum += qlp_coeff[1] * data[i-2]; 433 sum += qlp_coeff[0] * data[i-1]; 434 residual[i] = data[i] - (sum >> lp_quantization); 435 } 436 } 437 else { /* order == 5 */ 438 for(i = 0; i < (int)data_len; i++) { 439 sum = 0; 440 sum += qlp_coeff[4] * data[i-5]; 441 sum += qlp_coeff[3] * data[i-4]; 442 sum += qlp_coeff[2] * data[i-3]; 443 sum += qlp_coeff[1] * data[i-2]; 444 sum += qlp_coeff[0] * data[i-1]; 445 residual[i] = data[i] - (sum >> lp_quantization); 446 } 447 } 448 } 449 } 450 else { 451 if(order > 2) { 452 if(order == 4) { 453 for(i = 0; i < (int)data_len; i++) { 454 sum = 0; 455 sum += qlp_coeff[3] * data[i-4]; 456 sum += qlp_coeff[2] * data[i-3]; 457 sum += qlp_coeff[1] * data[i-2]; 458 sum += qlp_coeff[0] * data[i-1]; 459 residual[i] = data[i] - (sum >> lp_quantization); 460 } 461 } 462 else { /* order == 3 */ 463 for(i = 0; i < (int)data_len; i++) { 464 sum = 0; 465 sum += qlp_coeff[2] * data[i-3]; 466 sum += qlp_coeff[1] * data[i-2]; 467 sum += qlp_coeff[0] * data[i-1]; 468 residual[i] = data[i] - (sum >> lp_quantization); 469 } 470 } 471 } 472 else { 473 if(order == 2) { 474 for(i = 0; i < (int)data_len; i++) { 475 sum = 0; 476 sum += qlp_coeff[1] * data[i-2]; 477 sum += qlp_coeff[0] * data[i-1]; 478 residual[i] = data[i] - (sum >> lp_quantization); 479 } 480 } 481 else { /* order == 1 */ 482 for(i = 0; i < (int)data_len; i++) 483 residual[i] = data[i] - ((qlp_coeff[0] * data[i-1]) >> lp_quantization); 484 } 485 } 486 } 487 } 488 else { /* order > 12 */ 489 for(i = 0; i < (int)data_len; i++) { 490 sum = 0; 491 switch(order) { 492 case 32: sum += qlp_coeff[31] * data[i-32]; 493 case 31: sum += qlp_coeff[30] * data[i-31]; 494 case 30: sum += qlp_coeff[29] * data[i-30]; 495 case 29: sum += qlp_coeff[28] * data[i-29]; 496 case 28: sum += qlp_coeff[27] * data[i-28]; 497 case 27: sum += qlp_coeff[26] * data[i-27]; 498 case 26: sum += qlp_coeff[25] * data[i-26]; 499 case 25: sum += qlp_coeff[24] * data[i-25]; 500 case 24: sum += qlp_coeff[23] * data[i-24]; 501 case 23: sum += qlp_coeff[22] * data[i-23]; 502 case 22: sum += qlp_coeff[21] * data[i-22]; 503 case 21: sum += qlp_coeff[20] * data[i-21]; 504 case 20: sum += qlp_coeff[19] * data[i-20]; 505 case 19: sum += qlp_coeff[18] * data[i-19]; 506 case 18: sum += qlp_coeff[17] * data[i-18]; 507 case 17: sum += qlp_coeff[16] * data[i-17]; 508 case 16: sum += qlp_coeff[15] * data[i-16]; 509 case 15: sum += qlp_coeff[14] * data[i-15]; 510 case 14: sum += qlp_coeff[13] * data[i-14]; 511 case 13: sum += qlp_coeff[12] * data[i-13]; 512 sum += qlp_coeff[11] * data[i-12]; 513 sum += qlp_coeff[10] * data[i-11]; 514 sum += qlp_coeff[ 9] * data[i-10]; 515 sum += qlp_coeff[ 8] * data[i- 9]; 516 sum += qlp_coeff[ 7] * data[i- 8]; 517 sum += qlp_coeff[ 6] * data[i- 7]; 518 sum += qlp_coeff[ 5] * data[i- 6]; 519 sum += qlp_coeff[ 4] * data[i- 5]; 520 sum += qlp_coeff[ 3] * data[i- 4]; 521 sum += qlp_coeff[ 2] * data[i- 3]; 522 sum += qlp_coeff[ 1] * data[i- 2]; 523 sum += qlp_coeff[ 0] * data[i- 1]; 524 } 525 residual[i] = data[i] - (sum >> lp_quantization); 526 } 527 } 528} 529#endif 530 531void FLAC__lpc_compute_residual_from_qlp_coefficients_wide(const FLAC__int32 *data, unsigned data_len, const FLAC__int32 qlp_coeff[], unsigned order, int lp_quantization, FLAC__int32 residual[]) 532#if defined(FLAC__OVERFLOW_DETECT) || !defined(FLAC__LPC_UNROLLED_FILTER_LOOPS) 533{ 534 unsigned i, j; 535 FLAC__int64 sum; 536 const FLAC__int32 *history; 537 538#ifdef FLAC__OVERFLOW_DETECT_VERBOSE 539 fprintf(stderr,"FLAC__lpc_compute_residual_from_qlp_coefficients_wide: data_len=%d, order=%u, lpq=%d",data_len,order,lp_quantization); 540 for(i=0;i<order;i++) 541 fprintf(stderr,", q[%u]=%d",i,qlp_coeff[i]); 542 fprintf(stderr,"\n"); 543#endif 544 FLAC__ASSERT(order > 0); 545 546 for(i = 0; i < data_len; i++) { 547 sum = 0; 548 history = data; 549 for(j = 0; j < order; j++) 550 sum += (FLAC__int64)qlp_coeff[j] * (FLAC__int64)(*(--history)); 551 if(FLAC__bitmath_silog2_wide(sum >> lp_quantization) > 32) { 552#if defined _MSC_VER 553 fprintf(stderr,"FLAC__lpc_compute_residual_from_qlp_coefficients_wide: OVERFLOW, i=%u, sum=%I64d\n", i, sum >> lp_quantization); 554#else 555 fprintf(stderr,"FLAC__lpc_compute_residual_from_qlp_coefficients_wide: OVERFLOW, i=%u, sum=%lld\n", i, (long long)(sum >> lp_quantization)); 556#endif 557 break; 558 } 559 if(FLAC__bitmath_silog2_wide((FLAC__int64)(*data) - (sum >> lp_quantization)) > 32) { 560#if defined _MSC_VER 561 fprintf(stderr,"FLAC__lpc_compute_residual_from_qlp_coefficients_wide: OVERFLOW, i=%u, data=%d, sum=%I64d, residual=%I64d\n", i, *data, sum >> lp_quantization, (FLAC__int64)(*data) - (sum >> lp_quantization)); 562#else 563 fprintf(stderr,"FLAC__lpc_compute_residual_from_qlp_coefficients_wide: OVERFLOW, i=%u, data=%d, sum=%lld, residual=%lld\n", i, *data, (long long)(sum >> lp_quantization), (long long)((FLAC__int64)(*data) - (sum >> lp_quantization))); 564#endif 565 break; 566 } 567 *(residual++) = *(data++) - (FLAC__int32)(sum >> lp_quantization); 568 } 569} 570#else /* fully unrolled version for normal use */ 571{ 572 int i; 573 FLAC__int64 sum; 574 575 FLAC__ASSERT(order > 0); 576 FLAC__ASSERT(order <= 32); 577 578 /* 579 * We do unique versions up to 12th order since that's the subset limit. 580 * Also they are roughly ordered to match frequency of occurrence to 581 * minimize branching. 582 */ 583 if(order <= 12) { 584 if(order > 8) { 585 if(order > 10) { 586 if(order == 12) { 587 for(i = 0; i < (int)data_len; i++) { 588 sum = 0; 589 sum += qlp_coeff[11] * (FLAC__int64)data[i-12]; 590 sum += qlp_coeff[10] * (FLAC__int64)data[i-11]; 591 sum += qlp_coeff[9] * (FLAC__int64)data[i-10]; 592 sum += qlp_coeff[8] * (FLAC__int64)data[i-9]; 593 sum += qlp_coeff[7] * (FLAC__int64)data[i-8]; 594 sum += qlp_coeff[6] * (FLAC__int64)data[i-7]; 595 sum += qlp_coeff[5] * (FLAC__int64)data[i-6]; 596 sum += qlp_coeff[4] * (FLAC__int64)data[i-5]; 597 sum += qlp_coeff[3] * (FLAC__int64)data[i-4]; 598 sum += qlp_coeff[2] * (FLAC__int64)data[i-3]; 599 sum += qlp_coeff[1] * (FLAC__int64)data[i-2]; 600 sum += qlp_coeff[0] * (FLAC__int64)data[i-1]; 601 residual[i] = data[i] - (FLAC__int32)(sum >> lp_quantization); 602 } 603 } 604 else { /* order == 11 */ 605 for(i = 0; i < (int)data_len; i++) { 606 sum = 0; 607 sum += qlp_coeff[10] * (FLAC__int64)data[i-11]; 608 sum += qlp_coeff[9] * (FLAC__int64)data[i-10]; 609 sum += qlp_coeff[8] * (FLAC__int64)data[i-9]; 610 sum += qlp_coeff[7] * (FLAC__int64)data[i-8]; 611 sum += qlp_coeff[6] * (FLAC__int64)data[i-7]; 612 sum += qlp_coeff[5] * (FLAC__int64)data[i-6]; 613 sum += qlp_coeff[4] * (FLAC__int64)data[i-5]; 614 sum += qlp_coeff[3] * (FLAC__int64)data[i-4]; 615 sum += qlp_coeff[2] * (FLAC__int64)data[i-3]; 616 sum += qlp_coeff[1] * (FLAC__int64)data[i-2]; 617 sum += qlp_coeff[0] * (FLAC__int64)data[i-1]; 618 residual[i] = data[i] - (FLAC__int32)(sum >> lp_quantization); 619 } 620 } 621 } 622 else { 623 if(order == 10) { 624 for(i = 0; i < (int)data_len; i++) { 625 sum = 0; 626 sum += qlp_coeff[9] * (FLAC__int64)data[i-10]; 627 sum += qlp_coeff[8] * (FLAC__int64)data[i-9]; 628 sum += qlp_coeff[7] * (FLAC__int64)data[i-8]; 629 sum += qlp_coeff[6] * (FLAC__int64)data[i-7]; 630 sum += qlp_coeff[5] * (FLAC__int64)data[i-6]; 631 sum += qlp_coeff[4] * (FLAC__int64)data[i-5]; 632 sum += qlp_coeff[3] * (FLAC__int64)data[i-4]; 633 sum += qlp_coeff[2] * (FLAC__int64)data[i-3]; 634 sum += qlp_coeff[1] * (FLAC__int64)data[i-2]; 635 sum += qlp_coeff[0] * (FLAC__int64)data[i-1]; 636 residual[i] = data[i] - (FLAC__int32)(sum >> lp_quantization); 637 } 638 } 639 else { /* order == 9 */ 640 for(i = 0; i < (int)data_len; i++) { 641 sum = 0; 642 sum += qlp_coeff[8] * (FLAC__int64)data[i-9]; 643 sum += qlp_coeff[7] * (FLAC__int64)data[i-8]; 644 sum += qlp_coeff[6] * (FLAC__int64)data[i-7]; 645 sum += qlp_coeff[5] * (FLAC__int64)data[i-6]; 646 sum += qlp_coeff[4] * (FLAC__int64)data[i-5]; 647 sum += qlp_coeff[3] * (FLAC__int64)data[i-4]; 648 sum += qlp_coeff[2] * (FLAC__int64)data[i-3]; 649 sum += qlp_coeff[1] * (FLAC__int64)data[i-2]; 650 sum += qlp_coeff[0] * (FLAC__int64)data[i-1]; 651 residual[i] = data[i] - (FLAC__int32)(sum >> lp_quantization); 652 } 653 } 654 } 655 } 656 else if(order > 4) { 657 if(order > 6) { 658 if(order == 8) { 659 for(i = 0; i < (int)data_len; i++) { 660 sum = 0; 661 sum += qlp_coeff[7] * (FLAC__int64)data[i-8]; 662 sum += qlp_coeff[6] * (FLAC__int64)data[i-7]; 663 sum += qlp_coeff[5] * (FLAC__int64)data[i-6]; 664 sum += qlp_coeff[4] * (FLAC__int64)data[i-5]; 665 sum += qlp_coeff[3] * (FLAC__int64)data[i-4]; 666 sum += qlp_coeff[2] * (FLAC__int64)data[i-3]; 667 sum += qlp_coeff[1] * (FLAC__int64)data[i-2]; 668 sum += qlp_coeff[0] * (FLAC__int64)data[i-1]; 669 residual[i] = data[i] - (FLAC__int32)(sum >> lp_quantization); 670 } 671 } 672 else { /* order == 7 */ 673 for(i = 0; i < (int)data_len; i++) { 674 sum = 0; 675 sum += qlp_coeff[6] * (FLAC__int64)data[i-7]; 676 sum += qlp_coeff[5] * (FLAC__int64)data[i-6]; 677 sum += qlp_coeff[4] * (FLAC__int64)data[i-5]; 678 sum += qlp_coeff[3] * (FLAC__int64)data[i-4]; 679 sum += qlp_coeff[2] * (FLAC__int64)data[i-3]; 680 sum += qlp_coeff[1] * (FLAC__int64)data[i-2]; 681 sum += qlp_coeff[0] * (FLAC__int64)data[i-1]; 682 residual[i] = data[i] - (FLAC__int32)(sum >> lp_quantization); 683 } 684 } 685 } 686 else { 687 if(order == 6) { 688 for(i = 0; i < (int)data_len; i++) { 689 sum = 0; 690 sum += qlp_coeff[5] * (FLAC__int64)data[i-6]; 691 sum += qlp_coeff[4] * (FLAC__int64)data[i-5]; 692 sum += qlp_coeff[3] * (FLAC__int64)data[i-4]; 693 sum += qlp_coeff[2] * (FLAC__int64)data[i-3]; 694 sum += qlp_coeff[1] * (FLAC__int64)data[i-2]; 695 sum += qlp_coeff[0] * (FLAC__int64)data[i-1]; 696 residual[i] = data[i] - (FLAC__int32)(sum >> lp_quantization); 697 } 698 } 699 else { /* order == 5 */ 700 for(i = 0; i < (int)data_len; i++) { 701 sum = 0; 702 sum += qlp_coeff[4] * (FLAC__int64)data[i-5]; 703 sum += qlp_coeff[3] * (FLAC__int64)data[i-4]; 704 sum += qlp_coeff[2] * (FLAC__int64)data[i-3]; 705 sum += qlp_coeff[1] * (FLAC__int64)data[i-2]; 706 sum += qlp_coeff[0] * (FLAC__int64)data[i-1]; 707 residual[i] = data[i] - (FLAC__int32)(sum >> lp_quantization); 708 } 709 } 710 } 711 } 712 else { 713 if(order > 2) { 714 if(order == 4) { 715 for(i = 0; i < (int)data_len; i++) { 716 sum = 0; 717 sum += qlp_coeff[3] * (FLAC__int64)data[i-4]; 718 sum += qlp_coeff[2] * (FLAC__int64)data[i-3]; 719 sum += qlp_coeff[1] * (FLAC__int64)data[i-2]; 720 sum += qlp_coeff[0] * (FLAC__int64)data[i-1]; 721 residual[i] = data[i] - (FLAC__int32)(sum >> lp_quantization); 722 } 723 } 724 else { /* order == 3 */ 725 for(i = 0; i < (int)data_len; i++) { 726 sum = 0; 727 sum += qlp_coeff[2] * (FLAC__int64)data[i-3]; 728 sum += qlp_coeff[1] * (FLAC__int64)data[i-2]; 729 sum += qlp_coeff[0] * (FLAC__int64)data[i-1]; 730 residual[i] = data[i] - (FLAC__int32)(sum >> lp_quantization); 731 } 732 } 733 } 734 else { 735 if(order == 2) { 736 for(i = 0; i < (int)data_len; i++) { 737 sum = 0; 738 sum += qlp_coeff[1] * (FLAC__int64)data[i-2]; 739 sum += qlp_coeff[0] * (FLAC__int64)data[i-1]; 740 residual[i] = data[i] - (FLAC__int32)(sum >> lp_quantization); 741 } 742 } 743 else { /* order == 1 */ 744 for(i = 0; i < (int)data_len; i++) 745 residual[i] = data[i] - (FLAC__int32)((qlp_coeff[0] * (FLAC__int64)data[i-1]) >> lp_quantization); 746 } 747 } 748 } 749 } 750 else { /* order > 12 */ 751 for(i = 0; i < (int)data_len; i++) { 752 sum = 0; 753 switch(order) { 754 case 32: sum += qlp_coeff[31] * (FLAC__int64)data[i-32]; 755 case 31: sum += qlp_coeff[30] * (FLAC__int64)data[i-31]; 756 case 30: sum += qlp_coeff[29] * (FLAC__int64)data[i-30]; 757 case 29: sum += qlp_coeff[28] * (FLAC__int64)data[i-29]; 758 case 28: sum += qlp_coeff[27] * (FLAC__int64)data[i-28]; 759 case 27: sum += qlp_coeff[26] * (FLAC__int64)data[i-27]; 760 case 26: sum += qlp_coeff[25] * (FLAC__int64)data[i-26]; 761 case 25: sum += qlp_coeff[24] * (FLAC__int64)data[i-25]; 762 case 24: sum += qlp_coeff[23] * (FLAC__int64)data[i-24]; 763 case 23: sum += qlp_coeff[22] * (FLAC__int64)data[i-23]; 764 case 22: sum += qlp_coeff[21] * (FLAC__int64)data[i-22]; 765 case 21: sum += qlp_coeff[20] * (FLAC__int64)data[i-21]; 766 case 20: sum += qlp_coeff[19] * (FLAC__int64)data[i-20]; 767 case 19: sum += qlp_coeff[18] * (FLAC__int64)data[i-19]; 768 case 18: sum += qlp_coeff[17] * (FLAC__int64)data[i-18]; 769 case 17: sum += qlp_coeff[16] * (FLAC__int64)data[i-17]; 770 case 16: sum += qlp_coeff[15] * (FLAC__int64)data[i-16]; 771 case 15: sum += qlp_coeff[14] * (FLAC__int64)data[i-15]; 772 case 14: sum += qlp_coeff[13] * (FLAC__int64)data[i-14]; 773 case 13: sum += qlp_coeff[12] * (FLAC__int64)data[i-13]; 774 sum += qlp_coeff[11] * (FLAC__int64)data[i-12]; 775 sum += qlp_coeff[10] * (FLAC__int64)data[i-11]; 776 sum += qlp_coeff[ 9] * (FLAC__int64)data[i-10]; 777 sum += qlp_coeff[ 8] * (FLAC__int64)data[i- 9]; 778 sum += qlp_coeff[ 7] * (FLAC__int64)data[i- 8]; 779 sum += qlp_coeff[ 6] * (FLAC__int64)data[i- 7]; 780 sum += qlp_coeff[ 5] * (FLAC__int64)data[i- 6]; 781 sum += qlp_coeff[ 4] * (FLAC__int64)data[i- 5]; 782 sum += qlp_coeff[ 3] * (FLAC__int64)data[i- 4]; 783 sum += qlp_coeff[ 2] * (FLAC__int64)data[i- 3]; 784 sum += qlp_coeff[ 1] * (FLAC__int64)data[i- 2]; 785 sum += qlp_coeff[ 0] * (FLAC__int64)data[i- 1]; 786 } 787 residual[i] = data[i] - (FLAC__int32)(sum >> lp_quantization); 788 } 789 } 790} 791#endif 792 793#endif /* !defined FLAC__INTEGER_ONLY_LIBRARY */ 794 795void FLAC__lpc_restore_signal(const FLAC__int32 residual[], unsigned data_len, const FLAC__int32 qlp_coeff[], unsigned order, int lp_quantization, FLAC__int32 data[]) 796#if defined(FLAC__OVERFLOW_DETECT) || !defined(FLAC__LPC_UNROLLED_FILTER_LOOPS) 797{ 798 FLAC__int64 sumo; 799 unsigned i, j; 800 FLAC__int32 sum; 801 const FLAC__int32 *r = residual, *history; 802 803#ifdef FLAC__OVERFLOW_DETECT_VERBOSE 804 fprintf(stderr,"FLAC__lpc_restore_signal: data_len=%d, order=%u, lpq=%d",data_len,order,lp_quantization); 805 for(i=0;i<order;i++) 806 fprintf(stderr,", q[%u]=%d",i,qlp_coeff[i]); 807 fprintf(stderr,"\n"); 808#endif 809 FLAC__ASSERT(order > 0); 810 811 for(i = 0; i < data_len; i++) { 812 sumo = 0; 813 sum = 0; 814 history = data; 815 for(j = 0; j < order; j++) { 816 sum += qlp_coeff[j] * (*(--history)); 817 sumo += (FLAC__int64)qlp_coeff[j] * (FLAC__int64)(*history); 818#ifdef FLAC__OVERFLOW_DETECT_VERBOSE 819#if defined _MSC_VER 820 if(sumo > 2147483647I64 || sumo < -2147483648I64) 821 fprintf(stderr,"FLAC__lpc_restore_signal: OVERFLOW, i=%u, j=%u, c=%d, d=%d, sumo=%I64d\n",i,j,qlp_coeff[j],*history,sumo); 822#else 823 if(sumo > 2147483647ll || sumo < -2147483648ll) 824 fprintf(stderr,"FLAC__lpc_restore_signal: OVERFLOW, i=%u, j=%u, c=%d, d=%d, sumo=%lld\n",i,j,qlp_coeff[j],*history,(long long)sumo); 825#endif 826#endif 827 } 828 *(data++) = *(r++) + (sum >> lp_quantization); 829 } 830 831 /* Here's a slower but clearer version: 832 for(i = 0; i < data_len; i++) { 833 sum = 0; 834 for(j = 0; j < order; j++) 835 sum += qlp_coeff[j] * data[i-j-1]; 836 data[i] = residual[i] + (sum >> lp_quantization); 837 } 838 */ 839} 840#else /* fully unrolled version for normal use */ 841{ 842 int i; 843 FLAC__int32 sum; 844 845 FLAC__ASSERT(order > 0); 846 FLAC__ASSERT(order <= 32); 847 848 /* 849 * We do unique versions up to 12th order since that's the subset limit. 850 * Also they are roughly ordered to match frequency of occurrence to 851 * minimize branching. 852 */ 853 if(order <= 12) { 854 if(order > 8) { 855 if(order > 10) { 856 if(order == 12) { 857 for(i = 0; i < (int)data_len; i++) { 858 sum = 0; 859 sum += qlp_coeff[11] * data[i-12]; 860 sum += qlp_coeff[10] * data[i-11]; 861 sum += qlp_coeff[9] * data[i-10]; 862 sum += qlp_coeff[8] * data[i-9]; 863 sum += qlp_coeff[7] * data[i-8]; 864 sum += qlp_coeff[6] * data[i-7]; 865 sum += qlp_coeff[5] * data[i-6]; 866 sum += qlp_coeff[4] * data[i-5]; 867 sum += qlp_coeff[3] * data[i-4]; 868 sum += qlp_coeff[2] * data[i-3]; 869 sum += qlp_coeff[1] * data[i-2]; 870 sum += qlp_coeff[0] * data[i-1]; 871 data[i] = residual[i] + (sum >> lp_quantization); 872 } 873 } 874 else { /* order == 11 */ 875 for(i = 0; i < (int)data_len; i++) { 876 sum = 0; 877 sum += qlp_coeff[10] * data[i-11]; 878 sum += qlp_coeff[9] * data[i-10]; 879 sum += qlp_coeff[8] * data[i-9]; 880 sum += qlp_coeff[7] * data[i-8]; 881 sum += qlp_coeff[6] * data[i-7]; 882 sum += qlp_coeff[5] * data[i-6]; 883 sum += qlp_coeff[4] * data[i-5]; 884 sum += qlp_coeff[3] * data[i-4]; 885 sum += qlp_coeff[2] * data[i-3]; 886 sum += qlp_coeff[1] * data[i-2]; 887 sum += qlp_coeff[0] * data[i-1]; 888 data[i] = residual[i] + (sum >> lp_quantization); 889 } 890 } 891 } 892 else { 893 if(order == 10) { 894 for(i = 0; i < (int)data_len; i++) { 895 sum = 0; 896 sum += qlp_coeff[9] * data[i-10]; 897 sum += qlp_coeff[8] * data[i-9]; 898 sum += qlp_coeff[7] * data[i-8]; 899 sum += qlp_coeff[6] * data[i-7]; 900 sum += qlp_coeff[5] * data[i-6]; 901 sum += qlp_coeff[4] * data[i-5]; 902 sum += qlp_coeff[3] * data[i-4]; 903 sum += qlp_coeff[2] * data[i-3]; 904 sum += qlp_coeff[1] * data[i-2]; 905 sum += qlp_coeff[0] * data[i-1]; 906 data[i] = residual[i] + (sum >> lp_quantization); 907 } 908 } 909 else { /* order == 9 */ 910 for(i = 0; i < (int)data_len; i++) { 911 sum = 0; 912 sum += qlp_coeff[8] * data[i-9]; 913 sum += qlp_coeff[7] * data[i-8]; 914 sum += qlp_coeff[6] * data[i-7]; 915 sum += qlp_coeff[5] * data[i-6]; 916 sum += qlp_coeff[4] * data[i-5]; 917 sum += qlp_coeff[3] * data[i-4]; 918 sum += qlp_coeff[2] * data[i-3]; 919 sum += qlp_coeff[1] * data[i-2]; 920 sum += qlp_coeff[0] * data[i-1]; 921 data[i] = residual[i] + (sum >> lp_quantization); 922 } 923 } 924 } 925 } 926 else if(order > 4) { 927 if(order > 6) { 928 if(order == 8) { 929 for(i = 0; i < (int)data_len; i++) { 930 sum = 0; 931 sum += qlp_coeff[7] * data[i-8]; 932 sum += qlp_coeff[6] * data[i-7]; 933 sum += qlp_coeff[5] * data[i-6]; 934 sum += qlp_coeff[4] * data[i-5]; 935 sum += qlp_coeff[3] * data[i-4]; 936 sum += qlp_coeff[2] * data[i-3]; 937 sum += qlp_coeff[1] * data[i-2]; 938 sum += qlp_coeff[0] * data[i-1]; 939 data[i] = residual[i] + (sum >> lp_quantization); 940 } 941 } 942 else { /* order == 7 */ 943 for(i = 0; i < (int)data_len; i++) { 944 sum = 0; 945 sum += qlp_coeff[6] * data[i-7]; 946 sum += qlp_coeff[5] * data[i-6]; 947 sum += qlp_coeff[4] * data[i-5]; 948 sum += qlp_coeff[3] * data[i-4]; 949 sum += qlp_coeff[2] * data[i-3]; 950 sum += qlp_coeff[1] * data[i-2]; 951 sum += qlp_coeff[0] * data[i-1]; 952 data[i] = residual[i] + (sum >> lp_quantization); 953 } 954 } 955 } 956 else { 957 if(order == 6) { 958 for(i = 0; i < (int)data_len; i++) { 959 sum = 0; 960 sum += qlp_coeff[5] * data[i-6]; 961 sum += qlp_coeff[4] * data[i-5]; 962 sum += qlp_coeff[3] * data[i-4]; 963 sum += qlp_coeff[2] * data[i-3]; 964 sum += qlp_coeff[1] * data[i-2]; 965 sum += qlp_coeff[0] * data[i-1]; 966 data[i] = residual[i] + (sum >> lp_quantization); 967 } 968 } 969 else { /* order == 5 */ 970 for(i = 0; i < (int)data_len; i++) { 971 sum = 0; 972 sum += qlp_coeff[4] * data[i-5]; 973 sum += qlp_coeff[3] * data[i-4]; 974 sum += qlp_coeff[2] * data[i-3]; 975 sum += qlp_coeff[1] * data[i-2]; 976 sum += qlp_coeff[0] * data[i-1]; 977 data[i] = residual[i] + (sum >> lp_quantization); 978 } 979 } 980 } 981 } 982 else { 983 if(order > 2) { 984 if(order == 4) { 985 for(i = 0; i < (int)data_len; i++) { 986 sum = 0; 987 sum += qlp_coeff[3] * data[i-4]; 988 sum += qlp_coeff[2] * data[i-3]; 989 sum += qlp_coeff[1] * data[i-2]; 990 sum += qlp_coeff[0] * data[i-1]; 991 data[i] = residual[i] + (sum >> lp_quantization); 992 } 993 } 994 else { /* order == 3 */ 995 for(i = 0; i < (int)data_len; i++) { 996 sum = 0; 997 sum += qlp_coeff[2] * data[i-3]; 998 sum += qlp_coeff[1] * data[i-2]; 999 sum += qlp_coeff[0] * data[i-1]; 1000 data[i] = residual[i] + (sum >> lp_quantization); 1001 } 1002 } 1003 } 1004 else { 1005 if(order == 2) { 1006 for(i = 0; i < (int)data_len; i++) { 1007 sum = 0; 1008 sum += qlp_coeff[1] * data[i-2]; 1009 sum += qlp_coeff[0] * data[i-1]; 1010 data[i] = residual[i] + (sum >> lp_quantization); 1011 } 1012 } 1013 else { /* order == 1 */ 1014 for(i = 0; i < (int)data_len; i++) 1015 data[i] = residual[i] + ((qlp_coeff[0] * data[i-1]) >> lp_quantization); 1016 } 1017 } 1018 } 1019 } 1020 else { /* order > 12 */ 1021 for(i = 0; i < (int)data_len; i++) { 1022 sum = 0; 1023 switch(order) { 1024 case 32: sum += qlp_coeff[31] * data[i-32]; 1025 case 31: sum += qlp_coeff[30] * data[i-31]; 1026 case 30: sum += qlp_coeff[29] * data[i-30]; 1027 case 29: sum += qlp_coeff[28] * data[i-29]; 1028 case 28: sum += qlp_coeff[27] * data[i-28]; 1029 case 27: sum += qlp_coeff[26] * data[i-27]; 1030 case 26: sum += qlp_coeff[25] * data[i-26]; 1031 case 25: sum += qlp_coeff[24] * data[i-25]; 1032 case 24: sum += qlp_coeff[23] * data[i-24]; 1033 case 23: sum += qlp_coeff[22] * data[i-23]; 1034 case 22: sum += qlp_coeff[21] * data[i-22]; 1035 case 21: sum += qlp_coeff[20] * data[i-21]; 1036 case 20: sum += qlp_coeff[19] * data[i-20]; 1037 case 19: sum += qlp_coeff[18] * data[i-19]; 1038 case 18: sum += qlp_coeff[17] * data[i-18]; 1039 case 17: sum += qlp_coeff[16] * data[i-17]; 1040 case 16: sum += qlp_coeff[15] * data[i-16]; 1041 case 15: sum += qlp_coeff[14] * data[i-15]; 1042 case 14: sum += qlp_coeff[13] * data[i-14]; 1043 case 13: sum += qlp_coeff[12] * data[i-13]; 1044 sum += qlp_coeff[11] * data[i-12]; 1045 sum += qlp_coeff[10] * data[i-11]; 1046 sum += qlp_coeff[ 9] * data[i-10]; 1047 sum += qlp_coeff[ 8] * data[i- 9]; 1048 sum += qlp_coeff[ 7] * data[i- 8]; 1049 sum += qlp_coeff[ 6] * data[i- 7]; 1050 sum += qlp_coeff[ 5] * data[i- 6]; 1051 sum += qlp_coeff[ 4] * data[i- 5]; 1052 sum += qlp_coeff[ 3] * data[i- 4]; 1053 sum += qlp_coeff[ 2] * data[i- 3]; 1054 sum += qlp_coeff[ 1] * data[i- 2]; 1055 sum += qlp_coeff[ 0] * data[i- 1]; 1056 } 1057 data[i] = residual[i] + (sum >> lp_quantization); 1058 } 1059 } 1060} 1061#endif 1062 1063void FLAC__lpc_restore_signal_wide(const FLAC__int32 residual[], unsigned data_len, const FLAC__int32 qlp_coeff[], unsigned order, int lp_quantization, FLAC__int32 data[]) 1064#if defined(FLAC__OVERFLOW_DETECT) || !defined(FLAC__LPC_UNROLLED_FILTER_LOOPS) 1065{ 1066 unsigned i, j; 1067 FLAC__int64 sum; 1068 const FLAC__int32 *r = residual, *history; 1069 1070#ifdef FLAC__OVERFLOW_DETECT_VERBOSE 1071 fprintf(stderr,"FLAC__lpc_restore_signal_wide: data_len=%d, order=%u, lpq=%d",data_len,order,lp_quantization); 1072 for(i=0;i<order;i++) 1073 fprintf(stderr,", q[%u]=%d",i,qlp_coeff[i]); 1074 fprintf(stderr,"\n"); 1075#endif 1076 FLAC__ASSERT(order > 0); 1077 1078 for(i = 0; i < data_len; i++) { 1079 sum = 0; 1080 history = data; 1081 for(j = 0; j < order; j++) 1082 sum += (FLAC__int64)qlp_coeff[j] * (FLAC__int64)(*(--history)); 1083 if(FLAC__bitmath_silog2_wide(sum >> lp_quantization) > 32) { 1084#ifdef FLAC__OVERFLOW_DETECT_VERBOSE 1085#ifdef _MSC_VER 1086 fprintf(stderr,"FLAC__lpc_restore_signal_wide: OVERFLOW, i=%u, sum=%I64d\n", i, sum >> lp_quantization); 1087#else 1088 fprintf(stderr,"FLAC__lpc_restore_signal_wide: OVERFLOW, i=%u, sum=%lld\n", i, (long long)(sum >> lp_quantization)); 1089#endif 1090#endif 1091 break; 1092 } 1093 if(FLAC__bitmath_silog2_wide((FLAC__int64)(*r) + (sum >> lp_quantization)) > 32) { 1094#ifdef FLAC__OVERFLOW_DETECT_VERBOSE 1095#ifdef _MSC_VER 1096 fprintf(stderr,"FLAC__lpc_restore_signal_wide: OVERFLOW, i=%u, residual=%d, sum=%I64d, data=%I64d\n", i, *r, sum >> lp_quantization, (FLAC__int64)(*r) + (sum >> lp_quantization)); 1097#else 1098 fprintf(stderr,"FLAC__lpc_restore_signal_wide: OVERFLOW, i=%u, residual=%d, sum=%lld, data=%lld\n", i, *r, (long long)(sum >> lp_quantization), (long long)((FLAC__int64)(*r) + (sum >> lp_quantization))); 1099#endif 1100#endif 1101 break; 1102 } 1103 *(data++) = *(r++) + (FLAC__int32)(sum >> lp_quantization); 1104 } 1105} 1106#else /* fully unrolled version for normal use */ 1107{ 1108 int i; 1109 FLAC__int64 sum; 1110 1111 FLAC__ASSERT(order > 0); 1112 FLAC__ASSERT(order <= 32); 1113 1114 /* 1115 * We do unique versions up to 12th order since that's the subset limit. 1116 * Also they are roughly ordered to match frequency of occurrence to 1117 * minimize branching. 1118 */ 1119 if(order <= 12) { 1120 if(order > 8) { 1121 if(order > 10) { 1122 if(order == 12) { 1123 for(i = 0; i < (int)data_len; i++) { 1124 sum = 0; 1125 sum += qlp_coeff[11] * (FLAC__int64)data[i-12]; 1126 sum += qlp_coeff[10] * (FLAC__int64)data[i-11]; 1127 sum += qlp_coeff[9] * (FLAC__int64)data[i-10]; 1128 sum += qlp_coeff[8] * (FLAC__int64)data[i-9]; 1129 sum += qlp_coeff[7] * (FLAC__int64)data[i-8]; 1130 sum += qlp_coeff[6] * (FLAC__int64)data[i-7]; 1131 sum += qlp_coeff[5] * (FLAC__int64)data[i-6]; 1132 sum += qlp_coeff[4] * (FLAC__int64)data[i-5]; 1133 sum += qlp_coeff[3] * (FLAC__int64)data[i-4]; 1134 sum += qlp_coeff[2] * (FLAC__int64)data[i-3]; 1135 sum += qlp_coeff[1] * (FLAC__int64)data[i-2]; 1136 sum += qlp_coeff[0] * (FLAC__int64)data[i-1]; 1137 data[i] = residual[i] + (FLAC__int32)(sum >> lp_quantization); 1138 } 1139 } 1140 else { /* order == 11 */ 1141 for(i = 0; i < (int)data_len; i++) { 1142 sum = 0; 1143 sum += qlp_coeff[10] * (FLAC__int64)data[i-11]; 1144 sum += qlp_coeff[9] * (FLAC__int64)data[i-10]; 1145 sum += qlp_coeff[8] * (FLAC__int64)data[i-9]; 1146 sum += qlp_coeff[7] * (FLAC__int64)data[i-8]; 1147 sum += qlp_coeff[6] * (FLAC__int64)data[i-7]; 1148 sum += qlp_coeff[5] * (FLAC__int64)data[i-6]; 1149 sum += qlp_coeff[4] * (FLAC__int64)data[i-5]; 1150 sum += qlp_coeff[3] * (FLAC__int64)data[i-4]; 1151 sum += qlp_coeff[2] * (FLAC__int64)data[i-3]; 1152 sum += qlp_coeff[1] * (FLAC__int64)data[i-2]; 1153 sum += qlp_coeff[0] * (FLAC__int64)data[i-1]; 1154 data[i] = residual[i] + (FLAC__int32)(sum >> lp_quantization); 1155 } 1156 } 1157 } 1158 else { 1159 if(order == 10) { 1160 for(i = 0; i < (int)data_len; i++) { 1161 sum = 0; 1162 sum += qlp_coeff[9] * (FLAC__int64)data[i-10]; 1163 sum += qlp_coeff[8] * (FLAC__int64)data[i-9]; 1164 sum += qlp_coeff[7] * (FLAC__int64)data[i-8]; 1165 sum += qlp_coeff[6] * (FLAC__int64)data[i-7]; 1166 sum += qlp_coeff[5] * (FLAC__int64)data[i-6]; 1167 sum += qlp_coeff[4] * (FLAC__int64)data[i-5]; 1168 sum += qlp_coeff[3] * (FLAC__int64)data[i-4]; 1169 sum += qlp_coeff[2] * (FLAC__int64)data[i-3]; 1170 sum += qlp_coeff[1] * (FLAC__int64)data[i-2]; 1171 sum += qlp_coeff[0] * (FLAC__int64)data[i-1]; 1172 data[i] = residual[i] + (FLAC__int32)(sum >> lp_quantization); 1173 } 1174 } 1175 else { /* order == 9 */ 1176 for(i = 0; i < (int)data_len; i++) { 1177 sum = 0; 1178 sum += qlp_coeff[8] * (FLAC__int64)data[i-9]; 1179 sum += qlp_coeff[7] * (FLAC__int64)data[i-8]; 1180 sum += qlp_coeff[6] * (FLAC__int64)data[i-7]; 1181 sum += qlp_coeff[5] * (FLAC__int64)data[i-6]; 1182 sum += qlp_coeff[4] * (FLAC__int64)data[i-5]; 1183 sum += qlp_coeff[3] * (FLAC__int64)data[i-4]; 1184 sum += qlp_coeff[2] * (FLAC__int64)data[i-3]; 1185 sum += qlp_coeff[1] * (FLAC__int64)data[i-2]; 1186 sum += qlp_coeff[0] * (FLAC__int64)data[i-1]; 1187 data[i] = residual[i] + (FLAC__int32)(sum >> lp_quantization); 1188 } 1189 } 1190 } 1191 } 1192 else if(order > 4) { 1193 if(order > 6) { 1194 if(order == 8) { 1195 for(i = 0; i < (int)data_len; i++) { 1196 sum = 0; 1197 sum += qlp_coeff[7] * (FLAC__int64)data[i-8]; 1198 sum += qlp_coeff[6] * (FLAC__int64)data[i-7]; 1199 sum += qlp_coeff[5] * (FLAC__int64)data[i-6]; 1200 sum += qlp_coeff[4] * (FLAC__int64)data[i-5]; 1201 sum += qlp_coeff[3] * (FLAC__int64)data[i-4]; 1202 sum += qlp_coeff[2] * (FLAC__int64)data[i-3]; 1203 sum += qlp_coeff[1] * (FLAC__int64)data[i-2]; 1204 sum += qlp_coeff[0] * (FLAC__int64)data[i-1]; 1205 data[i] = residual[i] + (FLAC__int32)(sum >> lp_quantization); 1206 } 1207 } 1208 else { /* order == 7 */ 1209 for(i = 0; i < (int)data_len; i++) { 1210 sum = 0; 1211 sum += qlp_coeff[6] * (FLAC__int64)data[i-7]; 1212 sum += qlp_coeff[5] * (FLAC__int64)data[i-6]; 1213 sum += qlp_coeff[4] * (FLAC__int64)data[i-5]; 1214 sum += qlp_coeff[3] * (FLAC__int64)data[i-4]; 1215 sum += qlp_coeff[2] * (FLAC__int64)data[i-3]; 1216 sum += qlp_coeff[1] * (FLAC__int64)data[i-2]; 1217 sum += qlp_coeff[0] * (FLAC__int64)data[i-1]; 1218 data[i] = residual[i] + (FLAC__int32)(sum >> lp_quantization); 1219 } 1220 } 1221 } 1222 else { 1223 if(order == 6) { 1224 for(i = 0; i < (int)data_len; i++) { 1225 sum = 0; 1226 sum += qlp_coeff[5] * (FLAC__int64)data[i-6]; 1227 sum += qlp_coeff[4] * (FLAC__int64)data[i-5]; 1228 sum += qlp_coeff[3] * (FLAC__int64)data[i-4]; 1229 sum += qlp_coeff[2] * (FLAC__int64)data[i-3]; 1230 sum += qlp_coeff[1] * (FLAC__int64)data[i-2]; 1231 sum += qlp_coeff[0] * (FLAC__int64)data[i-1]; 1232 data[i] = residual[i] + (FLAC__int32)(sum >> lp_quantization); 1233 } 1234 } 1235 else { /* order == 5 */ 1236 for(i = 0; i < (int)data_len; i++) { 1237 sum = 0; 1238 sum += qlp_coeff[4] * (FLAC__int64)data[i-5]; 1239 sum += qlp_coeff[3] * (FLAC__int64)data[i-4]; 1240 sum += qlp_coeff[2] * (FLAC__int64)data[i-3]; 1241 sum += qlp_coeff[1] * (FLAC__int64)data[i-2]; 1242 sum += qlp_coeff[0] * (FLAC__int64)data[i-1]; 1243 data[i] = residual[i] + (FLAC__int32)(sum >> lp_quantization); 1244 } 1245 } 1246 } 1247 } 1248 else { 1249 if(order > 2) { 1250 if(order == 4) { 1251 for(i = 0; i < (int)data_len; i++) { 1252 sum = 0; 1253 sum += qlp_coeff[3] * (FLAC__int64)data[i-4]; 1254 sum += qlp_coeff[2] * (FLAC__int64)data[i-3]; 1255 sum += qlp_coeff[1] * (FLAC__int64)data[i-2]; 1256 sum += qlp_coeff[0] * (FLAC__int64)data[i-1]; 1257 data[i] = residual[i] + (FLAC__int32)(sum >> lp_quantization); 1258 } 1259 } 1260 else { /* order == 3 */ 1261 for(i = 0; i < (int)data_len; i++) { 1262 sum = 0; 1263 sum += qlp_coeff[2] * (FLAC__int64)data[i-3]; 1264 sum += qlp_coeff[1] * (FLAC__int64)data[i-2]; 1265 sum += qlp_coeff[0] * (FLAC__int64)data[i-1]; 1266 data[i] = residual[i] + (FLAC__int32)(sum >> lp_quantization); 1267 } 1268 } 1269 } 1270 else { 1271 if(order == 2) { 1272 for(i = 0; i < (int)data_len; i++) { 1273 sum = 0; 1274 sum += qlp_coeff[1] * (FLAC__int64)data[i-2]; 1275 sum += qlp_coeff[0] * (FLAC__int64)data[i-1]; 1276 data[i] = residual[i] + (FLAC__int32)(sum >> lp_quantization); 1277 } 1278 } 1279 else { /* order == 1 */ 1280 for(i = 0; i < (int)data_len; i++) 1281 data[i] = residual[i] + (FLAC__int32)((qlp_coeff[0] * (FLAC__int64)data[i-1]) >> lp_quantization); 1282 } 1283 } 1284 } 1285 } 1286 else { /* order > 12 */ 1287 for(i = 0; i < (int)data_len; i++) { 1288 sum = 0; 1289 switch(order) { 1290 case 32: sum += qlp_coeff[31] * (FLAC__int64)data[i-32]; 1291 case 31: sum += qlp_coeff[30] * (FLAC__int64)data[i-31]; 1292 case 30: sum += qlp_coeff[29] * (FLAC__int64)data[i-30]; 1293 case 29: sum += qlp_coeff[28] * (FLAC__int64)data[i-29]; 1294 case 28: sum += qlp_coeff[27] * (FLAC__int64)data[i-28]; 1295 case 27: sum += qlp_coeff[26] * (FLAC__int64)data[i-27]; 1296 case 26: sum += qlp_coeff[25] * (FLAC__int64)data[i-26]; 1297 case 25: sum += qlp_coeff[24] * (FLAC__int64)data[i-25]; 1298 case 24: sum += qlp_coeff[23] * (FLAC__int64)data[i-24]; 1299 case 23: sum += qlp_coeff[22] * (FLAC__int64)data[i-23]; 1300 case 22: sum += qlp_coeff[21] * (FLAC__int64)data[i-22]; 1301 case 21: sum += qlp_coeff[20] * (FLAC__int64)data[i-21]; 1302 case 20: sum += qlp_coeff[19] * (FLAC__int64)data[i-20]; 1303 case 19: sum += qlp_coeff[18] * (FLAC__int64)data[i-19]; 1304 case 18: sum += qlp_coeff[17] * (FLAC__int64)data[i-18]; 1305 case 17: sum += qlp_coeff[16] * (FLAC__int64)data[i-17]; 1306 case 16: sum += qlp_coeff[15] * (FLAC__int64)data[i-16]; 1307 case 15: sum += qlp_coeff[14] * (FLAC__int64)data[i-15]; 1308 case 14: sum += qlp_coeff[13] * (FLAC__int64)data[i-14]; 1309 case 13: sum += qlp_coeff[12] * (FLAC__int64)data[i-13]; 1310 sum += qlp_coeff[11] * (FLAC__int64)data[i-12]; 1311 sum += qlp_coeff[10] * (FLAC__int64)data[i-11]; 1312 sum += qlp_coeff[ 9] * (FLAC__int64)data[i-10]; 1313 sum += qlp_coeff[ 8] * (FLAC__int64)data[i- 9]; 1314 sum += qlp_coeff[ 7] * (FLAC__int64)data[i- 8]; 1315 sum += qlp_coeff[ 6] * (FLAC__int64)data[i- 7]; 1316 sum += qlp_coeff[ 5] * (FLAC__int64)data[i- 6]; 1317 sum += qlp_coeff[ 4] * (FLAC__int64)data[i- 5]; 1318 sum += qlp_coeff[ 3] * (FLAC__int64)data[i- 4]; 1319 sum += qlp_coeff[ 2] * (FLAC__int64)data[i- 3]; 1320 sum += qlp_coeff[ 1] * (FLAC__int64)data[i- 2]; 1321 sum += qlp_coeff[ 0] * (FLAC__int64)data[i- 1]; 1322 } 1323 data[i] = residual[i] + (FLAC__int32)(sum >> lp_quantization); 1324 } 1325 } 1326} 1327#endif 1328 1329#ifndef FLAC__INTEGER_ONLY_LIBRARY 1330 1331FLAC__double FLAC__lpc_compute_expected_bits_per_residual_sample(FLAC__double lpc_error, unsigned total_samples) 1332{ 1333 FLAC__double error_scale; 1334 1335 FLAC__ASSERT(total_samples > 0); 1336 1337 error_scale = 0.5 * M_LN2 * M_LN2 / (FLAC__double)total_samples; 1338 1339 return FLAC__lpc_compute_expected_bits_per_residual_sample_with_error_scale(lpc_error, error_scale); 1340} 1341 1342FLAC__double FLAC__lpc_compute_expected_bits_per_residual_sample_with_error_scale(FLAC__double lpc_error, FLAC__double error_scale) 1343{ 1344 if(lpc_error > 0.0) { 1345 FLAC__double bps = (FLAC__double)0.5 * log(error_scale * lpc_error) / M_LN2; 1346 if(bps >= 0.0) 1347 return bps; 1348 else 1349 return 0.0; 1350 } 1351 else if(lpc_error < 0.0) { /* error should not be negative but can happen due to inadequate floating-point resolution */ 1352 return 1e32; 1353 } 1354 else { 1355 return 0.0; 1356 } 1357} 1358 1359unsigned FLAC__lpc_compute_best_order(const FLAC__double lpc_error[], unsigned max_order, unsigned total_samples, unsigned overhead_bits_per_order) 1360{ 1361 unsigned order, index, best_index; /* 'index' the index into lpc_error; index==order-1 since lpc_error[0] is for order==1, lpc_error[1] is for order==2, etc */ 1362 FLAC__double bits, best_bits, error_scale; 1363 1364 FLAC__ASSERT(max_order > 0); 1365 FLAC__ASSERT(total_samples > 0); 1366 1367 error_scale = 0.5 * M_LN2 * M_LN2 / (FLAC__double)total_samples; 1368 1369 best_index = 0; 1370 best_bits = (unsigned)(-1); 1371 1372 for(index = 0, order = 1; index < max_order; index++, order++) { 1373 bits = FLAC__lpc_compute_expected_bits_per_residual_sample_with_error_scale(lpc_error[index], error_scale) * (FLAC__double)(total_samples - order) + (FLAC__double)(order * overhead_bits_per_order); 1374 if(bits < best_bits) { 1375 best_index = index; 1376 best_bits = bits; 1377 } 1378 } 1379 1380 return best_index+1; /* +1 since index of lpc_error[] is order-1 */ 1381} 1382 1383#endif /* !defined FLAC__INTEGER_ONLY_LIBRARY */ 1384