home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS - Coast to Coast / simteldosarchivecoasttocoast2.iso / c / lpc05b.zip / MATDURBN.C < prev    next >
C/C++ Source or Header  |  1992-05-22  |  3KB  |  141 lines
  1. /*
  2. *-----------------------------------------------------------------------------
  3. *    file:    matdurbn.c
  4. *    desc:    Levinson-Durbin algorithm
  5. *
  6. *        (specially for LPC, PARCOR calculation)
  7. *
  8. *    by:    ko shu pui, patrick
  9. *    date:
  10. *    revi:    21 may 92 v0.3
  11. *    ref:
  12. *
  13. *       [1] "Fundementals of Speech Signal Processing," Shuzo Saito,
  14. *       Kazuo Nakata, Academic Press, New York, 1985.
  15. *
  16. *-----------------------------------------------------------------------------
  17. */
  18.  
  19. #include <stdio.h>
  20.  
  21. #include "matrix.h"
  22.  
  23. /*
  24. *-----------------------------------------------------------------------------
  25. *    funct:    mat_durbin
  26. *    desct:    Levinson-Durbin algorithm
  27. *
  28. *        This function solve the linear eqns Ax = B:
  29. *
  30. *        |  v0   v1   v2  .. vn-1 | |  a1   |    |  v1   |
  31. *        |  v1   v0   v1  .. vn-2 | |  a2   |    |  v2   |
  32. *        |  v2   v1   v0  .. vn-3 | |  a3   |  = |  ..   |
  33. *        |  ...                   | |  ..   |    |  ..   |
  34. *        |  vn-1 vn-2 ..  .. v0   | |  an   |    |  vn   |
  35. *
  36. *        where A is a symmetric Toeplitz matrix and B
  37. *        in the above format (related to A)
  38. *
  39. *    given:    R = autocorrelated matrix (v0, v1, ... vn) (dim (n+1) x 1)
  40. *    retrn:    x (of Ax = B) - LPC coefficients
  41. *        P = PARCOR coefficients (dim n x 1 )
  42. *        E = residue power (dim (n+1) x 1 )
  43. *-----------------------------------------------------------------------------
  44. */
  45. MATRIX mat_durbin( R, P, E )
  46. MATRIX R, P, E;
  47. {
  48.     int    i, i1, j, ji, p, n;
  49.     MATRIX    W, A, X;
  50.  
  51.     p = MatRow(R) - 1;
  52.     W = mat_creat( p, 1, UNDEFINED );
  53.     A = mat_creat( p+1, p+1, UNDEFINED );
  54.     
  55.  
  56.     W[0][0] = R[1][0];
  57.     E[0][0] = R[0][0];
  58.  
  59.     for (i=1; i<=p; i++)
  60.         {
  61.         P[i-1][0] = W[i-1][0] / E[i-1][0];
  62.  
  63.         E[i][0] = E[i-1][0] * (1.0 - P[i-1][0] * P[i-1][0]);
  64.         if (E[i][0] <= 0.0)
  65.             E[i][0] = E[i-1][0];
  66.  
  67.         A[i-1][i-1] = -P[i-1][0];
  68.  
  69.         i1 = i-1;
  70.         if (i1 >= 1)
  71.             {
  72.             for (j=1; j<=i1; j++)
  73.             {
  74.             ji = i - j;
  75.             A[j-1][i-1] = A[j-1][i1-1] - P[i-1][0] * A[ji-1][i1-1];
  76.             }
  77.             }
  78.  
  79.         if (i != p)
  80.             {
  81.             W[i][0] = R[i+1][0];
  82.             for (j=1; j<=i; j++)
  83.                 W[i][0] += A[j-1][i-1] * R[i-j+1][0];
  84.             }
  85.         }
  86.  
  87.     X = mat_creat( p, 1, UNDEFINED );
  88.     for (i=0; i<p; i++)
  89.         {
  90.         X[i][0] = -A[i][p-1];
  91.         }
  92.  
  93.     mat_free( A );
  94.     mat_free( W );
  95.     return (X);
  96. }
  97.  
  98. /*
  99. *-----------------------------------------------------------------------------
  100. *    funct:    mat_lsolve_durbin
  101. *    desct:    Solve simultaneous linear eqns using
  102. *        Levinson-Durbin algorithm
  103. *
  104. *        This function solve the linear eqns Ax = B:
  105. *
  106. *        |  v0   v1   v2  .. vn-1 | |  a1   |    |  v1   |
  107. *        |  v1   v0   v1  .. vn-2 | |  a2   |    |  v2   |
  108. *        |  v2   v1   v0  .. vn-3 | |  a3   |  = |  ..   |
  109. *        |  ...                   | |  ..   |    |  ..   |
  110. *        |  vn-1 vn-2 ..  .. v0   | |  an   |    |  vn   |
  111. *
  112. *    domain:    where A is a symmetric Toeplitz matrix and B
  113. *        in the above format (related to A)
  114. *
  115. *    given:    A, B
  116. *    retrn:    x (of Ax = B)
  117. *
  118. *    WARNING: this version of Levinson-Durbin method may cause some lost
  119. *        of accuracy, it is only suitable for applications such as
  120. *        speech processing, etc where you bear this in mind
  121. *-----------------------------------------------------------------------------
  122. */
  123. MATRIX mat_lsolve_durbin( A, B, P, E )
  124. MATRIX A, B, P, E;
  125. {
  126.     MATRIX    R, X;
  127.     int    i, n;
  128.  
  129.     n = MatRow(A);
  130.     R = mat_creat(n+1, 1, UNDEFINED);
  131.     for (i=0; i<n; i++)
  132.         {
  133.         R[i][0] = A[i][0];
  134.         }
  135.     R[n][0] = B[n-1][0];
  136.  
  137.     X = mat_durbin( R, P, E );
  138.     mat_free( R );
  139.     return (X);
  140. }
  141.