lib/lsp.c

   1 /********************************************************************
   2  *                                                                  *
   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.                            *
   7  *                                                                  *
   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/                                             *
  11  *                                                                  *
  12  ********************************************************************
  13
  14   function: LSP (also called LSF) conversion routines
  15   author: Monty <monty@xiph.org>
  16   modifications by: Monty
  17   last modification date: Nov 16 1999
  18
  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:
  22
  23   http://www2.xtdl.com/~rothwlr/lsfpaper/lsfpage.html
  24
  25  ********************************************************************/
  26
  27 #include <math.h>
  28 #include <string.h>
  29 #include <stdlib.h>
  30 #include "os.h"
  31
  32 void vorbis_lsp_to_lpc(double *lsp,double *lpc,int m){
  33   int i,j,m2=m/2;
  34   double *O=alloca(sizeof(double)*m2);
  35   double *E=alloca(sizeof(double)*m2);
  36   double A;
  37   double *Ae=alloca(sizeof(double)*(m2+1));
  38   double *Ao=alloca(sizeof(double)*(m2+1));
  39   double B;
  40   double *Be=alloca(sizeof(double)*(m2));
  41   double *Bo=alloca(sizeof(double)*(m2));
  42   double temp;
  43
  44   /* even/odd roots setup */
  45   for(i=0;i<m2;i++){
  46     O[i]=-2.*cos(lsp[i*2]);
  47     E[i]=-2.*cos(lsp[i*2+1]);
  48   }
  49
  50   /* set up impulse response */
  51   for(j=0;j<m2;j++){
  52     Ae[j]=0.;
  53     Ao[j]=1.;
  54     Be[j]=0.;
  55     Bo[j]=1.;
  56   }
  57   Ao[j]=1.;
  58   Ae[j]=1.;
  59
  60   /* run impulse response */
  61   for(i=1;i<m+1;i++){
  62     A=B=0.;
  63     for(j=0;j<m2;j++){
  64       temp=O[j]*Ao[j]+Ae[j];
  65       Ae[j]=Ao[j];
  66       Ao[j]=A;
  67       A+=temp;
  68
  69       temp=E[j]*Bo[j]+Be[j];
  70       Be[j]=Bo[j];
  71       Bo[j]=B;
  72       B+=temp;
  73     }
  74     lpc[i-1]=(A+Ao[j]+B-Ae[j])/2;
  75     Ao[j]=A;
  76     Ae[j]=B;
  77   }
  78 }
  79
  80 static void kw(double *r,int n) {
  81   double *s=alloca(sizeof(double)*(n/2+1));
  82   double *c=alloca(sizeof(double)*(n+1));
  83   int i, j, k;
  84
  85   s[0] = 1.0;
  86   s[1] = -2.0;
  87   s[2] = 2.0;
  88   for(i=3;i<=n/2;i++) s[i] = s[i-2];
  89
  90   for(k=0;k<=n;k++) {
  91     c[k] = r[k];
  92     j = 1;
  93     for(i=k+2;i<=n;i+=2) {
  94       c[k] += s[j]*r[i];
  95       s[j] -= s[j-1];
  96       j++;
  97     }
  98   }
  99   for(k=0;k<=n;k++) r[k] = c[k];
 100 }
 101
 102
 103 static int comp(const void *a,const void *b){
 104   return(*(double *)a<*(double *)b);
 105 }
 106
 107 /* CACM algorithm 283. */
 108 static void cacm283(double *a,int ord,double *r){
 109   int i, k;
 110   double val, p, delta, error;
 111   double rooti;
 112
 113   for(i=0; i<ord;i++) r[i] = 2.0 * (i+0.5) / ord - 1.0;
 114
 115   for(error=1 ; error > 1.e-12; ) {
 116     error = 0;
 117     for( i=0; i<ord; i++) {  /* Update each point. */
 118       rooti = r[i];
 119       val = a[ord];
 120       p = a[ord];
 121       for(k=ord-1; k>= 0; k--) {
 122         val = val * rooti + a[k];
 123         if (k != i) p *= rooti - r[k];
 124       }
 125       delta = val/p;
 126       r[i] -= delta;
 127       error += delta*delta;
 128     }
 129   }
 130
 131   /* Replaced the original bubble sort with a real sort.  With your
 132      help, we can eliminate the bubble sort in our lifetime. --Monty */
 133
 134   qsort(r,ord,sizeof(double),comp);
 135
 136 }
 137
 138 /* Convert lpc coefficients to lsp coefficients */
 139 void vorbis_lpc_to_lsp(double *lpc,double *lsp,int m){
 140   int order2=m/2;
 141   double *g1=alloca(sizeof(double)*(order2+1));
 142   double *g2=alloca(sizeof(double)*(order2+1));
 143   double *g1r=alloca(sizeof(double)*(order2+1));
 144   double *g2r=alloca(sizeof(double)*(order2+1));
 145   int i;
 146
 147   /* Compute the lengths of the x polynomials. */
 148   /* Compute the first half of K & R F1 & F2 polynomials. */
 149   /* Compute half of the symmetric and antisymmetric polynomials. */
 150   /* Remove the roots at +1 and -1. */
 151
 152   g1[order2] = 1.0;
 153   for(i=0;i<order2;i++) g1[order2-i-1] = lpc[i]+lpc[m-i-1];
 154   g2[order2] = 1.0;
 155   for(i=0;i<order2;i++) g2[order2-i-1] = lpc[i]-lpc[m-i-1];
 156
 157   for(i=0; i<order2;i++) g1[order2-i-1] -= g1[order2-i];
 158   for(i=0; i<order2;i++) g2[order2-i-1] += g2[order2-i];
 159
 160   /* Convert into polynomials in cos(alpha) */
 161   kw(g1,order2);
 162   kw(g2,order2);
 163
 164   /* Find the roots of the 2 even polynomials.*/
 165
 166   cacm283(g1,order2,g1r);
 167   cacm283(g2,order2,g2r);
 168
 169   for(i=0;i<m;i+=2){
 170     lsp[i] = acos(g1r[i/2]*.5);
 171     lsp[i+1] = acos(g2r[i/2]*.5);
 172   }
 173 }