comparison libgsmfr2/lpc.c @ 268:0cfb7c95cce2

libgsmfr2: integrate lpc.c from libgsm
author Mychaela Falconia <falcon@freecalypso.org>
date Sun, 14 Apr 2024 01:50:48 +0000
parents
children
comparison
equal deleted inserted replaced
267:65d3304502bd 268:0cfb7c95cce2
1 /*
2 * This C source file has been adapted from TU-Berlin libgsm source,
3 * original notice follows:
4 *
5 * Copyright 1992 by Jutta Degener and Carsten Bormann, Technische
6 * Universitaet Berlin. See the accompanying file "COPYRIGHT" for
7 * details. THERE IS ABSOLUTELY NO WARRANTY FOR THIS SOFTWARE.
8 */
9
10 #include <stdint.h>
11 #include <assert.h>
12 #include "tw_gsmfr.h"
13 #include "typedef.h"
14 #include "ed_state.h"
15 #include "ed_internal.h"
16
17 /*
18 * 4.2.4 .. 4.2.7 LPC ANALYSIS SECTION
19 */
20
21 /* 4.2.4 */
22
23 static void Autocorrelation (
24 word * s, /* [0..159] IN/OUT */
25 longword * L_ACF) /* [0..8] OUT */
26 /*
27 * The goal is to compute the array L_ACF[k]. The signal s[i] must
28 * be scaled in order to avoid an overflow situation.
29 */
30 {
31 register int k, i;
32
33 word temp, smax, scalauto;
34
35 /* Dynamic scaling of the array s[0..159]
36 */
37
38 /* Search for the maximum.
39 */
40 smax = 0;
41 for (k = 0; k <= 159; k++) {
42 temp = GSM_ABS( s[k] );
43 if (temp > smax) smax = temp;
44 }
45
46 /* Computation of the scaling factor.
47 */
48 if (smax == 0) scalauto = 0;
49 else {
50 assert(smax > 0);
51 scalauto = 4 - gsm_norm( (longword)smax << 16 );/* sub(4,..) */
52 }
53
54 /* Scaling of the array s[0...159]
55 */
56
57 if (scalauto > 0) {
58
59 # define SCALE(n) \
60 case n: for (k = 0; k <= 159; k++) \
61 s[k] = GSM_MULT_R( s[k], 16384 >> (n-1) );\
62 break;
63
64 switch (scalauto) {
65 SCALE(1)
66 SCALE(2)
67 SCALE(3)
68 SCALE(4)
69 }
70 # undef SCALE
71 }
72
73 /* Compute the L_ACF[..].
74 */
75 {
76 word * sp = s;
77 word sl = *sp;
78
79 # define STEP(k) L_ACF[k] += ((longword)sl * sp[ -(k) ]);
80
81 # define NEXTI sl = *++sp
82
83
84 for (k = 9; k--; L_ACF[k] = 0) ;
85
86 STEP (0);
87 NEXTI;
88 STEP(0); STEP(1);
89 NEXTI;
90 STEP(0); STEP(1); STEP(2);
91 NEXTI;
92 STEP(0); STEP(1); STEP(2); STEP(3);
93 NEXTI;
94 STEP(0); STEP(1); STEP(2); STEP(3); STEP(4);
95 NEXTI;
96 STEP(0); STEP(1); STEP(2); STEP(3); STEP(4); STEP(5);
97 NEXTI;
98 STEP(0); STEP(1); STEP(2); STEP(3); STEP(4); STEP(5); STEP(6);
99 NEXTI;
100 STEP(0); STEP(1); STEP(2); STEP(3); STEP(4); STEP(5); STEP(6); STEP(7);
101
102 for (i = 8; i <= 159; i++) {
103
104 NEXTI;
105
106 STEP(0);
107 STEP(1); STEP(2); STEP(3); STEP(4);
108 STEP(5); STEP(6); STEP(7); STEP(8);
109 }
110
111 for (k = 9; k--; L_ACF[k] <<= 1) ;
112
113 }
114 /* Rescaling of the array s[0..159]
115 */
116 if (scalauto > 0) {
117 assert(scalauto <= 4);
118 for (k = 160; k--; *s++ <<= scalauto) ;
119 }
120 }
121
122 /* 4.2.5 */
123
124 static void Reflection_coefficients (
125 longword * L_ACF, /* 0...8 IN */
126 register word * r /* 0...7 OUT */
127 )
128 {
129 register int i, m, n;
130 register word temp;
131 register longword ltmp;
132 word ACF[9]; /* 0..8 */
133 word P[ 9]; /* 0..8 */
134 word K[ 9]; /* 2..8 */
135
136 /* Schur recursion with 16 bits arithmetic.
137 */
138
139 if (L_ACF[0] == 0) {
140 for (i = 8; i--; *r++ = 0) ;
141 return;
142 }
143
144 assert( L_ACF[0] != 0 );
145 temp = gsm_norm( L_ACF[0] );
146
147 assert(temp >= 0 && temp < 32);
148
149 /* ? overflow ? */
150 for (i = 0; i <= 8; i++) ACF[i] = SASR( L_ACF[i] << temp, 16 );
151
152 /* Initialize array P[..] and K[..] for the recursion.
153 */
154
155 for (i = 1; i <= 7; i++) K[ i ] = ACF[ i ];
156 for (i = 0; i <= 8; i++) P[ i ] = ACF[ i ];
157
158 /* Compute reflection coefficients
159 */
160 for (n = 1; n <= 8; n++, r++) {
161
162 temp = P[1];
163 temp = GSM_ABS(temp);
164 if (P[0] < temp) {
165 for (i = n; i <= 8; i++) *r++ = 0;
166 return;
167 }
168
169 *r = gsm_div( temp, P[0] );
170
171 assert(*r >= 0);
172 if (P[1] > 0) *r = -*r; /* r[n] = sub(0, r[n]) */
173 assert (*r != MIN_WORD);
174 if (n == 8) return;
175
176 /* Schur recursion
177 */
178 temp = GSM_MULT_R( P[1], *r );
179 P[0] = GSM_ADD( P[0], temp );
180
181 for (m = 1; m <= 8 - n; m++) {
182 temp = GSM_MULT_R( K[ m ], *r );
183 P[m] = GSM_ADD( P[ m+1 ], temp );
184
185 temp = GSM_MULT_R( P[ m+1 ], *r );
186 K[m] = GSM_ADD( K[ m ], temp );
187 }
188 }
189 }
190
191 /* 4.2.6 */
192
193 static void Transformation_to_Log_Area_Ratios (
194 register word * r /* 0..7 IN/OUT */
195 )
196 /*
197 * The following scaling for r[..] and LAR[..] has been used:
198 *
199 * r[..] = integer( real_r[..]*32768. ); -1 <= real_r < 1.
200 * LAR[..] = integer( real_LAR[..] * 16384 );
201 * with -1.625 <= real_LAR <= 1.625
202 */
203 {
204 register word temp;
205 register int i;
206
207 /* Computation of the LAR[0..7] from the r[0..7]
208 */
209 for (i = 1; i <= 8; i++, r++) {
210
211 temp = *r;
212 temp = GSM_ABS(temp);
213 assert(temp >= 0);
214
215 if (temp < 22118) {
216 temp >>= 1;
217 } else if (temp < 31130) {
218 assert( temp >= 11059 );
219 temp -= 11059;
220 } else {
221 assert( temp >= 26112 );
222 temp -= 26112;
223 temp <<= 2;
224 }
225
226 *r = *r < 0 ? -temp : temp;
227 assert( *r != MIN_WORD );
228 }
229 }
230
231 /* 4.2.7 */
232
233 static void Quantization_and_coding (
234 register word * LAR /* [0..7] IN/OUT */
235 )
236 {
237 register word temp;
238 longword ltmp;
239
240 /* This procedure needs four tables; the following equations
241 * give the optimum scaling for the constants:
242 *
243 * A[0..7] = integer( real_A[0..7] * 1024 )
244 * B[0..7] = integer( real_B[0..7] * 512 )
245 * MAC[0..7] = maximum of the LARc[0..7]
246 * MIC[0..7] = minimum of the LARc[0..7]
247 */
248
249 # undef STEP
250 # define STEP( A, B, MAC, MIC ) \
251 temp = GSM_MULT( A, *LAR ); \
252 temp = GSM_ADD( temp, B ); \
253 temp = GSM_ADD( temp, 256 ); \
254 temp = SASR( temp, 9 ); \
255 *LAR = temp>MAC ? MAC - MIC : (temp<MIC ? 0 : temp - MIC); \
256 LAR++;
257
258 STEP( 20480, 0, 31, -32 );
259 STEP( 20480, 0, 31, -32 );
260 STEP( 20480, 2048, 15, -16 );
261 STEP( 20480, -2560, 15, -16 );
262
263 STEP( 13964, 94, 7, -8 );
264 STEP( 15360, -1792, 7, -8 );
265 STEP( 8534, -341, 3, -4 );
266 STEP( 9036, -1144, 3, -4 );
267
268 # undef STEP
269 }
270
271 void Gsm_LPC_Analysis (
272 struct gsmfr_0610_state *S,
273 word * s, /* 0..159 signals IN/OUT */
274 word * LARc) /* 0..7 LARc's OUT */
275 {
276 longword L_ACF[9];
277
278 Autocorrelation (s, L_ACF );
279 Reflection_coefficients (L_ACF, LARc );
280 Transformation_to_Log_Area_Ratios (LARc);
281 Quantization_and_coding (LARc);
282 }