stop using lround() for now
[platform/upstream/flac.git] / src / libFLAC / lpc.c
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 void FLAC__lpc_window_data(const FLAC__real in[], const FLAC__real window[], FLAC__real out[], unsigned data_len)
53 {
54         unsigned i;
55         for(i = 0; i < data_len; i++)
56                 out[i] = in[i] * window[i];
57 }
58
59 void FLAC__lpc_compute_autocorrelation(const FLAC__real data[], unsigned data_len, unsigned lag, FLAC__real autoc[])
60 {
61         /* a readable, but slower, version */
62 #if 0
63         FLAC__real d;
64         unsigned i;
65
66         FLAC__ASSERT(lag > 0);
67         FLAC__ASSERT(lag <= data_len);
68
69         /*
70          * Technically we should subtract the mean first like so:
71          *   for(i = 0; i < data_len; i++)
72          *     data[i] -= mean;
73          * but it appears not to make enough of a difference to matter, and
74          * most signals are already closely centered around zero
75          */
76         while(lag--) {
77                 for(i = lag, d = 0.0; i < data_len; i++)
78                         d += data[i] * data[i - lag];
79                 autoc[lag] = d;
80         }
81 #endif
82
83         /*
84          * this version tends to run faster because of better data locality
85          * ('data_len' is usually much larger than 'lag')
86          */
87         FLAC__real d;
88         unsigned sample, coeff;
89         const unsigned limit = data_len - lag;
90
91         FLAC__ASSERT(lag > 0);
92         FLAC__ASSERT(lag <= data_len);
93
94         for(coeff = 0; coeff < lag; coeff++)
95                 autoc[coeff] = 0.0;
96         for(sample = 0; sample <= limit; sample++) {
97                 d = data[sample];
98                 for(coeff = 0; coeff < lag; coeff++)
99                         autoc[coeff] += d * data[sample+coeff];
100         }
101         for(; sample < data_len; sample++) {
102                 d = data[sample];
103                 for(coeff = 0; coeff < data_len - sample; coeff++)
104                         autoc[coeff] += d * data[sample+coeff];
105         }
106 }
107
108 void FLAC__lpc_compute_lp_coefficients(const FLAC__real autoc[], unsigned *max_order, FLAC__real lp_coeff[][FLAC__MAX_LPC_ORDER], FLAC__double error[])
109 {
110         unsigned i, j;
111         FLAC__double r, err, ref[FLAC__MAX_LPC_ORDER], lpc[FLAC__MAX_LPC_ORDER];
112
113         FLAC__ASSERT(0 != max_order);
114         FLAC__ASSERT(0 < *max_order);
115         FLAC__ASSERT(*max_order <= FLAC__MAX_LPC_ORDER);
116         FLAC__ASSERT(autoc[0] != 0.0);
117
118         err = autoc[0];
119
120         for(i = 0; i < *max_order; i++) {
121                 /* Sum up this iteration's reflection coefficient. */
122                 r = -autoc[i+1];
123                 for(j = 0; j < i; j++)
124                         r -= lpc[j] * autoc[i-j];
125                 ref[i] = (r/=err);
126
127                 /* Update LPC coefficients and total error. */
128                 lpc[i]=r;
129                 for(j = 0; j < (i>>1); j++) {
130                         FLAC__double tmp = lpc[j];
131                         lpc[j] += r * lpc[i-1-j];
132                         lpc[i-1-j] += r * tmp;
133                 }
134                 if(i & 1)
135                         lpc[j] += lpc[j] * r;
136
137                 err *= (1.0 - r * r);
138
139                 /* save this order */
140                 for(j = 0; j <= i; j++)
141                         lp_coeff[i][j] = (FLAC__real)(-lpc[j]); /* negate FIR filter coeff to get predictor coeff */
142                 error[i] = err;
143
144                 /* see SF bug #1601812 http://sourceforge.net/tracker/index.php?func=detail&aid=1601812&group_id=13478&atid=113478 */
145                 if(err == 0.0) {
146                         *max_order = i+1;
147                         return;
148                 }
149         }
150 }
151
152 int FLAC__lpc_quantize_coefficients(const FLAC__real lp_coeff[], unsigned order, unsigned precision, FLAC__int32 qlp_coeff[], int *shift)
153 {
154         unsigned i;
155         FLAC__double cmax;
156         FLAC__int32 qmax, qmin;
157
158         FLAC__ASSERT(precision > 0);
159         FLAC__ASSERT(precision >= FLAC__MIN_QLP_COEFF_PRECISION);
160
161         /* drop one bit for the sign; from here on out we consider only |lp_coeff[i]| */
162         precision--;
163         qmax = 1 << precision;
164         qmin = -qmax;
165         qmax--;
166
167         /* calc cmax = max( |lp_coeff[i]| ) */
168         cmax = 0.0;
169         for(i = 0; i < order; i++) {
170                 const FLAC__double d = fabs(lp_coeff[i]);
171                 if(d > cmax)
172                         cmax = d;
173         }
174
175         if(cmax <= 0.0) {
176                 /* => coefficients are all 0, which means our constant-detect didn't work */
177                 return 2;
178         }
179         else {
180                 const int max_shiftlimit = (1 << (FLAC__SUBFRAME_LPC_QLP_SHIFT_LEN-1)) - 1;
181                 const int min_shiftlimit = -max_shiftlimit - 1;
182                 int log2cmax;
183
184                 (void)frexp(cmax, &log2cmax);
185                 log2cmax--;
186                 *shift = (int)precision - log2cmax - 1;
187
188                 if(*shift > max_shiftlimit)
189                         *shift = max_shiftlimit;
190                 else if(*shift < min_shiftlimit)
191                         return 1;
192         }
193
194         if(*shift >= 0) {
195                 FLAC__double error = 0.0;
196                 FLAC__int32 q;
197                 for(i = 0; i < order; i++) {
198                         error += lp_coeff[i] * (1 << *shift);
199 #if 1 /* unfortunately lround() is C99 */
200                         if(error >= 0.0)
201                                 q = (FLAC__int32)(error + 0.5);
202                         else
203                                 q = (FLAC__int32)(error - 0.5);
204 #else
205                         q = lround(error);
206 #endif
207 #ifdef FLAC__OVERFLOW_DETECT
208                         if(q > qmax+1) /* we expect q==qmax+1 occasionally due to rounding */
209                                 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]);
210                         else if(q < qmin)
211                                 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]);
212 #endif
213                         if(q > qmax)
214                                 q = qmax;
215                         else if(q < qmin)
216                                 q = qmin;
217                         error -= q;
218                         qlp_coeff[i] = q;
219                 }
220         }
221         /* negative shift is very rare but due to design flaw, negative shift is
222          * a NOP in the decoder, so it must be handled specially by scaling down
223          * coeffs
224          */
225         else {
226                 const int nshift = -(*shift);
227                 FLAC__double error = 0.0;
228                 FLAC__int32 q;
229 #ifdef DEBUG
230                 fprintf(stderr,"FLAC__lpc_quantize_coefficients: negative shift=%d order=%u cmax=%f\n", *shift, order, cmax);
231 #endif
232                 for(i = 0; i < order; i++) {
233                         error += lp_coeff[i] / (1 << nshift);
234 #if 1 /* unfortunately lround() is C99 */
235                         if(error >= 0.0)
236                                 q = (FLAC__int32)(error + 0.5);
237                         else
238                                 q = (FLAC__int32)(error - 0.5);
239 #else
240                         q = lround(error);
241 #endif
242 #ifdef FLAC__OVERFLOW_DETECT
243                         if(q > qmax+1) /* we expect q==qmax+1 occasionally due to rounding */
244                                 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]);
245                         else if(q < qmin)
246                                 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]);
247 #endif
248                         if(q > qmax)
249                                 q = qmax;
250                         else if(q < qmin)
251                                 q = qmin;
252                         error -= q;
253                         qlp_coeff[i] = q;
254                 }
255                 *shift = 0;
256         }
257
258         return 0;
259 }
260
261 void 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[])
262 {
263 #ifdef FLAC__OVERFLOW_DETECT
264         FLAC__int64 sumo;
265 #endif
266         unsigned i, j;
267         FLAC__int32 sum;
268         const FLAC__int32 *history;
269
270 #ifdef FLAC__OVERFLOW_DETECT_VERBOSE
271         fprintf(stderr,"FLAC__lpc_compute_residual_from_qlp_coefficients: data_len=%d, order=%u, lpq=%d",data_len,order,lp_quantization);
272         for(i=0;i<order;i++)
273                 fprintf(stderr,", q[%u]=%d",i,qlp_coeff[i]);
274         fprintf(stderr,"\n");
275 #endif
276         FLAC__ASSERT(order > 0);
277
278         for(i = 0; i < data_len; i++) {
279 #ifdef FLAC__OVERFLOW_DETECT
280                 sumo = 0;
281 #endif
282                 sum = 0;
283                 history = data;
284                 for(j = 0; j < order; j++) {
285                         sum += qlp_coeff[j] * (*(--history));
286 #ifdef FLAC__OVERFLOW_DETECT
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 #endif
296                 }
297                 *(residual++) = *(data++) - (sum >> lp_quantization);
298         }
299
300         /* Here's a slower but clearer version:
301         for(i = 0; i < data_len; i++) {
302                 sum = 0;
303                 for(j = 0; j < order; j++)
304                         sum += qlp_coeff[j] * data[i-j-1];
305                 residual[i] = data[i] - (sum >> lp_quantization);
306         }
307         */
308 }
309
310 void 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[])
311 {
312         unsigned i, j;
313         FLAC__int64 sum;
314         const FLAC__int32 *history;
315
316 #ifdef FLAC__OVERFLOW_DETECT_VERBOSE
317         fprintf(stderr,"FLAC__lpc_compute_residual_from_qlp_coefficients_wide: data_len=%d, order=%u, lpq=%d",data_len,order,lp_quantization);
318         for(i=0;i<order;i++)
319                 fprintf(stderr,", q[%u]=%d",i,qlp_coeff[i]);
320         fprintf(stderr,"\n");
321 #endif
322         FLAC__ASSERT(order > 0);
323
324         for(i = 0; i < data_len; i++) {
325                 sum = 0;
326                 history = data;
327                 for(j = 0; j < order; j++)
328                         sum += (FLAC__int64)qlp_coeff[j] * (FLAC__int64)(*(--history));
329 #ifdef FLAC__OVERFLOW_DETECT
330                 if(FLAC__bitmath_silog2_wide(sum >> lp_quantization) > 32) {
331 #if defined _MSC_VER
332                         fprintf(stderr,"FLAC__lpc_compute_residual_from_qlp_coefficients_wide: OVERFLOW, i=%u, sum=%I64d\n", i, sum >> lp_quantization);
333 #else
334                         fprintf(stderr,"FLAC__lpc_compute_residual_from_qlp_coefficients_wide: OVERFLOW, i=%u, sum=%lld\n", i, (long long)(sum >> lp_quantization));
335 #endif
336                         break;
337                 }
338                 if(FLAC__bitmath_silog2_wide((FLAC__int64)(*data) - (sum >> lp_quantization)) > 32) {
339 #if defined _MSC_VER
340                         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));
341 #else
342                         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)));
343 #endif
344                         break;
345                 }
346 #endif
347                 *(residual++) = *(data++) - (FLAC__int32)(sum >> lp_quantization);
348         }
349 }
350
351 #endif /* !defined FLAC__INTEGER_ONLY_LIBRARY */
352
353 void FLAC__lpc_restore_signal(const FLAC__int32 residual[], unsigned data_len, const FLAC__int32 qlp_coeff[], unsigned order, int lp_quantization, FLAC__int32 data[])
354 {
355 #ifdef FLAC__OVERFLOW_DETECT
356         FLAC__int64 sumo;
357 #endif
358         unsigned i, j;
359         FLAC__int32 sum;
360         const FLAC__int32 *r = residual, *history;
361
362 #ifdef FLAC__OVERFLOW_DETECT_VERBOSE
363         fprintf(stderr,"FLAC__lpc_restore_signal: data_len=%d, order=%u, lpq=%d",data_len,order,lp_quantization);
364         for(i=0;i<order;i++)
365                 fprintf(stderr,", q[%u]=%d",i,qlp_coeff[i]);
366         fprintf(stderr,"\n");
367 #endif
368         FLAC__ASSERT(order > 0);
369
370         for(i = 0; i < data_len; i++) {
371 #ifdef FLAC__OVERFLOW_DETECT
372                 sumo = 0;
373 #endif
374                 sum = 0;
375                 history = data;
376                 for(j = 0; j < order; j++) {
377                         sum += qlp_coeff[j] * (*(--history));
378 #ifdef FLAC__OVERFLOW_DETECT
379                         sumo += (FLAC__int64)qlp_coeff[j] * (FLAC__int64)(*history);
380 #if defined _MSC_VER
381                         if(sumo > 2147483647I64 || sumo < -2147483648I64)
382                                 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);
383 #else
384                         if(sumo > 2147483647ll || sumo < -2147483648ll)
385                                 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);
386 #endif
387 #endif
388                 }
389                 *(data++) = *(r++) + (sum >> lp_quantization);
390         }
391
392         /* Here's a slower but clearer version:
393         for(i = 0; i < data_len; i++) {
394                 sum = 0;
395                 for(j = 0; j < order; j++)
396                         sum += qlp_coeff[j] * data[i-j-1];
397                 data[i] = residual[i] + (sum >> lp_quantization);
398         }
399         */
400 }
401
402 void 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[])
403 {
404         unsigned i, j;
405         FLAC__int64 sum;
406         const FLAC__int32 *r = residual, *history;
407
408 #ifdef FLAC__OVERFLOW_DETECT_VERBOSE
409         fprintf(stderr,"FLAC__lpc_restore_signal_wide: data_len=%d, order=%u, lpq=%d",data_len,order,lp_quantization);
410         for(i=0;i<order;i++)
411                 fprintf(stderr,", q[%u]=%d",i,qlp_coeff[i]);
412         fprintf(stderr,"\n");
413 #endif
414         FLAC__ASSERT(order > 0);
415
416         for(i = 0; i < data_len; i++) {
417                 sum = 0;
418                 history = data;
419                 for(j = 0; j < order; j++)
420                         sum += (FLAC__int64)qlp_coeff[j] * (FLAC__int64)(*(--history));
421 #ifdef FLAC__OVERFLOW_DETECT
422                 if(FLAC__bitmath_silog2_wide(sum >> lp_quantization) > 32) {
423 #ifdef _MSC_VER
424                         fprintf(stderr,"FLAC__lpc_restore_signal_wide: OVERFLOW, i=%u, sum=%I64d\n", i, sum >> lp_quantization);
425 #else
426                         fprintf(stderr,"FLAC__lpc_restore_signal_wide: OVERFLOW, i=%u, sum=%lld\n", i, (long long)(sum >> lp_quantization));
427 #endif
428                         break;
429                 }
430                 if(FLAC__bitmath_silog2_wide((FLAC__int64)(*r) + (sum >> lp_quantization)) > 32) {
431 #ifdef _MSC_VER
432                         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));
433 #else
434                         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)));
435 #endif
436                         break;
437                 }
438 #endif
439                 *(data++) = *(r++) + (FLAC__int32)(sum >> lp_quantization);
440         }
441 }
442
443 #ifndef FLAC__INTEGER_ONLY_LIBRARY
444
445 FLAC__double FLAC__lpc_compute_expected_bits_per_residual_sample(FLAC__double lpc_error, unsigned total_samples)
446 {
447         FLAC__double error_scale;
448
449         FLAC__ASSERT(total_samples > 0);
450
451         error_scale = 0.5 * M_LN2 * M_LN2 / (FLAC__double)total_samples;
452
453         return FLAC__lpc_compute_expected_bits_per_residual_sample_with_error_scale(lpc_error, error_scale);
454 }
455
456 FLAC__double FLAC__lpc_compute_expected_bits_per_residual_sample_with_error_scale(FLAC__double lpc_error, FLAC__double error_scale)
457 {
458         if(lpc_error > 0.0) {
459                 FLAC__double bps = (FLAC__double)0.5 * log(error_scale * lpc_error) / M_LN2;
460                 if(bps >= 0.0)
461                         return bps;
462                 else
463                         return 0.0;
464         }
465         else if(lpc_error < 0.0) { /* error should not be negative but can happen due to inadequate floating-point resolution */
466                 return 1e32;
467         }
468         else {
469                 return 0.0;
470         }
471 }
472
473 unsigned FLAC__lpc_compute_best_order(const FLAC__double lpc_error[], unsigned max_order, unsigned total_samples, unsigned overhead_bits_per_order)
474 {
475         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 */
476         FLAC__double bits, best_bits, error_scale;
477
478         FLAC__ASSERT(max_order > 0);
479         FLAC__ASSERT(total_samples > 0);
480
481         error_scale = 0.5 * M_LN2 * M_LN2 / (FLAC__double)total_samples;
482
483         best_index = 0;
484         best_bits = (unsigned)(-1);
485
486         for(index = 0, order = 1; index < max_order; index++, order++) {
487                 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);
488                 if(bits < best_bits) {
489                         best_index = index;
490                         best_bits = bits;
491                 }
492         }
493
494         return best_index+1; /* +1 since index of lpc_error[] is order-1 */
495 }
496
497 #endif /* !defined FLAC__INTEGER_ONLY_LIBRARY */