X-Git-Url: http://review.tizen.org/git/?a=blobdiff_plain;f=lib%2Flsp.c;h=de4ff73e3893d4f44addac8d5c8ee3be409bc240;hb=143caf4023a90c09a5eb685fdd46fb9b9c36b1ee;hp=fe95fc9a9da6781725fe048bffa161342ab40a3f;hpb=63bb426cb18452a6576d67abbaaa859df28d41b9;p=platform%2Fupstream%2Flibvorbis.git diff --git a/lib/lsp.c b/lib/lsp.c index fe95fc9..de4ff73 100644 --- a/lib/lsp.c +++ b/lib/lsp.c @@ -5,19 +5,20 @@ * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE * * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. * * * - * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2001 * - * by 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.20 2001/12/20 01:00:27 segher Exp $ + last mod: $Id$ The LSP generation code is taken (with minimal modification and a few bugfixes) from "On the Computation of the LSP Frequencies" by - Joseph Rothweiler , available at: - - http://www2.xtdl.com/~rothwlr/lsfpaper/lsfpage.html + Joseph Rothweiler (see http://www.rothweiler.us for contact info). + The paper is available at: + + http://www.myown1.com/joe/lsf ********************************************************************/ @@ -45,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 */ -#undef 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 @@ -61,11 +63,11 @@ /* 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; - do{ + while(c--){ q*=ftmp[0]-w; p*=ftmp[1]-w; ftmp+=2; - }while(--c); + } if(m&1){ /* odd order filter; slightly assymetric */ @@ -98,10 +100,10 @@ void vorbis_lsp_to_curve(float *curve,int *map,int n,int ln,float *lsp,int m, } q=frexp(p+q,&qexp); - q=vorbis_fromdBlook(amp* - vorbis_invsqlook(q)* - vorbis_invsq2explook(qexp+m)- - ampoffset); + q=vorbis_fromdBlook(amp* + vorbis_invsqlook(q)* + vorbis_invsq2explook(qexp+m)- + ampoffset); do{ curve[i++]*=q; @@ -117,26 +119,26 @@ void vorbis_lsp_to_curve(float *curve,int *map,int n,int ln,float *lsp,int m, compilers (like gcc) that can't inline across modules */ -static int MLOOP_1[64]={ +static const int MLOOP_1[64]={ 0,10,11,11, 12,12,12,12, 13,13,13,13, 13,13,13,13, 14,14,14,14, 14,14,14,14, 14,14,14,14, 14,14,14,14, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, }; -static int MLOOP_2[64]={ +static const int MLOOP_2[64]={ 0,4,5,5, 6,6,6,6, 7,7,7,7, 7,7,7,7, 8,8,8,8, 8,8,8,8, 8,8,8,8, 8,8,8,8, 9,9,9,9, 9,9,9,9, 9,9,9,9, 9,9,9,9, 9,9,9,9, 9,9,9,9, 9,9,9,9, 9,9,9,9, }; -static int MLOOP_3[8]={0,1,2,2,3,3,3,3}; +static const int MLOOP_3[8]={0,1,2,2,3,3,3,3}; /* 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){ /* 0 <= m < 256 */ @@ -160,15 +162,15 @@ void vorbis_lsp_to_curve(float *curve,int *map,int n,int ln,float *lsp,int m, for(j=3;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); qexp+=shift; } if(!(shift=MLOOP_1[(pi|qi)>>25])) if(!(shift=MLOOP_2[(pi|qi)>>19])) - shift=MLOOP_3[(pi|qi)>>16]; + shift=MLOOP_3[(pi|qi)>>16]; /* pi,qi normalized collectively, both tracked using qexp */ @@ -180,9 +182,9 @@ void vorbis_lsp_to_curve(float *curve,int *map,int n,int ln,float *lsp,int m, qexp+=shift; if(!(shift=MLOOP_1[(pi|qi)>>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]; + pi>>=shift; qi>>=shift; qexp+=shift-14*((m+1)>>1); @@ -198,8 +200,8 @@ void vorbis_lsp_to_curve(float *curve,int *map,int n,int ln,float *lsp,int m, /* 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 */ - + worth tracking step by step */ + pi>>=shift; qi>>=shift; qexp+=shift-7*m; @@ -207,36 +209,36 @@ void vorbis_lsp_to_curve(float *curve,int *map,int n,int ln,float *lsp,int 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; - + } - + /* 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; } } -#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 @@ -244,7 +246,7 @@ void vorbis_lsp_to_curve(float *curve,int *map,int n,int ln,float *lsp,int m, /* 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); } /* Newton-Raphson-Maehly actually functioned as a decent root finder, @@ -310,7 +309,6 @@ static int comp(const void *a,const void *b){ #define EPSILON 10e-7 static int Laguerre_With_Deflation(float *a,int ord,float *r){ int i,m; - double lastdelta=0.f; double *defl=alloca(sizeof(*defl)*(ord+1)); for(i=0;i<=ord;i++)defl[i]=a[i]; @@ -320,25 +318,25 @@ static int Laguerre_With_Deflation(float *a,int ord,float *r){ /* 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]; + 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 */ + return(-1); /* complex root! The LPC generator handed us a bad filter */ if(pp>0){ - denom = pp + sqrt(denom); - if(denom-(EPSILON))denom=-(EPSILON); + denom = pp - sqrt(denom); + if(denom>-(EPSILON))denom=-(EPSILON); } delta = m*p/denom; @@ -346,14 +344,13 @@ static int Laguerre_With_Deflation(float *a,int ord,float *r){ if(delta<0.f)delta*=-1; - if(fabs(delta/new)<10e-12)break; - lastdelta=delta; + 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++; @@ -370,27 +367,27 @@ static int Newton_Raphson(float *a,int ord,float *r){ 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]; + pp= pp* rooti + p; + p = p * rooti + a[k]; } delta = p/pp; root[i] -= delta; error+= delta*delta; } - + if(count>40)return(-1); - + count++; } @@ -420,12 +417,12 @@ int vorbis_lpc_to_lsp(float *lpc,float *lsp,int m){ /* 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[g1_order] = 1.f; for(i=1;i<=g1_order;i++) g1[g1_order-i] = lpc[i-1]+lpc[m-i]; g2[g2_order] = 1.f; for(i=1;i<=g2_order;i++) g2[g2_order-i] = lpc[i-1]-lpc[m-i]; - + if(g1_order>g2_order){ for(i=2; i<=g2_order;i++) g2[g2_order-i] += g2[g2_order-i+2]; }else{