~ubuntu-branches/ubuntu/saucy/iaxmodem/saucy

« back to all changes in this revision

Viewing changes to lib/spandsp/src/gsm0610_lpc.c

  • Committer: Bazaar Package Importer
  • Author(s): Julien BLACHE
  • Date: 2006-10-28 10:54:55 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20061028105455-qdnaih9tmq0uc29w
Tags: 0.1.15~dfsg-1
* New upstream release.
* debian/rules:
  + Use new ~dfsg versionning scheme.
  + Do not remove config.* in get-orig-source.
* debian/patches/11_build_configure-stamp.dpatch:
  + Updated.
* debian/iaxmodem.init:
  + Added LSB header.

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
/*
 
2
 * SpanDSP - a series of DSP components for telephony
 
3
 *
 
4
 * gsm0610_lpc.c - GSM 06.10 full rate speech codec.
 
5
 *
 
6
 * Written by Steve Underwood <steveu@coppice.org>
 
7
 *
 
8
 * Copyright (C) 2006 Steve Underwood
 
9
 *
 
10
 * All rights reserved.
 
11
 *
 
12
 * This program is free software; you can redistribute it and/or modify
 
13
 * it under the terms of the GNU General Public License as published by
 
14
 * the Free Software Foundation; either version 2 of the License, or
 
15
 * (at your option) any later version.
 
16
 *
 
17
 * This program is distributed in the hope that it will be useful,
 
18
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 
19
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
20
 * GNU General Public License for more details.
 
21
 *
 
22
 * You should have received a copy of the GNU General Public License
 
23
 * along with this program; if not, write to the Free Software
 
24
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 
25
 *
 
26
 * This code is based on the widely used GSM 06.10 code available from
 
27
 * http://kbs.cs.tu-berlin.de/~jutta/toast.html
 
28
 *
 
29
 * $Id: gsm0610_lpc.c,v 1.4 2006/06/13 13:29:46 steveu Exp $
 
30
 */
 
31
 
 
32
/*! \file */
 
33
 
 
34
#ifdef HAVE_CONFIG_H
 
35
#include <config.h>
 
36
#endif
 
37
 
 
38
#include <assert.h>
 
39
#include <inttypes.h>
 
40
#include <tgmath.h>
 
41
#include <stdlib.h>
 
42
 
 
43
#include "spandsp/telephony.h"
 
44
#include "spandsp/dc_restore.h"
 
45
#include "spandsp/gsm0610.h"
 
46
 
 
47
#include "gsm0610_local.h"
 
48
 
 
49
/* 4.2.4 .. 4.2.7 LPC ANALYSIS SECTION */
 
50
 
 
51
static const uint8_t bitoff[256] =
 
52
{
 
53
     8, 7, 6, 6, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4,
 
54
     3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
 
55
     2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
 
56
     2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
 
57
     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
 
58
     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
 
59
     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
 
60
     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
 
61
     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 
62
     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 
63
     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 
64
     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 
65
     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 
66
     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 
67
     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 
68
     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
 
69
};
 
70
 
 
71
/* The number of left shifts needed to normalize the 32 bit
 
72
   variable L_var1 for positive values on the interval
 
73
  
 
74
   with minimum of
 
75
   minimum of 1073741824  (01000000000000000000000000000000) and 
 
76
   maximum of 2147483647  (01111111111111111111111111111111)
 
77
 
 
78
   and for negative values on the interval with
 
79
   minimum of -2147483648 (-10000000000000000000000000000000) and
 
80
   maximum of -1073741824 ( -1000000000000000000000000000000).
 
81
  
 
82
   In order to normalize the result, the following
 
83
   operation must be done: L_norm_var1 = L_var1 << norm(L_var1);
 
84
  
 
85
   (That's 'ffs', only from the left, not the right..)
 
86
*/
 
87
 
 
88
int16_t gsm0610_norm(int32_t a)
 
89
{
 
90
    assert(a != 0);
 
91
 
 
92
    if (a < 0)
 
93
    {
 
94
        if (a <= -1073741824)
 
95
            return  0;
 
96
        /*endif*/
 
97
        a = ~a;
 
98
    }
 
99
    /*endif*/
 
100
 
 
101
    return  (a & 0xFFFF0000)
 
102
            ?
 
103
               ((a & 0xFF000000)
 
104
                ?
 
105
                -1 + bitoff[0xFF & (a >> 24)]
 
106
                :
 
107
                7 + bitoff[0xFF & (a >> 16)])
 
108
            :
 
109
               ((a & 0xFF00)
 
110
                ?
 
111
                15 + bitoff[0xFF & (a >> 8)]
 
112
                :
 
113
                23 + bitoff[0xFF & a]);
 
114
}
 
115
/*- End of function --------------------------------------------------------*/
 
116
 
 
117
/*
 
118
   (From p. 46, end of section 4.2.5)
 
119
  
 
120
   NOTE: The following lines gives [sic] one correct implementation
 
121
         of the div(num, denum) arithmetic operation.  Compute div
 
122
         which is the integer division of num by denom: with
 
123
         denom >= num > 0
 
124
*/
 
