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.9 2000/08/19 11:46:28 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_curve(double *curve,int n,double *lsp,int m,double amp,
  43                          double *w){
  44   int i,j,k;
  45   double *coslsp=alloca(m*sizeof(double));
  46   for(i=0;i<m;i++)coslsp[i]=2*cos(lsp[i]);
  47
  48   for(k=0;k<n;k++){
  49     double p=.70710678118654752440;
  50     double q=.70710678118654752440;
  51     for(j=0;j<m;){
  52       p*= *w-coslsp[j++];
  53       q*= *w-coslsp[j++];
  54     }
  55     curve[k]=amp/sqrt(p*p*(1.+ *w*.5)+q*q*(1.- *w*.5));
  56     w++;
  57   }
  58 }
  59
  60 static void cheby(double *g, int ord) {
  61   int i, j;
  62
  63   g[0] *= 0.5;
  64   for(i=2; i<= ord; i++) {
  65     for(j=ord; j >= i; j--) {
  66       g[j-2] -= g[j];
  67       g[j] += g[j];
  68     }
  69   }
  70 }
  71
  72 static int comp(const void *a,const void *b){
  73   if(*(double *)a<*(double *)b)
  74     return(1);
  75   else
  76     return(-1);
  77 }
  78
  79 /* CACM algorithm 283. */
  80 static void cacm283(double *a,int ord,double *r){
  81   int i, k;
  82   double val, p, delta, error;
  83   double rooti;
  84
  85   for(i=0; i<ord;i++) r[i] = 2.0 * (i+0.5) / ord - 1.0;
  86
  87   for(error=1 ; error > 1.e-12; ) {
  88     error = 0;
  89     for( i=0; i<ord; i++) {  /* Update each point. */
  90       rooti = r[i];
  91       val = a[ord];
  92       p = a[ord];
  93       for(k=ord-1; k>= 0; k--) {
  94         val = val * rooti + a[k];
  95         if (k != i) p *= rooti - r[k];
  96       }
  97       delta = val/p;
  98       r[i] -= delta;
  99
 100       error += delta*delta;
 101     }
 102   }
 103
 104   /* Replaced the original bubble sort with a real sort.  With your
 105      help, we can eliminate the bubble sort in our lifetime. --Monty */
 106
 107   qsort(r,ord,sizeof(double),comp);
 108
 109 }
 110
 111 /* Convert lpc coefficients to lsp coefficients */
 112 void vorbis_lpc_to_lsp(double *lpc,double *lsp,int m){
 113   int order2=m/2;
 114   double *g1=alloca(sizeof(double)*(order2+1));
 115   double *g2=alloca(sizeof(double)*(order2+1));
 116   double *g1r=alloca(sizeof(double)*(order2+1));
 117   double *g2r=alloca(sizeof(double)*(order2+1));
 118   int i;
 119
 120   /* Compute the lengths of the x polynomials. */
 121   /* Compute the first half of K & R F1 & F2 polynomials. */
 122   /* Compute half of the symmetric and antisymmetric polynomials. */
 123   /* Remove the roots at +1 and -1. */
 124
 125   g1[order2] = 1.0;
 126   for(i=0;i<order2;i++) g1[order2-i-1] = lpc[i]+lpc[m-i-1];
 127   g2[order2] = 1.0;
 128   for(i=0;i<order2;i++) g2[order2-i-1] = lpc[i]-lpc[m-i-1];
 129
 130   for(i=0; i<order2;i++) g1[order2-i-1] -= g1[order2-i];
 131   for(i=0; i<order2;i++) g2[order2-i-1] += g2[order2-i];
 132
 133   /* Convert into polynomials in cos(alpha) */
 134   cheby(g1,order2);
 135   cheby(g2,order2);
 136
 137   /* Find the roots of the 2 even polynomials.*/
 138
 139   cacm283(g1,order2,g1r);
 140   cacm283(g2,order2,g2r);
 141
 142   for(i=0;i<m;i+=2){
 143     lsp[i] = acos(g1r[i/2]);
 144     lsp[i+1] = acos(g2r[i/2]);
 145   }
 146 }