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