← Back to branch summary

~ubuntu-branches/ubuntu/saucy/flac/saucy

« back to all changes in this revision

Viewing changes to src/libFLAC/lpc.c

  • Committer: Bazaar Package Importer
  • Author(s): Matt Zimmerman
  • Date: 2001-12-10 03:09:22 UTC
  • Revision ID: james.westby@ubuntu.com-20011210030922-0vdtpz6a7mfwefo5
Tags: upstream-1.0.2
Import upstream version 1.0.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/libFLAC/lpc.c
 
1
/* libFLAC - Free Lossless Audio Codec library
 
2
 * Copyright (C) 2000,2001  Josh Coalson
 
3
 *
 
4
 * This library is free software; you can redistribute it and/or
 
5
 * modify it under the terms of the GNU Library General Public
 
6
 * License as published by the Free Software Foundation; either
 
7
 * version 2 of the License, or (at your option) any later version.
 
8
 *
 
9
 * This library is distributed in the hope that it will be useful,
 
10
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 
11
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 
12
 * Library General Public License for more details.
 
13
 *
 
14
 * You should have received a copy of the GNU Library General Public
 
15
 * License along with this library; if not, write to the
 
16
 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 
17
 * Boston, MA  02111-1307, USA.
 
18
 */
 
19
 
 
20
#include <math.h>
 
21
#include <stdio.h>
 
22
#include "FLAC/assert.h"
 
23
#include "FLAC/format.h"
 
24
#include "private/lpc.h"
 
25
 
 
26
#ifndef M_LN2
 
27
/* math.h in VC++ doesn't seem to have this (how Microsoft is that?) */
 
28
#define M_LN2 0.69314718055994530942
 
29
#endif
 
30
 
 
31
void FLAC__lpc_compute_autocorrelation(const FLAC__real data[], unsigned data_len, unsigned lag, FLAC__real autoc[])
 
32
{
 
33
        /* a readable, but slower, version */
 
34
#if 0
 
35
        FLAC__real d;
 
36
        unsigned i;
 
37
 
 
38
        FLAC__ASSERT(lag > 0);
 
39
        FLAC__ASSERT(lag <= data_len);
 
40
 
 
41
        while(lag--) {
 
42
                for(i = lag, d = 0.0; i < data_len; i++)
 
43
                        d += data[i] * data[i - lag];
 
44
                autoc[lag] = d;
 
45
        }
 
46
#endif
 
47
 
 
48
        /*
 
49
         * this version tends to run faster because of better data locality
 
50
         * ('data_len' is usually much larger than 'lag')
 
51
         */
 
52
        FLAC__real d;
 
53
        unsigned sample, coeff;
 
54
        const unsigned limit = data_len - lag;
 
55
 
 
56
        FLAC__ASSERT(lag > 0);
 
57
        FLAC__ASSERT(lag <= data_len);
 
58
 
 
59
        for(coeff = 0; coeff < lag; coeff++)
 
60
                autoc[coeff] = 0.0;
 
61
        for(sample = 0; sample <= limit; sample++) {
 
62
                d = data[sample];
 
63
                for(coeff = 0; coeff < lag; coeff++)
 
64
                        autoc[coeff] += d * data[sample+coeff];
 
65
        }
 
66
        for(; sample < data_len; sample++) {
 
67
                d = data[sample];
 
68
                for(coeff = 0; coeff < data_len - sample; coeff++)
 
69
                        autoc[coeff] += d * data[sample+coeff];
 
70
        }
 
71
}
 
72
 
 
73
void FLAC__lpc_compute_lp_coefficients(const FLAC__real autoc[], unsigned max_order, FLAC__real lp_coeff[][FLAC__MAX_LPC_ORDER], FLAC__real error[])
 
