summaryrefslogtreecommitdiff
path: root/libm/double/levnsn.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
committerEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
commit1077fa4d772832f77a677ce7fb7c2d513b959e3f (patch)
tree579bee13fb0b58d2800206366ec2caecbb15f3fc /libm/double/levnsn.c
parent22358dd7ce7bb49792204b698f01a6f69b9c8e08 (diff)
uClibc now has a math library. muahahahaha!
-Erik
Diffstat (limited to 'libm/double/levnsn.c')
-rw-r--r--libm/double/levnsn.c82
1 files changed, 82 insertions, 0 deletions
diff --git a/libm/double/levnsn.c b/libm/double/levnsn.c
new file mode 100644
index 000000000..3fda5d6bd
--- /dev/null
+++ b/libm/double/levnsn.c
@@ -0,0 +1,82 @@
+/* Levnsn.c */
+/* Levinson-Durbin LPC
+ *
+ * | R0 R1 R2 ... RN-1 | | A1 | | -R1 |
+ * | R1 R0 R1 ... RN-2 | | A2 | | -R2 |
+ * | R2 R1 R0 ... RN-3 | | A3 | = | -R3 |
+ * | ... | | ...| | ... |
+ * | RN-1 RN-2... R0 | | AN | | -RN |
+ *
+ * Ref: John Makhoul, "Linear Prediction, A Tutorial Review"
+ * Proc. IEEE Vol. 63, PP 561-580 April, 1975.
+ *
+ * R is the input autocorrelation function. R0 is the zero lag
+ * term. A is the output array of predictor coefficients. Note
+ * that a filter impulse response has a coefficient of 1.0 preceding
+ * A1. E is an array of mean square error for each prediction order
+ * 1 to N. REFL is an output array of the reflection coefficients.
+ */
+
+#define abs(x) ( (x) < 0 ? -(x) : (x) )
+
+int levnsn( n, r, a, e, refl )
+int n;
+double r[], a[], e[], refl[];
+{
+int k, km1, i, kmi, j;
+double ai, akk, err, err1, r0, t, akmi;
+double *pa, *pr;
+
+for( i=0; i<n; i++ )
+ {
+ a[i] = 0.0;
+ e[i] = 0.0;
+ refl[i] = 0.0;
+ }
+r0 = r[0];
+e[0] = r0;
+err = r0;
+
+akk = -r[1]/err;
+err = (1.0 - akk*akk) * err;
+e[1] = err;
+a[1] = akk;
+refl[1] = akk;
+
+if( err < 1.0e-2 )
+ return 0;
+
+for( k=2; k<n; k++ )
+ {
+ t = 0.0;
+ pa = &a[1];
+ pr = &r[k-1];
+ for( j=1; j<k; j++ )
+ t += *pa++ * *pr--;
+ akk = -( r[k] + t )/err;
+ refl[k] = akk;
+ km1 = k/2;
+ for( j=1; j<=km1; j++ )
+ {
+ kmi = k-j;
+ ai = a[j];
+ akmi = a[kmi];
+ a[j] = ai + akk*akmi;
+ if( i == kmi )
+ goto nxtk;
+ a[kmi] = akmi + akk*ai;
+ }
+nxtk:
+ a[k] = akk;
+ err1 = (1.0 - akk*akk)*err;
+ e[k] = err1;
+ if( err1 < 0 )
+ err1 = -err1;
+/* err1 = abs(err1);*/
+/* if( (err1 < 1.0e-2) || (err1 >= err) )*/
+ if( err1 < 1.0e-2 )
+ return 0;
+ err = err1;
+ }
+ return 0;
+}