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