1 /********************************************************************
3 * THIS FILE IS PART OF THE Ogg Vorbis SOFTWARE CODEC SOURCE CODE. *
4 * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
5 * THE GNU PUBLIC LICENSE 2, WHICH IS INCLUDED WITH THIS SOURCE. *
6 * PLEASE READ THESE TERMS DISTRIBUTING. *
8 * THE OggSQUISH SOURCE CODE IS (C) COPYRIGHT 1994-1999 *
9 * by 1999 Monty <monty@xiph.org> and The XIPHOPHORUS Company *
10 * http://www.xiph.org/ *
12 ********************************************************************
14 function: LSP (also called LSF) conversion routines
15 author: Monty <monty@xiph.org>
16 modifications by: Monty
17 last modification date: Aug 03 1999
19 The LSP generation code is taken (with minimal modification) from
20 "On the Computation of the LSP Frequencies" by Joseph Rothweiler
21 <rothwlr@altavista.net>, available at:
23 http://www2.xtdl.com/~rothwlr/lsfpaper/lsfpage.html
25 ********************************************************************/
31 void vorbis_lsp_to_lpc(double *lsp,double *lpc,int m){
34 double A,Ae[m2+1],Ao[m2+1];
35 double B,Be[m2],Bo[m2];
38 /* even/odd roots setup */
40 O[i]=-2.*cos(lsp[i*2]);
41 E[i]=-2.*cos(lsp[i*2+1]);
44 /* set up impulse response */
54 /* run impulse response */
58 temp=O[j]*Ao[j]+Ae[j];
63 temp=E[j]*Bo[j]+Be[j];
68 lpc[i-1]=(A+Ao[j]+B-Ae[j])/2;
74 static void kw(double *r,int n) {
82 for(i=3;i<=n/2;i++) s[i] = s[i-2];
87 for(i=k+2;i<=n;i+=2) {
93 for(k=0;k<=n;k++) r[k] = c[k];
97 static int comp(const void *a,const void *b){
98 return(*(double *)a<*(double *)b);
101 /* CACM algorithm 283. */
102 static void cacm283(double *a,int ord,double *r){
104 double val, p, delta, error;
107 for(i=0; i<ord;i++) r[i] = 2.0 * (i+0.5) / ord - 1.0;
109 for(error=1 ; error > 1.e-12; ) {
111 for( i=0; i<ord; i++) { /* Update each point. */
115 for(k=ord-1; k>= 0; k--) {
116 val = val * rooti + a[k];
117 if (k != i) p *= rooti - r[k];
121 error += delta*delta;
125 /* Replaced the original bubble sort with a real sort. With your
126 help, we can eliminate the bubble sort in our lifetime. --Monty */
128 qsort(r,ord,sizeof(double),comp);
132 /* Convert lpc coefficients to lsp coefficients */
133 void vorbis_lpc_to_lsp(double *lpc,double *lsp,int m){
135 double g1[order2+1], g2[order2+1];
136 double g1r[order2+1], g2r[order2+1];
139 /* Compute the lengths of the x polynomials. */
140 /* Compute the first half of K & R F1 & F2 polynomials. */
141 /* Compute half of the symmetric and antisymmetric polynomials. */
142 /* Remove the roots at +1 and -1. */
145 for(i=0;i<order2;i++) g1[order2-i-1] = lpc[i]+lpc[m-i-1];
147 for(i=0;i<order2;i++) g2[order2-i-1] = lpc[i]-lpc[m-i-1];
149 for(i=0; i<order2;i++) g1[order2-i-1] -= g1[order2-i];
150 for(i=0; i<order2;i++) g2[order2-i-1] += g2[order2-i];
152 /* Convert into polynomials in cos(alpha) */
156 /* Find the roots of the 2 even polynomials.*/
158 cacm283(g1,order2,g1r);
159 cacm283(g2,order2,g2r);
162 lsp[i] = acos(g1r[i/2]*.5);
163 lsp[i+1] = acos(g2r[i/2]*.5);