[go: nahoru, domu]

1/*
2 ** Copyright 2003-2010, VisualOn, Inc.
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 express or implied.
13 ** See the License for the specific language governing permissions and
14 ** limitations under the License.
15 */
16
17/***********************************************************************
18*  File: az_isp.c
19*
20*  Description:
21*-----------------------------------------------------------------------*
22* Compute the ISPs from  the LPC coefficients  (order=M)                *
23*-----------------------------------------------------------------------*
24*                                                                       *
25* The ISPs are the roots of the two polynomials F1(z) and F2(z)         *
26* defined as                                                            *
27*               F1(z) = A(z) + z^-m A(z^-1)                             *
28*  and          F2(z) = A(z) - z^-m A(z^-1)                             *
29*                                                                       *
30* For a even order m=2n, F1(z) has M/2 conjugate roots on the unit      *
31* circle and F2(z) has M/2-1 conjugate roots on the unit circle in      *
32* addition to two roots at 0 and pi.                                    *
33*                                                                       *
34* For a 16th order LP analysis, F1(z) and F2(z) can be written as       *
35*                                                                       *
36*   F1(z) = (1 + a[M])   PRODUCT  (1 - 2 cos(w_i) z^-1 + z^-2 )         *
37*                        i=0,2,4,6,8,10,12,14                           *
38*                                                                       *
39*   F2(z) = (1 - a[M]) (1 - z^-2) PRODUCT (1 - 2 cos(w_i) z^-1 + z^-2 ) *
40*                                 i=1,3,5,7,9,11,13                     *
41*                                                                       *
42* The ISPs are the M-1 frequencies w_i, i=0...M-2 plus the last         *
43* predictor coefficient a[M].                                           *
44*-----------------------------------------------------------------------*
45
46************************************************************************/
47
48#include "typedef.h"
49#include "basic_op.h"
50#include "oper_32b.h"
51#include "stdio.h"
52#include "grid100.tab"
53
54#define M   16
55#define NC  (M/2)
56
57/* local function */
58static __inline Word16 Chebps2(Word16 x, Word16 f[], Word32 n);
59
60void Az_isp(
61        Word16 a[],                           /* (i) Q12 : predictor coefficients                 */
62        Word16 isp[],                         /* (o) Q15 : Immittance spectral pairs              */
63        Word16 old_isp[]                      /* (i)     : old isp[] (in case not found M roots)  */
64       )
65{
66    Word32 i, j, nf, ip, order;
67    Word16 xlow, ylow, xhigh, yhigh, xmid, ymid, xint;
68    Word16 x, y, sign, exp;
69    Word16 *coef;
70    Word16 f1[NC + 1], f2[NC];
71    Word32 t0;
72    /*-------------------------------------------------------------*
73     * find the sum and diff polynomials F1(z) and F2(z)           *
74     *      F1(z) = [A(z) + z^M A(z^-1)]                           *
75     *      F2(z) = [A(z) - z^M A(z^-1)]/(1-z^-2)                  *
76     *                                                             *
77     * for (i=0; i<NC; i++)                                        *
78     * {                                                           *
79     *   f1[i] = a[i] + a[M-i];                                    *
80     *   f2[i] = a[i] - a[M-i];                                    *
81     * }                                                           *
82     * f1[NC] = 2.0*a[NC];                                         *
83     *                                                             *
84     * for (i=2; i<NC; i++)            Divide by (1-z^-2)          *
85     *   f2[i] += f2[i-2];                                         *
86     *-------------------------------------------------------------*/
87    for (i = 0; i < NC; i++)
88    {
89        t0 = a[i] << 15;
90        f1[i] = vo_round(t0 + (a[M - i] << 15));        /* =(a[i]+a[M-i])/2 */
91        f2[i] = vo_round(t0 - (a[M - i] << 15));        /* =(a[i]-a[M-i])/2 */
92    }
93    f1[NC] = a[NC];
94    for (i = 2; i < NC; i++)               /* Divide by (1-z^-2) */
95        f2[i] = add1(f2[i], f2[i - 2]);
96
97    /*---------------------------------------------------------------------*
98     * Find the ISPs (roots of F1(z) and F2(z) ) using the                 *
99     * Chebyshev polynomial evaluation.                                    *
100     * The roots of F1(z) and F2(z) are alternatively searched.            *
101     * We start by finding the first root of F1(z) then we switch          *
102     * to F2(z) then back to F1(z) and so on until all roots are found.    *
103     *                                                                     *
104     *  - Evaluate Chebyshev pol. at grid points and check for sign change.*
105     *  - If sign change track the root by subdividing the interval        *
106     *    2 times and ckecking sign change.                                *
107     *---------------------------------------------------------------------*/
108    nf = 0;                                  /* number of found frequencies */
109    ip = 0;                                  /* indicator for f1 or f2      */
110    coef = f1;
111    order = NC;
112    xlow = vogrid[0];
113    ylow = Chebps2(xlow, coef, order);
114    j = 0;
115    while ((nf < M - 1) && (j < GRID_POINTS))
116    {
117        j ++;
118        xhigh = xlow;
119        yhigh = ylow;
120        xlow = vogrid[j];
121        ylow = Chebps2(xlow, coef, order);
122        if ((ylow * yhigh) <= (Word32) 0)
123        {
124            /* divide 2 times the interval */
125            for (i = 0; i < 2; i++)
126            {
127                xmid = (xlow >> 1) + (xhigh >> 1);        /* xmid = (xlow + xhigh)/2 */
128                ymid = Chebps2(xmid, coef, order);
129                if ((ylow * ymid) <= (Word32) 0)
130                {
131                    yhigh = ymid;
132                    xhigh = xmid;
133                } else
134                {
135                    ylow = ymid;
136                    xlow = xmid;
137                }
138            }
139            /*-------------------------------------------------------------*
140             * Linear interpolation                                        *
141             *    xint = xlow - ylow*(xhigh-xlow)/(yhigh-ylow);            *
142             *-------------------------------------------------------------*/
143            x = xhigh - xlow;
144            y = yhigh - ylow;
145            if (y == 0)
146            {
147                xint = xlow;
148            } else
149            {
150                sign = y;
151                y = abs_s(y);
152                exp = norm_s(y);
153                y = y << exp;
154                y = div_s((Word16) 16383, y);
155                t0 = x * y;
156                t0 = (t0 >> (19 - exp));
157                y = vo_extract_l(t0);         /* y= (xhigh-xlow)/(yhigh-ylow) in Q11 */
158                if (sign < 0)
159                    y = -y;
160                t0 = ylow * y;      /* result in Q26 */
161                t0 = (t0 >> 10);        /* result in Q15 */
162                xint = vo_sub(xlow, vo_extract_l(t0));        /* xint = xlow - ylow*y */
163            }
164            isp[nf] = xint;
165            xlow = xint;
166            nf++;
167            if (ip == 0)
168            {
169                ip = 1;
170                coef = f2;
171                order = NC - 1;
172            } else
173            {
174                ip = 0;
175                coef = f1;
176                order = NC;
177            }
178            ylow = Chebps2(xlow, coef, order);
179        }
180    }
181    /* Check if M-1 roots found */
182    if(nf < M - 1)
183    {
184        for (i = 0; i < M; i++)
185        {
186            isp[i] = old_isp[i];
187        }
188    } else
189    {
190        isp[M - 1] = a[M] << 3;                      /* From Q12 to Q15 with saturation */
191    }
192    return;
193}
194
195/*--------------------------------------------------------------*
196* function  Chebps2:                                           *
197*           ~~~~~~~                                            *
198*    Evaluates the Chebishev polynomial series                 *
199*--------------------------------------------------------------*
200*                                                              *
201*  The polynomial order is                                     *
202*     n = M/2   (M is the prediction order)                    *
203*  The polynomial is given by                                  *
204*    C(x) = f(0)T_n(x) + f(1)T_n-1(x) + ... +f(n-1)T_1(x) + f(n)/2 *
205* Arguments:                                                   *
206*  x:     input value of evaluation; x = cos(frequency) in Q15 *
207*  f[]:   coefficients of the pol.                      in Q11 *
208*  n:     order of the pol.                                    *
209*                                                              *
210* The value of C(x) is returned. (Satured to +-1.99 in Q14)    *
211*                                                              *
212*--------------------------------------------------------------*/
213
214static __inline Word16 Chebps2(Word16 x, Word16 f[], Word32 n)
215{
216    Word32 i, cheb;
217    Word16 b0_h, b0_l, b1_h, b1_l, b2_h, b2_l;
218    Word32 t0;
219
220    /* Note: All computation are done in Q24. */
221
222    t0 = f[0] << 13;
223    b2_h = t0 >> 16;
224    b2_l = (t0 & 0xffff)>>1;
225
226    t0 = ((b2_h * x)<<1) + (((b2_l * x)>>15)<<1);
227    t0 <<= 1;
228    t0 += (f[1] << 13);                     /* + f[1] in Q24        */
229
230    b1_h = t0 >> 16;
231    b1_l = (t0 & 0xffff) >> 1;
232
233    for (i = 2; i < n; i++)
234    {
235        t0 = ((b1_h * x)<<1) + (((b1_l * x)>>15)<<1);
236
237        t0 += (b2_h * (-16384))<<1;
238        t0 += (f[i] << 12);
239        t0 <<= 1;
240        t0 -= (b2_l << 1);                  /* t0 = 2.0*x*b1 - b2 + f[i]; */
241
242        b0_h = t0 >> 16;
243        b0_l = (t0 & 0xffff) >> 1;
244
245        b2_l = b1_l;                         /* b2 = b1; */
246        b2_h = b1_h;
247        b1_l = b0_l;                         /* b1 = b0; */
248        b1_h = b0_h;
249    }
250
251    t0 = ((b1_h * x)<<1) + (((b1_l * x)>>15)<<1);
252    t0 += (b2_h * (-32768))<<1;             /* t0 = x*b1 - b2          */
253    t0 -= (b2_l << 1);
254    t0 += (f[n] << 12);                     /* t0 = x*b1 - b2 + f[i]/2 */
255
256    t0 = L_shl2(t0, 6);                     /* Q24 to Q30 with saturation */
257
258    cheb = extract_h(t0);                  /* Result in Q14              */
259
260    if (cheb == -32768)
261    {
262        cheb = -32767;                     /* to avoid saturation in Az_isp */
263    }
264    return (cheb);
265}
266
267
268
269