← Back to branch summary

~ubuntu-branches/ubuntu/precise/flac/precise-updates

« back to all changes in this revision

Viewing changes to src/libFLAC/lpc.c

  • Committer: Bazaar Package Importer
  • Author(s): Joshua Kwan
  • Date: 2007-05-29 22:56:36 UTC
  • mfrom: (1.1.4 upstream)
  • Revision ID: james.westby@ubuntu.com-20070529225636-ljeff8xxip09qaap
Tags: 1.1.4-1
http://bugs.debian.org/405167
http://bugs.debian.org/411311
* New upstream release. closes: #405167, #411311
  - libOggFLAC and libOggFLAC++ have been merged into libFLAC, so
    remove their corresponding packages.
  - Because of the API changes required to effect the above, there has
    been yet another soname bump. libflac7 -> libflac8 and
    libflac++5 -> libflac++6. Emails have been dispatched to the
    maintainers of dependent packages.
* Some notes on patches that were removed:
  - 02_stdin_stdout, 06_manpage_mention_utf8_convert: merged upstream
  - 08_manpage_warnings: Upstream has changed the manpage so it defintely
    can't fit in in 80 cols, so just forget about it. We'll live.
  - 05_eof_warnings_are_errors: Upstream decided to add a -w option to
    flac to treat all warnings as errors. I am going to defer to that
    for now, but if people think it's stupid let me know and I'll port
    the patch forward.
  - 04_stack_smasher: was a backport from 1.1.3, so it's obsolete.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/libFLAC/lpc.c
1
1
/* libFLAC - Free Lossless Audio Codec library
2
 
 * Copyright (C) 2000,2001,2002,2003,2004,2005  Josh Coalson
 
2
 * Copyright (C) 2000,2001,2002,2003,2004,2005,2006,2007  Josh Coalson
3
3
 *
4
4
 * Redistribution and use in source and binary forms, with or without
5
5
 * modification, are permitted provided that the following conditions
29
29
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
30
30
 */
31
31
 
 
32
#if HAVE_CONFIG_H
 
33
#  include <config.h>
 
34
#endif
 
35
 
32
36
#include <math.h>
33
37
#include "FLAC/assert.h"
34
38
#include "FLAC/format.h"
45
49
#define M_LN2 0.69314718055994530942
46
50
#endif
47
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
 
48
59
void FLAC__lpc_compute_autocorrelation(const FLAC__real data[], unsigned data_len, unsigned lag, FLAC__real autoc[])
49
60
{
50
61
        /* a readable, but slower, version */
55
66
        FLAC__ASSERT(lag > 0);
56
67
        FLAC__ASSERT(lag <= data_len);
57
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
         */
58
76
        while(lag--) {
59
77
                for(i = lag, d = 0.0; i < data_len; i++)
60
78
                        d += data[i] * data[i - lag];
87
105
        }
88
106
}
89
107
 
90
 
void FLAC__lpc_compute_lp_coefficients(const FLAC__real autoc[], unsigned max_order, FLAC__real lp_coeff[][FLAC__MAX_LPC_ORDER], FLAC__double error[])
 
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[])
91
109
{
92
110
        unsigned i, j;
93
111
        FLAC__double r, err, ref[FLAC__MAX_LPC_ORDER], lpc[FLAC__MAX_LPC_ORDER];
94
112
 
95
 
        FLAC__ASSERT(0 < max_order);
96
 
        FLAC__ASSERT(max_order <= FLAC__MAX_LPC_ORDER);
 
113
        FLAC__ASSERT(0 != max_order);
 
114
        FLAC__ASSERT(0 < *max_order);
 
115
        FLAC__ASSERT(*max_order <= FLAC__MAX_LPC_ORDER);
97
116
        FLAC__ASSERT(autoc[0] != 0.0);
98
117
 
99
118
        err = autoc[0];
100
119
 
101
 
        for(i = 0; i < max_order; i++) {
 
120
        for(i = 0; i < *max_order; i++) {
102
121
                /* Sum up this iteration's reflection coefficient. */
103
122
                r = -autoc[i+1];
104
123
                for(j = 0; j < i; j++)
121
140
                for(j = 0; j <= i; j++)
122
141
                        lp_coeff[i][j] = (FLAC__real)(-lpc[j]); /* negate FIR filter coeff to get predictor coeff */
123
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
                }
124
149
        }
125
150
}
126
151
 
