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