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: pitch_f4.c                                                *
19*                                                                      *
20*      Description: Find the closed loop pitch period with             *
21*	            1/4 subsample resolution.                          *
22*                                                                      *
23************************************************************************/
24
25#include "typedef.h"
26#include "basic_op.h"
27#include "math_op.h"
28#include "acelp.h"
29#include "cnst.h"
30
31#define UP_SAMP      4
32#define L_INTERPOL1  4
33
34#define UNUSED(x) (void)(x)
35
36/* Local functions */
37
38#ifdef ASM_OPT
39void Norm_corr_asm(
40		Word16 exc[],                         /* (i)     : excitation buffer                     */
41		Word16 xn[],                          /* (i)     : target vector                         */
42		Word16 h[],                           /* (i) Q15 : impulse response of synth/wgt filters */
43		Word16 L_subfr,
44		Word16 t_min,                         /* (i)     : minimum value of pitch lag.           */
45		Word16 t_max,                         /* (i)     : maximum value of pitch lag.           */
46		Word16 corr_norm[]                    /* (o) Q15 : normalized correlation                */
47		);
48#else
49static void Norm_Corr(
50		Word16 exc[],                         /* (i)     : excitation buffer                     */
51		Word16 xn[],                          /* (i)     : target vector                         */
52		Word16 h[],                           /* (i) Q15 : impulse response of synth/wgt filters */
53		Word16 L_subfr,
54		Word16 t_min,                         /* (i)     : minimum value of pitch lag.           */
55		Word16 t_max,                         /* (i)     : maximum value of pitch lag.           */
56		Word16 corr_norm[]                    /* (o) Q15 : normalized correlation                */
57		);
58#endif
59
60static Word16 Interpol_4(                  /* (o)  : interpolated value  */
61		Word16 * x,                           /* (i)  : input vector        */
62		Word32 frac                           /* (i)  : fraction (-4..+3)   */
63		);
64
65
66Word16 Pitch_fr4(                          /* (o)     : pitch period.                         */
67		Word16 exc[],                         /* (i)     : excitation buffer                     */
68		Word16 xn[],                          /* (i)     : target vector                         */
69		Word16 h[],                           /* (i) Q15 : impulse response of synth/wgt filters */
70		Word16 t0_min,                        /* (i)     : minimum value in the searched range.  */
71		Word16 t0_max,                        /* (i)     : maximum value in the searched range.  */
72		Word16 * pit_frac,                    /* (o)     : chosen fraction (0, 1, 2 or 3).       */
73		Word16 i_subfr,                       /* (i)     : indicator for first subframe.         */
74		Word16 t0_fr2,                        /* (i)     : minimum value for resolution 1/2      */
75		Word16 t0_fr1,                        /* (i)     : minimum value for resolution 1        */
76		Word16 L_subfr                        /* (i)     : Length of subframe                    */
77		)
78{
79	Word32 fraction, i;
80	Word16 t_min, t_max;
81	Word16 max, t0, step, temp;
82	Word16 *corr;
83	Word16 corr_v[40];                     /* Total length = t0_max-t0_min+1+2*L_inter */
84
85	/* Find interval to compute normalized correlation */
86
87	t_min = t0_min - L_INTERPOL1;
88	t_max = t0_max + L_INTERPOL1;
89	corr = &corr_v[-t_min];
90	/* Compute normalized correlation between target and filtered excitation */
91#ifdef ASM_OPT               /* asm optimization branch */
92    Norm_corr_asm(exc, xn, h, L_subfr, t_min, t_max, corr);
93#else
94	Norm_Corr(exc, xn, h, L_subfr, t_min, t_max, corr);
95#endif
96
97	/* Find integer pitch */
98
99	max = corr[t0_min];
100	t0 = t0_min;
101	for (i = t0_min + 1; i <= t0_max; i++)
102	{
103		if (corr[i] >= max)
104		{
105			max = corr[i];
106			t0 = i;
107		}
108	}
109	/* If first subframe and t0 >= t0_fr1, do not search fractionnal pitch */
110	if ((i_subfr == 0) && (t0 >= t0_fr1))
111	{
112		*pit_frac = 0;
113		return (t0);
114	}
115	/*------------------------------------------------------------------*
116	 * Search fractionnal pitch with 1/4 subsample resolution.          *
117	 * Test the fractions around t0 and choose the one which maximizes  *
118	 * the interpolated normalized correlation.                         *
119	 *------------------------------------------------------------------*/
120
121	step = 1;               /* 1/4 subsample resolution */
122	fraction = -3;
123	if ((t0_fr2 == PIT_MIN)||((i_subfr == 0) && (t0 >= t0_fr2)))
124	{
125		step = 2;              /* 1/2 subsample resolution */
126		fraction = -2;
127	}
128	if(t0 == t0_min)
129	{
130		fraction = 0;
131	}
132	max = Interpol_4(&corr[t0], fraction);
133
134	for (i = fraction + step; i <= 3; i += step)
135	{
136		temp = Interpol_4(&corr[t0], i);
137		if(temp > max)
138		{
139			max = temp;
140			fraction = i;
141		}
142	}
143	/* limit the fraction value in the interval [0,1,2,3] */
144	if (fraction < 0)
145	{
146		fraction += UP_SAMP;
147		t0 -= 1;
148	}
149	*pit_frac = fraction;
150	return (t0);
151}
152
153
154/***********************************************************************************
155* Function:  Norm_Corr()                                                            *
156*                                                                                   *
157* Description: Find the normalized correlation between the target vector and the    *
158* filtered past excitation.                                                         *
159* (correlation between target and filtered excitation divided by the                *
160*  square root of energy of target and filtered excitation).                        *
161************************************************************************************/
162#ifndef ASM_OPT
163static void Norm_Corr(
164		Word16 exc[],                         /* (i)     : excitation buffer                     */
165		Word16 xn[],                          /* (i)     : target vector                         */
166		Word16 h[],                           /* (i) Q15 : impulse response of synth/wgt filters */
167		Word16 L_subfr,
168		Word16 t_min,                         /* (i)     : minimum value of pitch lag.           */
169		Word16 t_max,                         /* (i)     : maximum value of pitch lag.           */
170		Word16 corr_norm[])                   /* (o) Q15 : normalized correlation                */
171{
172	Word32 i, k, t;
173	Word32 corr, exp_corr, norm, exp, scale;
174	Word16 exp_norm, excf[L_SUBFR], tmp;
175	Word32 L_tmp, L_tmp1, L_tmp2;
176        UNUSED(L_subfr);
177
178	/* compute the filtered excitation for the first delay t_min */
179	k = -t_min;
180
181#ifdef ASM_OPT              /* asm optimization branch */
182	Convolve_asm(&exc[k], h, excf, 64);
183#else
184	Convolve(&exc[k], h, excf, 64);
185#endif
186
187	/* Compute rounded down 1/sqrt(energy of xn[]) */
188	L_tmp = 0;
189	for (i = 0; i < 64; i+=4)
190	{
191		L_tmp += (xn[i] * xn[i]);
192		L_tmp += (xn[i+1] * xn[i+1]);
193		L_tmp += (xn[i+2] * xn[i+2]);
194		L_tmp += (xn[i+3] * xn[i+3]);
195	}
196
197	L_tmp = (L_tmp << 1) + 1;
198	exp = norm_l(L_tmp);
199	exp = (32 - exp);
200	//exp = exp + 2;                     /* energy of xn[] x 2 + rounded up     */
201	scale = -(exp >> 1);           /* (1<<scale) < 1/sqrt(energy rounded) */
202
203	/* loop for every possible period */
204
205	for (t = t_min; t <= t_max; t++)
206	{
207		/* Compute correlation between xn[] and excf[] */
208		L_tmp  = 0;
209		L_tmp1 = 0;
210		for (i = 0; i < 64; i+=4)
211		{
212			L_tmp  += (xn[i] * excf[i]);
213			L_tmp1 += (excf[i] * excf[i]);
214			L_tmp  += (xn[i+1] * excf[i+1]);
215			L_tmp1 += (excf[i+1] * excf[i+1]);
216			L_tmp  += (xn[i+2] * excf[i+2]);
217			L_tmp1 += (excf[i+2] * excf[i+2]);
218			L_tmp  += (xn[i+3] * excf[i+3]);
219			L_tmp1 += (excf[i+3] * excf[i+3]);
220		}
221
222		L_tmp = (L_tmp << 1) + 1;
223		L_tmp1 = (L_tmp1 << 1) + 1;
224
225		exp = norm_l(L_tmp);
226		L_tmp = (L_tmp << exp);
227		exp_corr = (30 - exp);
228		corr = extract_h(L_tmp);
229
230		exp = norm_l(L_tmp1);
231		L_tmp = (L_tmp1 << exp);
232		exp_norm = (30 - exp);
233
234		Isqrt_n(&L_tmp, &exp_norm);
235		norm = extract_h(L_tmp);
236
237		/* Normalize correlation = correlation * (1/sqrt(energy)) */
238
239		L_tmp = vo_L_mult(corr, norm);
240
241		L_tmp2 = exp_corr + exp_norm + scale;
242		if(L_tmp2 < 0)
243		{
244			L_tmp2 = -L_tmp2;
245			L_tmp = L_tmp >> L_tmp2;
246		}
247		else
248		{
249			L_tmp = L_tmp << L_tmp2;
250		}
251
252		corr_norm[t] = vo_round(L_tmp);
253		/* modify the filtered excitation excf[] for the next iteration */
254
255		if(t != t_max)
256		{
257			k = -(t + 1);
258			tmp = exc[k];
259			for (i = 63; i > 0; i--)
260			{
261				excf[i] = add1(vo_mult(tmp, h[i]), excf[i - 1]);
262			}
263			excf[0] = vo_mult(tmp, h[0]);
264		}
265	}
266	return;
267}
268
269#endif
270/************************************************************************************
271* Function: Interpol_4()                                                             *
272*                                                                                    *
273* Description: For interpolating the normalized correlation with 1/4 resolution.     *
274**************************************************************************************/
275
276/* 1/4 resolution interpolation filter (-3 dB at 0.791*fs/2) in Q14 */
277static Word16 inter4_1[4][8] =
278{
279	{-12, 420, -1732, 5429, 13418, -1242, 73, 32},
280	{-26, 455, -2142, 9910, 9910,  -2142, 455, -26},
281	{32,  73, -1242, 13418, 5429, -1732, 420, -12},
282	{206, -766, 1376, 14746, 1376, -766, 206, 0}
283};
284
285/*** Coefficients in floating point
286static float inter4_1[UP_SAMP*L_INTERPOL1+1] = {
2870.900000,
2880.818959,  0.604850,  0.331379,  0.083958,
289-0.075795, -0.130717, -0.105685, -0.046774,
2900.004467,  0.027789,  0.025642,  0.012571,
2910.001927, -0.001571, -0.000753,  0.000000};
292***/
293
294static Word16 Interpol_4(                  /* (o)  : interpolated value  */
295		Word16 * x,                           /* (i)  : input vector        */
296		Word32 frac                           /* (i)  : fraction (-4..+3)   */
297		)
298{
299	Word16 sum;
300	Word32  k, L_sum;
301	Word16 *ptr;
302
303	if (frac < 0)
304	{
305		frac += UP_SAMP;
306		x--;
307	}
308	x = x - L_INTERPOL1 + 1;
309	k = UP_SAMP - 1 - frac;
310	ptr = &(inter4_1[k][0]);
311
312	L_sum  = vo_mult32(x[0], (*ptr++));
313	L_sum += vo_mult32(x[1], (*ptr++));
314	L_sum += vo_mult32(x[2], (*ptr++));
315	L_sum += vo_mult32(x[3], (*ptr++));
316	L_sum += vo_mult32(x[4], (*ptr++));
317	L_sum += vo_mult32(x[5], (*ptr++));
318	L_sum += vo_mult32(x[6], (*ptr++));
319	L_sum += vo_mult32(x[7], (*ptr++));
320
321	sum = extract_h(L_add(L_shl2(L_sum, 2), 0x8000));
322	return (sum);
323}
324
325
326
327
328