# HG changeset patch
# User michael
# Date 1152902918 0
# Node ID 33af013504d506a28a4c351ace05793e976f2791
# Parent  af598a3e23407232b86624b29e06095112c4fe76
optionally (use_lpc=2) support Cholesky factorization for finding the lpc coeficients
  this will find the coefficients which minimize the sum of the squared errors,
  levinson-durbin recursion OTOH is only strictly correct if the autocorrelation matrix is a
  toeplitz matrix which it is only if the blocksize is infinite, this is also why applying
  a window (like the welch winodw we currently use) improves the lpc coefficients generated
  by levinson-durbin recursion ...

optionally (use_lpc>2) support iterative linear least abs() solver using cholesky
  factorization with adjusted weights in each iteration

compression gain for both is small, and multiple passes are of course dead slow

diff -r af598a3e2340 -r 33af013504d5 flacenc.c
--- a/flacenc.c	Fri Jul 14 18:38:23 2006 +0000
+++ b/flacenc.c	Fri Jul 14 18:48:38 2006 +0000
@@ -21,6 +21,7 @@
 #include "bitstream.h"
 #include "crc.h"
 #include "golomb.h"
+#include "lls.h"
 
 #define FLAC_MAX_CH  8
 #define FLAC_MIN_BLOCKSIZE  16
@@ -236,10 +237,12 @@
 
     /* set compression option overrides from AVCodecContext */
     if(avctx->use_lpc >= 0) {
-        s->options.use_lpc = !!avctx->use_lpc;
+        s->options.use_lpc = clip(avctx->use_lpc, 0, 11);
     }
-    av_log(avctx, AV_LOG_DEBUG, " use lpc: %s\n",
-           s->options.use_lpc? "yes" : "no");
+    if(s->options.use_lpc == 1)
+        av_log(avctx, AV_LOG_DEBUG, " use lpc: Levinson-Durbin recursion with Welch window\n");
+    else if(s->options.use_lpc > 1)
+        av_log(avctx, AV_LOG_DEBUG, " use lpc: Cholesky factorization\n");
 
     if(avctx->min_prediction_order >= 0) {
         if(s->options.use_lpc) {
@@ -725,21 +728,49 @@
  */
 static int lpc_calc_coefs(const int32_t *samples, int blocksize, int max_order,
                           int precision, int32_t coefs[][MAX_LPC_ORDER],
-                          int *shift)
+                          int *shift, int use_lpc)
 {
     double autoc[MAX_LPC_ORDER+1];
     double ref[MAX_LPC_ORDER];
     double lpc[MAX_LPC_ORDER][MAX_LPC_ORDER];
-    int i;
+    int i, j, pass;
     int opt_order;
 
     assert(max_order >= MIN_LPC_ORDER && max_order <= MAX_LPC_ORDER);
 
-    compute_autocorr(samples, blocksize, max_order+1, autoc);
+    if(use_lpc == 1){
+        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;
+
+        for(pass=0; pass<use_lpc-1; pass++){
+            av_init_lls(&m[pass&1], max_order/*3*/);
+
+            for(i=max_order; i<blocksize; i++){
+                for(j=0; j<=max_order; j++)
+                    var[j]= samples[i-j];
 
-    compute_lpc_coefs(autoc, max_order, lpc, ref);
+                if(pass){
+                    eval= av_evaluate_lls(&m[(pass-1)&1], var+1);
+                    eval= (512>>pass) + fabs(eval - var[0]);
+                    for(j=0; j<=max_order; j++)
+                        var[j]= samples[i-j] / sqrt(eval);
+                }
 
-    opt_order = estimate_best_order(ref, max_order);
+                av_update_lls(&m[pass&1], var, 1.0);
+            }
+            av_solve_lls(&m[pass&1], 0.001);
+            opt_order= max_order; //FIXME
+        }
+
+        for(i=0; i<opt_order; i++)
+            lpc[opt_order-1][i]= m[(pass-1)&1].coeff[i];
+    }
 
     i = opt_order-1;
     quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i]);
@@ -865,7 +896,7 @@
     }
 
     /* LPC */
-    sub->order = lpc_calc_coefs(smp, n, max_order, precision, coefs, shift);
+    sub->order = lpc_calc_coefs(smp, n, max_order, precision, coefs, shift, ctx->options.use_lpc);
     sub->type = FLAC_SUBFRAME_LPC;
     sub->type_code = sub->type | (sub->order-1);
     sub->shift = shift[sub->order-1];