125
static int16_t gsm_div(int16_t num, int16_t denom)
 
126
{
 
127
    int32_t L_num;
 
128
    int32_t L_denom;
 
129
    int16_t div;
 
130
    int k;
 
131
 
 
132
    /* The parameter num sometimes becomes zero.
 
133
       Although this is explicitly guarded against in 4.2.5,
 
134
       we assume that the result should then be zero as well. */
 
135
 
 
136
    assert(num >= 0  &&  denom >= num);
 
137
    if (num == 0)
 
138
        return  0;
 
139
    /*endif*/
 
140
    L_num = num;
 
141
    L_denom = denom;
 
142
    div = 0;
 
143
    k = 15;
 
144
    while (k--)
 
145
    {
 
146
        div <<= 1;
 
147
        L_num <<= 1;
 
148
 
 
149
        if (L_num >= L_denom)
 
150
        {
 
151
            L_num -= L_denom;
 
152
            div++;
 
153
        }
 
154
        /*endif*/
 
155
    }
 
156
    /*endwhile*/
 
157
 
 
158
    return  div;
 
159
}
 
160
/*- End of function --------------------------------------------------------*/
 
161
 
 
162
#if defined(__GNUC__)  &&  defined(__i386__)
 
163
 
 
164
#if defined(__CYGWIN32__)  ||  defined(WIN32)
 
165
const int32_t xxx_lower_bound = 0x80008000;
 
166
const int32_t xxx_upper_bound = 0x7fff7fff;
 
167
 
 
168
const int64_t xxx_ones = 0x0001000100010001LL;
 
169
#else
 
170
const int32_t _xxx_lower_bound = 0x80008000;
 
171
const int32_t _xxx_upper_bound = 0x7fff7fff;
 
172
 
 
173
const int64_t _xxx_ones = 0x0001000100010001LL;
 
174
#endif
 
175
 
 
176
void gsm0610_vec_vsraw(const int16_t *p, int n, int bits)
 
177
{
 
178
    if (n == 0)
 
179
        return;
 
180
    /*endif*/
 
181
    __asm__ __volatile__(
 
182
        " leal -16(%%esi,%%eax,2),%%edx;\n"         /* edx = top - 16 */
 
183
        " emms;\n"
 
184
        " movd %%ecx,%%mm3;\n"
 
185
        " movq _xxx_ones,%%mm2;\n"
 
186
        " psllw %%mm3,%%mm2;\n"
 
187
        " psrlw $1,%%mm2;\n"
 
188
        " cmpl %%edx,%%esi;"
 
189
        " ja 4f;\n"
 
190
 
 
191
         " .p2align 2;\n"
 
192
        /* 8 words per iteration */
 
193
        "6:\n"
 
194
        " movq (%%esi),%%mm0;\n"
 
195
        " movq 8(%%esi),%%mm1;\n"
 
196
        " paddsw %%mm2,%%mm0;\n"
 
197
        " psraw %%mm3,%%mm0;\n"
 
198
        " paddsw %%mm2,%%mm1;\n"
 
199
        " psraw %%mm3,%%mm1;\n"
 
200
        " movq %%mm0,(%%esi);\n"
 
201
        " movq %%mm1,8(%%esi);\n"
 
202
        " addl $16,%%esi;\n"
 
203
        " cmpl %%edx,%%esi;\n"
 
204
        " jbe 6b;\n"
 
205
 
 
206
        " .p2align 2;\n"
 
207
        "4:\n"
 
208
        " addl $12,%%edx;\n"                        /* now edx = top-4 */
 
209
        " cmpl %%edx,%%esi;\n"
 
210
        " ja 3f;\n"
 
211
 
 
212
        " .p2align 2;\n"
 
213
        /* do up to 6 words, two per iteration */
 
214
        "5:\n"
 
215
        " movd  (%%esi),%%mm0;\n"
 
216
        " paddsw %%mm2,%%mm0;\n"
 
217
        " psraw %%mm3,%%mm0;\n"
 
218
        " movd %%mm0,(%%esi);\n"
 
219
        " addl $4,%%esi;\n"
 
220
        " cmpl %%edx,%%esi;\n"
 
221
        " jbe 5b;\n"
 
222
 
 
223
        " .p2align 2;\n"
 
224
        "3:\n"
 
225
        " addl $2,%%edx;\n"                        /* now edx = top-2 */
 
226
        " cmpl %%edx,%%esi;\n"
 
227
        " ja 2f;\n"
 
228
        
 
229
        " movzwl (%%esi),%%eax;\n"
 
230
        " movd %%eax,%%mm0;\n"
 
231
        " paddsw %%mm2,%%mm0;\n"
 
232
        " psraw %%mm3,%%mm0;\n"
 
233
        " movd %%mm0,%%eax;\n"
 
234
        " movw %%ax,(%%esi);\n"
 
235
 
 
236
        " .p2align 2;\n"
 
237
        "2:\n"
 
238
        " emms;\n"
 
239
        :
 
240
        : "S" (p), "a" (n), "c" (bits)
 
241
        : "edx"
 
242
    );
 
243
}
 
