• Main Page
  • Related Pages
  • Modules
  • Data Structures
  • Files
  • File List
  • Globals

libavcodec/lpc.c

Go to the documentation of this file.
00001 
00022 #include "libavutil/lls.h"
00023 #include "dsputil.h"
00024 
00025 #define LPC_USE_DOUBLE
00026 #include "lpc.h"
00027 
00028 
00032 static void apply_welch_window(const int32_t *data, int len, double *w_data)
00033 {
00034     int i, n2;
00035     double w;
00036     double c;
00037 
00038     assert(!(len&1)); //the optimization in r11881 does not support odd len
00039                       //if someone wants odd len extend the change in r11881
00040 
00041     n2 = (len >> 1);
00042     c = 2.0 / (len - 1.0);
00043 
00044     w_data+=n2;
00045       data+=n2;
00046     for(i=0; i<n2; i++) {
00047         w = c - n2 + i;
00048         w = 1.0 - (w * w);
00049         w_data[-i-1] = data[-i-1] * w;
00050         w_data[+i  ] = data[+i  ] * w;
00051     }
00052 }
00053 
00058 void ff_lpc_compute_autocorr(const int32_t *data, int len, int lag,
00059                              double *autoc)
00060 {
00061     int i, j;
00062     double tmp[len + lag + 1];
00063     double *data1= tmp + lag;
00064 
00065     apply_welch_window(data, len, data1);
00066 
00067     for(j=0; j<lag; j++)
00068         data1[j-lag]= 0.0;
00069     data1[len] = 0.0;
00070 
00071     for(j=0; j<lag; j+=2){
00072         double sum0 = 1.0, sum1 = 1.0;
00073         for(i=j; i<len; i++){
00074             sum0 += data1[i] * data1[i-j];
00075             sum1 += data1[i] * data1[i-j-1];
00076         }
00077         autoc[j  ] = sum0;
00078         autoc[j+1] = sum1;
00079     }
00080 
00081     if(j==lag){
00082         double sum = 1.0;
00083         for(i=j-1; i<len; i+=2){
00084             sum += data1[i  ] * data1[i-j  ]
00085                  + data1[i+1] * data1[i-j+1];
00086         }
00087         autoc[j] = sum;
00088     }
00089 }
00090 
00094 static void quantize_lpc_coefs(double *lpc_in, int order, int precision,
00095                                int32_t *lpc_out, int *shift, int max_shift, int zero_shift)
00096 {
00097     int i;
00098     double cmax, error;
00099     int32_t qmax;
00100     int sh;
00101 
00102     /* define maximum levels */
00103     qmax = (1 << (precision - 1)) - 1;
00104 
00105     /* find maximum coefficient value */
00106     cmax = 0.0;
00107     for(i=0; i<order; i++) {
00108         cmax= FFMAX(cmax, fabs(lpc_in[i]));
00109     }
00110 
00111     /* if maximum value quantizes to zero, return all zeros */
00112     if(cmax * (1 << max_shift) < 1.0) {
00113         *shift = zero_shift;
00114         memset(lpc_out, 0, sizeof(int32_t) * order);
00115         return;
00116     }
00117 
00118     /* calculate level shift which scales max coeff to available bits */
00119     sh = max_shift;
00120     while((cmax * (1 << sh) > qmax) && (sh > 0)) {
00121         sh--;
00122     }
00123 
00124     /* since negative shift values are unsupported in decoder, scale down
00125        coefficients instead */
00126     if(sh == 0 && cmax > qmax) {
00127         double scale = ((double)qmax) / cmax;
00128         for(i=0; i<order; i++) {
00129             lpc_in[i] *= scale;
00130         }
00131     }
00132 
00133     /* output quantized coefficients and level shift */
00134     error=0;
00135     for(i=0; i<order; i++) {
00136         error -= lpc_in[i] * (1 << sh);
00137         lpc_out[i] = av_clip(lrintf(error), -qmax, qmax);
00138         error -= lpc_out[i];
00139     }
00140     *shift = sh;
00141 }
00142 
00143 static int estimate_best_order(double *ref, int min_order, int max_order)
00144 {
00145     int i, est;
00146 
00147     est = min_order;
00148     for(i=max_order-1; i>=min_order-1; i--) {
00149         if(ref[i] > 0.10) {
00150             est = i+1;
00151             break;
00152         }
00153     }
00154     return est;
00155 }
00156 
00165 int ff_lpc_calc_coefs(DSPContext *s,
00166                       const int32_t *samples, int blocksize, int min_order,
00167                       int max_order, int precision,
00168                       int32_t coefs[][MAX_LPC_ORDER], int *shift, int use_lpc,
00169                       int omethod, int max_shift, int zero_shift)
00170 {
00171     double autoc[MAX_LPC_ORDER+1];
00172     double ref[MAX_LPC_ORDER];
00173     double lpc[MAX_LPC_ORDER][MAX_LPC_ORDER];
00174     int i, j, pass;
00175     int opt_order;
00176 
00177     assert(max_order >= MIN_LPC_ORDER && max_order <= MAX_LPC_ORDER && use_lpc > 0);
00178 
00179     if(use_lpc == 1){
00180         s->lpc_compute_autocorr(samples, blocksize, max_order, autoc);
00181 
00182         compute_lpc_coefs(autoc, max_order, &lpc[0][0], MAX_LPC_ORDER, 0, 1);
00183 
00184         for(i=0; i<max_order; i++)
00185             ref[i] = fabs(lpc[i][i]);
00186     }else{
00187         LLSModel m[2];
00188         double var[MAX_LPC_ORDER+1], av_uninit(weight);
00189 
00190         for(pass=0; pass<use_lpc-1; pass++){
00191             av_init_lls(&m[pass&1], max_order);
00192 
00193             weight=0;
00194             for(i=max_order; i<blocksize; i++){
00195                 for(j=0; j<=max_order; j++)
00196                     var[j]= samples[i-j];
00197 
00198                 if(pass){
00199                     double eval, inv, rinv;
00200                     eval= av_evaluate_lls(&m[(pass-1)&1], var+1, max_order-1);
00201                     eval= (512>>pass) + fabs(eval - var[0]);
00202                     inv = 1/eval;
00203                     rinv = sqrt(inv);
00204                     for(j=0; j<=max_order; j++)
00205                         var[j] *= rinv;
00206                     weight += inv;
00207                 }else
00208                     weight++;
00209 
00210                 av_update_lls(&m[pass&1], var, 1.0);
00211             }
00212             av_solve_lls(&m[pass&1], 0.001, 0);
00213         }
00214 
00215         for(i=0; i<max_order; i++){
00216             for(j=0; j<max_order; j++)
00217                 lpc[i][j]=-m[(pass-1)&1].coeff[i][j];
00218             ref[i]= sqrt(m[(pass-1)&1].variance[i] / weight) * (blocksize - max_order) / 4000;
00219         }
00220         for(i=max_order-1; i>0; i--)
00221             ref[i] = ref[i-1] - ref[i];
00222     }
00223     opt_order = max_order;
00224 
00225     if(omethod == ORDER_METHOD_EST) {
00226         opt_order = estimate_best_order(ref, min_order, max_order);
00227         i = opt_order-1;
00228         quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i], max_shift, zero_shift);
00229     } else {
00230         for(i=min_order-1; i<max_order; i++) {
00231             quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i], max_shift, zero_shift);
00232         }
00233     }
00234 
00235     return opt_order;
00236 }

Generated on Fri Sep 16 2011 17:17:39 for FFmpeg by  doxygen 1.7.1