74
{
 
75
        unsigned i, j;
 
76
        double r, err, ref[FLAC__MAX_LPC_ORDER], lpc[FLAC__MAX_LPC_ORDER];
 
77
 
 
78
        FLAC__ASSERT(0 < max_order);
 
79
        FLAC__ASSERT(max_order <= FLAC__MAX_LPC_ORDER);
 
80
        FLAC__ASSERT(autoc[0] != 0.0);
 
81
 
 
82
        err = autoc[0];
 
83
 
 
84
        for(i = 0; i < max_order; i++) {
 
85
                /* Sum up this iteration's reflection coefficient. */
 
86
                r = -autoc[i+1];
 
87
                for(j = 0; j < i; j++)
 
88
                        r -= lpc[j] * autoc[i-j];
 
89
                ref[i] = (r/=err);
 
90
 
 
91
                /* Update LPC coefficients and total error. */
 
92
                lpc[i]=r;
 
93
                for(j = 0; j < (i>>1); j++) {
 
94
                        double tmp = lpc[j];
 
95
                        lpc[j] += r * lpc[i-1-j];
 
96
                        lpc[i-1-j] += r * tmp;
 
97
                }
 
98
                if(i & 1)
 
99
                        lpc[j] += lpc[j] * r;
 
100
 
 
101
                err *= (1.0 - r * r);
 
102
 
 
103
                /* save this order */
 
104
                for(j = 0; j <= i; j++)
 
105
                        lp_coeff[i][j] = (FLAC__real)(-lpc[j]); /* negate FIR filter coeff to get predictor coeff */
 
106
                error[i] = (FLAC__real)err;
 
107
        }
 
108
}
 
109
 
 
110
int FLAC__lpc_quantize_coefficients(const FLAC__real lp_coeff[], unsigned order, unsigned precision, unsigned bits_per_sample, FLAC__int32 qlp_coeff[], int *shift)
 
111
{
 
112
        unsigned i;
 
113
        double d, cmax = -1e32;
 
114
        FLAC__int32 qmax, qmin;
 
115
        const int max_shiftlimit = (1 << (FLAC__SUBFRAME_LPC_QLP_SHIFT_LEN-1)) - 1;
 
116
        const int min_shiftlimit = -max_shiftlimit - 1;
 
117
 
 
118
        FLAC__ASSERT(bits_per_sample > 0);
 
119
        FLAC__ASSERT(bits_per_sample <= sizeof(FLAC__int32)*8);
 
120
        FLAC__ASSERT(precision > 0);
 
121
        FLAC__ASSERT(precision >= FLAC__MIN_QLP_COEFF_PRECISION);
 
122
        FLAC__ASSERT(precision + bits_per_sample < sizeof(FLAC__int32)*8);
 
123
#ifdef NDEBUG
 
124
        (void)bits_per_sample; /* silence compiler warning about unused parameter */
 
125
#endif
 
126
 
 
127
        /* drop one bit for the sign; from here on out we consider only |lp_coeff[i]| */
 
128
        precision--;
 
129
        qmax = 1 << precision;
 
130
        qmin = -qmax;
 
131
        qmax--;
 
132
 
 
133
        for(i = 0; i < order; i++) {
 
134
                if(lp_coeff[i] == 0.0)
 
135
                        continue;
 
136
                d = fabs(lp_coeff[i]);
 
137
                if(d > cmax)
 
138
                        cmax = d;
 
139
        }
 
140
redo_it:
 
141
        if(cmax <= 0.0) {
 
142
                /* => coefficients are all 0, which means our constant-detect didn't work */
 
143
                return 2;
 
144
        }
 
145
        else {
 
146
                int log2cmax;
 
147
 
 
148
                (void)frexp(cmax, &log2cmax);   
 
149
                log2cmax--;
 
150
                *shift = (int)precision - log2cmax - 1;
 
151
 
 
152
                if(*shift < min_shiftlimit || *shift > max_shiftlimit) {
 
153
                        return 1;
 
154
                }
 
155
        }
 
156
 
 
157
        if(*shift >= 0) {
 
158
                for(i = 0; i < order; i++) {
 
159
                        qlp_coeff[i] = (FLAC__int32)floor((double)lp_coeff[i] * (double)(1 << *shift));
 
160
 
 
161
                        /* double-check the result */
 
162
                        if(qlp_coeff[i] > qmax || qlp_coeff[i] < qmin) {
 
163
#ifdef FLAC__OVERFLOW_DETECT
 
164
                                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, (double)lp_coeff[i] * (double)(1 << *shift), floor((double)lp_coeff[i] * (double)(1 << *shift)));
 
165
#endif
 
166
                                cmax *= 2.0;
 
167
                                goto redo_it;
 
168
                        }
 
169
                }
 
170
        }
 
