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