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.8 2000/05/08 20:49:49 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 /* Note that the lpc-lsp conversion finds the roots of polynomial with
  26    an iterative root polisher (CACM algorithm 283).  It *is* possible
  27    to confuse this algorithm into not converging; that should only
  28    happen with absurdly closely spaced roots (very sharp peaks in the
  29    LPC f response) which in turn should be impossible in our use of
  30    the code.  If this *does* happen anyway, it's a bug in the floor
  31    finder; find the cause of the confusion (probably a single bin
  32    spike or accidental near-double-limit resolution problems) and
  33    correct it. */
  34
  35 #include <math.h>
  36 #include <string.h>
  37 #include <stdlib.h>
  38 #include "lsp.h"
  39 #include "os.h"
  40 #include "misc.h"
  41
  42 void vorbis_lsp_to_lpc(double *lsp,double *lpc,int m){
  43   int i,j,m2=m/2;
  44   double *O=alloca(sizeof(double)*m2);
  45   double *E=alloca(sizeof(double)*m2);
  46   double A;
  47   double *Ae=alloca(sizeof(double)*(m2+1));
  48   double *Ao=alloca(sizeof(double)*(m2+1));
  49   double B;
  50   double *Be=alloca(sizeof(double)*(m2));
  51   double *Bo=alloca(sizeof(double)*(m2));
  52   double temp;
  53
  54   /* even/odd roots setup */
  55   for(i=0;i<m2;i++){
  56     O[i]=-2.*cos(lsp[i*2]);
  57     E[i]=-2.*cos(lsp[i*2+1]);
  58   }
  59
  60   /* set up impulse response */
  61   for(j=0;j<m2;j++){
  62     Ae[j]=0.;
  63     Ao[j]=1.;
  64     Be[j]=0.;
  65     Bo[j]=1.;
  66   }
  67   Ao[j]=1.;
  68   Ae[j]=1.;
  69
  70   /* run impulse response */
  71   for(i=1;i<m+1;i++){
  72     A=B=0.;
  73     for(j=0;j<m2;j++){
  74       temp=O[j]*Ao[j]+Ae[j];
  75       Ae[j]=Ao[j];
  76       Ao[j]=A;
  77       A+=temp;
  78
  79       temp=E[j]*Bo[j]+Be[j];
  80       Be[j]=Bo[j];
  81       Bo[j]=B;
  82       B+=temp;
  83     }
  84     lpc[i-1]=(A+Ao[j]+B-Ae[j])/2;
  85     Ao[j]=A;
  86     Ae[j]=B;
  87   }
  88 }
  89
  90 static void cheby(double *g, int ord) {
  91   int i, j;
  92
  93   g[0] *= 0.5;
  94   for(i=2; i<= ord; i++) {
  95     for(j=ord; j >= i; j--) {
  96       g[j-2] -= g[j];
  97       g[j] += g[j];
  98     }
  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
 130       error += delta*delta;
 131     }
 132   }
 133
 134   /* Replaced the original bubble sort with a real sort.  With your
 135      help, we can eliminate the bubble sort in our lifetime. --Monty */
 136
 137   qsort(r,ord,sizeof(double),comp);
 138
 139 }
 140
 141 /* Convert lpc coefficients to lsp coefficients */
 142 void vorbis_lpc_to_lsp(double *lpc,double *lsp,int m){
 143   int order2=m/2;
 144   double *g1=alloca(sizeof(double)*(order2+1));
 145   double *g2=alloca(sizeof(double)*(order2+1));
 146   double *g1r=alloca(sizeof(double)*(order2+1));
 147   double *g2r=alloca(sizeof(double)*(order2+1));
 148   int i;
 149
 150   /* Compute the lengths of the x polynomials. */
 151   /* Compute the first half of K & R F1 & F2 polynomials. */
 152   /* Compute half of the symmetric and antisymmetric polynomials. */
 153   /* Remove the roots at +1 and -1. */
 154
 155   g1[order2] = 1.0;
 156   for(i=0;i<order2;i++) g1[order2-i-1] = lpc[i]+lpc[m-i-1];
 157   g2[order2] = 1.0;
 158   for(i=0;i<order2;i++) g2[order2-i-1] = lpc[i]-lpc[m-i-1];
 159
 160   for(i=0; i<order2;i++) g1[order2-i-1] -= g1[order2-i];
 161   for(i=0; i<order2;i++) g2[order2-i-1] += g2[order2-i];
 162
 163   /* Convert into polynomials in cos(alpha) */
 164   cheby(g1,order2);
 165   cheby(g2,order2);
 166
 167   /* Find the roots of the 2 even polynomials.*/
 168
 169   cacm283(g1,order2,g1r);
 170   cacm283(g2,order2,g2r);
 171
 172   for(i=0;i<m;i+=2){
 173     lsp[i] = acos(g1r[i/2]);
 174     lsp[i+1] = acos(g2r[i/2]);
 175   }
 176 }