127
152
int FLAC__lpc_quantize_coefficients(const FLAC__real lp_coeff[], unsigned order, unsigned precision, FLAC__int32 qlp_coeff[], int *shift)
128
153
{
129
154
        unsigned i;
130
 
        FLAC__double d, cmax = -1e32;
 
155
        FLAC__double cmax;
131
156
        FLAC__int32 qmax, qmin;
132
 
        const int max_shiftlimit = (1 << (FLAC__SUBFRAME_LPC_QLP_SHIFT_LEN-1)) - 1;
133
 
        const int min_shiftlimit = -max_shiftlimit - 1;
134
157
 
135
158
        FLAC__ASSERT(precision > 0);
136
159
        FLAC__ASSERT(precision >= FLAC__MIN_QLP_COEFF_PRECISION);
141
164
        qmin = -qmax;
142
165
        qmax--;
143
166
 
 
167
        /* calc cmax = max( |lp_coeff[i]| ) */
 
168
        cmax = 0.0;
144
169
        for(i = 0; i < order; i++) {
145
 
                if(lp_coeff[i] == 0.0)
146
 
                        continue;
147
 
                d = fabs(lp_coeff[i]);
 
170
                const FLAC__double d = fabs(lp_coeff[i]);
148
171
                if(d > cmax)
149
172
                        cmax = d;
150
173
        }
151
 
redo_it:
 
174
 
152
175
        if(cmax <= 0.0) {
153
176
                /* => coefficients are all 0, which means our constant-detect didn't work */
154
177
                return 2;
155
178
        }
156
179
        else {
 
180
                const int max_shiftlimit = (1 << (FLAC__SUBFRAME_LPC_QLP_SHIFT_LEN-1)) - 1;
 
181
                const int min_shiftlimit = -max_shiftlimit - 1;
157
182
                int log2cmax;
158
183
 
159
184
                (void)frexp(cmax, &log2cmax);
160
185
                log2cmax--;
161
186
                *shift = (int)precision - log2cmax - 1;
162
187
 
163
 
                if(*shift < min_shiftlimit || *shift > max_shiftlimit) {
164
 
#if 0
165
 
                        /*@@@ this does not seem to help at all, but was not extensively tested either: */
166
 
                        if(*shift > max_shiftlimit)
167
 
                                *shift = max_shiftlimit;
168
 
                        else
169
 
#endif
170
 
                                return 1;
171
 
                }
 
188
                if(*shift > max_shiftlimit)
 
189
                        *shift = max_shiftlimit;
 
190
                else if(*shift < min_shiftlimit)
 
191
                        return 1;
172
192
        }
173
193
 
174
194
        if(*shift >= 0) {
 
195
                FLAC__double error = 0.0;
 
196
                FLAC__int32 q;
175
197
                for(i = 0; i < order; i++) {
176
 
                        qlp_coeff[i] = (FLAC__int32)floor((FLAC__double)lp_coeff[i] * (FLAC__double)(1 << *shift));
177
 
 
178
 
                        /* double-check the result */
179
 
                        if(qlp_coeff[i] > qmax || qlp_coeff[i] < qmin) {
 
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
180
207
#ifdef FLAC__OVERFLOW_DETECT
181
 
                                fprintf(stderr,"FLAC__lpc_quantize_coefficients: compensating for overflow, qlp_coeff[%u]=%d, lp_coeff[%u]=%f, cmax=%f, precision=%u, shift=%d, q=%f, f(q)=%f\n", i, qlp_coeff[i], i, lp_coeff[i], cmax, precision, *shift, (FLAC__double)lp_coeff[i] * (FLAC__double)(1 << *shift), floor((FLAC__double)lp_coeff[i] * (FLAC__double)(1 << *shift)));
 
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]);
182
212
#endif
183
 
                                cmax *= 2.0;
184
 
                                goto redo_it;
185
 
                        }
 
213
                        if(q > qmax)
 
214
                                q = qmax;
 
215
                        else if(q < qmin)
 
216
                                q = qmin;
 
217
                        error -= q;
 
218
                        qlp_coeff[i] = q;
186
219
                }
