1/* ------------------------------------------------------------------
2 * Copyright (C) 1998-2009 PacketVideo
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 *      http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either
13 * express or implied.
14 * See the License for the specific language governing permissions
15 * and limitations under the License.
16 * -------------------------------------------------------------------
17 */
18/****************************************************************************************
19Portions of this file are derived from the following 3GPP standard:
20
21    3GPP TS 26.173
22    ANSI-C code for the Adaptive Multi-Rate - Wideband (AMR-WB) speech codec
23    Available from http://www.3gpp.org
24
25(C) 2007, 3GPP Organizational Partners (ARIB, ATIS, CCSA, ETSI, TTA, TTC)
26Permission to distribute, modify and use this file under the standard license
27terms listed above has been obtained from the copyright holder.
28****************************************************************************************/
29/*
30------------------------------------------------------------------------------
31
32
33
34 Filename: isp_az.cpp
35
36     Date: 05/08/2004
37
38------------------------------------------------------------------------------
39 REVISION HISTORY
40
41
42 Description:
43
44------------------------------------------------------------------------------
45 INPUT AND OUTPUT DEFINITIONS
46
47     int16 isp[],              (i) Q15 : Immittance spectral pairs
48     int16 a[],                (o) Q12 : predictor coefficients (order=M)
49     int16 m,                  (i)     : order
50     int16 adaptive_scaling    (i) 0   : adaptive scaling disabled
51                                   1   : adaptive scaling enabled
52
53
54------------------------------------------------------------------------------
55 FUNCTION DESCRIPTION
56
57    Compute the LPC coefficients from isp (order=M)
58------------------------------------------------------------------------------
59 REQUIREMENTS
60
61
62------------------------------------------------------------------------------
63 REFERENCES
64
65------------------------------------------------------------------------------
66 PSEUDO-CODE
67
68------------------------------------------------------------------------------
69*/
70
71
72/*----------------------------------------------------------------------------
73; INCLUDES
74----------------------------------------------------------------------------*/
75
76#include "pv_amr_wb_type_defs.h"
77#include "pvamrwbdecoder_basic_op.h"
78#include "pvamrwbdecoder_cnst.h"
79#include "pvamrwbdecoder_acelp.h"
80#include "pvamrwb_math_op.h"
81
82/*----------------------------------------------------------------------------
83; MACROS
84; Define module specific macros here
85----------------------------------------------------------------------------*/
86
87
88/*----------------------------------------------------------------------------
89; DEFINES
90; Include all pre-processor statements here. Include conditional
91; compile variables also.
92----------------------------------------------------------------------------*/
93#define NC (M/2)
94#define NC16k (M16k/2)
95
96/*----------------------------------------------------------------------------
97; LOCAL FUNCTION DEFINITIONS
98; Function Prototype declaration
99----------------------------------------------------------------------------*/
100
101
102#ifdef __cplusplus
103extern "C"
104{
105#endif
106
107    void Get_isp_pol(int16 * isp, int32 * f, int16 n);
108    void Get_isp_pol_16kHz(int16 * isp, int32 * f, int16 n);
109
110#ifdef __cplusplus
111}
112#endif
113
114/*----------------------------------------------------------------------------
115; LOCAL STORE/BUFFER/POINTER DEFINITIONS
116; Variable declaration - defined here and used outside this module
117----------------------------------------------------------------------------*/
118
119/*----------------------------------------------------------------------------
120; EXTERNAL FUNCTION REFERENCES
121; Declare functions defined elsewhere and referenced in this module
122----------------------------------------------------------------------------*/
123
124/*----------------------------------------------------------------------------
125; EXTERNAL GLOBAL STORE/BUFFER/POINTER REFERENCES
126; Declare variables used in this module but defined elsewhere
127----------------------------------------------------------------------------*/
128
129/*----------------------------------------------------------------------------
130; FUNCTION CODE
131----------------------------------------------------------------------------*/
132
133void Isp_Az(
134    int16 isp[],            /* (i) Q15 : Immittance spectral pairs         */
135    int16 a[],              /* (o) Q12 : predictor coefficients (order=M)  */
136    int16 m,                /* (i)     : order                     */
137    int16 adaptive_scaling  /* (i) 0   : adaptive scaling disabled */
138    /*     1   : adaptive scaling enabled  */
139)
140{
141    int16 i, j;
142    int32 f1[NC16k + 1], f2[NC16k];
143    int16 nc;
144    int32 t0;
145    int32 t1;
146    int16 q, q_sug;
147    int32 tmax;
148
149    nc = m >> 1;
150
151
152    if (nc > 8)
153    {
154        Get_isp_pol_16kHz(&isp[0], f1, nc);
155        for (i = 0; i <= nc; i++)
156        {
157            f1[i] = shl_int32(f1[i], 2);
158        }
159        Get_isp_pol_16kHz(&isp[1], f2, nc - 1);
160        for (i = 0; i <= nc - 1; i++)
161        {
162            f2[i] = shl_int32(f2[i], 2);
163        }
164    }
165    else
166    {
167        Get_isp_pol(&isp[0], f1, nc);
168        Get_isp_pol(&isp[1], f2, nc - 1);
169    }
170
171    /*
172     *  Multiply F2(z) by (1 - z^-2)
173     */
174
175    for (i = nc - 1; i > 1; i--)
176    {
177        f2[i] -= f2[i - 2];      /* f2[i] -= f2[i-2]; */
178    }
179
180    /*
181     *  Scale F1(z) by (1+isp[m-1])  and  F2(z) by (1-isp[m-1])
182     */
183
184    for (i = 0; i < nc; i++)
185    {
186        /* f1[i] *= (1.0 + isp[M-1]); */
187
188        /* f2[i] *= (1.0 - isp[M-1]); */
189        t0 = f1[i];
190        t1 = f2[i];
191        t0 = fxp_mul32_by_16b(t0, isp[m - 1]) << 1;
192        t1 = fxp_mul32_by_16b(t1, isp[m - 1]) << 1;
193        f1[i] += t0;
194        f2[i] -= t1;
195
196    }
197
198    /*
199     *  A(z) = (F1(z)+F2(z))/2
200     *  F1(z) is symmetric and F2(z) is antisymmetric
201     */
202
203    /* a[0] = 1.0; */
204    a[0] = 4096;
205    tmax = 1;
206    j = m - 1;
207    for (i = 1;  i < nc; i++)
208    {
209        /* a[i] = 0.5*(f1[i] + f2[i]); */
210
211        t0 = add_int32(f1[i], f2[i]);          /* f1[i] + f2[i]             */
212        /* compute t1 = abs(t0) */
213        t1 = t0 - (t0 < 0);
214        t1 = t1 ^(t1 >> 31);  /* t1 = t1 ^sign(t1) */
215
216        tmax |= t1;
217        /* from Q23 to Q12 and * 0.5 */
218        a[i] = (int16)((t0 >> 12) + ((t0 >> 11) & 1));
219
220
221        /* a[j] = 0.5*(f1[i] - f2[i]); */
222
223        t0 = sub_int32(f1[i], f2[i]);          /* f1[i] - f2[i]             */
224        /* compute t1 = abs(t0) */
225        t1 = t0 - (t0 < 0);
226        t1 = t1 ^(t1 >> 31);  /* t1 = t1 ^sign(t1) */
227
228        tmax |= t1;
229
230        /* from Q23 to Q12 and * 0.5 */
231        a[j--] = (int16)((t0 >> 12) + ((t0 >> 11) & 1));
232
233    }
234
235    /* rescale data if overflow has occured and reprocess the loop */
236
237
238    if (adaptive_scaling == 1)
239    {
240        q = 4 - normalize_amr_wb(tmax);        /* adaptive scaling enabled */
241    }
242    else
243    {
244        q = 0;                   /* adaptive scaling disabled */
245    }
246
247
248    if (q > 0)
249    {
250        q_sug = 12 + q;
251        for (i = 1, j = m - 1; i < nc; i++, j--)
252        {
253            /* a[i] = 0.5*(f1[i] + f2[i]); */
254
255            t0 = add_int32(f1[i], f2[i]);          /* f1[i] + f2[i]             */
256            /* from Q23 to Q12 and * 0.5 */
257            a[i] = (int16)((t0 >> q_sug) + ((t0 >> (q_sug - 1)) & 1));
258
259
260            /* a[j] = 0.5*(f1[i] - f2[i]); */
261
262            t0 = sub_int32(f1[i], f2[i]);          /* f1[i] - f2[i]             */
263            /* from Q23 to Q12 and * 0.5 */
264            a[j] = (int16)((t0 >> q_sug) + ((t0 >> (q_sug - 1)) & 1));
265
266        }
267        a[0] >>=  q;
268    }
269    else
270    {
271        q_sug = 12;
272        q     = 0;
273    }
274
275    /* a[NC] = 0.5*f1[NC]*(1.0 + isp[M-1]); */
276
277
278    t0 = (int32)(((int64)f1[nc] * isp[m - 1]) >> 16) << 1;
279
280
281    t0 = add_int32(f1[nc], t0);
282
283    /* from Q23 to Q12 and * 0.5 */
284    a[nc] = (int16)((t0 >> q_sug) + ((t0 >> (q_sug - 1)) & 1));
285    a[m] = shr_rnd(isp[m - 1], (3 + q));           /* from Q15 to Q12          */
286
287    /* a[m] = isp[m-1]; */
288
289
290    return;
291}
292
293
294
295/*
296Get_isp_pol
297------------------------------------------------------------------------------
298 INPUT AND OUTPUT DEFINITIONS
299
300   isp[]   : isp vector (cosine domaine)         in Q15
301   f[]     : the coefficients of F1 or F2        in Q23
302   n       : == NC for F1(z); == NC-1 for F2(z)
303
304
305------------------------------------------------------------------------------
306 FUNCTION DESCRIPTION
307
308    Find the polynomial F1(z) or F2(z) from the ISPs.
309  This is performed by expanding the product polynomials:
310
311  F1(z) =   product   ( 1 - 2 isp_i z^-1 + z^-2 )
312          i=0,2,4,6,8
313  F2(z) =   product   ( 1 - 2 isp_i z^-1 + z^-2 )
314          i=1,3,5,7
315
316  where isp_i are the ISPs in the cosine domain.
317------------------------------------------------------------------------------
318 REQUIREMENTS
319
320
321------------------------------------------------------------------------------
322 REFERENCES
323
324------------------------------------------------------------------------------
325 PSEUDO-CODE
326
327----------------------------------------------------------------------------*/
328
329/*----------------------------------------------------------------------------
330; FUNCTION CODE
331----------------------------------------------------------------------------*/
332
333
334void Get_isp_pol(int16 * isp, int32 * f, int16 n)
335{
336    int16 i, j;
337    int32 t0;
338
339
340    /* All computation in Q23 */
341
342    f[0] = 0x00800000;                        /* f[0] = 1.0;        in Q23  */
343    f[1] = -isp[0] << 9;                      /* f[1] = -2.0*isp[0] in Q23  */
344
345    f += 2;                                   /* Advance f pointer          */
346    isp += 2;                                 /* Advance isp pointer        */
347
348    for (i = 2; i <= n; i++)
349    {
350        *f = f[-2];
351
352        for (j = 1; j < i; j++)
353        {
354
355            t0 = fxp_mul32_by_16b(f[-1], *isp);
356            t0 = shl_int32(t0, 2);
357
358            *f -= t0;                      /* *f -= t0            */
359            *(f) += f[-2];                 /* *f += f[-2]         */
360            f--;
361
362
363        }
364        *f -= *isp << 9;
365
366        f += i;                            /* Advance f pointer   */
367        isp += 2;                          /* Advance isp pointer */
368    }
369}
370
371void Get_isp_pol_16kHz(int16 * isp, int32 * f, int16 n)
372{
373    int16 i, j;
374    int32 t0;
375
376    /* All computation in Q23 */
377
378    f[0] = 0x00200000;                        /* f[0] = 0.25;        in Q23  */
379
380    f[1] = -isp[0] << 7;                      /* f[1] = -0.5*isp[0] in Q23  */
381
382    f += 2;                                   /* Advance f pointer          */
383    isp += 2;                                 /* Advance isp pointer        */
384
385    for (i = 2; i <= n; i++)
386    {
387        *f = f[-2];
388
389        for (j = 1; j < i; j++, f--)
390        {
391            t0 = fxp_mul32_by_16b(f[-1], *isp);
392            t0 = shl_int32(t0, 2);
393
394            *f -= t0;                      /* *f -= t0            */
395            *f += f[-2];                   /* *f += f[-2]         */
396        }
397        *f -= *isp << 7;
398        f += i;                            /* Advance f pointer   */
399        isp += 2;                          /* Advance isp pointer */
400    }
401    return;
402}
403
404