3 * Copyright (c) 2006 Justin Ruggles <jruggle@earthlink.net>
5 * This file is part of FFmpeg.
7 * FFmpeg is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU Lesser General Public
9 * License as published by the Free Software Foundation; either
10 * version 2.1 of the License, or (at your option) any later version.
12 * FFmpeg is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Lesser General Public License for more details.
17 * You should have received a copy of the GNU Lesser General Public
18 * License along with FFmpeg; if not, write to the Free Software
19 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22 #include "libavutil/lls.h"
28 * Levinson-Durbin recursion.
29 * Produces LPC coefficients from autocorrelation data.
31 static void compute_lpc_coefs(const double *autoc, int max_order,
32 double lpc[][MAX_LPC_ORDER], double *ref)
36 double lpc_tmp[MAX_LPC_ORDER];
38 for(i=0; i<max_order; i++) lpc_tmp[i] = 0;
41 for(i=0; i<max_order; i++) {
44 r -= lpc_tmp[j] * autoc[i-j];
55 lpc_tmp[j] += r * lpc_tmp[i-1-j];
56 lpc_tmp[i-1-j] += r * tmp;
59 lpc_tmp[j] += lpc_tmp[j] * r;
63 lpc[i][j] = -lpc_tmp[j];
69 * Quantize LPC coefficients
71 static void quantize_lpc_coefs(double *lpc_in, int order, int precision,
72 int32_t *lpc_out, int *shift, int max_shift, int zero_shift)
79 /* define maximum levels */
80 qmax = (1 << (precision - 1)) - 1;
82 /* find maximum coefficient value */
84 for(i=0; i<order; i++) {
85 cmax= FFMAX(cmax, fabs(lpc_in[i]));
88 /* if maximum value quantizes to zero, return all zeros */
89 if(cmax * (1 << max_shift) < 1.0) {
91 memset(lpc_out, 0, sizeof(int32_t) * order);
95 /* calculate level shift which scales max coeff to available bits */
97 while((cmax * (1 << sh) > qmax) && (sh > 0)) {
101 /* since negative shift values are unsupported in decoder, scale down
102 coefficients instead */
103 if(sh == 0 && cmax > qmax) {
104 double scale = ((double)qmax) / cmax;
105 for(i=0; i<order; i++) {
110 /* output quantized coefficients and level shift */
112 for(i=0; i<order; i++) {
113 error += lpc_in[i] * (1 << sh);
114 lpc_out[i] = av_clip(lrintf(error), -qmax, qmax);
120 static int estimate_best_order(double *ref, int max_order)
125 for(i=max_order-1; i>=0; i--) {
135 * Calculate LPC coefficients for multiple orders
137 int ff_lpc_calc_coefs(DSPContext *s,
138 const int32_t *samples, int blocksize, int max_order,
139 int precision, int32_t coefs[][MAX_LPC_ORDER],
140 int *shift, int use_lpc, int omethod, int max_shift, int zero_shift)
142 double autoc[MAX_LPC_ORDER+1];
143 double ref[MAX_LPC_ORDER];
144 double lpc[MAX_LPC_ORDER][MAX_LPC_ORDER];
148 assert(max_order >= MIN_LPC_ORDER && max_order <= MAX_LPC_ORDER);
151 s->flac_compute_autocorr(samples, blocksize, max_order, autoc);
153 compute_lpc_coefs(autoc, max_order, lpc, ref);
156 double var[MAX_LPC_ORDER+1], weight;
158 for(pass=0; pass<use_lpc-1; pass++){
159 av_init_lls(&m[pass&1], max_order);
162 for(i=max_order; i<blocksize; i++){
163 for(j=0; j<=max_order; j++)
164 var[j]= samples[i-j];
167 double eval, inv, rinv;
168 eval= av_evaluate_lls(&m[(pass-1)&1], var+1, max_order-1);
169 eval= (512>>pass) + fabs(eval - var[0]);
172 for(j=0; j<=max_order; j++)
178 av_update_lls(&m[pass&1], var, 1.0);
180 av_solve_lls(&m[pass&1], 0.001, 0);
183 for(i=0; i<max_order; i++){
184 for(j=0; j<max_order; j++)
185 lpc[i][j]= m[(pass-1)&1].coeff[i][j];
186 ref[i]= sqrt(m[(pass-1)&1].variance[i] / weight) * (blocksize - max_order) / 4000;
188 for(i=max_order-1; i>0; i--)
189 ref[i] = ref[i-1] - ref[i];
191 opt_order = max_order;
193 if(omethod == ORDER_METHOD_EST) {
194 opt_order = estimate_best_order(ref, max_order);
196 quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i], max_shift, zero_shift);
198 for(i=0; i<max_order; i++) {
199 quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i], max_shift, zero_shift);