00001 00022 #include "libavutil/lls.h" 00023 00024 #define LPC_USE_DOUBLE 00025 #include "lpc.h" 00026 00027 00031 static void lpc_apply_welch_window_c(const int32_t *data, int len, 00032 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 static void lpc_compute_autocorr_c(const double *data, int len, int lag, 00059 double *autoc) 00060 { 00061 int i, j; 00062 00063 for(j=0; j<lag; j+=2){ 00064 double sum0 = 1.0, sum1 = 1.0; 00065 for(i=j; i<len; i++){ 00066 sum0 += data[i] * data[i-j]; 00067 sum1 += data[i] * data[i-j-1]; 00068 } 00069 autoc[j ] = sum0; 00070 autoc[j+1] = sum1; 00071 } 00072 00073 if(j==lag){ 00074 double sum = 1.0; 00075 for(i=j-1; i<len; i+=2){ 00076 sum += data[i ] * data[i-j ] 00077 + data[i+1] * data[i-j+1]; 00078 } 00079 autoc[j] = sum; 00080 } 00081 } 00082 00086 static void quantize_lpc_coefs(double *lpc_in, int order, int precision, 00087 int32_t *lpc_out, int *shift, int max_shift, int zero_shift) 00088 { 00089 int i; 00090 double cmax, error; 00091 int32_t qmax; 00092 int sh; 00093 00094 /* define maximum levels */ 00095 qmax = (1 << (precision - 1)) - 1; 00096 00097 /* find maximum coefficient value */ 00098 cmax = 0.0; 00099 for(i=0; i<order; i++) { 00100 cmax= FFMAX(cmax, fabs(lpc_in[i])); 00101 } 00102 00103 /* if maximum value quantizes to zero, return all zeros */ 00104 if(cmax * (1 << max_shift) < 1.0) { 00105 *shift = zero_shift; 00106 memset(lpc_out, 0, sizeof(int32_t) * order); 00107 return; 00108 } 00109 00110 /* calculate level shift which scales max coeff to available bits */ 00111 sh = max_shift; 00112 while((cmax * (1 << sh) > qmax) && (sh > 0)) { 00113 sh--; 00114 } 00115 00116 /* since negative shift values are unsupported in decoder, scale down 00117 coefficients instead */ 00118 if(sh == 0 && cmax > qmax) { 00119 double scale = ((double)qmax) / cmax; 00120 for(i=0; i<order; i++) { 00121 lpc_in[i] *= scale; 00122 } 00123 } 00124 00125 /* output quantized coefficients and level shift */ 00126 error=0; 00127 for(i=0; i<order; i++) { 00128 error -= lpc_in[i] * (1 << sh); 00129 lpc_out[i] = av_clip(lrintf(error), -qmax, qmax); 00130 error -= lpc_out[i]; 00131 } 00132 *shift = sh; 00133 } 00134 00135 static int estimate_best_order(double *ref, int min_order, int max_order) 00136 { 00137 int i, est; 00138 00139 est = min_order; 00140 for(i=max_order-1; i>=min_order-1; i--) { 00141 if(ref[i] > 0.10) { 00142 est = i+1; 00143 break; 00144 } 00145 } 00146 return est; 00147 } 00148 00155 int ff_lpc_calc_coefs(LPCContext *s, 00156 const int32_t *samples, int blocksize, int min_order, 00157 int max_order, int precision, 00158 int32_t coefs[][MAX_LPC_ORDER], int *shift, 00159 enum FFLPCType lpc_type, int lpc_passes, 00160 int omethod, int max_shift, int zero_shift) 00161 { 00162 double autoc[MAX_LPC_ORDER+1]; 00163 double ref[MAX_LPC_ORDER]; 00164 double lpc[MAX_LPC_ORDER][MAX_LPC_ORDER]; 00165 int i, j, pass; 00166 int opt_order; 00167 00168 assert(max_order >= MIN_LPC_ORDER && max_order <= MAX_LPC_ORDER && 00169 lpc_type > FF_LPC_TYPE_FIXED); 00170 00171 /* reinit LPC context if parameters have changed */ 00172 if (blocksize != s->blocksize || max_order != s->max_order || 00173 lpc_type != s->lpc_type) { 00174 ff_lpc_end(s); 00175 ff_lpc_init(s, blocksize, max_order, lpc_type); 00176 } 00177 00178 if (lpc_type == FF_LPC_TYPE_LEVINSON) { 00179 double *windowed_samples = s->windowed_samples + max_order; 00180 00181 s->lpc_apply_welch_window(samples, blocksize, windowed_samples); 00182 00183 s->lpc_compute_autocorr(windowed_samples, blocksize, max_order, autoc); 00184 00185 compute_lpc_coefs(autoc, max_order, &lpc[0][0], MAX_LPC_ORDER, 0, 1); 00186 00187 for(i=0; i<max_order; i++) 00188 ref[i] = fabs(lpc[i][i]); 00189 } else if (lpc_type == FF_LPC_TYPE_CHOLESKY) { 00190 LLSModel m[2]; 00191 double var[MAX_LPC_ORDER+1], av_uninit(weight); 00192 00193 for(pass=0; pass<lpc_passes; pass++){ 00194 av_init_lls(&m[pass&1], max_order); 00195 00196 weight=0; 00197 for(i=max_order; i<blocksize; i++){ 00198 for(j=0; j<=max_order; j++) 00199 var[j]= samples[i-j]; 00200 00201 if(pass){ 00202 double eval, inv, rinv; 00203 eval= av_evaluate_lls(&m[(pass-1)&1], var+1, max_order-1); 00204 eval= (512>>pass) + fabs(eval - var[0]); 00205 inv = 1/eval; 00206 rinv = sqrt(inv); 00207 for(j=0; j<=max_order; j++) 00208 var[j] *= rinv; 00209 weight += inv; 00210 }else 00211 weight++; 00212 00213 av_update_lls(&m[pass&1], var, 1.0); 00214 } 00215 av_solve_lls(&m[pass&1], 0.001, 0); 00216 } 00217 00218 for(i=0; i<max_order; i++){ 00219 for(j=0; j<max_order; j++) 00220 lpc[i][j]=-m[(pass-1)&1].coeff[i][j]; 00221 ref[i]= sqrt(m[(pass-1)&1].variance[i] / weight) * (blocksize - max_order) / 4000; 00222 } 00223 for(i=max_order-1; i>0; i--) 00224 ref[i] = ref[i-1] - ref[i]; 00225 } 00226 opt_order = max_order; 00227 00228 if(omethod == ORDER_METHOD_EST) { 00229 opt_order = estimate_best_order(ref, min_order, max_order); 00230 i = opt_order-1; 00231 quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i], max_shift, zero_shift); 00232 } else { 00233 for(i=min_order-1; i<max_order; i++) { 00234 quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i], max_shift, zero_shift); 00235 } 00236 } 00237 00238 return opt_order; 00239 } 00240 00241 av_cold int ff_lpc_init(LPCContext *s, int blocksize, int max_order, 00242 enum FFLPCType lpc_type) 00243 { 00244 s->blocksize = blocksize; 00245 s->max_order = max_order; 00246 s->lpc_type = lpc_type; 00247 00248 if (lpc_type == FF_LPC_TYPE_LEVINSON) { 00249 s->windowed_samples = av_mallocz((blocksize + max_order + 2) * 00250 sizeof(*s->windowed_samples)); 00251 if (!s->windowed_samples) 00252 return AVERROR(ENOMEM); 00253 } else { 00254 s->windowed_samples = NULL; 00255 } 00256 00257 s->lpc_apply_welch_window = lpc_apply_welch_window_c; 00258 s->lpc_compute_autocorr = lpc_compute_autocorr_c; 00259 00260 if (HAVE_MMX) 00261 ff_lpc_init_x86(s); 00262 00263 return 0; 00264 } 00265 00266 av_cold void ff_lpc_end(LPCContext *s) 00267 { 00268 av_freep(&s->windowed_samples); 00269 }