171
        else { /* (*shift < 0) */
 
172
                const int nshift = -(*shift);
 
173
#ifdef DEBUG
 
174
                fprintf(stderr,"FLAC__lpc_quantize_coefficients: negative shift = %d\n", *shift);
 
175
#endif
 
176
                for(i = 0; i < order; i++) {
 
177
                        qlp_coeff[i] = (FLAC__int32)floor((double)lp_coeff[i] / (double)(1 << nshift));
 
178
 
 
179
                        /* double-check the result */
 
180
                        if(qlp_coeff[i] > qmax || qlp_coeff[i] < qmin) {
 
181
#ifdef FLAC__OVERFLOW_DETECT
 
182
                                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, (double)lp_coeff[i] / (double)(1 << nshift), floor((double)lp_coeff[i] / (double)(1 << nshift)));
 
183
#endif
 
184
                                cmax *= 2.0;
 
185
                                goto redo_it;
 
186
                        }
 
187
                }
 
188
        }
 
189
 
 
190
        return 0;
 
191
}
 
192
 
 
193
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[])
 
194
{
 
195
#ifdef FLAC__OVERFLOW_DETECT
 
196
        FLAC__int64 sumo;
 
197
#endif
 
198
        unsigned i, j;
 
199
        FLAC__int32 sum;
 
200
        const FLAC__int32 *history;
 
201
 
 
202
#ifdef FLAC__OVERFLOW_DETECT_VERBOSE
 
203
        fprintf(stderr,"FLAC__lpc_compute_residual_from_qlp_coefficients: data_len=%d, order=%u, lpq=%d",data_len,order,lp_quantization);
 
204
        for(i=0;i<order;i++)
 
205
                fprintf(stderr,", q[%u]=%d",i,qlp_coeff[i]);
 
206
        fprintf(stderr,"\n");
 
207
#endif
 
208
        FLAC__ASSERT(order > 0);
 
209
 
 
210
        for(i = 0; i < data_len; i++) {
 
211
#ifdef FLAC__OVERFLOW_DETECT
 
212
                sumo = 0;
 
213
#endif
 
214
                sum = 0;
 
215
                history = data;
 
216
                for(j = 0; j < order; j++) {
 
217
                        sum += qlp_coeff[j] * (*(--history));
 
218
#ifdef FLAC__OVERFLOW_DETECT
 
219
                        sumo += (FLAC__int64)qlp_coeff[j] * (FLAC__int64)(*history);
 
220
#if defined _MSC_VER || defined __MINGW32__ /* don't know how to do 64-bit literals in VC++ */
 
221
                        if(sumo < 0) sumo = -sumo;
 
222
                        if(sumo > 2147483647)
 
223
#else
 
224
                        if(sumo > 2147483647ll || sumo < -2147483648ll)
 
225
#endif
 
226
                        {
 
227
                                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);
 
228
                        }
 
229
#endif
 
230
                }
 
231
                *(residual++) = *(data++) - (sum >> lp_quantization);
 
232
        }
 
233
 
 
234
        /* Here's a slower but clearer version:
 
235
        for(i = 0; i < data_len; i++) {
 
236
                sum = 0;
 
237
                for(j = 0; j < order; j++)
 
238
                        sum += qlp_coeff[j] * data[i-j-1];
 
239
                residual[i] = data[i] - (sum >> lp_quantization);
 
240
        }
 
241
        */
 
242
}
 
243
 
 
244
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[])
 
245
{
 
246
#ifdef FLAC__OVERFLOW_DETECT
 
247
        FLAC__int64 sumo;
 
248
#endif
 
249
        unsigned i, j;
 
250
        FLAC__int32 sum;
 
251
        const FLAC__int32 *history;
 
252
 
 
253
#ifdef FLAC__OVERFLOW_DETECT_VERBOSE
 
254
        fprintf(stderr,"FLAC__lpc_restore_signal: data_len=%d, order=%u, lpq=%d",data_len,order,lp_quantization);
 
255
        for(i=0;i<order;i++)
 
256
                fprintf(stderr,", q[%u]=%d",i,qlp_coeff[i]);
 
257
        fprintf(stderr,"\n");
 
258
#endif
 
259
        FLAC__ASSERT(order > 0);
 
260
 
 
261
        for(i = 0; i < data_len; i++) {
 
262
#ifdef FLAC__OVERFLOW_DETECT
 
263
                sumo = 0;
 
264
#endif
 
265
                sum = 0;
 
266
                history = data;
 
267
                for(j = 0; j < order; j++) {
 
268
                        sum += qlp_coeff[j] * (*(--history));
 
269
#ifdef FLAC__OVERFLOW_DETECT
 
270
                        sumo += (FLAC__int64)qlp_coeff[j] * (FLAC__int64)(*history);
 
271
#if defined _MSC_VER || defined __MINGW32__ /* don't know how to do 64-bit literals in VC++ */
 
272
                        if(sumo < 0) sumo = -sumo;
 
273
                        if(sumo > 2147483647)
 
274
#else
 
275
                        if(sumo > 2147483647ll || sumo < -2147483648ll)
 
276
#endif
 
277
                        {
 
278
                                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);
 
279
                        }
 
280
#endif
 
281
                }
 
282
                *(data++) = *(residual++) + (sum >> lp_quantization);
 
283
        }
 
