lpc10/invert.c

   1 /*
   2
   3  * Revision 1.1  1996/08/19  22:32:00  jaf
   4  * Initial revision
   5  *
   6
   7 */
   8
   9 /*  -- translated by f2c (version 19951025).
  10    You must link the resulting object file with the libraries:
  11         -lf2c -lm   (in that order)
  12 */
  13
  14 #include "f2c.h"
  15
  16 extern int invert_(integer *order, real *phi, real *psi, real *rc);
  17
  18 /* **************************************************************** */
  19
  20 /*      INVERT Version 45G */
  21
  22 /*
  23  * Revision 1.1  1996/08/19  22:32:00  jaf
  24  * Initial revision
  25  * */
  26 /* Revision 1.3  1996/03/18  20:52:47  jaf */
  27 /* Just added a few comments about which array indices of the arguments */
  28 /* are used, and mentioning that this subroutine has no local state. */
  29
  30 /* Revision 1.2  1996/03/13  16:51:32  jaf */
  31 /* Comments added explaining that none of the local variables of this */
  32 /* subroutine need to be saved from one invocation to the next. */
  33
  34 /* Eliminated a comment from the original, describing a local array X */
  35 /* that appeared nowhere in the code. */
  36
  37 /* Revision 1.1  1996/02/07 14:47:20  jaf */
  38 /* Initial revision */
  39
  40
  41 /* **************************************************************** */
  42
  43 /*  Invert a covariance matrix using Choleski decomposition method. */
  44
  45 /* Input: */
  46 /*  ORDER            - Analysis order */
  47 /*  PHI(ORDER,ORDER) - Covariance matrix */
  48 /*                    Indices (I,J) read, where ORDER .GE. I .GE. J .GE. 1.*/
  49 /*                     All other indices untouched. */
  50 /*  PSI(ORDER)       - Column vector to be predicted */
  51 /*                     Indices 1 through ORDER read. */
  52 /* Output: */
  53 /*  RC(ORDER)        - Pseudo reflection coefficients */
  54 /*                    Indices 1 through ORDER written, and then possibly read.
  55 */
  56 /* Internal: */
  57 /*  V(ORDER,ORDER)   - Temporary matrix */
  58 /*                     Same indices written as read from PHI. */
  59 /*                     Many indices may be read and written again after */
  60 /*                     initially being copied from PHI, but all indices */
  61 /*                     are written before being read. */
  62
  63 /*  NOTE: Temporary matrix V is not needed and may be replaced */
  64 /*    by PHI if the original PHI values do not need to be preserved. */
  65
  66 /* Subroutine */ int invert_(integer *order, real *phi, real *psi, real *rc)
  67 {
  68     /* System generated locals */
  69     integer phi_dim1, phi_offset, i__1, i__2, i__3;
  70     real r__1, r__2;
  71
  72     /* Local variables */
  73     real save;
  74     integer i__, j, k;
  75     real v[100] /* was [10][10] */;
  76
  77 /*       Arguments */
  78
  79 /*   LPC Configuration parameters: */
  80 /* Frame size, Prediction order, Pitch period */
  81 /*      Parameters/constants */
  82 /*       Local variables that need not be saved */
  83 /*  Decompose PHI into V * D * V' where V is a triangular matrix whose */
  84 /*   main diagonal elements are all 1, V' is the transpose of V, and */
  85 /*   D is a vector.  Here D(n) is stored in location V(n,n). */
  86     /* Parameter adjustments */
  87     --rc;
  88     --psi;
  89     phi_dim1 = *order;
  90     phi_offset = phi_dim1 + 1;
  91     phi -= phi_offset;
  92
  93     /* Function Body */
  94     i__1 = *order;
  95     for (j = 1; j <= i__1; ++j) {
  96         i__2 = *order;
  97         for (i__ = j; i__ <= i__2; ++i__) {
  98             v[i__ + j * 10 - 11] = phi[i__ + j * phi_dim1];
  99         }
 100         i__2 = j - 1;
 101         for (k = 1; k <= i__2; ++k) {
 102             save = v[j + k * 10 - 11] * v[k + k * 10 - 11];
 103             i__3 = *order;
 104             for (i__ = j; i__ <= i__3; ++i__) {
 105                 v[i__ + j * 10 - 11] -= v[i__ + k * 10 - 11] * save;
 106             }
 107         }
 108 /*  Compute intermediate results, which are similar to RC's */
 109         if ((r__1 = v[j + j * 10 - 11], abs(r__1)) < 1e-10f) {
 110             goto L100;
 111         }
 112         rc[j] = psi[j];
 113         i__2 = j - 1;
 114         for (k = 1; k <= i__2; ++k) {
 115             rc[j] -= rc[k] * v[j + k * 10 - 11];
 116         }
 117         v[j + j * 10 - 11] = 1.f / v[j + j * 10 - 11];
 118         rc[j] *= v[j + j * 10 - 11];
 119 /* Computing MAX */
 120 /* Computing MIN */
 121         r__2 = rc[j];
 122         r__1 = min(r__2,.999f);
 123         rc[j] = max(r__1,-.999f);
 124     }
 125     return 0;
 126 /*  Zero out higher order RC's if algorithm terminated early */
 127 L100:
 128     i__1 = *order;
 129     for (i__ = j; i__ <= i__1; ++i__) {
 130         rc[i__] = 0.f;
 131     }
 132 /*  Back substitute for PC's (if needed) */
 133 /* 110  DO J = ORDER,1,-1 */
 134 /*         PC(J) = RC(J) */
 135 /*         DO I = 1,J-1 */
 136 /*            PC(J) = PC(J) - PC(I)*V(J,I) */
 137 /*         END DO */
 138 /*      END DO */
 139     return 0;
 140 } /* invert_ */
 141