-/**
- * Apply Welch window function to audio block
- */
-static void apply_welch_window(const int32_t *data, int len, double *w_data)
-{
- int i, n2;
- double w;
- double c;
-
- n2 = (len >> 1);
- c = 2.0 / (len - 1.0);
- for(i=0; i<n2; i++) {
- w = c - i - 1.0;
- w = 1.0 - (w * w);
- w_data[i] = data[i] * w;
- w_data[len-1-i] = data[len-1-i] * w;
- }
-}
-
-/**
- * Calculates autocorrelation data from audio samples
- * A Welch window function is applied before calculation.
- */
-void ff_flac_compute_autocorr(const int32_t *data, int len, int lag,
- double *autoc)
-{
- int i, j;
- double tmp[len + lag + 1];
- double *data1= tmp + lag;
-
- apply_welch_window(data, len, data1);
-
- for(j=0; j<lag; j++)
- data1[j-lag]= 0.0;
- data1[len] = 0.0;
-
- for(j=0; j<lag; j+=2){
- double sum0 = 1.0, sum1 = 1.0;
- for(i=0; i<len; i++){
- sum0 += data1[i] * data1[i-j];
- sum1 += data1[i] * data1[i-j-1];
- }
- autoc[j ] = sum0;
- autoc[j+1] = sum1;
- }
-
- if(j==lag){
- double sum = 1.0;
- for(i=0; i<len; i+=2){
- sum += data1[i ] * data1[i-j ]
- + data1[i+1] * data1[i-j+1];
- }
- autoc[j] = sum;
- }
-}
-
-/**
- * Levinson-Durbin recursion.
- * Produces LPC coefficients from autocorrelation data.
- */
-static void compute_lpc_coefs(const double *autoc, int max_order,
- double lpc[][MAX_LPC_ORDER], double *ref)
-{
- int i, j, i2;
- double r, err, tmp;
- double lpc_tmp[MAX_LPC_ORDER];
-
- for(i=0; i<max_order; i++) lpc_tmp[i] = 0;
- err = autoc[0];
-
- for(i=0; i<max_order; i++) {
- r = -autoc[i+1];
- for(j=0; j<i; j++) {
- r -= lpc_tmp[j] * autoc[i-j];
- }
- r /= err;
- ref[i] = fabs(r);
-
- err *= 1.0 - (r * r);
-
- i2 = (i >> 1);
- lpc_tmp[i] = r;
- for(j=0; j<i2; j++) {
- tmp = lpc_tmp[j];
- lpc_tmp[j] += r * lpc_tmp[i-1-j];
- lpc_tmp[i-1-j] += r * tmp;
- }
- if(i & 1) {
- lpc_tmp[j] += lpc_tmp[j] * r;
- }
-
- for(j=0; j<=i; j++) {
- lpc[i][j] = -lpc_tmp[j];
- }
- }
-}
-
-/**
- * Quantize LPC coefficients
- */
-static void quantize_lpc_coefs(double *lpc_in, int order, int precision,
- int32_t *lpc_out, int *shift)
-{
- int i;
- double cmax, error;
- int32_t qmax;
- int sh;
-
- /* define maximum levels */
- qmax = (1 << (precision - 1)) - 1;
-
- /* find maximum coefficient value */
- cmax = 0.0;
- for(i=0; i<order; i++) {
- cmax= FFMAX(cmax, fabs(lpc_in[i]));
- }
-
- /* if maximum value quantizes to zero, return all zeros */
- if(cmax * (1 << MAX_LPC_SHIFT) < 1.0) {
- *shift = 0;
- memset(lpc_out, 0, sizeof(int32_t) * order);
- return;
- }
-
- /* calculate level shift which scales max coeff to available bits */
- sh = MAX_LPC_SHIFT;
- while((cmax * (1 << sh) > qmax) && (sh > 0)) {
- sh--;
- }
-
- /* since negative shift values are unsupported in decoder, scale down
- coefficients instead */
- if(sh == 0 && cmax > qmax) {
- double scale = ((double)qmax) / cmax;
- for(i=0; i<order; i++) {
- lpc_in[i] *= scale;
- }
- }
-
- /* output quantized coefficients and level shift */
- error=0;
- for(i=0; i<order; i++) {
- error += lpc_in[i] * (1 << sh);
- lpc_out[i] = av_clip(lrintf(error), -qmax, qmax);
- error -= lpc_out[i];
- }
- *shift = sh;
-}
-
-static int estimate_best_order(double *ref, int max_order)
-{
- int i, est;
-
- est = 1;
- for(i=max_order-1; i>=0; i--) {
- if(ref[i] > 0.10) {
- est = i+1;
- break;
- }
- }
- return est;
-}
-
-/**
- * Calculate LPC coefficients for multiple orders
- */
-static int lpc_calc_coefs(FlacEncodeContext *s,
- const int32_t *samples, int blocksize, int max_order,
- int precision, int32_t coefs[][MAX_LPC_ORDER],
- int *shift, int use_lpc, int omethod)
-{
- double autoc[MAX_LPC_ORDER+1];
- double ref[MAX_LPC_ORDER];
- double lpc[MAX_LPC_ORDER][MAX_LPC_ORDER];
- int i, j, pass;
- int opt_order;
-
- assert(max_order >= MIN_LPC_ORDER && max_order <= MAX_LPC_ORDER);
-
- if(use_lpc == 1){
- s->dsp.flac_compute_autocorr(samples, blocksize, max_order, autoc);
-
- compute_lpc_coefs(autoc, max_order, lpc, ref);
- }else{
- LLSModel m[2];
- double var[MAX_LPC_ORDER+1], eval, weight;
-
- for(pass=0; pass<use_lpc-1; pass++){
- av_init_lls(&m[pass&1], max_order);
-
- weight=0;
- for(i=max_order; i<blocksize; i++){
- for(j=0; j<=max_order; j++)
- var[j]= samples[i-j];
-
- if(pass){
- eval= av_evaluate_lls(&m[(pass-1)&1], var+1, max_order-1);
- eval= (512>>pass) + fabs(eval - var[0]);
- for(j=0; j<=max_order; j++)
- var[j]/= sqrt(eval);
- weight += 1/eval;
- }else
- weight++;
-
- av_update_lls(&m[pass&1], var, 1.0);
- }
- av_solve_lls(&m[pass&1], 0.001, 0);
- }
-
- for(i=0; i<max_order; i++){
- for(j=0; j<max_order; j++)
- lpc[i][j]= m[(pass-1)&1].coeff[i][j];
- ref[i]= sqrt(m[(pass-1)&1].variance[i] / weight) * (blocksize - max_order) / 4000;
- }
- for(i=max_order-1; i>0; i--)
- ref[i] = ref[i-1] - ref[i];
- }
- opt_order = max_order;
-
- if(omethod == ORDER_METHOD_EST) {
- opt_order = estimate_best_order(ref, max_order);
- i = opt_order-1;
- quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i]);
- } else {
- for(i=0; i<max_order; i++) {
- quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i]);
- }
- }
-
- return opt_order;
-}
-
-
-static void encode_residual_verbatim(int32_t *res, int32_t *smp, int n)
-{
- assert(n > 0);
- memcpy(res, smp, n * sizeof(int32_t));
-}