18e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/********************************************************************
28e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels *                                                                  *
38e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
48e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
58e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
68e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
78e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels *                                                                  *
88e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2009             *
98e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels * by the Xiph.Org Foundation http://www.xiph.org/                  *
108e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels *                                                                  *
118e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels ********************************************************************
128e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
138e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  function: LSP (also called LSF) conversion routines
148e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  last mod: $Id: lsp.c 16227 2009-07-08 06:58:46Z xiphmont $
158e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
168e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  The LSP generation code is taken (with minimal modification and a
178e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  few bugfixes) from "On the Computation of the LSP Frequencies" by
188e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  Joseph Rothweiler (see http://www.rothweiler.us for contact info).
198e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  The paper is available at:
208e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
218e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  http://www.myown1.com/joe/lsf
228e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
238e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels ********************************************************************/
248e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
258e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/* Note that the lpc-lsp conversion finds the roots of polynomial with
268e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels   an iterative root polisher (CACM algorithm 283).  It *is* possible
278e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels   to confuse this algorithm into not converging; that should only
288e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels   happen with absurdly closely spaced roots (very sharp peaks in the
298e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels   LPC f response) which in turn should be impossible in our use of
308e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels   the code.  If this *does* happen anyway, it's a bug in the floor
318e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels   finder; find the cause of the confusion (probably a single bin
328e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels   spike or accidental near-float-limit resolution problems) and
338e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels   correct it. */
348e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
358e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#include <math.h>
368e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#include <string.h>
378e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#include <stdlib.h>
388e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#include "lsp.h"
398e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#include "os.h"
408e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#include "misc.h"
418e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#include "lookup.h"
428e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#include "scales.h"
438e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
448e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/* three possible LSP to f curve functions; the exact computation
458e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels   (float), a lookup based float implementation, and an integer
468e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels   implementation.  The float lookup is likely the optimal choice on
478e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels   any machine with an FPU.  The integer implementation is *not* fixed
488e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels   point (due to the need for a large dynamic range and thus a
498e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels   seperately tracked exponent) and thus much more complex than the
508e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels   relatively simple float implementations. It's mostly for future
518e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels   work on a fully fixed point implementation for processors like the
528e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels   ARM family. */
538e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
548e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/* define either of these (preferably FLOAT_LOOKUP) to have faster
558e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels   but less precise implementation. */
568e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#undef FLOAT_LOOKUP
578e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#undef INT_LOOKUP
588e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
598e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#ifdef FLOAT_LOOKUP
608e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#include "lookup.c" /* catch this in the build system; we #include for
618e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels                       compilers (like gcc) that can't inline across
628e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels                       modules */
638e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
648e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/* side effect: changes *lsp to cosines of lsp */
658e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsvoid vorbis_lsp_to_curve(float *curve,int *map,int n,int ln,float *lsp,int m,
668e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels                            float amp,float ampoffset){
678e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  int i;
688e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  float wdel=M_PI/ln;
698e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  vorbis_fpu_control fpu;
708e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
718e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  vorbis_fpu_setround(&fpu);
728e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  for(i=0;i<m;i++)lsp[i]=vorbis_coslook(lsp[i]);
738e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
748e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  i=0;
758e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  while(i<n){
768e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    int k=map[i];
778e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    int qexp;
788e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    float p=.7071067812f;
798e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    float q=.7071067812f;
808e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    float w=vorbis_coslook(wdel*k);
818e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    float *ftmp=lsp;
828e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    int c=m>>1;
838e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
848e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    do{
858e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      q*=ftmp[0]-w;
868e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      p*=ftmp[1]-w;
878e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      ftmp+=2;
888e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    }while(--c);
898e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
908e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    if(m&1){
918e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      /* odd order filter; slightly assymetric */
928e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      /* the last coefficient */
938e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      q*=ftmp[0]-w;
948e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      q*=q;
958e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      p*=p*(1.f-w*w);
968e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    }else{
978e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      /* even order filter; still symmetric */
988e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      q*=q*(1.f+w);
998e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      p*=p*(1.f-w);
1008e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    }
1018e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
1028e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    q=frexp(p+q,&qexp);
1038e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    q=vorbis_fromdBlook(amp*
1048e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels                        vorbis_invsqlook(q)*
1058e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels                        vorbis_invsq2explook(qexp+m)-
1068e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels                        ampoffset);
1078e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
1088e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    do{
1098e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      curve[i++]*=q;
1108e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    }while(map[i]==k);
1118e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  }
1128e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  vorbis_fpu_restore(fpu);
1138e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels}
1148e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
1158e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#else
1168e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
1178e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#ifdef INT_LOOKUP
1188e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#include "lookup.c" /* catch this in the build system; we #include for
1198e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels                       compilers (like gcc) that can't inline across
1208e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels                       modules */
1218e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
1228e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsstatic const int MLOOP_1[64]={
1238e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels   0,10,11,11, 12,12,12,12, 13,13,13,13, 13,13,13,13,
1248e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  14,14,14,14, 14,14,14,14, 14,14,14,14, 14,14,14,14,
1258e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15,
1268e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15,
1278e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels};
1288e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
1298e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsstatic const int MLOOP_2[64]={
1308e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  0,4,5,5, 6,6,6,6, 7,7,7,7, 7,7,7,7,
1318e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  8,8,8,8, 8,8,8,8, 8,8,8,8, 8,8,8,8,
1328e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  9,9,9,9, 9,9,9,9, 9,9,9,9, 9,9,9,9,
1338e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  9,9,9,9, 9,9,9,9, 9,9,9,9, 9,9,9,9,
1348e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels};
1358e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
1368e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsstatic const int MLOOP_3[8]={0,1,2,2,3,3,3,3};
1378e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
1388e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
1398e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/* side effect: changes *lsp to cosines of lsp */
1408e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsvoid vorbis_lsp_to_curve(float *curve,int *map,int n,int ln,float *lsp,int m,
1418e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels                            float amp,float ampoffset){
1428e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
1438e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  /* 0 <= m < 256 */
1448e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
1458e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  /* set up for using all int later */
1468e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  int i;
1478e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  int ampoffseti=rint(ampoffset*4096.f);
1488e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  int ampi=rint(amp*16.f);
1498e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  long *ilsp=alloca(m*sizeof(*ilsp));
1508e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  for(i=0;i<m;i++)ilsp[i]=vorbis_coslook_i(lsp[i]/M_PI*65536.f+.5f);
1518e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
1528e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  i=0;
1538e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  while(i<n){
1548e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    int j,k=map[i];
1558e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    unsigned long pi=46341; /* 2**-.5 in 0.16 */
1568e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    unsigned long qi=46341;
1578e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    int qexp=0,shift;
1588e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    long wi=vorbis_coslook_i(k*65536/ln);
1598e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
1608e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    qi*=labs(ilsp[0]-wi);
1618e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    pi*=labs(ilsp[1]-wi);
1628e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
1638e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    for(j=3;j<m;j+=2){
1648e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      if(!(shift=MLOOP_1[(pi|qi)>>25]))
1658e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels        if(!(shift=MLOOP_2[(pi|qi)>>19]))
1668e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels          shift=MLOOP_3[(pi|qi)>>16];
1678e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      qi=(qi>>shift)*labs(ilsp[j-1]-wi);
1688e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      pi=(pi>>shift)*labs(ilsp[j]-wi);
1698e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      qexp+=shift;
1708e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    }
1718e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    if(!(shift=MLOOP_1[(pi|qi)>>25]))
1728e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      if(!(shift=MLOOP_2[(pi|qi)>>19]))
1738e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels        shift=MLOOP_3[(pi|qi)>>16];
1748e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
1758e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    /* pi,qi normalized collectively, both tracked using qexp */
1768e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
1778e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    if(m&1){
1788e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      /* odd order filter; slightly assymetric */
1798e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      /* the last coefficient */
1808e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      qi=(qi>>shift)*labs(ilsp[j-1]-wi);
1818e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      pi=(pi>>shift)<<14;
1828e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      qexp+=shift;
1838e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
1848e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      if(!(shift=MLOOP_1[(pi|qi)>>25]))
1858e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels        if(!(shift=MLOOP_2[(pi|qi)>>19]))
1868e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels          shift=MLOOP_3[(pi|qi)>>16];
1878e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
1888e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      pi>>=shift;
1898e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      qi>>=shift;
1908e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      qexp+=shift-14*((m+1)>>1);
1918e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
1928e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      pi=((pi*pi)>>16);
1938e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      qi=((qi*qi)>>16);
1948e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      qexp=qexp*2+m;
1958e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
1968e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      pi*=(1<<14)-((wi*wi)>>14);
1978e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      qi+=pi>>14;
1988e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
1998e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    }else{
2008e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      /* even order filter; still symmetric */
2018e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
2028e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      /* p*=p(1-w), q*=q(1+w), let normalization drift because it isn't
2038e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels         worth tracking step by step */
2048e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
2058e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      pi>>=shift;
2068e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      qi>>=shift;
2078e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      qexp+=shift-7*m;
2088e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
2098e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      pi=((pi*pi)>>16);
2108e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      qi=((qi*qi)>>16);
2118e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      qexp=qexp*2+m;
2128e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
2138e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      pi*=(1<<14)-wi;
2148e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      qi*=(1<<14)+wi;
2158e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      qi=(qi+pi)>>14;
2168e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
2178e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    }
2188e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
2198e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
2208e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    /* we've let the normalization drift because it wasn't important;
2218e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels       however, for the lookup, things must be normalized again.  We
2228e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels       need at most one right shift or a number of left shifts */
2238e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
2248e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    if(qi&0xffff0000){ /* checks for 1.xxxxxxxxxxxxxxxx */
2258e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      qi>>=1; qexp++;
2268e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    }else
2278e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      while(qi && !(qi&0x8000)){ /* checks for 0.0xxxxxxxxxxxxxxx or less*/
2288e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels        qi<<=1; qexp--;
2298e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      }
2308e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
2318e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    amp=vorbis_fromdBlook_i(ampi*                     /*  n.4         */
2328e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels                            vorbis_invsqlook_i(qi,qexp)-
2338e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels                                                      /*  m.8, m+n<=8 */
2348e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels                            ampoffseti);              /*  8.12[0]     */
2358e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
2368e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    curve[i]*=amp;
2378e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    while(map[++i]==k)curve[i]*=amp;
2388e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  }
2398e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels}
2408e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
2418e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#else
2428e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
2438e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/* old, nonoptimized but simple version for any poor sap who needs to
2448e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels   figure out what the hell this code does, or wants the other
2458e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels   fraction of a dB precision */
2468e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
2478e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/* side effect: changes *lsp to cosines of lsp */
2488e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsvoid vorbis_lsp_to_curve(float *curve,int *map,int n,int ln,float *lsp,int m,
2498e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels                            float amp,float ampoffset){
2508e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  int i;
2518e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  float wdel=M_PI/ln;
2528e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  for(i=0;i<m;i++)lsp[i]=2.f*cos(lsp[i]);
2538e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
2548e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  i=0;
2558e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  while(i<n){
2568e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    int j,k=map[i];
2578e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    float p=.5f;
2588e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    float q=.5f;
2598e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    float w=2.f*cos(wdel*k);
2608e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    for(j=1;j<m;j+=2){
2618e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      q *= w-lsp[j-1];
2628e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      p *= w-lsp[j];
2638e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    }
2648e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    if(j==m){
2658e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      /* odd order filter; slightly assymetric */
2668e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      /* the last coefficient */
2678e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      q*=w-lsp[j-1];
2688e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      p*=p*(4.f-w*w);
2698e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      q*=q;
2708e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    }else{
2718e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      /* even order filter; still symmetric */
2728e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      p*=p*(2.f-w);
2738e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      q*=q*(2.f+w);
2748e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    }
2758e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
2768e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    q=fromdB(amp/sqrt(p+q)-ampoffset);
2778e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
2788e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    curve[i]*=q;
2798e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    while(map[++i]==k)curve[i]*=q;
2808e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  }
2818e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels}
2828e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
2838e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#endif
2848e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#endif
2858e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
2868e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsstatic void cheby(float *g, int ord) {
2878e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  int i, j;
2888e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
2898e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  g[0] *= .5f;
2908e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  for(i=2; i<= ord; i++) {
2918e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    for(j=ord; j >= i; j--) {
2928e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      g[j-2] -= g[j];
2938e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      g[j] += g[j];
2948e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    }
2958e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  }
2968e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels}
2978e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
2988e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsstatic int comp(const void *a,const void *b){
2998e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  return (*(float *)a<*(float *)b)-(*(float *)a>*(float *)b);
3008e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels}
3018e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
3028e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/* Newton-Raphson-Maehly actually functioned as a decent root finder,
3038e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels   but there are root sets for which it gets into limit cycles
3048e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels   (exacerbated by zero suppression) and fails.  We can't afford to
3058e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels   fail, even if the failure is 1 in 100,000,000, so we now use
3068e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels   Laguerre and later polish with Newton-Raphson (which can then
3078e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels   afford to fail) */
3088e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
3098e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels#define EPSILON 10e-7
3108e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsstatic int Laguerre_With_Deflation(float *a,int ord,float *r){
3118e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  int i,m;
3128e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  double lastdelta=0.f;
3138e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  double *defl=alloca(sizeof(*defl)*(ord+1));
3148e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  for(i=0;i<=ord;i++)defl[i]=a[i];
3158e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
3168e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  for(m=ord;m>0;m--){
3178e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    double new=0.f,delta;
3188e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
3198e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    /* iterate a root */
3208e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    while(1){
3218e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      double p=defl[m],pp=0.f,ppp=0.f,denom;
3228e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
3238e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      /* eval the polynomial and its first two derivatives */
3248e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      for(i=m;i>0;i--){
3258e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels        ppp = new*ppp + pp;
3268e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels        pp  = new*pp  + p;
3278e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels        p   = new*p   + defl[i-1];
3288e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      }
3298e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
3308e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      /* Laguerre's method */
3318e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      denom=(m-1) * ((m-1)*pp*pp - m*p*ppp);
3328e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      if(denom<0)
3338e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels        return(-1);  /* complex root!  The LPC generator handed us a bad filter */
3348e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
3358e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      if(pp>0){
3368e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels        denom = pp + sqrt(denom);
3378e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels        if(denom<EPSILON)denom=EPSILON;
3388e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      }else{
3398e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels        denom = pp - sqrt(denom);
3408e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels        if(denom>-(EPSILON))denom=-(EPSILON);
3418e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      }
3428e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
3438e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      delta  = m*p/denom;
3448e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      new   -= delta;
3458e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
3468e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      if(delta<0.f)delta*=-1;
3478e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
3488e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      if(fabs(delta/new)<10e-12)break;
3498e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      lastdelta=delta;
3508e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    }
3518e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
3528e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    r[m-1]=new;
3538e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
3548e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    /* forward deflation */
3558e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
3568e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    for(i=m;i>0;i--)
3578e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      defl[i-1]+=new*defl[i];
3588e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    defl++;
3598e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
3608e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  }
3618e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  return(0);
3628e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels}
3638e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
3648e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
3658e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/* for spit-and-polish only */
3668e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsstatic int Newton_Raphson(float *a,int ord,float *r){
3678e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  int i, k, count=0;
3688e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  double error=1.f;
3698e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  double *root=alloca(ord*sizeof(*root));
3708e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
3718e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  for(i=0; i<ord;i++) root[i] = r[i];
3728e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
3738e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  while(error>1e-20){
3748e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    error=0;
3758e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
3768e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    for(i=0; i<ord; i++) { /* Update each point. */
3778e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      double pp=0.,delta;
3788e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      double rooti=root[i];
3798e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      double p=a[ord];
3808e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      for(k=ord-1; k>= 0; k--) {
3818e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
3828e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels        pp= pp* rooti + p;
3838e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels        p = p * rooti + a[k];
3848e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      }
3858e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
3868e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      delta = p/pp;
3878e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      root[i] -= delta;
3888e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels      error+= delta*delta;
3898e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    }
3908e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
3918e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    if(count>40)return(-1);
3928e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
3938e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    count++;
3948e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  }
3958e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
3968e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  /* Replaced the original bubble sort with a real sort.  With your
3978e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels     help, we can eliminate the bubble sort in our lifetime. --Monty */
3988e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
3998e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  for(i=0; i<ord;i++) r[i] = root[i];
4008e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  return(0);
4018e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels}
4028e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
4038e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
4048e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels/* Convert lpc coefficients to lsp coefficients */
4058e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckelsint vorbis_lpc_to_lsp(float *lpc,float *lsp,int m){
4068e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  int order2=(m+1)>>1;
4078e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  int g1_order,g2_order;
4088e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  float *g1=alloca(sizeof(*g1)*(order2+1));
4098e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  float *g2=alloca(sizeof(*g2)*(order2+1));
4108e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  float *g1r=alloca(sizeof(*g1r)*(order2+1));
4118e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  float *g2r=alloca(sizeof(*g2r)*(order2+1));
4128e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  int i;
4138e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
4148e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  /* even and odd are slightly different base cases */
4158e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  g1_order=(m+1)>>1;
4168e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  g2_order=(m)  >>1;
4178e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
4188e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  /* Compute the lengths of the x polynomials. */
4198e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  /* Compute the first half of K & R F1 & F2 polynomials. */
4208e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  /* Compute half of the symmetric and antisymmetric polynomials. */
4218e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  /* Remove the roots at +1 and -1. */
4228e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
4238e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  g1[g1_order] = 1.f;
4248e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  for(i=1;i<=g1_order;i++) g1[g1_order-i] = lpc[i-1]+lpc[m-i];
4258e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  g2[g2_order] = 1.f;
4268e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  for(i=1;i<=g2_order;i++) g2[g2_order-i] = lpc[i-1]-lpc[m-i];
4278e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
4288e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  if(g1_order>g2_order){
4298e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    for(i=2; i<=g2_order;i++) g2[g2_order-i] += g2[g2_order-i+2];
4308e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  }else{
4318e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    for(i=1; i<=g1_order;i++) g1[g1_order-i] -= g1[g1_order-i+1];
4328e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    for(i=1; i<=g2_order;i++) g2[g2_order-i] += g2[g2_order-i+1];
4338e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  }
4348e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
4358e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  /* Convert into polynomials in cos(alpha) */
4368e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  cheby(g1,g1_order);
4378e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  cheby(g2,g2_order);
4388e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
4398e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  /* Find the roots of the 2 even polynomials.*/
4408e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  if(Laguerre_With_Deflation(g1,g1_order,g1r) ||
4418e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels     Laguerre_With_Deflation(g2,g2_order,g2r))
4428e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    return(-1);
4438e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
4448e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  Newton_Raphson(g1,g1_order,g1r); /* if it fails, it leaves g1r alone */
4458e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  Newton_Raphson(g2,g2_order,g2r); /* if it fails, it leaves g2r alone */
4468e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
4478e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  qsort(g1r,g1_order,sizeof(*g1r),comp);
4488e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  qsort(g2r,g2_order,sizeof(*g2r),comp);
4498e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
4508e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  for(i=0;i<g1_order;i++)
4518e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    lsp[i*2] = acos(g1r[i]);
4528e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels
4538e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  for(i=0;i<g2_order;i++)
4548e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels    lsp[i*2+1] = acos(g2r[i]);
4558e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels  return(0);
4568e01cdce135d5d816f92d7bb83f9a930aa1b45aeLucas Eckels}
457