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 25#define LPC_USE_DOUBLE 26#include "lpc.h" 27 28 29/** 30 * Apply Welch window function to audio block 31 */ 32static void apply_welch_window(const int32_t *data, int len, double *w_data) 33{ 34 int i, n2; 35 double w; 36 double c; 37 38 assert(!(len&1)); //the optimization in r11881 does not support odd len 39 //if someone wants odd len extend the change in r11881 40 41 n2 = (len >> 1); 42 c = 2.0 / (len - 1.0); 43 44 w_data+=n2; 45 data+=n2; 46 for(i=0; i<n2; i++) { 47 w = c - n2 + i; 48 w = 1.0 - (w * w); 49 w_data[-i-1] = data[-i-1] * w; 50 w_data[+i ] = data[+i ] * w; 51 } 52} 53 54/** 55 * Calculates autocorrelation data from audio samples 56 * A Welch window function is applied before calculation. 57 */ 58void ff_lpc_compute_autocorr(const int32_t *data, int len, int lag, 59 double *autoc) 60{ 61 int i, j; 62 double tmp[len + lag + 1]; 63 double *data1= tmp + lag; 64 65 apply_welch_window(data, len, data1); 66 67 for(j=0; j<lag; j++) 68 data1[j-lag]= 0.0; 69 data1[len] = 0.0; 70 71 for(j=0; j<lag; j+=2){ 72 double sum0 = 1.0, sum1 = 1.0; 73 for(i=j; i<len; i++){ 74 sum0 += data1[i] * data1[i-j]; 75 sum1 += data1[i] * data1[i-j-1]; 76 } 77 autoc[j ] = sum0; 78 autoc[j+1] = sum1; 79 } 80 81 if(j==lag){ 82 double sum = 1.0; 83 for(i=j-1; i<len; i+=2){ 84 sum += data1[i ] * data1[i-j ] 85 + data1[i+1] * data1[i-j+1]; 86 } 87 autoc[j] = sum; 88 } 89} 90 91/** 92 * Quantize LPC coefficients 93 */ 94static void quantize_lpc_coefs(double *lpc_in, int order, int precision, 95 int32_t *lpc_out, int *shift, int max_shift, int zero_shift) 96{ 97 int i; 98 double cmax, error; 99 int32_t qmax; 100 int sh; 101 102 /* define maximum levels */ 103 qmax = (1 << (precision - 1)) - 1; 104 105 /* find maximum coefficient value */ 106 cmax = 0.0; 107 for(i=0; i<order; i++) { 108 cmax= FFMAX(cmax, fabs(lpc_in[i])); 109 } 110 111 /* if maximum value quantizes to zero, return all zeros */ 112 if(cmax * (1 << max_shift) < 1.0) { 113 *shift = zero_shift; 114 memset(lpc_out, 0, sizeof(int32_t) * order); 115 return; 116 } 117 118 /* calculate level shift which scales max coeff to available bits */ 119 sh = max_shift; 120 while((cmax * (1 << sh) > qmax) && (sh > 0)) { 121 sh--; 122 } 123 124 /* since negative shift values are unsupported in decoder, scale down 125 coefficients instead */ 126 if(sh == 0 && cmax > qmax) { 127 double scale = ((double)qmax) / cmax; 128 for(i=0; i<order; i++) { 129 lpc_in[i] *= scale; 130 } 131 } 132 133 /* output quantized coefficients and level shift */ 134 error=0; 135 for(i=0; i<order; i++) { 136 error -= lpc_in[i] * (1 << sh); 137 lpc_out[i] = av_clip(lrintf(error), -qmax, qmax); 138 error -= lpc_out[i]; 139 } 140 *shift = sh; 141} 142 143static int estimate_best_order(double *ref, int min_order, int max_order) 144{ 145 int i, est; 146 147 est = min_order; 148 for(i=max_order-1; i>=min_order-1; i--) { 149 if(ref[i] > 0.10) { 150 est = i+1; 151 break; 152 } 153 } 154 return est; 155} 156 157/** 158 * Calculate LPC coefficients for multiple orders 159 * 160 * @param use_lpc LPC method for determining coefficients 161 * 0 = LPC with fixed pre-defined coeffs 162 * 1 = LPC with coeffs determined by Levinson-Durbin recursion 163 * 2+ = LPC with coeffs determined by Cholesky factorization using (use_lpc-1) passes. 164 */ 165int ff_lpc_calc_coefs(DSPContext *s, 166 const int32_t *samples, int blocksize, int min_order, 167 int max_order, int precision, 168 int32_t coefs[][MAX_LPC_ORDER], int *shift, int use_lpc, 169 int omethod, int max_shift, int zero_shift) 170{ 171 double autoc[MAX_LPC_ORDER+1]; 172 double ref[MAX_LPC_ORDER]; 173 double lpc[MAX_LPC_ORDER][MAX_LPC_ORDER]; 174 int i, j, pass; 175 int opt_order; 176 177 assert(max_order >= MIN_LPC_ORDER && max_order <= MAX_LPC_ORDER && use_lpc > 0); 178 179 if(use_lpc == 1){ 180 s->lpc_compute_autocorr(samples, blocksize, max_order, autoc); 181 182 compute_lpc_coefs(autoc, max_order, &lpc[0][0], MAX_LPC_ORDER, 0, 1); 183 184 for(i=0; i<max_order; i++) 185 ref[i] = fabs(lpc[i][i]); 186 }else{ 187 LLSModel m[2]; 188 double var[MAX_LPC_ORDER+1], av_uninit(weight); 189 190 for(pass=0; pass<use_lpc-1; pass++){ 191 av_init_lls(&m[pass&1], max_order); 192 193 weight=0; 194 for(i=max_order; i<blocksize; i++){ 195 for(j=0; j<=max_order; j++) 196 var[j]= samples[i-j]; 197 198 if(pass){ 199 double eval, inv, rinv; 200 eval= av_evaluate_lls(&m[(pass-1)&1], var+1, max_order-1); 201 eval= (512>>pass) + fabs(eval - var[0]); 202 inv = 1/eval; 203 rinv = sqrt(inv); 204 for(j=0; j<=max_order; j++) 205 var[j] *= rinv; 206 weight += inv; 207 }else 208 weight++; 209 210 av_update_lls(&m[pass&1], var, 1.0); 211 } 212 av_solve_lls(&m[pass&1], 0.001, 0); 213 } 214 215 for(i=0; i<max_order; i++){ 216 for(j=0; j<max_order; j++) 217 lpc[i][j]=-m[(pass-1)&1].coeff[i][j]; 218 ref[i]= sqrt(m[(pass-1)&1].variance[i] / weight) * (blocksize - max_order) / 4000; 219 } 220 for(i=max_order-1; i>0; i--) 221 ref[i] = ref[i-1] - ref[i]; 222 } 223 opt_order = max_order; 224 225 if(omethod == ORDER_METHOD_EST) { 226 opt_order = estimate_best_order(ref, min_order, max_order); 227 i = opt_order-1; 228 quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i], max_shift, zero_shift); 229 } else { 230 for(i=min_order-1; i<max_order; i++) { 231 quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i], max_shift, zero_shift); 232 } 233 } 234 235 return opt_order; 236} 237