187
220
        }
188
 
        else { /* (*shift < 0) */
 
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 {
189
226
                const int nshift = -(*shift);
 
227
                FLAC__double error = 0.0;
 
228
                FLAC__int32 q;
190
229
#ifdef DEBUG
191
 
                fprintf(stderr,"FLAC__lpc_quantize_coefficients: negative shift = %d\n", *shift);
 
230
                fprintf(stderr,"FLAC__lpc_quantize_coefficients: negative shift=%d order=%u cmax=%f\n", *shift, order, cmax);
192
231
#endif
193
232
                for(i = 0; i < order; i++) {
194
 
                        qlp_coeff[i] = (FLAC__int32)floor((FLAC__double)lp_coeff[i] / (FLAC__double)(1 << nshift));
195
 
 
196
 
                        /* double-check the result */
197
 
                        if(qlp_coeff[i] > qmax || qlp_coeff[i] < qmin) {
 
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
198
242
#ifdef FLAC__OVERFLOW_DETECT
199
 
                                fprintf(stderr,"FLAC__lpc_quantize_coefficients: compensating for overflow, qlp_coeff[%u]=%d, lp_coeff[%u]=%f, cmax=%f, precision=%u, shift=%d, q=%f, f(q)=%f\n", i, qlp_coeff[i], i, lp_coeff[i], cmax, precision, *shift, (FLAC__double)lp_coeff[i] / (FLAC__double)(1 << nshift), floor((FLAC__double)lp_coeff[i] / (FLAC__double)(1 << nshift)));
 
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]);
200
247
#endif
201
 
                                cmax *= 2.0;
202
 
                                goto redo_it;
203
 
                        }
 
248
                        if(q > qmax)
 
249
                                q = qmax;
 
250
                        else if(q < qmin)
 
251
                                q = qmin;
 
252
                        error -= q;
 
253
                        qlp_coeff[i] = q;
204
254
                }
 
255
                *shift = 0;
205
256
        }
206
257
 
207
258
        return 0;
239
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);
240
291
#else
241
292
                        if(sumo > 2147483647ll || sumo < -2147483648ll)
242
 
                                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,sumo);
 
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);
243
294
#endif
244
295
#endif
245
296
                }
277
328
                        sum += (FLAC__int64)qlp_coeff[j] * (FLAC__int64)(*(--history));
278
329
#ifdef FLAC__OVERFLOW_DETECT
279
330
                if(FLAC__bitmath_silog2_wide(sum >> lp_quantization) > 32) {
280
 
                        fprintf(stderr,"FLAC__lpc_compute_residual_from_qlp_coefficients_wide: OVERFLOW, i=%u, sum=%lld\n", i, sum >> lp_quantization);
 
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
281
336
                        break;
282
337
                }
283
338
                if(FLAC__bitmath_silog2_wide((FLAC__int64)(*data) - (sum >> lp_quantization)) > 32) {
284
 
                        fprintf(stderr,"FLAC__lpc_compute_residual_from_qlp_coefficients_wide: OVERFLOW, i=%u, data=%d, sum=%lld, residual=%lld\n", i, *data, sum >> lp_quantization, (FLAC__int64)(*data) - (sum >> lp_quantization));
 
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
285
344
                        break;
286
345
                }
287
346
#endif
298
357
#endif
299
358
        unsigned i, j;
300
359
        FLAC__int32 sum;
301
 
        const FLAC__int32 *history;
 
360
        const FLAC__int32 *r = residual, *history;
302
361
 
303
362
#ifdef FLAC__OVERFLOW_DETECT_VERBOSE
304
363
        fprintf(stderr,"FLAC__lpc_restore_signal: data_len=%d, order=%u, lpq=%d",data_len,order,lp_quantization);
323
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);
324
383
#else
325
384
                        if(sumo > 2147483647ll || sumo < -2147483648ll)
