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