/**
* LPC utility code
- * Copyright (c) 2006 Justin Ruggles <jruggle@earthlink.net>
+ * Copyright (c) 2006 Justin Ruggles <justin.ruggles@gmail.com>
*
* This file is part of FFmpeg.
*
#include "libavutil/lls.h"
#include "dsputil.h"
+
+#define LPC_USE_DOUBLE
#include "lpc.h"
/**
- * Levinson-Durbin recursion.
- * Produces LPC coefficients from autocorrelation data.
+ * Apply Welch window function to audio block
*/
-static void compute_lpc_coefs(const double *autoc, int max_order,
- double lpc[][MAX_LPC_ORDER], double *ref)
+static void apply_welch_window(const int32_t *data, int len, double *w_data)
{
- 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];
- }
- }
+ int i, n2;
+ double w;
+ double c;
+
+ assert(!(len&1)); //the optimization in r11881 does not support odd len
+ //if someone wants odd len extend the change in r11881
+
+ n2 = (len >> 1);
+ c = 2.0 / (len - 1.0);
+
+ w_data+=n2;
+ data+=n2;
+ for(i=0; i<n2; i++) {
+ w = c - n2 + i;
+ w = 1.0 - (w * w);
+ w_data[-i-1] = data[-i-1] * w;
+ w_data[+i ] = data[+i ] * w;
+ }
+}
+
+/**
+ * Calculates autocorrelation data from audio samples
+ * A Welch window function is applied before calculation.
+ */
+void ff_lpc_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=j; 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=j-1; i<len; i+=2){
+ sum += data1[i ] * data1[i-j ]
+ + data1[i+1] * data1[i-j+1];
+ }
+ autoc[j] = sum;
+ }
}
/**
/* output quantized coefficients and level shift */
error=0;
for(i=0; i<order; i++) {
- error += lpc_in[i] * (1 << sh);
+ 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)
+static int estimate_best_order(double *ref, int min_order, int max_order)
{
int i, est;
- est = 1;
- for(i=max_order-1; i>=0; i--) {
+ est = min_order;
+ for(i=max_order-1; i>=min_order-1; i--) {
if(ref[i] > 0.10) {
est = i+1;
break;
/**
* Calculate LPC coefficients for multiple orders
+ *
+ * @param use_lpc LPC method for determining coefficients
+ * 0 = LPC with fixed pre-defined coeffs
+ * 1 = LPC with coeffs determined by Levinson-Durbin recursion
+ * 2+ = LPC with coeffs determined by Cholesky factorization using (use_lpc-1) passes.
*/
int ff_lpc_calc_coefs(DSPContext *s,
- const int32_t *samples, int blocksize, int min_order, int max_order,
- int precision, int32_t coefs[][MAX_LPC_ORDER],
- int *shift, int use_lpc, int omethod, int max_shift, int zero_shift)
+ const int32_t *samples, int blocksize, int min_order,
+ int max_order, int precision,
+ int32_t coefs[][MAX_LPC_ORDER], int *shift, int use_lpc,
+ int omethod, int max_shift, int zero_shift)
{
double autoc[MAX_LPC_ORDER+1];
double ref[MAX_LPC_ORDER];
int i, j, pass;
int opt_order;
- assert(max_order >= MIN_LPC_ORDER && max_order <= MAX_LPC_ORDER);
+ assert(max_order >= MIN_LPC_ORDER && max_order <= MAX_LPC_ORDER && use_lpc > 0);
if(use_lpc == 1){
- s->flac_compute_autocorr(samples, blocksize, max_order, autoc);
+ s->lpc_compute_autocorr(samples, blocksize, max_order, autoc);
+
+ compute_lpc_coefs(autoc, max_order, &lpc[0][0], MAX_LPC_ORDER, 0, 1);
- compute_lpc_coefs(autoc, max_order, lpc, ref);
+ for(i=0; i<max_order; i++)
+ ref[i] = fabs(lpc[i][i]);
}else{
LLSModel m[2];
- double var[MAX_LPC_ORDER+1], weight;
+ double var[MAX_LPC_ORDER+1], av_uninit(weight);
for(pass=0; pass<use_lpc-1; pass++){
av_init_lls(&m[pass&1], max_order);
for(i=0; i<max_order; i++){
for(j=0; j<max_order; j++)
- lpc[i][j]= m[(pass-1)&1].coeff[i][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--)
opt_order = max_order;
if(omethod == ORDER_METHOD_EST) {
- opt_order = estimate_best_order(ref, max_order);
+ opt_order = estimate_best_order(ref, min_order, max_order);
i = opt_order-1;
quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i], max_shift, zero_shift);
} else {