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