244
/*- End of function --------------------------------------------------------*/
 
245
 
 
246
int32_t gsm0610_vec_iprod(const int16_t *p, const int16_t *q, int n)
 
247
{
 
248
    int32_t iprod;
 
249
 
 
250
    __asm__ __volatile__(
 
251
        " emms;\n"
 
252
        " pxor %%mm0,%%mm0;\n"
 
253
        " leal -32(%%esi,%%eax,2),%%edx;\n"     /* edx = top - 32 */
 
254
 
 
255
        " cmpl %%edx,%%esi;\n"
 
256
        " ja 1f;\n"
 
257
 
 
258
        /* Work in blocks of 16 int16_t's until we are near the end */
 
259
        " .p2align 2;\n"
 
260
        "2:\n"
 
261
        " movq (%%edi),%%mm1;\n"
 
262
        " movq (%%esi),%%mm2;\n"
 
263
        " pmaddwd %%mm2,%%mm1;\n"
 
264
        " paddd %%mm1,%%mm0;\n"
 
265
        " movq 8(%%edi),%%mm1;\n"
 
266
        " movq 8(%%esi),%%mm2;\n"
 
267
        " pmaddwd %%mm2,%%mm1;\n"
 
268
        " paddd %%mm1,%%mm0;\n"
 
269
        " movq 16(%%edi),%%mm1;\n"
 
270
        " movq 16(%%esi),%%mm2;\n"
 
271
        " pmaddwd %%mm2,%%mm1;\n"
 
272
        " paddd %%mm1,%%mm0;\n"
 
273
        " movq 24(%%edi),%%mm1;\n"
 
274
        " movq 24(%%esi),%%mm2;\n"
 
275
        " pmaddwd %%mm2,%%mm1;\n"
 
276
        " paddd %%mm1,%%mm0;\n"
 
277
 
 
278
        " addl $32,%%esi;\n"
 
279
        " addl $32,%%edi;\n"
 
280
        " cmpl %%edx,%%esi;\n"
 
281
        " jbe 2b;\n"
 
282
 
 
283
        " .p2align 2;\n"
 
284
        "1:\n"
 
285
        " addl $24,%%edx;\n"                  /* now edx = top - 8 */
 
286
        " cmpl %%edx,%%esi;\n"
 
287
        " ja 3f;\n"
 
288
 
 
289
        /* Work in blocks of 4 int16_t's until we are near the end */
 
290
        " .p2align 2;\n"
 
291
        "4:\n"
 
292
        " movq (%%edi),%%mm1;\n"
 
293
        " movq (%%esi),%%mm2;\n"
 
294
        " pmaddwd %%mm2,%%mm1;\n"
 
295
        " paddd %%mm1,%%mm0;\n"
 
296
 
 
297
        " addl $8,%%esi;\n"
 
298
        " addl $8,%%edi;\n"
 
299
        " cmpl %%edx,%%esi;"
 
300
        " jbe 4b;\n"
 
301
 
 
302
        " .p2align 2;\n"
 
303
        "3:\n"
 
304
        " addl $4,%%edx;\n"                  /* now edx = top - 4 */
 
305
        " cmpl %%edx,%%esi;\n"
 
306
        " ja 5f;\n"
 
307
 
 
308
        /* Work in a block of 2 int16_t's */
 
309
        " movd (%%edi),%%mm1;\n"
 
310
        " movd (%%esi),%%mm2;\n"
 
311
        " pmaddwd %%mm2,%%mm1;\n"
 
312
        " paddd %%mm1,%%mm0;\n"
 
313
 
 
314
        " addl $4,%%esi;\n"
 
315
        " addl $4,%%edi;\n"
 
316
 
 
317
        " .p2align 2;\n"
 
318
        "5:\n"
 
319
        " addl $2,%%edx;\n"                  /* now edx = top - 2 */
 
320
        " cmpl %%edx,%%esi;\n"
 
321
        " ja 6f;\n"
 
322
 
 
323
        /* Deal with the very last int16_t, when n is odd */
 
324
        " movswl (%%edi),%%eax;\n"
 
325
        " andl $65535,%%eax;\n"
 
326
        " movd %%eax,%%mm1;\n"
 
327
        " movswl (%%esi),%%eax;\n"
 
328
        " andl $65535,%%eax;\n"
 
329
        " movd %%eax,%%mm2;\n"
 
330
        " pmaddwd %%mm2,%%mm1;\n"
 
331
        " paddd %%mm1,%%mm0;\n"
 
332
 
 
333
        " .p2align 2;\n"
 
334
        "6:\n"
 
335
        /* Merge the pieces of the answer */
 
336
        " movq %%mm0,%%mm1;\n"
 
337
        " punpckhdq %%mm0,%%mm1;\n"
 
338
        " paddd %%mm1,%%mm0;\n"
 
339
        /* et voila, eax has the final result */
 
340
        " movd %%mm0,%%eax;\n"
 
341
 
 
342
        " emms;\n"
 
343
        : "=a" (iprod)
 
344
        : "S" (p), "D" (q), "a" (n)
 
345
        : "cc"
 
346
    );
 
347
 
 
348
    return  iprod;
 
349
}
 
