X-Git-Url: http://review.tizen.org/git/?a=blobdiff_plain;f=lib%2Flsp.c;h=de4ff73e3893d4f44addac8d5c8ee3be409bc240;hb=143caf4023a90c09a5eb685fdd46fb9b9c36b1ee;hp=1adab140d82515df70144830eb25cbca8b063cfb;hpb=b1ab109c9a700c85aafde936a4ac9cdcaff045bf;p=platform%2Fupstream%2Flibvorbis.git diff --git a/lib/lsp.c b/lib/lsp.c index 1adab14..de4ff73 100644 --- a/lib/lsp.c +++ b/lib/lsp.c @@ -1,24 +1,24 @@ /******************************************************************** * * - * 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. * + * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. * + * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS * + * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE * + * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. * * * - * THE OggSQUISH SOURCE CODE IS (C) COPYRIGHT 1994-2000 * - * by Monty and The XIPHOPHORUS Company * - * http://www.xiph.org/ * + * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2009 * + * by the Xiph.Org Foundation http://www.xiph.org/ * * * ******************************************************************** function: LSP (also called LSF) conversion routines - last mod: $Id: lsp.c,v 1.10 2000/10/12 03:12:53 xiphmont Exp $ + last mod: $Id$ - The LSP generation code is taken (with minimal modification) from - "On the Computation of the LSP Frequencies" by Joseph Rothweiler - , available at: - - http://www2.xtdl.com/~rothwlr/lsfpaper/lsfpage.html + The LSP generation code is taken (with minimal modification and a + few bugfixes) from "On the Computation of the LSP Frequencies" by + Joseph Rothweiler (see http://www.rothweiler.us for contact info). + The paper is available at: + + http://www.myown1.com/joe/lsf ********************************************************************/ @@ -46,14 +46,15 @@ implementation. The float lookup is likely the optimal choice on any machine with an FPU. The integer implementation is *not* fixed point (due to the need for a large dynamic range and thus a - seperately tracked exponent) and thus much more complex than the + separately tracked exponent) and thus much more complex than the relatively simple float implementations. It's mostly for future work on a fully fixed point implementation for processors like the ARM family. */ -/* undefine both for the 'old' but more precise implementation */ -#define FLOAT_LOOKUP -#undef INT_LOOKUP +/* define either of these (preferably FLOAT_LOOKUP) to have faster + but less precise implementation. */ +#undef FLOAT_LOOKUP +#undef INT_LOOKUP #ifdef FLOAT_LOOKUP #include "lookup.c" /* catch this in the build system; we #include for @@ -62,31 +63,53 @@ /* side effect: changes *lsp to cosines of lsp */ void vorbis_lsp_to_curve(float *curve,int *map,int n,int ln,float *lsp,int m, - float amp,float ampoffset){ + float amp,float ampoffset){ int i; float wdel=M_PI/ln; + vorbis_fpu_control fpu; + + vorbis_fpu_setround(&fpu); for(i=0;i>1; + + while(c--){ + q*=ftmp[0]-w; + p*=ftmp[1]-w; + ftmp+=2; + } - for(j=0;j>25])) - if(!(shift=MLOOP_2[(pi|qi)>>19])) - shift=MLOOP_3[(pi|qi)>>16]; + if(!(shift=MLOOP_2[(pi|qi)>>19])) + shift=MLOOP_3[(pi|qi)>>16]; + qi=(qi>>shift)*labs(ilsp[j-1]-wi); pi=(pi>>shift)*labs(ilsp[j]-wi); - qi=(qi>>shift)*labs(ilsp[j+1]-wi); qexp+=shift; } if(!(shift=MLOOP_1[(pi|qi)>>25])) if(!(shift=MLOOP_2[(pi|qi)>>19])) - shift=MLOOP_3[(pi|qi)>>16]; - pi>>=shift; - qi>>=shift; - qexp+=shift-7*m; + shift=MLOOP_3[(pi|qi)>>16]; /* pi,qi normalized collectively, both tracked using qexp */ - /* p*=p(1-w), q*=q(1+w), let normalization drift because it isn't - worth tracking step by step */ + if(m&1){ + /* odd order filter; slightly assymetric */ + /* the last coefficient */ + qi=(qi>>shift)*labs(ilsp[j-1]-wi); + pi=(pi>>shift)<<14; + qexp+=shift; + + if(!(shift=MLOOP_1[(pi|qi)>>25])) + if(!(shift=MLOOP_2[(pi|qi)>>19])) + shift=MLOOP_3[(pi|qi)>>16]; + + pi>>=shift; + qi>>=shift; + qexp+=shift-14*((m+1)>>1); + + pi=((pi*pi)>>16); + qi=((qi*qi)>>16); + qexp=qexp*2+m; + + pi*=(1<<14)-((wi*wi)>>14); + qi+=pi>>14; + + }else{ + /* even order filter; still symmetric */ + + /* p*=p(1-w), q*=q(1+w), let normalization drift because it isn't + worth tracking step by step */ + + pi>>=shift; + qi>>=shift; + qexp+=shift-7*m; - pi=((pi*pi)>>16); - qi=((qi*qi)>>16); - qexp=qexp*2+m; + pi=((pi*pi)>>16); + qi=((qi*qi)>>16); + qexp=qexp*2+m; + + pi*=(1<<14)-wi; + qi*=(1<<14)+wi; + qi=(qi+pi)>>14; + + } - qi*=(1<<14)-wi; - pi*=(1<<14)+wi; - - qi=(qi+pi)>>14; /* we've let the normalization drift because it wasn't important; however, for the lookup, things must be normalized again. We need at most one right shift or a number of left shifts */ if(qi&0xffff0000){ /* checks for 1.xxxxxxxxxxxxxxxx */ - qi>>=1; qexp++; + qi>>=1; qexp++; }else while(qi && !(qi&0x8000)){ /* checks for 0.0xxxxxxxxxxxxxxx or less*/ - qi<<=1; qexp--; + qi<<=1; qexp--; } amp=vorbis_fromdBlook_i(ampi* /* n.4 */ - vorbis_invsqlook_i(qi,qexp)- - /* m.8, m+n<=8 */ - ampoffseti); /* 8.12[0] */ + vorbis_invsqlook_i(qi,qexp)- + /* m.8, m+n<=8 */ + ampoffseti); /* 8.12[0] */ - curve[i]=amp; - while(map[++i]==k)curve[i]=amp; + curve[i]*=amp; + while(map[++i]==k)curve[i]*=amp; } } -#else +#else /* old, nonoptimized but simple version for any poor sap who needs to - figure out what the hell this code does, or wants the other tiny + figure out what the hell this code does, or wants the other fraction of a dB precision */ /* side effect: changes *lsp to cosines of lsp */ void vorbis_lsp_to_curve(float *curve,int *map,int n,int ln,float *lsp,int m, - float amp,float ampoffset){ + float amp,float ampoffset){ int i; float wdel=M_PI/ln; - for(i=0;i= i; j--) { g[j-2] -= g[j]; - g[j] += g[j]; + g[j] += g[j]; } } } static int comp(const void *a,const void *b){ - if(*(float *)a<*(float *)b) - return(1); - else - return(-1); + return (*(float *)a<*(float *)b)-(*(float *)a>*(float *)b); } -/* 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; i1.e-20){ +/* Newton-Raphson-Maehly actually functioned as a decent root finder, + but there are root sets for which it gets into limit cycles + (exacerbated by zero suppression) and fails. We can't afford to + fail, even if the failure is 1 in 100,000,000, so we now use + Laguerre and later polish with Newton-Raphson (which can then + afford to fail) */ + +#define EPSILON 10e-7 +static int Laguerre_With_Deflation(float *a,int ord,float *r){ + int i,m; + double *defl=alloca(sizeof(*defl)*(ord+1)); + for(i=0;i<=ord;i++)defl[i]=a[i]; + + for(m=ord;m>0;m--){ + double new=0.f,delta; + + /* iterate a root */ + while(1){ + double p=defl[m],pp=0.f,ppp=0.f,denom; + + /* eval the polynomial and its first two derivatives */ + for(i=m;i>0;i--){ + ppp = new*ppp + pp; + pp = new*pp + p; + p = new*p + defl[i-1]; + } + + /* Laguerre's method */ + denom=(m-1) * ((m-1)*pp*pp - m*p*ppp); + if(denom<0) + return(-1); /* complex root! The LPC generator handed us a bad filter */ + + if(pp>0){ + denom = pp + sqrt(denom); + if(denom-(EPSILON))denom=-(EPSILON); + } + + delta = m*p/denom; + new -= delta; + + if(delta<0.f)delta*=-1; + + if(fabs(delta/new)<10e-12)break; + } + + r[m-1]=new; + + /* forward deflation */ + + for(i=m;i>0;i--) + defl[i-1]+=new*defl[i]; + defl++; + + } + return(0); +} + + +/* for spit-and-polish only */ +static int Newton_Raphson(float *a,int ord,float *r){ + int i, k, count=0; + double error=1.f; + double *root=alloca(ord*sizeof(*root)); + + for(i=0; i1e-20){ error=0; - + for(i=0; i= 0; k--) { - pp= pp* rooti + p; - p = p * rooti+ a[k]; - if (k != i) ac += 1./(rooti - root[k]); + pp= pp* rooti + p; + p = p * rooti + a[k]; } - 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. */ + delta = p/pp; root[i] -= delta; - error += delta*delta; - } - - if(maxiter && count>maxiter && error>=besterror)break; - - /* anything to help out the polisher; converge using doubles */ - if(!count || error100)maxiter=100; - } + error+= delta*delta; } + if(count>40)return(-1); + count++; } /* Replaced the original bubble sort with a real sort. With your help, we can eliminate the bubble sort in our lifetime. --Monty */ - - qsort(r,ord,sizeof(float),comp); + for(i=0; i>1; + int g1_order,g2_order; + float *g1=alloca(sizeof(*g1)*(order2+1)); + float *g2=alloca(sizeof(*g2)*(order2+1)); + float *g1r=alloca(sizeof(*g1r)*(order2+1)); + float *g2r=alloca(sizeof(*g2r)*(order2+1)); int i; + /* even and odd are slightly different base cases */ + g1_order=(m+1)>>1; + g2_order=(m) >>1; + /* Compute the lengths of the x polynomials. */ /* Compute the first half of K & R F1 & F2 polynomials. */ /* Compute half of the symmetric and antisymmetric polynomials. */ /* Remove the roots at +1 and -1. */ - - g1[order2] = 1.0; - for(i=0;ig2_order){ + for(i=2; i<=g2_order;i++) g2[g2_order-i] += g2[g2_order-i+2]; + }else{ + for(i=1; i<=g1_order;i++) g1[g1_order-i] -= g1[g1_order-i+1]; + for(i=1; i<=g2_order;i++) g2[g2_order-i] += g2[g2_order-i+1]; + } /* Convert into polynomials in cos(alpha) */ - cheby(g1,order2); - cheby(g2,order2); + cheby(g1,g1_order); + cheby(g2,g2_order); /* Find the roots of the 2 even polynomials.*/ - - Newton_Raphson_Maehly(g1,order2,g1r); - Newton_Raphson_Maehly(g2,order2,g2r); - - for(i=0;i