Libav 0.7.1
|
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 00157 int ff_lpc_calc_coefs(LPCContext *s, 00158 const int32_t *samples, int blocksize, int min_order, 00159 int max_order, int precision, 00160 int32_t coefs[][MAX_LPC_ORDER], int *shift, 00161 enum FFLPCType lpc_type, int lpc_passes, 00162 int omethod, int max_shift, int zero_shift) 00163 { 00164 double autoc[MAX_LPC_ORDER+1]; 00165 double ref[MAX_LPC_ORDER]; 00166 double lpc[MAX_LPC_ORDER][MAX_LPC_ORDER]; 00167 int i, j, pass; 00168 int opt_order; 00169 00170 assert(max_order >= MIN_LPC_ORDER && max_order <= MAX_LPC_ORDER && 00171 lpc_type > FF_LPC_TYPE_FIXED); 00172 00173 /* reinit LPC context if parameters have changed */ 00174 if (blocksize != s->blocksize || max_order != s->max_order || 00175 lpc_type != s->lpc_type) { 00176 ff_lpc_end(s); 00177 ff_lpc_init(s, blocksize, max_order, lpc_type); 00178 } 00179 00180 if (lpc_type == FF_LPC_TYPE_LEVINSON) { 00181 double *windowed_samples = s->windowed_samples + max_order; 00182 00183 s->lpc_apply_welch_window(samples, blocksize, windowed_samples); 00184 00185 s->lpc_compute_autocorr(windowed_samples, blocksize, max_order, autoc); 00186 00187 compute_lpc_coefs(autoc, max_order, &lpc[0][0], MAX_LPC_ORDER, 0, 1); 00188 00189 for(i=0; i<max_order; i++) 00190 ref[i] = fabs(lpc[i][i]); 00191 } else if (lpc_type == FF_LPC_TYPE_CHOLESKY) { 00192 LLSModel m[2]; 00193 double var[MAX_LPC_ORDER+1], av_uninit(weight); 00194 00195 for(pass=0; pass<lpc_passes; pass++){ 00196 av_init_lls(&m[pass&1], max_order); 00197 00198 weight=0; 00199 for(i=max_order; i<blocksize; i++){ 00200 for(j=0; j<=max_order; j++) 00201 var[j]= samples[i-j]; 00202 00203 if(pass){ 00204 double eval, inv, rinv; 00205 eval= av_evaluate_lls(&m[(pass-1)&1], var+1, max_order-1); 00206 eval= (512>>pass) + fabs(eval - var[0]); 00207 inv = 1/eval; 00208 rinv = sqrt(inv); 00209 for(j=0; j<=max_order; j++) 00210 var[j] *= rinv; 00211 weight += inv; 00212 }else 00213 weight++; 00214 00215 av_update_lls(&m[pass&1], var, 1.0); 00216 } 00217 av_solve_lls(&m[pass&1], 0.001, 0); 00218 } 00219 00220 for(i=0; i<max_order; i++){ 00221 for(j=0; j<max_order; j++) 00222 lpc[i][j]=-m[(pass-1)&1].coeff[i][j]; 00223 ref[i]= sqrt(m[(pass-1)&1].variance[i] / weight) * (blocksize - max_order) / 4000; 00224 } 00225 for(i=max_order-1; i>0; i--) 00226 ref[i] = ref[i-1] - ref[i]; 00227 } 00228 opt_order = max_order; 00229 00230 if(omethod == ORDER_METHOD_EST) { 00231 opt_order = estimate_best_order(ref, min_order, max_order); 00232 i = opt_order-1; 00233 quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i], max_shift, zero_shift); 00234 } else { 00235 for(i=min_order-1; i<max_order; i++) { 00236 quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i], max_shift, zero_shift); 00237 } 00238 } 00239 00240 return opt_order; 00241 } 00242 00243 av_cold int ff_lpc_init(LPCContext *s, int blocksize, int max_order, 00244 enum FFLPCType lpc_type) 00245 { 00246 s->blocksize = blocksize; 00247 s->max_order = max_order; 00248 s->lpc_type = lpc_type; 00249 00250 if (lpc_type == FF_LPC_TYPE_LEVINSON) { 00251 s->windowed_samples = av_mallocz((blocksize + max_order + 2) * 00252 sizeof(*s->windowed_samples)); 00253 if (!s->windowed_samples) 00254 return AVERROR(ENOMEM); 00255 } else { 00256 s->windowed_samples = NULL; 00257 } 00258 00259 s->lpc_apply_welch_window = lpc_apply_welch_window_c; 00260 s->lpc_compute_autocorr = lpc_compute_autocorr_c; 00261 00262 if (HAVE_MMX) 00263 ff_lpc_init_x86(s); 00264 00265 return 0; 00266 } 00267 00268 av_cold void ff_lpc_end(LPCContext *s) 00269 { 00270 av_freep(&s->windowed_samples); 00271 }