350
/*- End of function --------------------------------------------------------*/
 
351
 
 
352
int32_t gsm0610_vec_maxmin(const int16_t *p, int n, int16_t *out)
 
353
{
 
354
    int32_t lmax;
 
355
 
 
356
    __asm__ __volatile__(
 
357
        " emms;\n"
 
358
        " pushl %%edx;\n"
 
359
        " leal -8(%%esi,%%eax,2),%%edx;\n"
 
360
 
 
361
        " cmpl %%edx,%%esi;\n"
 
362
        " jbe 2f;\n"
 
363
        " movd _xxx_lower_bound,%%mm0;\n"
 
364
        " movd _xxx_upper_bound,%%mm1;\n"
 
365
        " jmp 1f;\n"
 
366
 
 
367
        " .p2align 2;\n"
 
368
        "2:\n"
 
369
        " movq (%%esi),%%mm0;\n"   /* mm0 will be max's */
 
370
        " movq %%mm0,%%mm1;\n"     /* mm1 will be min's */
 
371
        " addl $8,%%esi;\n"
 
372
        " cmpl %%edx,%%esi;\n"
 
373
        " ja 4f;\n"
 
374
 
 
375
        " .p2align 2;\n"
 
376
        "3:\n"
 
377
        " movq (%%esi),%%mm2;\n"
 
378
 
 
379
        " movq %%mm2,%%mm3;\n"
 
380
        " pcmpgtw %%mm0,%%mm3;\n"  /* mm3 is bitmask for words where mm2 > mm0 */ 
 
381
        " movq %%mm3,%%mm4;\n"
 
382
        " pand %%mm2,%%mm3;\n"     /* mm3 is mm2 masked to new max's */
 
383
        " pandn %%mm0,%%mm4;\n"    /* mm4 is mm0 masked to its max's */
 
384
        " por %%mm3,%%mm4;\n"
 
385
        " movq %%mm4,%%mm0;\n"     /* Now mm0 is updated max's */
 
386
        
 
387
        " movq %%mm1,%%mm3;\n"
 
388
        " pcmpgtw %%mm2,%%mm3;\n"  /* mm3 is bitmask for words where mm2 < mm1 */ 
 
389
        " pand %%mm3,%%mm2;\n"     /* mm2 is mm2 masked to new min's */
 
390
        " pandn %%mm1,%%mm3;\n"    /* mm3 is mm1 masked to its min's */
 
391
        " por %%mm3,%%mm2;\n"
 
392
        " movq %%mm2,%%mm1;\n"     /* now mm1 is updated min's */
 
393
 
 
394
        " addl $8,%%esi;\n"
 
395
        " cmpl %%edx,%%esi;\n"
 
396
        " jbe 3b;\n"
 
397
 
 
398
        " .p2align 2;\n"
 
399
        "4:\n"
 
400
        /* Merge down the 4-word max/mins to lower 2 words */
 
401
        " movq %%mm0,%%mm2;\n"
 
402
        " psrlq $32,%%mm2;\n"
 
403
        " movq %%mm2,%%mm3;\n"
 
404
        " pcmpgtw %%mm0,%%mm3;\n"  /* mm3 is bitmask for words where mm2 > mm0 */ 
 
405
        " pand %%mm3,%%mm2;\n"     /* mm2 is mm2 masked to new max's */
 
406
        " pandn %%mm0,%%mm3;\n"    /* mm3 is mm0 masked to its max's */
 
407
        " por %%mm3,%%mm2;\n"
 
408
        " movq %%mm2,%%mm0;\n"     /* now mm0 is updated max's */
 
409
 
 
410
        " movq %%mm1,%%mm2;\n"
 
411
        " psrlq $32,%%mm2;\n"
 
412
        " movq %%mm1,%%mm3;\n"
 
413
        " pcmpgtw %%mm2,%%mm3;\n"  /* mm3 is bitmask for words where mm2 < mm1 */ 
 
414
        " pand %%mm3,%%mm2;\n"     /* mm2 is mm2 masked to new min's */
 
415
        " pandn %%mm1,%%mm3;\n"    /* mm3 is mm1 masked to its min's */
 
416
        " por %%mm3,%%mm2;\n"
 
417
        " movq %%mm2,%%mm1;\n"     /* now mm1 is updated min's */
 
418
 
 
419
        " .p2align 2;\n"
 
420
        "1:\n"
 
421
        " addl $4,%%edx;\n"        /* now dx = top-4 */
 
422
        " cmpl %%edx,%%esi;\n"
 
423
        " ja 5f;\n"
 
424
        /* Here, there are >= 2 words of input remaining */
 
425
        " movd (%%esi),%%mm2;\n"
 
426
 
 
427
        " movq %%mm2,%%mm3;\n"
 
428
        " pcmpgtw %%mm0,%%mm3;\n"  /* mm3 is bitmask for words where mm2 > mm0 */ 
 
429
        " movq %%mm3,%%mm4;\n"
 
430
        " pand %%mm2,%%mm3;\n"     /* mm3 is mm2 masked to new max's */
 
431
        " pandn %%mm0,%%mm4;\n"    /* mm4 is mm0 masked to its max's */
 
432
        " por %%mm3,%%mm4;\n"
 
433
        " movq %%mm4,%%mm0;\n"     /* now mm0 is updated max's */
 
434
 
 
435
        " movq %%mm1,%%mm3;\n"
 
436
        " pcmpgtw %%mm2,%%mm3;\n"  /* mm3 is bitmask for words where mm2 < mm1 */ 
 
437
        " pand %%mm3,%%mm2;\n"     /* mm2 is mm2 masked to new min's */
 
438
        " pandn %%mm1,%%mm3;\n"    /* mm3 is mm1 masked to its min's */
 
439
        " por %%mm3,%%mm2;\n"
 
440
        " movq %%mm2,%%mm1;\n"     /* now mm1 is updated min's */
 
441
 
 
442
        " addl $4,%%esi;\n"
 
443
 
 
444
        " .p2align 2;\n"
 
445
        "5:\n"
 
446
        /* Merge down the 2-word max/mins to 1 word */
 
447
        " movq %%mm0,%%mm2;\n"
 
448
        " psrlq $16,%%mm2;\n"
 
449
        " movq %%mm2,%%mm3;\n"
 
450
        " pcmpgtw %%mm0,%%mm3;\n"  /* mm3 is bitmask for words where mm2 > mm0 */ 
 
451
        " pand %%mm3,%%mm2;\n"     /* mm2 is mm2 masked to new max's */
 
452
        " pandn %%mm0,%%mm3;\n"    /* mm3 is mm0 masked to its max's */
 
453
        " por %%mm3,%%mm2;\n"
 
454
        " movd %%mm2,%%ecx;\n"     /* cx is max so far */
 
455
 
 
456
        " movq %%mm1,%%mm2;\n"
 
457
        " psrlq $16,%%mm2;\n"
 
458
        " movq %%mm1,%%mm3;\n"
 
459
        " pcmpgtw %%mm2,%%mm3;\n"  /* mm3 is bitmask for words where mm2 < mm1 */ 
 
460
        " pand %%mm3,%%mm2;\n"     /* mm2 is mm2 masked to new min's */
 
461
        " pandn %%mm1,%%mm3;\n"    /* mm3 is mm1 masked to its min's */
 
462
        " por %%mm3,%%mm2;\n"
 
463
        " movd %%mm2,%%eax;\n"     /* ax is min so far */
 
464
        
 
465
        " addl $2,%%edx;\n"        /* now dx = top-2 */
 
466
        " cmpl %%edx,%%esi;\n"
 
467
        " ja 6f;\n"
 
468
 
 
469
        /* Here, there is one word of input left */
 
470
        " cmpw (%%esi),%%cx;\n"
 
471
        " jge 9f;\n"
 
472
        " movw (%%esi),%%cx;\n"
 
473
        " .p2align 2;\n"
 
474
        "9:\n"
 
475
        " cmpw (%%esi),%%ax;\n"
 
476
        " jle 6f;\n"
 
477
        " movw (%%esi),%%ax;\n"
 
478
 
 
479
        " .p2align 2;\n"
 
480
        "6:\n"
 
481
        /* (finally!) cx is the max, ax the min */
 
482
        " movswl %%cx,%%ecx;\n"
 
483
        " movswl %%ax,%%eax;\n"
 
484
 
 
485
        " popl %%edx;\n"            /* ptr to output max,min vals */
 
486
        " andl %%edx,%%edx;\n"
 
487
        " jz 7f;\n"
 
488
        " movw %%cx,(%%edx);\n"    /* max */
 
489
        " movw %%ax,2(%%edx);\n"   /* min */
 
490
        " .p2align 2;\n"
 
491
        "7:\n"
 
492
        /* Now calculate max absolute value */
 
493
        " negl %%eax;\n"
 
494
        " cmpl %%ecx,%%eax;\n"
 
495
        " jge 8f;\n"
 
496
        " movl %%ecx,%%eax;\n"
 
497
        " .p2align 2;\n"
 
498
        "8:\n"
 
499
        " emms;\n"
 
500
        : "=a" (lmax)
 
501
        : "S" (p), "a" (n), "d" (out)
 
502
        : "ecx"
 
503
    );
 
504
    return  lmax;
 
505
}
 
