← Back to branch summary

~ubuntu-branches/ubuntu/gutsy/audacity/gutsy-backports

« back to all changes in this revision

Viewing changes to lib-src/libsndfile/src/GSM610/lpc.c

  • Committer: Bazaar Package Importer
  • Author(s): John Dong
  • Date: 2008-02-18 21:58:19 UTC
  • mfrom: (13.1.2 hardy)
  • Revision ID: james.westby@ubuntu.com-20080218215819-tmbcf1rx238r8gdv
Tags: 1.3.4-1.1ubuntu1~gutsy1
Automated backport upload; no source changes.

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib-src/libsndfile/src/GSM610/lpc.c
1
 
/*
2
 
 * Copyright 1992 by Jutta Degener and Carsten Bormann, Technische
3
 
 * Universitaet Berlin.  See the accompanying file "COPYRIGHT" for
4
 
 * details.  THERE IS ABSOLUTELY NO WARRANTY FOR THIS SOFTWARE.
5
 
 */
6
 
 
7
 
#include <stdio.h>
8
 
#include <assert.h>
9
 
 
10
 
#include "gsm610_priv.h"
11
 
 
12
 
#include "gsm.h"
13
 
 
14
 
/*
15
 
 *  4.2.4 .. 4.2.7 LPC ANALYSIS SECTION
16
 
 */
17
 
 
18
 
/* 4.2.4 */
19
 
 
20
 
 
21
 
static void Autocorrelation (
22
 
        word     * s,           /* [0..159]     IN/OUT  */
23
 
        longword * L_ACF)       /* [0..8]       OUT     */
24
 
/*
25
 
 *  The goal is to compute the array L_ACF[k].  The signal s[i] must
26
 
 *  be scaled in order to avoid an overflow situation.
27
 
 */
28
 
{
29
 
        register int    k, i;
30
 
 
31
 
        word            temp, smax, scalauto;
32
 
 
33
 
#ifdef  USE_FLOAT_MUL
34
 
        float           float_s[160];
35
 
#endif
36
 
 
37
 
        /*  Dynamic scaling of the array  s[0..159]
38
 
         */
39
 
 
40
 
        /*  Search for the maximum.
41
 
         */
42
 
        smax = 0;
43
 
        for (k = 0; k <= 159; k++) {
44
 
                temp = GSM_ABS( s[k] );
45
 
                if (temp > smax) smax = temp;
46
 
        }
47
 
 
48
 
        /*  Computation of the scaling factor.
49
 
         */
50
 
        if (smax == 0) scalauto = 0;
51
 
        else {
52
 
                assert(smax > 0);
53
 
                scalauto = 4 - gsm_norm( (longword)smax << 16 );/* sub(4,..) */
54
 
        }
55
 
 
56
 
        /*  Scaling of the array s[0...159]
57
 
         */
58
 
 
59
 
        if (scalauto > 0) {
60
 
 
61
 
# ifdef USE_FLOAT_MUL
62
 
#   define SCALE(n)     \
63
 
        case n: for (k = 0; k <= 159; k++) \
64
 
                        float_s[k] = (float)    \
65
 
                                (s[k] = GSM_MULT_R(s[k], 16384 >> (n-1)));\
66
 
                break;
67
 
# else 
68
 
#   define SCALE(n)     \
69
 
        case n: for (k = 0; k <= 159; k++) \
70
 
                        s[k] = GSM_MULT_R( s[k], 16384 >> (n-1) );\
71
 
                break;
72
 
# endif /* USE_FLOAT_MUL */
73
 
 
74
 
                switch (scalauto) {
75
 
                SCALE(1)
76
 
                SCALE(2)
77
 
                SCALE(3)
78
 
                SCALE(4)
79
 
                }
80
 
# undef SCALE
81
 
        }
82
 
# ifdef USE_FLOAT_MUL
83
 
        else for (k = 0; k <= 159; k++) float_s[k] = (float) s[k];
84
 
# endif
85
 
 
86
 
        /*  Compute the L_ACF[..].
87
 
         */
88
 
        {
89
 
# ifdef USE_FLOAT_MUL
90
 
                register float * sp = float_s;
91
 
                register float   sl = *sp;
92
 
 
93
 
#               define STEP(k)   L_ACF[k] += (longword)(sl * sp[ -(k) ]);
94
 
# else
95
 
                word  * sp = s;
96
 
                word    sl = *sp;
97
 
 
98
 
#               define STEP(k)   L_ACF[k] += ((longword)sl * sp[ -(k) ]);
99
 
# endif
100
 
 
101
 
#       define NEXTI     sl = *++sp
102
 
 
103
 
 
104
 
        for (k = 9; k--; L_ACF[k] = 0) ;
105
 
 
106
 
        STEP (0);
107
 
        NEXTI;
108
 
        STEP(0); STEP(1);
109
 
        NEXTI;
110
 
        STEP(0); STEP(1); STEP(2);
111
 
        NEXTI;
112
 
        STEP(0); STEP(1); STEP(2); STEP(3);
113
 
        NEXTI;
114
 
        STEP(0); STEP(1); STEP(2); STEP(3); STEP(4);
115
 
        NEXTI;
116
 
        STEP(0); STEP(1); STEP(2); STEP(3); STEP(4); STEP(5);
117
 
        NEXTI;
118
 
        STEP(0); STEP(1); STEP(2); STEP(3); STEP(4); STEP(5); STEP(6);
119
 
        NEXTI;
120
 
        STEP(0); STEP(1); STEP(2); STEP(3); STEP(4); STEP(5); STEP(6); STEP(7);
121
 
 
122
 
        for (i = 8; i <= 159; i++) {
123
 
 
124
 
                NEXTI;
125
 
 
126
 
                STEP(0);
127
 
                STEP(1); STEP(2); STEP(3); STEP(4);
128
 
                STEP(5); STEP(6); STEP(7); STEP(8);
129
 
        }
130
 
 
131
 
        for (k = 9; k--; L_ACF[k] <<= 1) ; 
132
 
 
133
 
        }
134
 
        /*   Rescaling of the array s[0..159]
135
 
         */
136
 
        if (scalauto > 0) {
137
 
                assert(scalauto <= 4); 
138
 
                for (k = 160; k--; *s++ <<= scalauto) ;
139
 
        }
140
 
}
141
 
 
142
 
