-/**
- * Levinson-Durbin recursion.
- * Produces LPC coefficients from autocorrelation data.
- */
-static void compute_lpc_coefs(const double *autoc, int max_order,
- double lpc[][MAX_LPC_ORDER], double *ref)
-{
- int i, j, i2;
- double r, err, tmp;
- double lpc_tmp[MAX_LPC_ORDER];
-
- for(i=0; i<max_order; i++) lpc_tmp[i] = 0;
- err = autoc[0];
-
- for(i=0; i<max_order; i++) {
- r = -autoc[i+1];
- for(j=0; j<i; j++) {
- r -= lpc_tmp[j] * autoc[i-j];
- }
- r /= err;
- ref[i] = fabs(r);
-
- err *= 1.0 - (r * r);
-
- i2 = (i >> 1);
- lpc_tmp[i] = r;
- for(j=0; j<i2; j++) {
- tmp = lpc_tmp[j];
- lpc_tmp[j] += r * lpc_tmp[i-1-j];
- lpc_tmp[i-1-j] += r * tmp;
- }
- if(i & 1) {
- lpc_tmp[j] += lpc_tmp[j] * r;
- }
-
- for(j=0; j<=i; j++) {
- lpc[i][j] = -lpc_tmp[j];
- }
- }
-}
-
-/**
- * Quantize LPC coefficients
- */
-static void quantize_lpc_coefs(double *lpc_in, int order, int precision,
- int32_t *lpc_out, int *shift)
-{
- int i;
- double cmax, error;
- int32_t qmax;
- int sh;
-
- /* define maximum levels */
- qmax = (1 << (precision - 1)) - 1;
-
- /* find maximum coefficient value */
- cmax = 0.0;
- for(i=0; i<order; i++) {
- cmax= FFMAX(cmax, fabs(lpc_in[i]));
- }
-
- /* if maximum value quantizes to zero, return all zeros */
- if(cmax * (1 << MAX_LPC_SHIFT) < 1.0) {
- *shift = 0;
- memset(lpc_out, 0, sizeof(int32_t) * order);
- return;
- }
-
- /* calculate level shift which scales max coeff to available bits */
- sh = MAX_LPC_SHIFT;
- while((cmax * (1 << sh) > qmax) && (sh > 0)) {
- sh--;
- }
-
- /* since negative shift values are unsupported in decoder, scale down
- coefficients instead */
- if(sh == 0 && cmax > qmax) {
- double scale = ((double)qmax) / cmax;
- for(i=0; i<order; i++) {
- lpc_in[i] *= scale;
- }
- }
-
- /* output quantized coefficients and level shift */
- error=0;
- for(i=0; i<order; i++) {
- error += lpc_in[i] * (1 << sh);
- lpc_out[i] = av_clip(lrintf(error), -qmax, qmax);
- error -= lpc_out[i];
- }
- *shift = sh;
-}
-
-static int estimate_best_order(double *ref, int max_order)
-{
- int i, est;
-
- est = 1;
- for(i=max_order-1; i>=0; i--) {
- if(ref[i] > 0.10) {
- est = i+1;
- break;
- }
- }
- return est;
-}
-
-/**
- * Calculate LPC coefficients for multiple orders
- */
-static int lpc_calc_coefs(FlacEncodeContext *s,
- const int32_t *samples, int blocksize, int max_order,
- int precision, int32_t coefs[][MAX_LPC_ORDER],
- int *shift, int use_lpc, int omethod)
-{
- double autoc[MAX_LPC_ORDER+1];
- double ref[MAX_LPC_ORDER];
- double lpc[MAX_LPC_ORDER][MAX_LPC_ORDER];
- int i, j, pass;
- int opt_order;
-
- assert(max_order >= MIN_LPC_ORDER && max_order <= MAX_LPC_ORDER);
-
- if(use_lpc == 1){
- s->dsp.flac_compute_autocorr(samples, blocksize, max_order, autoc);
-
- compute_lpc_coefs(autoc, max_order, lpc, ref);
- }else{
- LLSModel m[2];
- double var[MAX_LPC_ORDER+1], weight;
-
- for(pass=0; pass<use_lpc-1; pass++){
- av_init_lls(&m[pass&1], max_order);
-
- weight=0;
- for(i=max_order; i<blocksize; i++){
- for(j=0; j<=max_order; j++)
- var[j]= samples[i-j];
-
- if(pass){
- double eval, inv, rinv;
- eval= av_evaluate_lls(&m[(pass-1)&1], var+1, max_order-1);
- eval= (512>>pass) + fabs(eval - var[0]);
- inv = 1/eval;
- rinv = sqrt(inv);
- for(j=0; j<=max_order; j++)
- var[j] *= rinv;
- weight += inv;
- }else
- weight++;
-
- av_update_lls(&m[pass&1], var, 1.0);
- }
- av_solve_lls(&m[pass&1], 0.001, 0);
- }
-
- for(i=0; i<max_order; i++){
- for(j=0; j<max_order; j++)
- lpc[i][j]= m[(pass-1)&1].coeff[i][j];
- ref[i]= sqrt(m[(pass-1)&1].variance[i] / weight) * (blocksize - max_order) / 4000;
- }
- for(i=max_order-1; i>0; i--)
- ref[i] = ref[i-1] - ref[i];
- }
- opt_order = max_order;
-
- if(omethod == ORDER_METHOD_EST) {
- opt_order = estimate_best_order(ref, max_order);
- i = opt_order-1;
- quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i]);
- } else {
- for(i=0; i<max_order; i++) {
- quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i]);
- }
- }
-
- return opt_order;
-}
-