506
/*- End of function --------------------------------------------------------*/
 
507
#endif
 
508
 
 
509
/* 4.2.4 */
 
510
static void autocorrelation(int16_t amp[GSM0610_FRAME_LEN], int32_t L_ACF[9])
 
511
{
 
512
    int k;
 
513
    int16_t smax;
 
514
    int16_t scalauto;
 
515
#if !(defined(__GNUC__)  &&  defined(__i386__))
 
516
    int i;
 
517
    int temp;
 
518
    int16_t *sp;
 
519
    int16_t sl;
 
520
#endif
 
521
    
 
522
    /* The goal is to compute the array L_ACF[k].  The signal s[i] must
 
523
       be scaled in order to avoid an overflow situation. */
 
524
 
 
525
    /* Dynamic scaling of the array  s[0..159] */
 
526
    /* Search for the maximum. */
 
527
#if defined(__GNUC__)  &&  defined(__i386__)
 
528
    smax = saturate(gsm0610_vec_maxmin(amp, GSM0610_FRAME_LEN, NULL));
 
529
#else
 
530
    for (smax = 0, k = 0;  k < GSM0610_FRAME_LEN;  k++)
 
531
    {
 
532
        temp = gsm_abs(amp[k]);
 
533
        if (temp > smax)
 
534
            smax = temp;
 
535
        /*endif*/
 
536
    }
 
537
    /*endfor*/
 
538
#endif
 
539
 
 
540
    /* Computation of the scaling factor. */
 
541
    if (smax == 0)
 
542
    {
 
543
        scalauto = 0;
 
544
    }
 
545
    else
 
546
    {
 
547
        assert(smax > 0);
 
548
        scalauto = (int16_t) (4 - gsm0610_norm((int32_t) smax << 16));
 
549
    }
 
550
    /*endif*/
 
551
 
 
552
    /* Scaling of the array s[0...159] */
 
553
#if defined(__GNUC__)  &&  defined(__i386__)
 
554
    if (scalauto > 0)
 
555
        gsm0610_vec_vsraw(amp, GSM0610_FRAME_LEN, scalauto);
 
556
    /*endif*/
 
557
#else
 
558
    if (scalauto > 0)
 
559
    {
 
560
        for (k = 0;  k < GSM0610_FRAME_LEN;  k++)
 
561
            amp[k] = gsm_mult_r(amp[k], 16384 >> (scalauto - 1));
 
562
        /*endfor*/
 
563
    }
 
564
    /*endif*/
 
565
#endif
 
566
 
 
567
    /* Compute the L_ACF[..]. */
 
568
#if defined(__GNUC__)  &&  defined(__i386__)
 
569
    for (k = 0;  k < 9;  k++)
 
570
        L_ACF[k] = gsm0610_vec_iprod(amp, amp + k, GSM0610_FRAME_LEN - k) << 1;
 
571
    /*endfor*/
 
572
#else
 
573
    sp = amp;
 
574
    sl = *sp;
 
575
    L_ACF[0] = ((int32_t) sl*sp[0]);
 
576
    sl = *++sp;
 
577
    L_ACF[0] += ((int32_t) sl*sp[0]);
 
578
    L_ACF[1] = ((int32_t) sl*sp[-1]);
 
579
    sl = *++sp;
 
580
    L_ACF[0] += ((int32_t) sl*sp[0]);
 
581
    L_ACF[1] += ((int32_t) sl*sp[-1]);
 
582
    L_ACF[2] = ((int32_t) sl*sp[-2]);
 
583
    sl = *++sp;
 
584
    L_ACF[0] += ((int32_t) sl*sp[0]);
 
585
    L_ACF[1] += ((int32_t) sl*sp[-1]);
 
586
    L_ACF[2] += ((int32_t) sl*sp[-2]);
 
587
    L_ACF[3] = ((int32_t) sl*sp[-3]);
 
588
    sl = *++sp;
 
589
    L_ACF[0] += ((int32_t) sl*sp[0]);
 
590
    L_ACF[1] += ((int32_t) sl*sp[-1]);
 
591
    L_ACF[2] += ((int32_t) sl*sp[-2]);
 
592
    L_ACF[3] += ((int32_t) sl*sp[-3]);
 
593
    L_ACF[4] = ((int32_t) sl*sp[-4]);
 
594
    sl = *++sp;
 
595
    L_ACF[0] += ((int32_t) sl*sp[0]);
 
596
    L_ACF[1] += ((int32_t) sl*sp[-1]);
 
597
    L_ACF[2] += ((int32_t) sl*sp[-2]);
 
598
    L_ACF[3] += ((int32_t) sl*sp[-3]);
 
599
    L_ACF[4] += ((int32_t) sl*sp[-4]);
 
600
    L_ACF[5] = ((int32_t) sl*sp[-5]);
 
601
    sl = *++sp;
 
602
    L_ACF[0] += ((int32_t) sl*sp[0]);
 
603
    L_ACF[1] += ((int32_t) sl*sp[-1]);
 
604
    L_ACF[2] += ((int32_t) sl*sp[-2]);
 
605
    L_ACF[3] += ((int32_t) sl*sp[-3]);
 
606
    L_ACF[4] += ((int32_t) sl*sp[-4]);
 
607
    L_ACF[5] += ((int32_t) sl*sp[-5]);
 
608
    L_ACF[6] = ((int32_t) sl*sp[-6]);
 
609
    sl = *++sp;
 
610
    L_ACF[0] += ((int32_t) sl*sp[0]);
 
611
    L_ACF[1] += ((int32_t) sl*sp[-1]);
 
612
    L_ACF[2] += ((int32_t) sl*sp[-2]);
 
613
    L_ACF[3] += ((int32_t) sl*sp[-3]);
 
614
    L_ACF[4] += ((int32_t) sl*sp[-4]);
 
615
    L_ACF[5] += ((int32_t) sl*sp[-5]);
 
616
    L_ACF[6] += ((int32_t) sl*sp[-6]);
 
617
    L_ACF[7] = ((int32_t) sl*sp[-7]);
 
618
    sl = *++sp;
 
619
    L_ACF[0] += ((int32_t) sl*sp[0]);
 
620
    L_ACF[1] += ((int32_t) sl*sp[-1]);
 
621
    L_ACF[2] += ((int32_t) sl*sp[-2]);
 
622
    L_ACF[3] += ((int32_t) sl*sp[-3]);
 
623
    L_ACF[4] += ((int32_t) sl*sp[-4]);
 
624
    L_ACF[5] += ((int32_t) sl*sp[-5]);
 
625
    L_ACF[6] += ((int32_t) sl*sp[-6]);
 
626
    L_ACF[7] += ((int32_t) sl*sp[-7]);
 
627
    L_ACF[8] = ((int32_t) sl*sp[-8]);
 
628
    for (i = 9;  i < GSM0610_FRAME_LEN;  i++)
 
629
    {
 
630
        sl = *++sp;
 
631
        L_ACF[0] += ((int32_t) sl*sp[0]);
 
632
        L_ACF[1] += ((int32_t) sl*sp[-1]);
 
633
        L_ACF[2] += ((int32_t) sl*sp[-2]);
 
634
        L_ACF[3] += ((int32_t) sl*sp[-3]);
 
635
        L_ACF[4] += ((int32_t) sl*sp[-4]);
 
636
        L_ACF[5] += ((int32_t) sl*sp[-5]);
 
637
        L_ACF[6] += ((int32_t) sl*sp[-6]);
 
638
        L_ACF[7] += ((int32_t) sl*sp[-7]);
 
639
        L_ACF[8] += ((int32_t) sl*sp[-8]);
 
640
    }
 
641
    /*endfor*/
 
642
    for (k = 0;  k < 9;  k++)
 
643
        L_ACF[k] <<= 1;
 
644
    /*endfor*/
 
645
#endif
 
646
    /* Rescaling of the array s[0..159] */
 
647
    if (scalauto > 0)
 
648
    {
 
649
        assert(scalauto <= 4); 
 
650
        for (k = 0;  k < GSM0610_FRAME_LEN;  k++)
 
651
            amp[k] <<= scalauto;
 
652
        /*endfor*/
 
653
    }
 
654
    /*endif*/
 
655
}
 