#if defined(USE_FLOAT_MUL) && defined(FAST)
143
 
 
144
 
static void Fast_Autocorrelation (
145
 
        word * s,               /* [0..159]     IN/OUT  */
146
 
        longword * L_ACF)       /* [0..8]       OUT     */
147
 
{
148
 
        register int    k, i;
149
 
        float f_L_ACF[9];
150
 
        float scale;
151
 
 
152
 
        float          s_f[160];
153
 
        register float *sf = s_f;
154
 
 
155
 
        for (i = 0; i < 160; ++i) sf[i] = s[i];
156
 
        for (k = 0; k <= 8; k++) {
157
 
                register float L_temp2 = 0;
158
 
                register float *sfl = sf - k;
159
 
                for (i = k; i < 160; ++i) L_temp2 += sf[i] * sfl[i];
160
 
                f_L_ACF[k] = L_temp2;
161
 
        }
162
 
        scale = MAX_LONGWORD / f_L_ACF[0];
163
 
 
164
 
        for (k = 0; k <= 8; k++) {
165
 
                L_ACF[k] = f_L_ACF[k] * scale;
166
 
        }
167
 
}
168
 
#endif  /* defined (USE_FLOAT_MUL) && defined (FAST) */
169
 
 
170
 
/* 4.2.5 */
171
 
 
172
 
static void Reflection_coefficients (
173
 
        longword        * L_ACF,                /* 0...8        IN      */
174
 
        register word   * r                     /* 0...7        OUT     */
175
 
)
176
 
{
177
 
        register int    i, m, n;
178
 
        register word   temp;
179
 
        word            ACF[9]; /* 0..8 */
180
 
        word            P[  9]; /* 0..8 */
181
 
        word            K[  9]; /* 2..8 */
182
 
 
183
 
        /*  Schur recursion with 16 bits arithmetic.
184
 
         */
185
 
 
186
 
        if (L_ACF[0] == 0) {
187
 
                for (i = 8; i--; *r++ = 0) ;
188
 
                return;
189
 
        }
190
 
 
191
 
        assert( L_ACF[0] != 0 );
192
 
        temp = gsm_norm( L_ACF[0] );
193
 
 
194
 
        assert(temp >= 0 && temp < 32);
195
 
 
196
 
        /* ? overflow ? */
197
 
        for (i = 0; i <= 8; i++) ACF[i] = SASR_L( L_ACF[i] << temp, 16 );
198
 
 
199
 
        /*   Initialize array P[..] and K[..] for the recursion.
200
 
         */
201
 
 
202
 
        for (i = 1; i <= 7; i++) K[ i ] = ACF[ i ];
203
 
        for (i = 0; i <= 8; i++) P[ i ] = ACF[ i ];
204
 
 
205
 
        /*   Compute reflection coefficients
206
 
         */
207
 
        for (n = 1; n <= 8; n++, r++) {
208
 
 
209
 
                temp = P[1];
210
 
                temp = GSM_ABS(temp);
211
 
                if (P[0] < temp) {
212
 
                        for (i = n; i <= 8; i++) *r++ = 0;
213
 
                        return;
214
 
                }
215
 
 
216
 
                *r = gsm_div( temp, P[0] );
217
 
 
218
 
                assert(*r >= 0);
219
 
                if (P[1] > 0) *r = -*r;         /* r[n] = sub(0, r[n]) */
220
 
                assert (*r != MIN_WORD);
221
 
                if (n == 8) return; 
222
 
 
223
 
                /*  Schur recursion
224
 
                 */
225
 
                temp = GSM_MULT_R( P[1], *r );
226
 
                P[0] = GSM_ADD( P[0], temp );
227
 
 
228
 
                for (m = 1; m <= 8 - n; m++) {
229
 
                        temp     = GSM_MULT_R( K[ m   ],    *r );
230
 
                        P[m]     = GSM_ADD(    P[ m+1 ],  temp );
231
 
 
232
 
                        temp     = GSM_MULT_R( P[ m+1 ],    *r );
233
 
                        K[m]     = GSM_ADD(    K[ m   ],  temp );
234
 
                }
235
 
        }
236
 
}
237
 
 
238
 
