[xiph-cvs] cvs commit: vorbis/lib analysis.c lsp.c
Monty
xiphmont at xiph.org
Tue Oct 10 15:20:15 PDT 2000
xiphmont 00/10/10 15:20:15
Modified: . Tag: branch_postbeta2 README
lib Tag: branch_postbeta2 analysis.c lsp.c
Log:
A new root polisher in lsp.c for Jack to test with his woogie EGCS.
Monty
Revision Changes Path
No revision
No revision
1.4.8.2 +3 -3 vorbis/README
Index: README
===================================================================
RCS file: /usr/local/cvsroot/vorbis/README,v
retrieving revision 1.4.8.1
retrieving revision 1.4.8.2
diff -u -r1.4.8.1 -r1.4.8.2
--- README 2000/10/05 22:41:25 1.4.8.1
+++ README 2000/10/10 22:20:14 1.4.8.2
@@ -2,8 +2,8 @@
* *
* THIS FILE IS PART OF THE Ogg Vorbis SOFTWARE CODEC SOURCE CODE. *
* USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
-* THE GNU PUBLIC LICENSE 2, WHICH IS INCLUDED WITH THIS SOURCE. *
-* PLEASE READ THESE TERMS DISTRIBUTING. *
+* THE GNU LESSER/LIBRARY PUBLIC LICENSE 2, WHICH IS INCLUDED WITH *
+* THIS SOURCE. PLEASE READ THESE TERMS DISTRIBUTING. *
* *
* THE OggSQUISH SOURCE CODE IS (C) COPYRIGHT 1994-2000 *
* by Monty <monty at xiph.org> and The XIPHOPHORUS Company *
@@ -104,4 +104,4 @@
Monty <monty at xiph.org>
-$Id: README,v 1.4.8.1 2000/10/05 22:41:25 xiphmont Exp $
+$Id: README,v 1.4.8.2 2000/10/10 22:20:14 xiphmont Exp $
No revision
No revision
1.33.2.4 +4 -4 vorbis/lib/analysis.c
Index: analysis.c
===================================================================
RCS file: /usr/local/cvsroot/vorbis/lib/analysis.c,v
retrieving revision 1.33.2.3
retrieving revision 1.33.2.4
diff -u -r1.33.2.3 -r1.33.2.4
--- analysis.c 2000/09/27 06:20:58 1.33.2.3
+++ analysis.c 2000/10/10 22:20:14 1.33.2.4
@@ -12,7 +12,7 @@
********************************************************************
function: single-block PCM analysis mode dispatch
- last mod: $Id: analysis.c,v 1.33.2.3 2000/09/27 06:20:58 jack Exp $
+ last mod: $Id: analysis.c,v 1.33.2.4 2000/10/10 22:20:14 xiphmont Exp $
********************************************************************/
@@ -53,10 +53,10 @@
if(vb->W){
oggpack_write(&vb->opb,vb->lW,1);
oggpack_write(&vb->opb,vb->nW,1);
- fprintf(stderr,"*");
- }else{
+ /*fprintf(stderr,"*");*/
+ }/*else{
fprintf(stderr,".");
- }
+ }*/
if(_mapping_P[type]->forward(vb,vd->mode[mode]))
return(-1);
1.9.2.6 +76 -20 vorbis/lib/lsp.c
Index: lsp.c
===================================================================
RCS file: /usr/local/cvsroot/vorbis/lib/lsp.c,v
retrieving revision 1.9.2.5
retrieving revision 1.9.2.6
diff -u -r1.9.2.5 -r1.9.2.6
--- lsp.c 2000/10/05 21:35:40 1.9.2.5
+++ lsp.c 2000/10/10 22:20:14 1.9.2.6
@@ -12,7 +12,7 @@
********************************************************************
function: LSP (also called LSF) conversion routines
- last mod: $Id: lsp.c,v 1.9.2.5 2000/10/05 21:35:40 xiphmont Exp $
+ last mod: $Id: lsp.c,v 1.9.2.6 2000/10/10 22:20:14 xiphmont Exp $
The LSP generation code is taken (with minimal modification) from
"On the Computation of the LSP Frequencies" by Joseph Rothweiler
@@ -241,31 +241,87 @@
return(-1);
}
-/* CACM algorithm 283. */
-static void cacm283(float *a,int ord,float *r){
- int i, k;
- double val, p, delta, error;
- double rooti;
+/* This is one of those 'mathemeticians should not write code' kind of
+ cases. Newton's method of polishing roots is straightforward
+ enough... except in those cases where it just fails in the real
+ world. In our case below, we're worried about a local mini/maxima
+ shooting a root estimation off to infinity, or the new estimation
+ chaotically oscillating about convergence (shouldn't actually be a
+ problem in our usage.
+
+ Maehly's modification (zero suppression, to prevent two tenative
+ roots from collapsing to the same actual root) similarly can
+ temporarily shoot a root off toward infinity. It would come
+ back... if it were not for the fact that machine representation has
+ limited dynamic range and resolution. This too is guarded by
+ limiting delta.
+
+ Last problem is convergence criteria; we don't know what a 'double'
+ is on our hardware/compiler, and the convergence limit is bounded
+ by roundoff noise. So, we hack convergence:
+
+ Require at most 1e-6 mean squared error for all zeroes. When
+ converging, start the clock ticking at 1e-6; limit our polishing to
+ as many more iterations as took us to get this far, 100 max.
+
+ Past max iters, quit when MSE is no longer decreasing *or* we go
+ below ~1e-20 MSE, whichever happens first. */
+
+static void Newton_Raphson_Maehly(float *a,int ord,float *r){
+ int i, k, count=0, maxiter=0;
+ double error=1.,besterror=1.;
+ double *root=alloca(ord*sizeof(double));
- for(i=0; i<ord;i++) r[i] = 2.0 * (i+0.5) / ord - 1.0;
+ for(i=0; i<ord;i++) root[i] = 2.0 * (i+0.5) / ord - 1.0;
- for(error=1 ; error > 1.e-12; ) {
- error = 0;
- for( i=0; i<ord; i++) { /* Update each point. */
- rooti = r[i];
- val = a[ord];
- p = a[ord];
+ while(error>1.e-20){
+ error=0;
+
+ for(i=0; i<ord; i++) { /* Update each point. */
+ double ac=0.,pp=0.,delta;
+ double rooti=root[i];
+ double p=a[ord];
for(k=ord-1; k>= 0; k--) {
- val = val * rooti + a[k];
- if (k != i) p *= rooti - r[k];
+
+ pp= pp* rooti + p;
+ p = p * rooti+ a[k];
+ if (k != i) ac += 1./(rooti - root[k]);
}
- delta = val/p;
- r[i] -= delta;
+ ac=p*ac;
+
+ delta = p/(pp-ac);
+
+ /* don't allow the correction to scream off into infinity if we
+ happened to polish right at a local mini/maximum */
+ if(delta<-3)delta=-3;
+ if(delta>3.)delta=3.; /* 3 is not a random choice; it's large
+ enough to make sure the first pass
+ can't accidentally limit two poles to
+ the same value in a fatal nonelastic
+ collision. */
+
+ root[i] -= delta;
error += delta*delta;
}
+
+ if(maxiter && count>maxiter && error>=besterror)break;
+
+ /* anything to help out the polisher; converge using doubles */
+ if(!count || error<besterror){
+ for(i=0; i<ord; i++) r[i]=root[i];
+ besterror=error;
+ if(error<1.e-6){ /* rough minimum criteria */
+ maxiter=count*2+10;
+ if(maxiter>100)maxiter=100;
+ }
+ }
+
+ count++;
}
-
+
+ fprintf(stderr,"***** error=%g, count=%d\n",error,count);
+
/* Replaced the original bubble sort with a real sort. With your
help, we can eliminate the bubble sort in our lifetime. --Monty */
@@ -301,8 +357,8 @@
/* Find the roots of the 2 even polynomials.*/
- cacm283(g1,order2,g1r);
- cacm283(g2,order2,g2r);
+ Newton_Raphson_Maehly(g1,order2,g1r);
+ Newton_Raphson_Maehly(g2,order2,g2r);
for(i=0;i<m;i+=2){
lsp[i] = acos(g1r[i/2]);
--- >8 ----
List archives: http://www.xiph.org/archives/
Ogg project homepage: http://www.xiph.org/ogg/
To unsubscribe from this list, send a message to 'cvs-request at xiph.org'
containing only the word 'unsubscribe' in the body. No subject is needed.
Unsubscribe messages sent to the list will be ignored/filtered.
More information about the commits
mailing list