326
 
                                fprintf(stderr,"FLAC__lpc_restore_signal: OVERFLOW, i=%u, j=%u, c=%d, d=%d, sumo=%lld\n",i,j,qlp_coeff[j],*history,sumo);
 
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);
327
386
#endif
328
387
#endif
329
388
                }
330
 
                *(data++) = *(residual++) + (sum >> lp_quantization);
 
389
                *(data++) = *(r++) + (sum >> lp_quantization);
331
390
        }
332
391
 
333
392
        /* Here's a slower but clearer version:
344
403
{
345
404
        unsigned i, j;
346
405
        FLAC__int64 sum;
347
 
        const FLAC__int32 *history;
 
406
        const FLAC__int32 *r = residual, *history;
348
407
 
349
408
#ifdef FLAC__OVERFLOW_DETECT_VERBOSE
350
409
        fprintf(stderr,"FLAC__lpc_restore_signal_wide: data_len=%d, order=%u, lpq=%d",data_len,order,lp_quantization);
361
420
                        sum += (FLAC__int64)qlp_coeff[j] * (FLAC__int64)(*(--history));
362
421
#ifdef FLAC__OVERFLOW_DETECT
363
422
                if(FLAC__bitmath_silog2_wide(sum >> lp_quantization) > 32) {
364
 
                        fprintf(stderr,"FLAC__lpc_restore_signal_wide: OVERFLOW, i=%u, sum=%lld\n", i, sum >> lp_quantization);
365
 
                        break;
366
 
                }
367
 
                if(FLAC__bitmath_silog2_wide((FLAC__int64)(*residual) + (sum >> lp_quantization)) > 32) {
368
 
                        fprintf(stderr,"FLAC__lpc_restore_signal_wide: OVERFLOW, i=%u, residual=%d, sum=%lld, data=%lld\n", i, *residual, sum >> lp_quantization, (FLAC__int64)(*residual) + (sum >> lp_quantization));
369
 
                        break;
370
 
                }
371
 
#endif
372
 
                *(data++) = *(residual++) + (FLAC__int32)(sum >> lp_quantization);
 
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);
373
440
        }
374
441
}
375
442
 
403
470
        }
404
471
}
405
472
 
406
 
unsigned FLAC__lpc_compute_best_order(const FLAC__double lpc_error[], unsigned max_order, unsigned total_samples, unsigned bits_per_signal_sample)
 
473
unsigned FLAC__lpc_compute_best_order(const FLAC__double lpc_error[], unsigned max_order, unsigned total_samples, unsigned overhead_bits_per_order)
407
474
{
408
 
        unsigned order, best_order;
409
 
        FLAC__double best_bits, tmp_bits, error_scale;
 
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;
410
477
 
411
478
        FLAC__ASSERT(max_order > 0);
412
479
        FLAC__ASSERT(total_samples > 0);
413
480
 
414
481
        error_scale = 0.5 * M_LN2 * M_LN2 / (FLAC__double)total_samples;
415
482
 
416
 
        best_order = 0;
417
 
        best_bits = FLAC__lpc_compute_expected_bits_per_residual_sample_with_error_scale(lpc_error[0], error_scale) * (FLAC__double)total_samples;
 
483
        best_index = 0;
 
484
        best_bits = (unsigned)(-1);
418
485
 
419
 
        for(order = 1; order < max_order; order++) {
420
 
                tmp_bits = FLAC__lpc_compute_expected_bits_per_residual_sample_with_error_scale(lpc_error[order], error_scale) * (FLAC__double)(total_samples - order) + (FLAC__double)(order * bits_per_signal_sample);
421
 
                if(tmp_bits < best_bits) {
422
 
                        best_order = order;
423
 
                        best_bits = tmp_bits;
 
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;
424
491
                }
425
492
        }
426
493
 
427
 
        return best_order+1; /* +1 since index of lpc_error[] is order-1 */
 
494
        return best_index+1; /* +1 since index of lpc_error[] is order-1 */
428
495
}
429
496
 
430
497
#endif /* !defined FLAC__INTEGER_ONLY_LIBRARY */