656
/*- End of function --------------------------------------------------------*/
 
657
 
 
658
/* 4.2.5 */
 
659
static void reflection_coefficients(int32_t L_ACF[9], int16_t r[8])
 
660
{
 
661
    int i;
 
662
    int m;
 
663
    int n;
 
664
    int16_t temp;
 
665
    int16_t ACF[9]; // 0..8
 
666
    int16_t P[9];   // 0..8
 
667
    int16_t K[9];   // 2..8
 
668
 
 
669
    /* Schur recursion with 16 bits arithmetic. */
 
670
    if (L_ACF[0] == 0)
 
671
    {
 
672
        for (i = 8;  i--;  *r++ = 0)
 
673
            ;
 
674
        /*endfor*/
 
675
        return;
 
676
    }
 
677
    /*endif*/
 
678
 
 
679
    assert(L_ACF[0] != 0);
 
680
    temp = gsm0610_norm(L_ACF[0]);
 
681
 
 
682
    assert(temp >= 0  &&  temp < 32);
 
683
 
 
684
    /* ? overflow ? */
 
685
    for (i = 0;  i <= 8;  i++)
 
686
        ACF[i] = (int16_t) ((L_ACF[i] << temp) >> 16);
 
687
    /*endfor*/
 
688
 
 
689
    /* Initialize array P[..] and K[..] for the recursion. */
 
690
    for (i = 1;  i <= 7;  i++)
 
691
        K[i] = ACF[i];
 
692
    /*endfor*/
 
693
    for (i = 0;  i <= 8;  i++)
 
694
        P[i] = ACF[i];
 
695
    /*endfor*/
 
696
    /* Compute reflection coefficients */
 
697
    for (n = 1;  n <= 8;  n++, r++)
 
698
    {
 
699
        temp = P[1];
 
700
        temp = gsm_abs (temp);
 
701
        if (P[0] < temp)
 
702
        {
 
703
            for (i = n;  i <= 8;  i++)
 
704
                *r++ = 0;
 
705
            /*endfor*/
 
706
            return;
 
707
        }
 
708
        /*endif*/
 
709
 
 
710
        *r = gsm_div(temp, P[0]);
 
711
 
 
712
        assert(*r >= 0);
 
713
        if (P[1] > 0)
 
714
            *r = -*r;     /* r[n] = sub(0, r[n]) */
 
715
        /*endif*/
 
716
        assert(*r != INT16_MIN);
 
717
        if (n == 8)
 
718
            return; 
 
719
        /*endif*/
 
720
 
 
721
        /* Schur recursion */
 
722
        temp = gsm_mult_r(P[1], *r);
 
723
        P[0] = gsm_add(P[0], temp);
 
724
 
 
725
        for (m = 1;  m <= 8 - n;  m++)
 
726
        {
 
727
            temp = gsm_mult_r(K[m], *r);
 
728
            P[m] = gsm_add(P[m + 1], temp);
 
729
 
 
730
            temp = gsm_mult_r(P[m + 1], *r);
 
731
            K[m] = gsm_add(K[m], temp);
 
732
        }
 
733
        /*endfor*/
 
734
    }
 
735
    /*endfor*/
 
736
}
 
