1 /********************************************************************
3 * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. *
4 * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS *
5 * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
6 * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. *
8 * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2002 *
9 * by the XIPHOPHORUS Company http://www.xiph.org/ *
11 ********************************************************************
13 function: LSP (also called LSF) conversion routines
14 last mod: $Id: lsp.c,v 1.23 2002/10/11 07:44:28 xiphmont Exp $
16 The LSP generation code is taken (with minimal modification and a
17 few bugfixes) from "On the Computation of the LSP Frequencies" by
18 Joseph Rothweiler <rothwlr@altavista.net>, available at:
20 http://www2.xtdl.com/~rothwlr/lsfpaper/lsfpage.html
22 ********************************************************************/
24 /* Note that the lpc-lsp conversion finds the roots of polynomial with
25 an iterative root polisher (CACM algorithm 283). It *is* possible
26 to confuse this algorithm into not converging; that should only
27 happen with absurdly closely spaced roots (very sharp peaks in the
28 LPC f response) which in turn should be impossible in our use of
29 the code. If this *does* happen anyway, it's a bug in the floor
30 finder; find the cause of the confusion (probably a single bin
31 spike or accidental near-float-limit resolution problems) and
43 /* three possible LSP to f curve functions; the exact computation
44 (float), a lookup based float implementation, and an integer
45 implementation. The float lookup is likely the optimal choice on
46 any machine with an FPU. The integer implementation is *not* fixed
47 point (due to the need for a large dynamic range and thus a
48 seperately tracked exponent) and thus much more complex than the
49 relatively simple float implementations. It's mostly for future
50 work on a fully fixed point implementation for processors like the
53 /* undefine both for the 'old' but more precise implementation */
58 #include "lookup.c" /* catch this in the build system; we #include for
59 compilers (like gcc) that can't inline across
62 /* side effect: changes *lsp to cosines of lsp */
63 void vorbis_lsp_to_curve(float *curve,int *map,int n,int ln,float *lsp,int m,
64 float amp,float ampoffset){
67 vorbis_fpu_control fpu;
69 vorbis_fpu_setround(&fpu);
70 for(i=0;i<m;i++)lsp[i]=vorbis_coslook(lsp[i]);
78 float w=vorbis_coslook(wdel*k);
89 /* odd order filter; slightly assymetric */
90 /* the last coefficient */
95 /* even order filter; still symmetric */
101 q=vorbis_fromdBlook(amp*
103 vorbis_invsq2explook(qexp+m)-
110 vorbis_fpu_restore(fpu);
116 #include "lookup.c" /* catch this in the build system; we #include for
117 compilers (like gcc) that can't inline across
120 static int MLOOP_1[64]={
121 0,10,11,11, 12,12,12,12, 13,13,13,13, 13,13,13,13,
122 14,14,14,14, 14,14,14,14, 14,14,14,14, 14,14,14,14,
123 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15,
124 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15,
127 static int MLOOP_2[64]={
128 0,4,5,5, 6,6,6,6, 7,7,7,7, 7,7,7,7,
129 8,8,8,8, 8,8,8,8, 8,8,8,8, 8,8,8,8,
130 9,9,9,9, 9,9,9,9, 9,9,9,9, 9,9,9,9,
131 9,9,9,9, 9,9,9,9, 9,9,9,9, 9,9,9,9,
134 static int MLOOP_3[8]={0,1,2,2,3,3,3,3};
137 /* side effect: changes *lsp to cosines of lsp */
138 void vorbis_lsp_to_curve(float *curve,int *map,int n,int ln,float *lsp,int m,
139 float amp,float ampoffset){
143 /* set up for using all int later */
145 int ampoffseti=rint(ampoffset*4096.f);
146 int ampi=rint(amp*16.f);
147 long *ilsp=alloca(m*sizeof(*ilsp));
148 for(i=0;i<m;i++)ilsp[i]=vorbis_coslook_i(lsp[i]/M_PI*65536.f+.5f);
153 unsigned long pi=46341; /* 2**-.5 in 0.16 */
154 unsigned long qi=46341;
156 long wi=vorbis_coslook_i(k*65536/ln);
158 qi*=labs(ilsp[0]-wi);
159 pi*=labs(ilsp[1]-wi);
162 if(!(shift=MLOOP_1[(pi|qi)>>25]))
163 if(!(shift=MLOOP_2[(pi|qi)>>19]))
164 shift=MLOOP_3[(pi|qi)>>16];
165 qi=(qi>>shift)*labs(ilsp[j-1]-wi);
166 pi=(pi>>shift)*labs(ilsp[j]-wi);
169 if(!(shift=MLOOP_1[(pi|qi)>>25]))
170 if(!(shift=MLOOP_2[(pi|qi)>>19]))
171 shift=MLOOP_3[(pi|qi)>>16];
173 /* pi,qi normalized collectively, both tracked using qexp */
176 /* odd order filter; slightly assymetric */
177 /* the last coefficient */
178 qi=(qi>>shift)*labs(ilsp[j-1]-wi);
182 if(!(shift=MLOOP_1[(pi|qi)>>25]))
183 if(!(shift=MLOOP_2[(pi|qi)>>19]))
184 shift=MLOOP_3[(pi|qi)>>16];
188 qexp+=shift-14*((m+1)>>1);
194 pi*=(1<<14)-((wi*wi)>>14);
198 /* even order filter; still symmetric */
200 /* p*=p(1-w), q*=q(1+w), let normalization drift because it isn't
201 worth tracking step by step */
218 /* we've let the normalization drift because it wasn't important;
219 however, for the lookup, things must be normalized again. We
220 need at most one right shift or a number of left shifts */
222 if(qi&0xffff0000){ /* checks for 1.xxxxxxxxxxxxxxxx */
225 while(qi && !(qi&0x8000)){ /* checks for 0.0xxxxxxxxxxxxxxx or less*/
229 amp=vorbis_fromdBlook_i(ampi* /* n.4 */
230 vorbis_invsqlook_i(qi,qexp)-
232 ampoffseti); /* 8.12[0] */
235 while(map[++i]==k)curve[i]*=amp;
241 /* old, nonoptimized but simple version for any poor sap who needs to
242 figure out what the hell this code does, or wants the other
243 fraction of a dB precision */
245 /* side effect: changes *lsp to cosines of lsp */
246 void vorbis_lsp_to_curve(float *curve,int *map,int n,int ln,float *lsp,int m,
247 float amp,float ampoffset){
250 for(i=0;i<m;i++)lsp[i]=2.f*cos(lsp[i]);
257 float w=2.f*cos(wdel*k);
263 /* odd order filter; slightly assymetric */
264 /* the last coefficient */
269 /* even order filter; still symmetric */
274 q=fromdB(amp/sqrt(p+q)-ampoffset);
277 while(map[++i]==k)curve[i]*=q;
284 static void cheby(float *g, int ord) {
288 for(i=2; i<= ord; i++) {
289 for(j=ord; j >= i; j--) {
296 static int comp(const void *a,const void *b){
297 return (*(float *)a<*(float *)b)-(*(float *)a>*(float *)b);
300 /* Newton-Raphson-Maehly actually functioned as a decent root finder,
301 but there are root sets for which it gets into limit cycles
302 (exacerbated by zero suppression) and fails. We can't afford to
303 fail, even if the failure is 1 in 100,000,000, so we now use
304 Laguerre and later polish with Newton-Raphson (which can then
307 #define EPSILON 10e-7
308 static int Laguerre_With_Deflation(float *a,int ord,float *r){
310 double lastdelta=0.f;
311 double *defl=alloca(sizeof(*defl)*(ord+1));
312 for(i=0;i<=ord;i++)defl[i]=a[i];
315 double new=0.f,delta;
319 double p=defl[m],pp=0.f,ppp=0.f,denom;
321 /* eval the polynomial and its first two derivatives */
325 p = new*p + defl[i-1];
328 /* Laguerre's method */
329 denom=(m-1) * ((m-1)*pp*pp - m*p*ppp);
331 return(-1); /* complex root! The LPC generator handed us a bad filter */
334 denom = pp + sqrt(denom);
335 if(denom<EPSILON)denom=EPSILON;
337 denom = pp - sqrt(denom);
338 if(denom>-(EPSILON))denom=-(EPSILON);
344 if(delta<0.f)delta*=-1;
346 if(fabs(delta/new)<10e-12)break;
352 /* forward deflation */
355 defl[i-1]+=new*defl[i];
363 /* for spit-and-polish only */
364 static int Newton_Raphson(float *a,int ord,float *r){
367 double *root=alloca(ord*sizeof(*root));
369 for(i=0; i<ord;i++) root[i] = r[i];
374 for(i=0; i<ord; i++) { /* Update each point. */
376 double rooti=root[i];
378 for(k=ord-1; k>= 0; k--) {
381 p = p * rooti + a[k];
389 if(count>40)return(-1);
394 /* Replaced the original bubble sort with a real sort. With your
395 help, we can eliminate the bubble sort in our lifetime. --Monty */
397 for(i=0; i<ord;i++) r[i] = root[i];
402 /* Convert lpc coefficients to lsp coefficients */
403 int vorbis_lpc_to_lsp(float *lpc,float *lsp,int m){
405 int g1_order,g2_order;
406 float *g1=alloca(sizeof(*g1)*(order2+1));
407 float *g2=alloca(sizeof(*g2)*(order2+1));
408 float *g1r=alloca(sizeof(*g1r)*(order2+1));
409 float *g2r=alloca(sizeof(*g2r)*(order2+1));
412 /* even and odd are slightly different base cases */
416 /* Compute the lengths of the x polynomials. */
417 /* Compute the first half of K & R F1 & F2 polynomials. */
418 /* Compute half of the symmetric and antisymmetric polynomials. */
419 /* Remove the roots at +1 and -1. */
422 for(i=1;i<=g1_order;i++) g1[g1_order-i] = lpc[i-1]+lpc[m-i];
424 for(i=1;i<=g2_order;i++) g2[g2_order-i] = lpc[i-1]-lpc[m-i];
426 if(g1_order>g2_order){
427 for(i=2; i<=g2_order;i++) g2[g2_order-i] += g2[g2_order-i+2];
429 for(i=1; i<=g1_order;i++) g1[g1_order-i] -= g1[g1_order-i+1];
430 for(i=1; i<=g2_order;i++) g2[g2_order-i] += g2[g2_order-i+1];
433 /* Convert into polynomials in cos(alpha) */
437 /* Find the roots of the 2 even polynomials.*/
438 if(Laguerre_With_Deflation(g1,g1_order,g1r) ||
439 Laguerre_With_Deflation(g2,g2_order,g2r))
442 Newton_Raphson(g1,g1_order,g1r); /* if it fails, it leaves g1r alone */
443 Newton_Raphson(g2,g2_order,g2r); /* if it fails, it leaves g2r alone */
445 qsort(g1r,g1_order,sizeof(*g1r),comp);
446 qsort(g2r,g2_order,sizeof(*g2r),comp);
448 for(i=0;i<g1_order;i++)
449 lsp[i*2] = acos(g1r[i]);
451 for(i=0;i<g2_order;i++)
452 lsp[i*2+1] = acos(g2r[i]);