Libav
|
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 }