284
 
 
285
        /* Here's a slower but clearer version:
 
286
        for(i = 0; i < data_len; i++) {
 
287
                sum = 0;
 
288
                for(j = 0; j < order; j++)
 
289
                        sum += qlp_coeff[j] * data[i-j-1];
 
290
                data[i] = residual[i] + (sum >> lp_quantization);
 
291
        }
 
292
        */
 
293
}
 
294
 
 
295
FLAC__real FLAC__lpc_compute_expected_bits_per_residual_sample(FLAC__real lpc_error, unsigned total_samples)
 
296
{
 
297
        double error_scale;
 
298
 
 
299
        FLAC__ASSERT(total_samples > 0);
 
300
 
 
301
        error_scale = 0.5 * M_LN2 * M_LN2 / (FLAC__real)total_samples;
 
302
 
 
303
        if(lpc_error > 0.0) {
 
304
                FLAC__real bps = (FLAC__real)((double)0.5 * log(error_scale * lpc_error) / M_LN2);
 
305
                if(bps >= 0.0)
 
306
                        return bps;
 
307
                else
 
308
                        return 0.0;
 
309
        }
 
310
        else if(lpc_error < 0.0) { /* error should not be negative but can happen due to inadequate float resolution */
 
311
                return (FLAC__real)1e32;
 
312
        }
 
313
        else {
 
314
                return 0.0;
 
315
        }
 
316
}
 
317
 
 
318
FLAC__real FLAC__lpc_compute_expected_bits_per_residual_sample_with_error_scale(FLAC__real lpc_error, double error_scale)
 
319
{
 
320
        if(lpc_error > 0.0) {
 
321
                FLAC__real bps = (FLAC__real)((double)0.5 * log(error_scale * lpc_error) / M_LN2);
 
322
                if(bps >= 0.0)
 
323
                        return bps;
 
324
                else
 
325
                        return 0.0;
 
326
        }
 
327
        else if(lpc_error < 0.0) { /* error should not be negative but can happen due to inadequate float resolution */
 
328
                return (FLAC__real)1e32;
 
329
        }
 
330
        else {
 
331
                return 0.0;
 
332
        }
 
333
}
 
334
 
 
335
unsigned FLAC__lpc_compute_best_order(const FLAC__real lpc_error[], unsigned max_order, unsigned total_samples, unsigned bits_per_signal_sample)
 
336
{
 
337
        unsigned order, best_order;
 
338
        FLAC__real best_bits, tmp_bits;
 
339
        double error_scale;
 
340
 
 
341
        FLAC__ASSERT(max_order > 0);
 
342
        FLAC__ASSERT(total_samples > 0);
 
343
 
 
344
        error_scale = 0.5 * M_LN2 * M_LN2 / (FLAC__real)total_samples;
 
345
 
 
346
        best_order = 0;
 
347
        best_bits = FLAC__lpc_compute_expected_bits_per_residual_sample_with_error_scale(lpc_error[0], error_scale) * (FLAC__real)total_samples;
 
348
 
 
349
        for(order = 1; order < max_order; order++) {
 
350
                tmp_bits = FLAC__lpc_compute_expected_bits_per_residual_sample_with_error_scale(lpc_error[order], error_scale) * (FLAC__real)(total_samples - order) + (FLAC__real)(order * bits_per_signal_sample);
 
351
                if(tmp_bits < best_bits) {
 
352
                        best_order = order;
 
353
                        best_bits = tmp_bits;
 
354
                }
 
355
        }
 
356
 
 
357
        return best_order+1; /* +1 since index of lpc_error[] is order-1 */
 
358
}