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: Aug 03 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
  31 void vorbis_lsp_to_lpc(double *lsp,double *lpc,int m){
  32   int i,j,m2=m/2;
  33   double O[m2],E[m2];
  34   double A,Ae[m2+1],Ao[m2+1];
  35   double B,Be[m2],Bo[m2];
  36   double temp;
  37
  38   /* even/odd roots setup */
  39   for(i=0;i<m2;i++){
  40     O[i]=-2.*cos(lsp[i*2]);
  41     E[i]=-2.*cos(lsp[i*2+1]);
  42   }
  43
  44   /* set up impulse response */
  45   for(j=0;j<m2;j++){
  46     Ae[j]=0.;
  47     Ao[j]=1.;
  48     Be[j]=0.;
  49     Bo[j]=1.;
  50   }
  51   Ao[j]=1.;
  52   Ae[j]=1.;
  53
  54   /* run impulse response */
  55   for(i=1;i<m+1;i++){
  56     A=B=0.;
  57     for(j=0;j<m2;j++){
  58       temp=O[j]*Ao[j]+Ae[j];
  59       Ae[j]=Ao[j];
  60       Ao[j]=A;
  61       A+=temp;
  62
  63       temp=E[j]*Bo[j]+Be[j];
  64       Be[j]=Bo[j];
  65       Bo[j]=B;
  66       B+=temp;
  67     }
  68     lpc[i-1]=(A+Ao[j]+B-Ae[j])/2;
  69     Ao[j]=A;
  70     Ae[j]=B;
  71   }
  72 }
  73
  74 static void kw(double *r,int n) {
  75   double s[n/2+1];
  76   double c[n+1];
  77   int i, j, k;
  78
  79   s[0] = 1.0;
  80   s[1] = -2.0;
  81   s[2] = 2.0;
  82   for(i=3;i<=n/2;i++) s[i] = s[i-2];
  83
  84   for(k=0;k<=n;k++) {
  85     c[k] = r[k];
  86     j = 1;
  87     for(i=k+2;i<=n;i+=2) {
  88       c[k] += s[j]*r[i];
  89       s[j] -= s[j-1];
  90       j++;
  91     }
  92   }
  93   for(k=0;k<=n;k++) r[k] = c[k];
  94 }
  95
  96
  97 static int comp(const void *a,const void *b){
  98   return(*(double *)a<*(double *)b);
  99 }
 100
 101 /* CACM algorithm 283. */
 102 static void cacm283(double *a,int ord,double *r){
 103   int i, k;
 104   double val, p, delta, error;
 105   double rooti;
 106
 107   for(i=0; i<ord;i++) r[i] = 2.0 * (i+0.5) / ord - 1.0;
 108
 109   for(error=1 ; error > 1.e-12; ) {
 110     error = 0;
 111     for( i=0; i<ord; i++) {  /* Update each point. */
 112       rooti = r[i];
 113       val = a[ord];
 114       p = a[ord];
 115       for(k=ord-1; k>= 0; k--) {
 116         val = val * rooti + a[k];
 117         if (k != i) p *= rooti - r[k];
 118       }
 119       delta = val/p;
 120       r[i] -= delta;
 121       error += delta*delta;
 122     }
 123   }
 124
 125   /* Replaced the original bubble sort with a real sort.  With your
 126      help, we can eliminate the bubble sort in our lifetime. --Monty */
 127
 128   qsort(r,ord,sizeof(double),comp);
 129
 130 }
 131
 132 /* Convert lpc coefficients to lsp coefficients */
 133 void vorbis_lpc_to_lsp(double *lpc,double *lsp,int m){
 134   int order2=m/2;
 135   double g1[order2+1], g2[order2+1];
 136   double g1r[order2+1], g2r[order2+1];
 137   int i;
 138
 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. */
 143
 144   g1[order2] = 1.0;
 145   for(i=0;i<order2;i++) g1[order2-i-1] = lpc[i]+lpc[m-i-1];
 146   g2[order2] = 1.0;
 147   for(i=0;i<order2;i++) g2[order2-i-1] = lpc[i]-lpc[m-i-1];
 148
 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];
 151
 152   /* Convert into polynomials in cos(alpha) */
 153   kw(g1,order2);
 154   kw(g2,order2);
 155
 156   /* Find the roots of the 2 even polynomials.*/
 157
 158   cacm283(g1,order2,g1r);
 159   cacm283(g2,order2,g2r);
 160
 161   for(i=0;i<m;i+=2){
 162     lsp[i] = acos(g1r[i/2]*.5);
 163     lsp[i+1] = acos(g2r[i/2]*.5);
 164   }
 165 }
 166