Commit 408ec4e2 authored by Michael Niedermayer's avatar Michael Niedermayer

calculate all coefficients for several orders during cholesky factorization,...

calculate all coefficients for several orders during cholesky factorization, the resulting coefficients are not strictly optimal though as there is a small difference in the autocorrelation matrixes which is ignored for the smaller orders

Originally committed as revision 5758 to svn://svn.ffmpeg.org/ffmpeg/trunk
parent 6ce704bb
......@@ -742,35 +742,41 @@ static int lpc_calc_coefs(const int32_t *samples, int blocksize, int max_order,
compute_autocorr(samples, blocksize, max_order+1, autoc);
compute_lpc_coefs(autoc, max_order, lpc, ref);
opt_order = estimate_best_order(ref, max_order);
}else{
LLSModel m[2];
double var[MAX_LPC_ORDER+1], eval;
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);
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);
opt_order= max_order; //FIXME
av_solve_lls(&m[pass&1], 0.001, 0);
}
for(i=0; i<opt_order; i++)
lpc[opt_order-1][i]= m[(pass-1)&1].coeff[i];
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 = estimate_best_order(ref, max_order);
i = opt_order-1;
quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i]);
......
......@@ -49,12 +49,11 @@ void av_update_lls(LLSModel *m, double *var, double decay){
}
}
double av_solve_lls(LLSModel *m, double threshold){
void av_solve_lls(LLSModel *m, double threshold, int min_order){
int i,j,k;
double (*factor)[MAX_VARS+1]= &m->covariance[1][0];
double (*covar )[MAX_VARS+1]= &m->covariance[1][1];
double *covar_y = m->covariance[0];
double variance;
int count= m->indep_count;
for(i=0; i<count; i++){
......@@ -75,33 +74,34 @@ double av_solve_lls(LLSModel *m, double threshold){
for(i=0; i<count; i++){
double sum= covar_y[i+1];
for(k=i-1; k>=0; k--)
sum -= factor[i][k]*m->coeff[k];
m->coeff[i]= sum / factor[i][i];
sum -= factor[i][k]*m->coeff[0][k];
m->coeff[0][i]= sum / factor[i][i];
}
for(i=count-1; i>=0; i--){
double sum= m->coeff[i];
for(k=i+1; k<count; k++)
sum -= factor[k][i]*m->coeff[k];
m->coeff[i]= sum / factor[i][i];
}
for(j=count-1; j>=min_order; j--){
for(i=j; i>=0; i--){
double sum= m->coeff[0][i];
for(k=i+1; k<=j; k++)
sum -= factor[k][i]*m->coeff[j][k];
m->coeff[j][i]= sum / factor[i][i];
}
variance= covar_y[0];
for(i=0; i<count; i++){
double sum= m->coeff[i]*covar[i][i] - 2*covar_y[i+1];
for(j=0; j<i; j++)
sum += 2*m->coeff[j]*covar[j][i];
variance += m->coeff[i]*sum;
m->variance[j]= covar_y[0];
for(i=0; i<=j; i++){
double sum= m->coeff[j][i]*covar[i][i] - 2*covar_y[i+1];
for(k=0; k<i; k++)
sum += 2*m->coeff[j][k]*covar[k][i];
m->variance[j] += m->coeff[j][i]*sum;
}
}
return variance;
}
double av_evaluate_lls(LLSModel *m, double *param){
double av_evaluate_lls(LLSModel *m, double *param, int order){
int i;
double out= 0;
for(i=0; i<m->indep_count; i++)
out+= param[i]*m->coeff[i];
for(i=0; i<=order; i++)
out+= param[i]*m->coeff[order][i];
return out;
}
......@@ -113,27 +113,35 @@ double av_evaluate_lls(LLSModel *m, double *param){
int main(){
LLSModel m;
int i;
int i, order;
av_init_lls(&m, 3);
for(i=0; i<100; i++){
double var[4];
double eval, variance;
#if 0
var[1] = rand() / (double)RAND_MAX;
var[2] = rand() / (double)RAND_MAX;
var[3] = rand() / (double)RAND_MAX;
var[2]= var[1] + var[3];
var[2]= var[1] + var[3]/2;
var[0] = var[1] + var[2] + var[3] + var[1]*var[2]/100;
eval= av_evaluate_lls(&m, var+1);
#else
var[0] = (rand() / (double)RAND_MAX - 0.5)*2;
var[1] = var[0] + rand() / (double)RAND_MAX - 0.5;
var[2] = var[1] + rand() / (double)RAND_MAX - 0.5;
var[3] = var[2] + rand() / (double)RAND_MAX - 0.5;
#endif
av_update_lls(&m, var, 0.99);
variance= av_solve_lls(&m, 0.001);
av_log(NULL, AV_LOG_DEBUG, "real:%f pred:%f var:%f coeffs:%f %f %f\n",
var[0], eval, sqrt(variance / (i+1)),
m.coeff[0], m.coeff[1], m.coeff[2]);
av_solve_lls(&m, 0.001, 0);
for(order=0; order<3; order++){
eval= av_evaluate_lls(&m, var+1, order);
av_log(NULL, AV_LOG_DEBUG, "real:%f order:%d pred:%f var:%f coeffs:%f %f %f\n",
var[0], order, eval, sqrt(m.variance[order] / (i+1)),
m.coeff[order][0], m.coeff[order][1], m.coeff[order][2]);
}
}
return 0;
}
......
......@@ -30,13 +30,14 @@
*/
typedef struct LLSModel{
double covariance[MAX_VARS+1][MAX_VARS+1];
double coeff[MAX_VARS];
double coeff[MAX_VARS][MAX_VARS];
double variance[MAX_VARS];
int indep_count;
}LLSModel;
void av_init_lls(LLSModel *m, int indep_count);
void av_update_lls(LLSModel *m, double *param, double decay);
double av_solve_lls(LLSModel *m, double threshold);
double av_evaluate_lls(LLSModel *m, double *param);
void av_solve_lls(LLSModel *m, double threshold, int min_order);
double av_evaluate_lls(LLSModel *m, double *param, int order);
#endif
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment