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
invert_(integer * order,real * phi,real * psi,real * rc)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
142