libs/Music/LPC.cpp

   1 #include "LPC.h"
   2
   3 // http://kbs.cs.tu-berlin.de/~jutta/gsm/lpc.html
   4
   5 #define P_MAX   8       /* order p of LPC analysis, typically 8..14 */
   6
   7 /* Invented by N. Levinson in 1947, modified by J. Durbin in 1959.
   8  */
   9 double                      /* returns minimum mean square error    */
  10                         levinson_durbin(
  11                         double const * ac,  /*  in: [0...p] autocorrelation values  */
  12         double       * ref, /* out: [0...p-1] reflection coef's     */
  13         double       * lpc) /*      [0...p-1] LPC coefficients      */
  14 {
  15         int i, j;  double r, error = ac[0];
  16
  17         if (ac[0] == 0) {
  18                 for (i = 0; i < P_MAX; i++) ref[i] = 0; return 0; }
  19
  20                 for (i = 0; i < P_MAX; i++) {
  21
  22                 /* Sum up this iteration's reflection coefficient.
  23                 */
  24                         r = -ac[i + 1];
  25                         for (j = 0; j < i; j++) r -= lpc[j] * ac[i - j];
  26                         ref[i] = r /= error;
  27
  28                 /*  Update LPC coefficients and total error.
  29                 */
  30                         lpc[i] = r;
  31                         for (j = 0; j < i/2; j++) {
  32                                 double tmp  = lpc[j];
  33                                 lpc[j]     += r * lpc[i-1-j];
  34                                 lpc[i-1-j] += r * tmp;
  35                         }
  36                         if (i % 2) lpc[j] += lpc[j] * r;
  37
  38                         error *= 1.0 - r * r;
  39                 }
  40                 return error;
  41 }
  42
  43 /* I. Sch"ur's recursion from 1917 is related to the Levinson-Durbin
  44  * method, but faster on parallel architectures; where Levinson-Durbin
  45  * would take time proportional to p * log(p), Sch"ur only requires
  46  * time proportional to p.  The GSM coder uses an integer version of
  47  * the Sch"ur recursion.
  48  */
  49 double                     /* returns the minimum mean square error */
  50                         schur(  double const * ac,      /*  in: [0...p] autocorrelation c's */
  51                                         double       * ref)     /* out: [0...p-1] reflection c's    */
  52 {
  53         int i, m;  double r, error = ac[0], G[2][P_MAX];
  54
  55         if (ac[0] == 0.0) {
  56                 for (i = 0; i < P_MAX; i++) ref[i] = 0;
  57                 return 0;
  58         }
  59
  60         /*  Initialize rows of the generator matrix G to ac[1...P]
  61         */
  62         for (i = 0; i < P_MAX; i++) G[0][i] = G[1][i] = ac[i+1];
  63
  64         for (i = 0;;) {
  65
  66                 /*  Calculate this iteration's reflection
  67                 *  coefficient and error.
  68                 */
  69                 ref[i] = r = -G[1][0] / error;
  70                 error += G[1][0] * r;
  71
  72                 if (++i >= P_MAX) return error;
  73
  74                 /*  Update the generator matrix.
  75                 *
  76                 *  Unlike the Levinson-Durbin summing of reflection
  77                 *  coefficients, this loop could be distributed to
  78                 *  p processors who each take only constant time.
  79                 */
  80                 for (m = 0; m < P_MAX - i; m++) {
  81                         G[1][m] = G[1][m+1] + r * G[0][m];
  82                         G[0][m] = G[1][m+1] * r + G[0][m];
  83                 }
  84         }
  85 }