737
/*- End of function --------------------------------------------------------*/
 
738
 
 
739
/* 4.2.6 */
 
740
static void transform_to_log_area_ratios(int16_t r[8])
 
741
{
 
742
    int16_t temp;
 
743
    int i;
 
744
 
 
745
    /* The following scaling for r[..] and LAR[..] has been used:
 
746
      
 
747
       r[..]   = integer (real_r[..]*32768.); -1 <= real_r < 1.
 
748
       LAR[..] = integer (real_LAR[..] * 16384);
 
749
       with -1.625 <= real_LAR <= 1.625
 
750
    */
 
751
 
 
752
    /* Computation of the LAR[0..7] from the r[0..7] */
 
753
    for (i = 1;  i <= 8;  i++, r++)
 
754
    {
 
755
        temp = *r;
 
756
        temp = gsm_abs(temp);
 
757
        assert(temp >= 0);
 
758
 
 
759
        if (temp < 22118)
 
760
        {
 
761
            temp >>= 1;
 
762
        }
 
763
        else if (temp < 31130)
 
764
        {
 
765
            assert(temp >= 11059);
 
766
            temp -= 11059;
 
767
        }
 
768
        else
 
769
        {
 
770
            assert(temp >= 26112);
 
771
            temp -= 26112;
 
772
            temp <<= 2;
 
773
        }
 
774
        /*endif*/
 
775
 
 
776
        *r = (*r < 0)  ?  -temp  :  temp;
 
777
        assert(*r != INT16_MIN);
 
778
    }
 
779
    /*endfor*/
 
780
}
 