/* 4.2.6 */
239
 
 
240
 
static void Transformation_to_Log_Area_Ratios (
241
 
        register word   * r                     /* 0..7    IN/OUT */
242
 
)
243
 
/*
244
 
 *  The following scaling for r[..] and LAR[..] has been used:
245
 
 *
246
 
 *  r[..]   = integer( real_r[..]*32768. ); -1 <= real_r < 1.
247
 
 *  LAR[..] = integer( real_LAR[..] * 16384 );
248
 
 *  with -1.625 <= real_LAR <= 1.625
249
 
 */
250
 
{
251
 
        register word   temp;
252
 
        register int    i;
253
 
 
254
 
 
255
 
        /* Computation of the LAR[0..7] from the r[0..7]
256
 
         */
257
 
        for (i = 1; i <= 8; i++, r++) {
258
 
 
259
 
                temp = *r;
260
 
                temp = GSM_ABS(temp);
261
 
                assert(temp >= 0);
262
 
 
263
 
                if (temp < 22118) {
264
 
                        temp >>= 1;
265
 
                } else if (temp < 31130) {
266
 
                        assert( temp >= 11059 );
267
 
                        temp -= 11059;
268
 
                } else {
269
 
                        assert( temp >= 26112 );
270
 
                        temp -= 26112;
271
 
                        temp <<= 2;
272
 
                }
273
 
 
274
 
                *r = *r < 0 ? -temp : temp;
275
 
                assert( *r != MIN_WORD );
276
 
        }
277
 
}
278
 
 
279
 
/* 4.2.7 */
280
 
 
281
 
static void Quantization_and_coding (
282
 
        register word * LAR     /* [0..7]       IN/OUT  */
283
 
)
284
 
{
285
 
        register word   temp;
286
 
 
287
 
        /*  This procedure needs four tables; the following equations
288
 
         *  give the optimum scaling for the constants:
289
 
         *  
290
 
         *  A[0..7] = integer( real_A[0..7] * 1024 )
291
 
         *  B[0..7] = integer( real_B[0..7] *  512 )
292
 
         *  MAC[0..7] = maximum of the LARc[0..7]
293
 
         *  MIC[0..7] = minimum of the LARc[0..7]
294
 
         */
295
 
 
296
 
#       undef STEP
297
 
#       define  STEP( A, B, MAC, MIC )          \
298
 
                temp = GSM_MULT( A,   *LAR );   \
299
 
                temp = GSM_ADD(  temp,   B );   \
300
 
                temp = GSM_ADD(  temp, 256 );   \
301
 
                temp = SASR_W(     temp,   9 ); \
302
 
                *LAR  =  temp>MAC ? MAC - MIC : (temp<MIC ? 0 : temp - MIC); \
303
 
                LAR++;
304
 
 
305
 
        STEP(  20480,     0,  31, -32 );
306
 
        STEP(  20480,     0,  31, -32 );
307
 
        STEP(  20480,  2048,  15, -16 );
308
 
        STEP(  20480, -2560,  15, -16 );
309
 
 
310
 
        STEP(  13964,    94,   7,  -8 );
311
 
        STEP(  15360, -1792,   7,  -8 );
312
 
        STEP(   8534,  -341,   3,  -4 );
313
 
        STEP(   9036, -1144,   3,  -4 );
314
 
 
315
 
#       undef   STEP
316
 
}
317
 
 
318
 
void Gsm_LPC_Analysis (
319
 
        struct gsm_state *S,
320
 
        word             * s,           /* 0..159 signals       IN/OUT  */
321
 
        word             * LARc)        /* 0..7   LARc's        OUT     */
322
 
{
323
 
        longword        L_ACF[9];
324
 
 
325
 
#if defined(USE_FLOAT_MUL) && defined(FAST)
326
 
        if (S->fast) Fast_Autocorrelation (s,     L_ACF );
327
 
        else
328
 
#endif
329
 
        Autocorrelation                   (s,     L_ACF );
330
 
        Reflection_coefficients           (L_ACF, LARc  );
331
 
        Transformation_to_Log_Area_Ratios (LARc);
332
 
        Quantization_and_coding           (LARc);
333
 
}
334
 
/*
335
 
** Do not edit or modify anything in this comment block.
336
 
** The arch-tag line is a file identity tag for the GNU Arch 
337
 
** revision control system.
338
 
**
339
 
** arch-tag: 63146664-a002-4e1e-8b7b-f0cc8a6a53da
340
 
*/
341