Home | History | Annotate | Download | only in src
      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: isp_az.c                                                  *
     19 *                                                                      *
     20 *      Description:Compute the LPC coefficients from isp (order=M)     *
     21 *                                                                      *
     22 ************************************************************************/
     23 
     24 #include "typedef.h"
     25 #include "basic_op.h"
     26 #include "oper_32b.h"
     27 #include "cnst.h"
     28 
     29 #define NC (M/2)
     30 #define NC16k (M16k/2)
     31 
     32 /* local function */
     33 
     34 static void Get_isp_pol(Word16 * isp, Word32 * f, Word16 n);
     35 static void Get_isp_pol_16kHz(Word16 * isp, Word32 * f, Word16 n);
     36 
     37 void Isp_Az(
     38 		Word16 isp[],                         /* (i) Q15 : Immittance spectral pairs            */
     39 		Word16 a[],                           /* (o) Q12 : predictor coefficients (order = M)   */
     40 		Word16 m,
     41 		Word16 adaptive_scaling               /* (i) 0   : adaptive scaling disabled */
     42 		                                      /*     1   : adaptive scaling enabled  */
     43 	   )
     44 {
     45 	Word32 i, j;
     46 	Word16 hi, lo;
     47 	Word32 f1[NC16k + 1], f2[NC16k];
     48 	Word16 nc;
     49 	Word32 t0;
     50 	Word16 q, q_sug;
     51 	Word32 tmax;
     52 
     53 	nc = (m >> 1);
     54 	if(nc > 8)
     55 	{
     56 		Get_isp_pol_16kHz(&isp[0], f1, nc);
     57 		for (i = 0; i <= nc; i++)
     58 		{
     59 			f1[i] = f1[i] << 2;
     60 		}
     61 	} else
     62 		Get_isp_pol(&isp[0], f1, nc);
     63 
     64 	if (nc > 8)
     65 	{
     66 		Get_isp_pol_16kHz(&isp[1], f2, (nc - 1));
     67 		for (i = 0; i <= nc - 1; i++)
     68 		{
     69 			f2[i] = f2[i] << 2;
     70 		}
     71 	} else
     72 		Get_isp_pol(&isp[1], f2, (nc - 1));
     73 
     74 	/*-----------------------------------------------------*
     75 	 *  Multiply F2(z) by (1 - z^-2)                       *
     76 	 *-----------------------------------------------------*/
     77 
     78 	for (i = (nc - 1); i > 1; i--)
     79 	{
     80 		f2[i] = vo_L_sub(f2[i], f2[i - 2]);          /* f2[i] -= f2[i-2]; */
     81 	}
     82 
     83 	/*----------------------------------------------------------*
     84 	 *  Scale F1(z) by (1+isp[m-1])  and  F2(z) by (1-isp[m-1]) *
     85 	 *----------------------------------------------------------*/
     86 
     87 	for (i = 0; i < nc; i++)
     88 	{
     89 		/* f1[i] *= (1.0 + isp[M-1]); */
     90 
     91 		hi = f1[i] >> 16;
     92 		lo = (f1[i] & 0xffff)>>1;
     93 
     94 		t0 = Mpy_32_16(hi, lo, isp[m - 1]);
     95 		f1[i] = vo_L_add(f1[i], t0);
     96 
     97 		/* f2[i] *= (1.0 - isp[M-1]); */
     98 
     99 		hi = f2[i] >> 16;
    100 		lo = (f2[i] & 0xffff)>>1;
    101 		t0 = Mpy_32_16(hi, lo, isp[m - 1]);
    102 		f2[i] = vo_L_sub(f2[i], t0);
    103 	}
    104 
    105 	/*-----------------------------------------------------*
    106 	 *  A(z) = (F1(z)+F2(z))/2                             *
    107 	 *  F1(z) is symmetric and F2(z) is antisymmetric      *
    108 	 *-----------------------------------------------------*/
    109 
    110 	/* a[0] = 1.0; */
    111 	a[0] = 4096;
    112 	tmax = 1;
    113 	for (i = 1, j = m - 1; i < nc; i++, j--)
    114 	{
    115 		/* a[i] = 0.5*(f1[i] + f2[i]); */
    116 
    117 		t0 = vo_L_add(f1[i], f2[i]);          /* f1[i] + f2[i]             */
    118 		tmax |= L_abs(t0);
    119 		a[i] = (Word16)(vo_L_shr_r(t0, 12)); /* from Q23 to Q12 and * 0.5 */
    120 
    121 		/* a[j] = 0.5*(f1[i] - f2[i]); */
    122 
    123 		t0 = vo_L_sub(f1[i], f2[i]);          /* f1[i] - f2[i]             */
    124 		tmax |= L_abs(t0);
    125 		a[j] = (Word16)(vo_L_shr_r(t0, 12)); /* from Q23 to Q12 and * 0.5 */
    126 	}
    127 
    128 	/* rescale data if overflow has occured and reprocess the loop */
    129 	if(adaptive_scaling == 1)
    130 		q = 4 - norm_l(tmax);        /* adaptive scaling enabled */
    131 	else
    132 		q = 0;                           /* adaptive scaling disabled */
    133 
    134 	if (q > 0)
    135 	{
    136 		q_sug = (12 + q);
    137 		for (i = 1, j = m - 1; i < nc; i++, j--)
    138 		{
    139 			/* a[i] = 0.5*(f1[i] + f2[i]); */
    140 			t0 = vo_L_add(f1[i], f2[i]);          /* f1[i] + f2[i]             */
    141 			a[i] = (Word16)(vo_L_shr_r(t0, q_sug)); /* from Q23 to Q12 and * 0.5 */
    142 
    143 			/* a[j] = 0.5*(f1[i] - f2[i]); */
    144 			t0 = vo_L_sub(f1[i], f2[i]);          /* f1[i] - f2[i]             */
    145 			a[j] = (Word16)(vo_L_shr_r(t0, q_sug)); /* from Q23 to Q12 and * 0.5 */
    146 		}
    147 		a[0] = shr(a[0], q);
    148 	}
    149 	else
    150 	{
    151 		q_sug = 12;
    152 		q     = 0;
    153 	}
    154 	/* a[NC] = 0.5*f1[NC]*(1.0 + isp[M-1]); */
    155 	hi = f1[nc] >> 16;
    156 	lo = (f1[nc] & 0xffff)>>1;
    157 	t0 = Mpy_32_16(hi, lo, isp[m - 1]);
    158 	t0 = vo_L_add(f1[nc], t0);
    159 	a[nc] = (Word16)(L_shr_r(t0, q_sug));    /* from Q23 to Q12 and * 0.5 */
    160 	/* a[m] = isp[m-1]; */
    161 
    162 	a[m] = vo_shr_r(isp[m - 1], (3 + q));           /* from Q15 to Q12          */
    163 	return;
    164 }
    165 
    166 /*-----------------------------------------------------------*
    167 * procedure Get_isp_pol:                                    *
    168 *           ~~~~~~~~~~~                                     *
    169 *   Find the polynomial F1(z) or F2(z) from the ISPs.       *
    170 * This is performed by expanding the product polynomials:   *
    171 *                                                           *
    172 * F1(z) =   product   ( 1 - 2 isp_i z^-1 + z^-2 )           *
    173 *         i=0,2,4,6,8                                       *
    174 * F2(z) =   product   ( 1 - 2 isp_i z^-1 + z^-2 )           *
    175 *         i=1,3,5,7                                         *
    176 *                                                           *
    177 * where isp_i are the ISPs in the cosine domain.            *
    178 *-----------------------------------------------------------*
    179 *                                                           *
    180 * Parameters:                                               *
    181 *  isp[]   : isp vector (cosine domaine)         in Q15     *
    182 *  f[]     : the coefficients of F1 or F2        in Q23     *
    183 *  n       : == NC for F1(z); == NC-1 for F2(z)             *
    184 *-----------------------------------------------------------*/
    185 
    186 static void Get_isp_pol(Word16 * isp, Word32 * f, Word16 n)
    187 {
    188 	Word16 hi, lo;
    189 	Word32 i, j, t0;
    190 	/* All computation in Q23 */
    191 
    192 	f[0] = vo_L_mult(4096, 1024);               /* f[0] = 1.0;        in Q23  */
    193 	f[1] = vo_L_mult(isp[0], -256);             /* f[1] = -2.0*isp[0] in Q23  */
    194 
    195 	f += 2;                                  /* Advance f pointer          */
    196 	isp += 2;                                /* Advance isp pointer        */
    197 	for (i = 2; i <= n; i++)
    198 	{
    199 		*f = f[-2];
    200 		for (j = 1; j < i; j++, f--)
    201 		{
    202 			hi = f[-1]>>16;
    203 			lo = (f[-1] & 0xffff)>>1;
    204 
    205 			t0 = Mpy_32_16(hi, lo, *isp);  /* t0 = f[-1] * isp    */
    206 			t0 = t0 << 1;
    207 			*f = vo_L_sub(*f, t0);              /* *f -= t0            */
    208 			*f = vo_L_add(*f, f[-2]);           /* *f += f[-2]         */
    209 		}
    210 		*f -= (*isp << 9);           /* *f -= isp<<8        */
    211 		f += i;                            /* Advance f pointer   */
    212 		isp += 2;                          /* Advance isp pointer */
    213 	}
    214 	return;
    215 }
    216 
    217 static void Get_isp_pol_16kHz(Word16 * isp, Word32 * f, Word16 n)
    218 {
    219 	Word16 hi, lo;
    220 	Word32 i, j, t0;
    221 
    222 	/* All computation in Q23 */
    223 	f[0] = L_mult(4096, 256);                /* f[0] = 1.0;        in Q23  */
    224 	f[1] = L_mult(isp[0], -64);              /* f[1] = -2.0*isp[0] in Q23  */
    225 
    226 	f += 2;                                  /* Advance f pointer          */
    227 	isp += 2;                                /* Advance isp pointer        */
    228 
    229 	for (i = 2; i <= n; i++)
    230 	{
    231 		*f = f[-2];
    232 		for (j = 1; j < i; j++, f--)
    233 		{
    234 			VO_L_Extract(f[-1], &hi, &lo);
    235 			t0 = Mpy_32_16(hi, lo, *isp);  /* t0 = f[-1] * isp    */
    236 			t0 = L_shl2(t0, 1);
    237 			*f = L_sub(*f, t0);              /* *f -= t0            */
    238 			*f = L_add(*f, f[-2]);           /* *f += f[-2]         */
    239 		}
    240 		*f = L_msu(*f, *isp, 64);            /* *f -= isp<<8        */
    241 		f += i;                            /* Advance f pointer   */
    242 		isp += 2;                          /* Advance isp pointer */
    243 	}
    244 	return;
    245 }
    246 
    247 
    248