781
/*- End of function --------------------------------------------------------*/
 
782
 
 
783
/* 4.2.7 */
 
784
static void quantization_and_coding(int16_t LAR[8])
 
785
{
 
786
    int16_t temp;
 
787
 
 
788
    /* This procedure needs four tables; the following equations
 
789
       give the optimum scaling for the constants:
 
790
        
 
791
       A[0..7] = integer(real_A[0..7] * 1024)
 
792
       B[0..7] = integer(real_B[0..7] *  512)
 
793
       MAC[0..7] = maximum of the LARc[0..7]
 
794
       MIC[0..7] = minimum of the LARc[0..7] */
 
795
 
 
796
#undef STEP
 
797
#define STEP(A,B,MAC,MIC)                                       \
 
798
        temp = gsm_mult(A, *LAR);                               \
 
799
        temp = gsm_add(temp, B);                                \
 
800
        temp = gsm_add(temp, 256);                              \
 
801
        temp >>= 9;                                             \
 
802
        *LAR  = (int16_t) ((temp > MAC)                         \
 
803
                         ?                                      \
 
804
                         MAC - MIC                              \
 
805
                         :                                      \
 
806
                         ((temp < MIC)  ?  0  :  temp - MIC));  \
 
807
        LAR++;
 
808
 
 
809
    STEP(20480,     0,  31, -32);
 
810
    STEP(20480,     0,  31, -32);
 
811
    STEP(20480,  2048,  15, -16);
 
812
    STEP(20480, -2560,  15, -16);
 
813
 
 
814
    STEP(13964,    94,   7,  -8);
 
815
    STEP(15360, -1792,   7,  -8);
 
816
    STEP( 8534,  -341,   3,  -4);
 
817
    STEP( 9036, -1144,   3,  -4);
 
818
#undef STEP
 
819
}
 
820
/*- End of function --------------------------------------------------------*/
 
821
 
 
822
void gsm0610_lpc_analysis(gsm0610_state_t *s,
 
823
                          int16_t amp[GSM0610_FRAME_LEN],
 
824
                          int16_t LARc[8])
 
825
{
 
826
    int32_t L_ACF[9];
 
827
 
 
828
    autocorrelation(amp, L_ACF);
 
829
    reflection_coefficients(L_ACF, LARc);
 
830
    transform_to_log_area_ratios(LARc);
 
831
    quantization_and_coding(LARc);
 
832
}
 
833
/*- End of function --------------------------------------------------------*/
 
834
/*- End of file ------------------------------------------------------------*/