git.sesse.net Git - ffmpeg/blob - libavcodec/lpc.c

   1 /**
   2  * LPC utility code
   3  * Copyright (c) 2006  Justin Ruggles <justin.ruggles@gmail.com>
   4  *
   5  * This file is part of FFmpeg.
   6  *
   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.
  11  *
  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.
  16  *
  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
  20  */
  21
  22 #include "libavutil/lls.h"
  23 #include "dsputil.h"
  24 #include "lpc.h"
  25
  26
  27 /**
  28  * Levinson-Durbin recursion.
  29  * Produces LPC coefficients from autocorrelation data.
  30  */
  31 static void compute_lpc_coefs(const double *autoc, int max_order,
  32                               double lpc[][MAX_LPC_ORDER], double *ref)
  33 {
  34     int i, j;
  35     double err = autoc[0];
  36     double lpc_tmp[MAX_LPC_ORDER];
  37
  38     for(i=0; i<max_order; i++) {
  39         double r = -autoc[i+1];
  40
  41         for(j=0; j<i; j++)
  42             r -= lpc_tmp[j] * autoc[i-j];
  43
  44         r /= err;
  45         ref[i] = fabs(r);
  46
  47         err *= 1.0 - (r * r);
  48
  49         lpc_tmp[i] = r;
  50         for(j=0; j < i>>1; j++) {
  51             double tmp = lpc_tmp[j];
  52             lpc_tmp[j] += r * lpc_tmp[i-1-j];
  53             lpc_tmp[i-1-j] += r * tmp;
  54         }
  55
  56         if(i & 1)
  57             lpc_tmp[j] += lpc_tmp[j] * r;
  58
  59         for(j=0; j<=i; j++)
  60             lpc[i][j] = -lpc_tmp[j];
  61     }
  62 }
  63
  64 /**
  65  * Quantize LPC coefficients
  66  */
  67 static void quantize_lpc_coefs(double *lpc_in, int order, int precision,
  68                                int32_t *lpc_out, int *shift, int max_shift, int zero_shift)
  69 {
  70     int i;
  71     double cmax, error;
  72     int32_t qmax;
  73     int sh;
  74
  75     /* define maximum levels */
  76     qmax = (1 << (precision - 1)) - 1;
  77
  78     /* find maximum coefficient value */
  79     cmax = 0.0;
  80     for(i=0; i<order; i++) {
  81         cmax= FFMAX(cmax, fabs(lpc_in[i]));
  82     }
  83
  84     /* if maximum value quantizes to zero, return all zeros */
  85     if(cmax * (1 << max_shift) < 1.0) {
  86         *shift = zero_shift;
  87         memset(lpc_out, 0, sizeof(int32_t) * order);
  88         return;
  89     }
  90
  91     /* calculate level shift which scales max coeff to available bits */
  92     sh = max_shift;
  93     while((cmax * (1 << sh) > qmax) && (sh > 0)) {
  94         sh--;
  95     }
  96
  97     /* since negative shift values are unsupported in decoder, scale down
  98        coefficients instead */
  99     if(sh == 0 && cmax > qmax) {
 100         double scale = ((double)qmax) / cmax;
 101         for(i=0; i<order; i++) {
 102             lpc_in[i] *= scale;
 103         }
 104     }
 105
 106     /* output quantized coefficients and level shift */
 107     error=0;
 108     for(i=0; i<order; i++) {
 109         error += lpc_in[i] * (1 << sh);
 110         lpc_out[i] = av_clip(lrintf(error), -qmax, qmax);
 111         error -= lpc_out[i];
 112     }
 113     *shift = sh;
 114 }
 115
 116 static int estimate_best_order(double *ref, int min_order, int max_order)
 117 {
 118     int i, est;
 119
 120     est = min_order;
 121     for(i=max_order-1; i>=min_order-1; i--) {
 122         if(ref[i] > 0.10) {
 123             est = i+1;
 124             break;
 125         }
 126     }
 127     return est;
 128 }
 129
 130 /**
 131  * Calculate LPC coefficients for multiple orders
 132  */
 133 int ff_lpc_calc_coefs(DSPContext *s,
 134                       const int32_t *samples, int blocksize, int min_order,
 135                       int max_order, int precision,
 136                       int32_t coefs[][MAX_LPC_ORDER], int *shift, int use_lpc,
 137                       int omethod, int max_shift, int zero_shift)
 138 {
 139     double autoc[MAX_LPC_ORDER+1];
 140     double ref[MAX_LPC_ORDER];
 141     double lpc[MAX_LPC_ORDER][MAX_LPC_ORDER];
 142     int i, j, pass;
 143     int opt_order;
 144
 145     assert(max_order >= MIN_LPC_ORDER && max_order <= MAX_LPC_ORDER);
 146
 147     if(use_lpc == 1){
 148         s->flac_compute_autocorr(samples, blocksize, max_order, autoc);
 149
 150         compute_lpc_coefs(autoc, max_order, lpc, ref);
 151     }else{
 152         LLSModel m[2];
 153         double var[MAX_LPC_ORDER+1], weight;
 154
 155         for(pass=0; pass<use_lpc-1; pass++){
 156             av_init_lls(&m[pass&1], max_order);
 157
 158             weight=0;
 159             for(i=max_order; i<blocksize; i++){
 160                 for(j=0; j<=max_order; j++)
 161                     var[j]= samples[i-j];
 162
 163                 if(pass){
 164                     double eval, inv, rinv;
 165                     eval= av_evaluate_lls(&m[(pass-1)&1], var+1, max_order-1);
 166                     eval= (512>>pass) + fabs(eval - var[0]);
 167                     inv = 1/eval;
 168                     rinv = sqrt(inv);
 169                     for(j=0; j<=max_order; j++)
 170                         var[j] *= rinv;
 171                     weight += inv;
 172                 }else
 173                     weight++;
 174
 175                 av_update_lls(&m[pass&1], var, 1.0);
 176             }
 177             av_solve_lls(&m[pass&1], 0.001, 0);
 178         }
 179
 180         for(i=0; i<max_order; i++){
 181             for(j=0; j<max_order; j++)
 182                 lpc[i][j]= m[(pass-1)&1].coeff[i][j];
 183             ref[i]= sqrt(m[(pass-1)&1].variance[i] / weight) * (blocksize - max_order) / 4000;
 184         }
 185         for(i=max_order-1; i>0; i--)
 186             ref[i] = ref[i-1] - ref[i];
 187     }
 188     opt_order = max_order;
 189
 190     if(omethod == ORDER_METHOD_EST) {
 191         opt_order = estimate_best_order(ref, min_order, max_order);
 192         i = opt_order-1;
 193         quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i], max_shift, zero_shift);
 194     } else {
 195         for(i=min_order-1; i<max_order; i++) {
 196             quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i], max_shift, zero_shift);
 197         }
 198     }
 199
 200     return opt_order;
 201 }