Coverage Report

Created: 2026-01-10 06:21

next uncovered line (L), next uncovered region (R), next uncovered branch (B)
/src/speex/libspeex/ltp.c
Line
Count
Source
1
/* Copyright (C) 2002-2006 Jean-Marc Valin
2
   File: ltp.c
3
   Long-Term Prediction functions
4
5
   Redistribution and use in source and binary forms, with or without
6
   modification, are permitted provided that the following conditions
7
   are met:
8
9
   - Redistributions of source code must retain the above copyright
10
   notice, this list of conditions and the following disclaimer.
11
12
   - Redistributions in binary form must reproduce the above copyright
13
   notice, this list of conditions and the following disclaimer in the
14
   documentation and/or other materials provided with the distribution.
15
16
   - Neither the name of the Xiph.org Foundation nor the names of its
17
   contributors may be used to endorse or promote products derived from
18
   this software without specific prior written permission.
19
20
   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21
   ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22
   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23
   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
24
   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
25
   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
26
   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
27
   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
28
   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
29
   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
30
   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31
*/
32
33
#ifdef HAVE_CONFIG_H
34
#include "config.h"
35
#endif
36
37
#include <math.h>
38
#include "ltp.h"
39
#include "stack_alloc.h"
40
#include "filters.h"
41
#include "math_approx.h"
42
#include "os_support.h"
43
44
#ifndef NULL
45
#define NULL 0
46
#endif
47
48
49
#ifdef _USE_SSE
50
#include "ltp_sse.h"
51
#elif defined (ARM4_ASM) || defined(ARM5E_ASM)
52
#include "ltp_arm4.h"
53
#elif defined (BFIN_ASM)
54
#include "ltp_bfin.h"
55
#endif
56
57
#ifndef OVERRIDE_INNER_PROD
58
spx_word32_t inner_prod(const spx_word16_t *x, const spx_word16_t *y, int len)
59
6.24M
{
60
6.24M
   spx_word32_t sum=0;
61
6.24M
   len >>= 2;
62
139M
   while(len--)
63
133M
   {
64
133M
      spx_word32_t part=0;
65
133M
      part = MAC16_16(part,*x++,*y++);
66
133M
      part = MAC16_16(part,*x++,*y++);
67
133M
      part = MAC16_16(part,*x++,*y++);
68
133M
      part = MAC16_16(part,*x++,*y++);
69
      /* HINT: If you had a 40-bit accumulator, you could shift only at the end */
70
133M
      sum = ADD32(sum,SHR32(part,6));
71
133M
   }
72
6.24M
   return sum;
73
6.24M
}
74
#endif
75
76
#ifndef DISABLE_ENCODER
77
78
#ifndef OVERRIDE_PITCH_XCORR
79
#if 0 /* HINT: Enable this for machines with enough registers (i.e. not x86) */
80
static void pitch_xcorr(const spx_word16_t *_x, const spx_word16_t *_y, spx_word32_t *corr, int len, int nb_pitch, char *stack)
81
{
82
   int i,j;
83
   for (i=0;i<nb_pitch;i+=4)
84
   {
85
      /* Compute correlation*/
86
      /*corr[nb_pitch-1-i]=inner_prod(x, _y+i, len);*/
87
      spx_word32_t sum1=0;
88
      spx_word32_t sum2=0;
89
      spx_word32_t sum3=0;
90
      spx_word32_t sum4=0;
91
      const spx_word16_t *y = _y+i;
92
      const spx_word16_t *x = _x;
93
      spx_word16_t y0, y1, y2, y3;
94
      /*y0=y[0];y1=y[1];y2=y[2];y3=y[3];*/
95
      y0=*y++;
96
      y1=*y++;
97
      y2=*y++;
98
      y3=*y++;
99
      for (j=0;j<len;j+=4)
100
      {
101
         spx_word32_t part1;
102
         spx_word32_t part2;
103
         spx_word32_t part3;
104
         spx_word32_t part4;
105
         part1 = MULT16_16(*x,y0);
106
         part2 = MULT16_16(*x,y1);
107
         part3 = MULT16_16(*x,y2);
108
         part4 = MULT16_16(*x,y3);
109
         x++;
110
         y0=*y++;
111
         part1 = MAC16_16(part1,*x,y1);
112
         part2 = MAC16_16(part2,*x,y2);
113
         part3 = MAC16_16(part3,*x,y3);
114
         part4 = MAC16_16(part4,*x,y0);
115
         x++;
116
         y1=*y++;
117
         part1 = MAC16_16(part1,*x,y2);
118
         part2 = MAC16_16(part2,*x,y3);
119
         part3 = MAC16_16(part3,*x,y0);
120
         part4 = MAC16_16(part4,*x,y1);
121
         x++;
122
         y2=*y++;
123
         part1 = MAC16_16(part1,*x,y3);
124
         part2 = MAC16_16(part2,*x,y0);
125
         part3 = MAC16_16(part3,*x,y1);
126
         part4 = MAC16_16(part4,*x,y2);
127
         x++;
128
         y3=*y++;
129
130
         sum1 = ADD32(sum1,SHR32(part1,6));
131
         sum2 = ADD32(sum2,SHR32(part2,6));
132
         sum3 = ADD32(sum3,SHR32(part3,6));
133
         sum4 = ADD32(sum4,SHR32(part4,6));
134
      }
135
      corr[nb_pitch-1-i]=sum1;
136
      corr[nb_pitch-2-i]=sum2;
137
      corr[nb_pitch-3-i]=sum3;
138
      corr[nb_pitch-4-i]=sum4;
139
   }
140
141
}
142
#else
143
static void pitch_xcorr(const spx_word16_t *_x, const spx_word16_t *_y, spx_word32_t *corr, int len, int nb_pitch, char *stack)
144
40.8k
{
145
40.8k
   int i;
146
5.00M
   for (i=0;i<nb_pitch;i++)
147
4.96M
   {
148
      /* Compute correlation*/
149
4.96M
      corr[nb_pitch-1-i]=inner_prod(_x, _y+i, len);
150
4.96M
   }
151
152
40.8k
}
153
#endif
154
#endif
155
156
#ifndef OVERRIDE_COMPUTE_PITCH_ERROR
157
static inline spx_word32_t compute_pitch_error(spx_word16_t *C, spx_word16_t *g, spx_word16_t pitch_control)
158
10.1M
{
159
10.1M
   spx_word32_t sum = 0;
160
10.1M
   sum = ADD32(sum,MULT16_16(MULT16_16_16(g[0],pitch_control),C[0]));
161
10.1M
   sum = ADD32(sum,MULT16_16(MULT16_16_16(g[1],pitch_control),C[1]));
162
10.1M
   sum = ADD32(sum,MULT16_16(MULT16_16_16(g[2],pitch_control),C[2]));
163
10.1M
   sum = SUB32(sum,MULT16_16(MULT16_16_16(g[0],g[1]),C[3]));
164
10.1M
   sum = SUB32(sum,MULT16_16(MULT16_16_16(g[2],g[1]),C[4]));
165
10.1M
   sum = SUB32(sum,MULT16_16(MULT16_16_16(g[2],g[0]),C[5]));
166
10.1M
   sum = SUB32(sum,MULT16_16(MULT16_16_16(g[0],g[0]),C[6]));
167
10.1M
   sum = SUB32(sum,MULT16_16(MULT16_16_16(g[1],g[1]),C[7]));
168
10.1M
   sum = SUB32(sum,MULT16_16(MULT16_16_16(g[2],g[2]),C[8]));
169
10.1M
   return sum;
170
10.1M
}
171
#endif
172
173
#ifndef OVERRIDE_OPEN_LOOP_NBEST_PITCH
174
void open_loop_nbest_pitch(spx_word16_t *sw, int start, int end, int len, int *pitch, spx_word16_t *gain, int N, char *stack)
175
40.8k
{
176
40.8k
   int i,j,k;
177
40.8k
   VARDECL(spx_word32_t *best_score);
178
40.8k
   VARDECL(spx_word32_t *best_ener);
179
40.8k
   spx_word32_t e0;
180
40.8k
   VARDECL(spx_word32_t *corr);
181
40.8k
#ifdef FIXED_POINT
182
   /* In fixed-point, we need only one (temporary) array of 32-bit values and two (corr16, ener16)
183
      arrays for (normalized) 16-bit values */
184
40.8k
   VARDECL(spx_word16_t *corr16);
185
40.8k
   VARDECL(spx_word16_t *ener16);
186
40.8k
   spx_word32_t *energy;
187
40.8k
   int cshift=0, eshift=0;
188
40.8k
   int scaledown = 0;
189
40.8k
   ALLOC(corr16, end-start+1, spx_word16_t);
190
40.8k
   ALLOC(ener16, end-start+1, spx_word16_t);
191
40.8k
   ALLOC(corr, end-start+1, spx_word32_t);
192
40.8k
   energy = corr;
193
#else
194
   /* In floating-point, we need to float arrays and no normalized copies */
195
   VARDECL(spx_word32_t *energy);
196
   spx_word16_t *corr16;
197
   spx_word16_t *ener16;
198
   ALLOC(energy, end-start+2, spx_word32_t);
199
   ALLOC(corr, end-start+1, spx_word32_t);
200
   corr16 = corr;
201
   ener16 = energy;
202
#endif
203
204
40.8k
   ALLOC(best_score, N, spx_word32_t);
205
40.8k
   ALLOC(best_ener, N, spx_word32_t);
206
264k
   for (i=0;i<N;i++)
207
223k
   {
208
223k
        best_score[i]=-1;
209
223k
        best_ener[i]=0;
210
223k
        pitch[i]=start;
211
223k
   }
212
213
40.8k
#ifdef FIXED_POINT
214
7.71M
   for (i=-end;i<len;i++)
215
7.68M
   {
216
7.68M
      if (ABS16(sw[i])>16383)
217
11.1k
      {
218
11.1k
         scaledown=1;
219
11.1k
         break;
220
11.1k
      }
221
7.68M
   }
222
   /* If the weighted input is close to saturation, then we scale it down */
223
40.8k
   if (scaledown)
224
11.1k
   {
225
2.60M
      for (i=-end;i<len;i++)
226
2.59M
      {
227
2.59M
         sw[i]=SHR16(sw[i],1);
228
2.59M
      }
229
11.1k
   }
230
40.8k
#endif
231
40.8k
   energy[0]=inner_prod(sw-start, sw-start, len);
232
40.8k
   e0=inner_prod(sw, sw, len);
233
4.96M
   for (i=start;i<end;i++)
234
4.92M
   {
235
      /* Update energy for next pitch*/
236
4.92M
      energy[i-start+1] = SUB32(ADD32(energy[i-start],SHR32(MULT16_16(sw[-i-1],sw[-i-1]),6)), SHR32(MULT16_16(sw[-i+len-1],sw[-i+len-1]),6));
237
4.92M
      if (energy[i-start+1] < 0)
238
0
         energy[i-start+1] = 0;
239
4.92M
   }
240
241
40.8k
#ifdef FIXED_POINT
242
40.8k
   eshift = normalize16(energy, ener16, 32766, end-start+1);
243
40.8k
#endif
244
245
   /* In fixed-point, this actually overrites the energy array (aliased to corr) */
246
40.8k
   pitch_xcorr(sw, sw-end, corr, len, end-start+1, stack);
247
248
40.8k
#ifdef FIXED_POINT
249
   /* Normalize to 180 so we can square it and it still fits in 16 bits */
250
40.8k
   cshift = normalize16(corr, corr16, 180, end-start+1);
251
   /* If we scaled weighted input down, we need to scale it up again (OK, so we've just lost the LSB, who cares?) */
252
40.8k
   if (scaledown)
253
11.1k
   {
254
2.60M
      for (i=-end;i<len;i++)
255
2.59M
      {
256
2.59M
         sw[i]=SHL16(sw[i],1);
257
2.59M
      }
258
11.1k
   }
259
40.8k
#endif
260
261
   /* Search for the best pitch prediction gain */
262
5.00M
   for (i=start;i<=end;i++)
263
4.96M
   {
264
4.96M
      spx_word16_t tmp = MULT16_16_16(corr16[i-start],corr16[i-start]);
265
      /* Instead of dividing the tmp by the energy, we multiply on the other side */
266
4.96M
      if (MULT16_16(tmp,best_ener[N-1])>MULT16_16(best_score[N-1],ADD16(1,ener16[i-start])))
267
752k
      {
268
         /* We can safely put it last and then check */
269
752k
         best_score[N-1]=tmp;
270
752k
         best_ener[N-1]=ener16[i-start]+1;
271
752k
         pitch[N-1]=i;
272
         /* Check if it comes in front of others */
273
2.43M
         for (j=0;j<N-1;j++)
274
2.31M
         {
275
2.31M
            if (MULT16_16(tmp,best_ener[j])>MULT16_16(best_score[j],ADD16(1,ener16[i-start])))
276
627k
            {
277
2.88M
               for (k=N-1;k>j;k--)
278
2.25M
               {
279
2.25M
                  best_score[k]=best_score[k-1];
280
2.25M
                  best_ener[k]=best_ener[k-1];
281
2.25M
                  pitch[k]=pitch[k-1];
282
2.25M
               }
283
627k
               best_score[j]=tmp;
284
627k
               best_ener[j]=ener16[i-start]+1;
285
627k
               pitch[j]=i;
286
627k
               break;
287
627k
            }
288
2.31M
         }
289
752k
      }
290
4.96M
   }
291
292
   /* Compute open-loop gain if necessary */
293
40.8k
   if (gain)
294
18.1k
   {
295
126k
      for (j=0;j<N;j++)
296
108k
      {
297
108k
         spx_word16_t g;
298
108k
         i=pitch[j];
299
108k
         g = DIV32(SHL32(EXTEND32(corr16[i-start]),cshift), 10+SHR32(MULT16_16(spx_sqrt(e0),spx_sqrt(SHL32(EXTEND32(ener16[i-start]),eshift))),6));
300
         /* FIXME: g = max(g,corr/energy) */
301
108k
         if (g<0)
302
44.7k
            g = 0;
303
108k
         gain[j]=g;
304
108k
      }
305
18.1k
   }
306
307
308
40.8k
}
309
#endif
310
311
#ifndef OVERRIDE_PITCH_GAIN_SEARCH_3TAP_VQ
312
static int pitch_gain_search_3tap_vq(
313
  const signed char *gain_cdbk,
314
  int                gain_cdbk_size,
315
  spx_word16_t      *C16,
316
  spx_word16_t       max_gain
317
)
318
119k
{
319
119k
  const signed char *ptr=gain_cdbk;
320
119k
  int                best_cdbk=0;
321
119k
  spx_word32_t       best_sum=-VERY_LARGE32;
322
119k
  spx_word32_t       sum=0;
323
119k
  spx_word16_t       g[3];
324
119k
  spx_word16_t       pitch_control=64;
325
119k
  spx_word16_t       gain_sum;
326
119k
  int                i;
327
328
10.2M
  for (i=0;i<gain_cdbk_size;i++) {
329
330
10.1M
    ptr = gain_cdbk+4*i;
331
10.1M
    g[0]=ADD16((spx_word16_t)ptr[0],32);
332
10.1M
    g[1]=ADD16((spx_word16_t)ptr[1],32);
333
10.1M
    g[2]=ADD16((spx_word16_t)ptr[2],32);
334
10.1M
    gain_sum = (spx_word16_t)ptr[3];
335
336
10.1M
    sum = compute_pitch_error(C16, g, pitch_control);
337
338
10.1M
    if (sum>best_sum && gain_sum<=max_gain) {
339
417k
      best_sum=sum;
340
417k
      best_cdbk=i;
341
417k
    }
342
10.1M
  }
343
344
119k
  return best_cdbk;
345
119k
}
346
#endif
347
348
/** Finds the best quantized 3-tap pitch predictor by analysis by synthesis */
349
static spx_word32_t pitch_gain_search_3tap(
350
const spx_word16_t target[],       /* Target vector */
351
const spx_coef_t ak[],          /* LPCs for this subframe */
352
const spx_coef_t awk1[],        /* Weighted LPCs #1 for this subframe */
353
const spx_coef_t awk2[],        /* Weighted LPCs #2 for this subframe */
354
spx_sig_t exc[],                /* Excitation */
355
const signed char *gain_cdbk,
356
int gain_cdbk_size,
357
int   pitch,                    /* Pitch value */
358
int   p,                        /* Number of LPC coeffs */
359
int   nsf,                      /* Number of samples in subframe */
360
SpeexBits *bits,
361
char *stack,
362
const spx_word16_t *exc2,
363
const spx_word16_t *r,
364
spx_word16_t *new_target,
365
int  *cdbk_index,
366
int plc_tuning,
367
spx_word32_t cumul_gain,
368
int scaledown
369
)
370
119k
{
371
119k
   int i,j;
372
119k
   VARDECL(spx_word16_t *tmp1);
373
119k
   VARDECL(spx_word16_t *e);
374
119k
   spx_word16_t *x[3];
375
119k
   spx_word32_t corr[3];
376
119k
   spx_word32_t A[3][3];
377
119k
   spx_word16_t gain[3];
378
119k
   spx_word32_t err;
379
119k
   spx_word16_t max_gain=128;
380
119k
   int          best_cdbk=0;
381
382
119k
   ALLOC(tmp1, 3*nsf, spx_word16_t);
383
119k
   ALLOC(e, nsf, spx_word16_t);
384
385
119k
   if (cumul_gain > 262144)
386
1.02k
      max_gain = 31;
387
388
119k
   x[0]=tmp1;
389
119k
   x[1]=tmp1+nsf;
390
119k
   x[2]=tmp1+2*nsf;
391
392
4.89M
   for (j=0;j<nsf;j++)
393
4.77M
      new_target[j] = target[j];
394
395
119k
   {
396
119k
      int bound;
397
119k
      VARDECL(spx_mem_t *mm);
398
119k
      int pp=pitch-1;
399
119k
      ALLOC(mm, p, spx_mem_t);
400
119k
      bound = nsf;
401
119k
      if (nsf-pp>0)
402
46.2k
         bound = pp;
403
4.17M
      for (j=0;j<bound;j++)
404
4.05M
         e[j]=exc2[j-pp];
405
119k
      bound = nsf;
406
119k
      if (nsf-pp-pitch>0)
407
16.3k
         bound = pp+pitch;
408
778k
      for (;j<bound;j++)
409
659k
         e[j]=exc2[j-pp-pitch];
410
188k
      for (;j<nsf;j++)
411
68.7k
         e[j]=0;
412
119k
#ifdef FIXED_POINT
413
      /* Scale target and excitation down if needed (avoiding overflow) */
414
119k
      if (scaledown)
415
26.6k
      {
416
1.09M
         for (j=0;j<nsf;j++)
417
1.06M
            e[j] = SHR16(e[j],1);
418
1.09M
         for (j=0;j<nsf;j++)
419
1.06M
            new_target[j] = SHR16(new_target[j],1);
420
26.6k
      }
421
119k
#endif
422
1.31M
      for (j=0;j<p;j++)
423
1.19M
         mm[j] = 0;
424
119k
      iir_mem16(e, ak, e, nsf, p, mm, stack);
425
1.31M
      for (j=0;j<p;j++)
426
1.19M
         mm[j] = 0;
427
119k
      filter10(e, awk1, awk2, e, nsf, mm, stack);
428
4.89M
      for (j=0;j<nsf;j++)
429
4.77M
         x[2][j] = e[j];
430
119k
   }
431
358k
   for (i=1;i>=0;i--)
432
238k
   {
433
238k
      spx_word16_t e0=exc2[-pitch-1+i];
434
238k
#ifdef FIXED_POINT
435
      /* Scale excitation down if needed (avoiding overflow) */
436
238k
      if (scaledown)
437
53.3k
         e0 = SHR16(e0,1);
438
238k
#endif
439
238k
      x[i][0]=MULT16_16_Q14(r[0], e0);
440
9.55M
      for (j=0;j<nsf-1;j++)
441
9.31M
         x[i][j+1]=ADD32(x[i+1][j],MULT16_16_P14(r[j+1], e0));
442
238k
   }
443
444
477k
   for (i=0;i<3;i++)
445
358k
      corr[i]=inner_prod(x[i],new_target,nsf);
446
477k
   for (i=0;i<3;i++)
447
1.07M
      for (j=0;j<=i;j++)
448
716k
         A[i][j]=A[j][i]=inner_prod(x[i],x[j],nsf);
449
450
119k
   {
451
119k
      spx_word32_t C[9];
452
119k
#ifdef FIXED_POINT
453
119k
      spx_word16_t C16[9];
454
#else
455
      spx_word16_t *C16=C;
456
#endif
457
119k
      C[0]=corr[2];
458
119k
      C[1]=corr[1];
459
119k
      C[2]=corr[0];
460
119k
      C[3]=A[1][2];
461
119k
      C[4]=A[0][1];
462
119k
      C[5]=A[0][2];
463
119k
      C[6]=A[2][2];
464
119k
      C[7]=A[1][1];
465
119k
      C[8]=A[0][0];
466
467
      /*plc_tuning *= 2;*/
468
119k
      if (plc_tuning<2)
469
0
         plc_tuning=2;
470
119k
      if (plc_tuning>30)
471
0
         plc_tuning=30;
472
119k
#ifdef FIXED_POINT
473
119k
      C[0] = SHL32(C[0],1);
474
119k
      C[1] = SHL32(C[1],1);
475
119k
      C[2] = SHL32(C[2],1);
476
119k
      C[3] = SHL32(C[3],1);
477
119k
      C[4] = SHL32(C[4],1);
478
119k
      C[5] = SHL32(C[5],1);
479
119k
      C[6] = MAC16_32_Q15(C[6],MULT16_16_16(plc_tuning,655),C[6]);
480
119k
      C[7] = MAC16_32_Q15(C[7],MULT16_16_16(plc_tuning,655),C[7]);
481
119k
      C[8] = MAC16_32_Q15(C[8],MULT16_16_16(plc_tuning,655),C[8]);
482
119k
      normalize16(C, C16, 32767, 9);
483
#else
484
      C[6]*=.5*(1+.02*plc_tuning);
485
      C[7]*=.5*(1+.02*plc_tuning);
486
      C[8]*=.5*(1+.02*plc_tuning);
487
#endif
488
489
119k
      best_cdbk = pitch_gain_search_3tap_vq(gain_cdbk, gain_cdbk_size, C16, max_gain);
490
491
119k
#ifdef FIXED_POINT
492
119k
      gain[0] = ADD16(32,(spx_word16_t)gain_cdbk[best_cdbk*4]);
493
119k
      gain[1] = ADD16(32,(spx_word16_t)gain_cdbk[best_cdbk*4+1]);
494
119k
      gain[2] = ADD16(32,(spx_word16_t)gain_cdbk[best_cdbk*4+2]);
495
      /*printf ("%d %d %d %d\n",gain[0],gain[1],gain[2], best_cdbk);*/
496
#else
497
      gain[0] = 0.015625*gain_cdbk[best_cdbk*4]  + .5;
498
      gain[1] = 0.015625*gain_cdbk[best_cdbk*4+1]+ .5;
499
      gain[2] = 0.015625*gain_cdbk[best_cdbk*4+2]+ .5;
500
#endif
501
119k
      *cdbk_index=best_cdbk;
502
119k
   }
503
504
119k
   SPEEX_MEMSET(exc, 0, nsf);
505
477k
   for (i=0;i<3;i++)
506
358k
   {
507
358k
      int j;
508
358k
      int tmp1, tmp3;
509
358k
      int pp=pitch+1-i;
510
358k
      tmp1=nsf;
511
358k
      if (tmp1>pp)
512
133k
         tmp1=pp;
513
12.6M
      for (j=0;j<tmp1;j++)
514
12.2M
         exc[j]=MAC16_16(exc[j],SHL16(gain[2-i],7),exc2[j-pp]);
515
358k
      tmp3=nsf;
516
358k
      if (tmp3>pp+pitch)
517
41.1k
         tmp3=pp+pitch;
518
2.24M
      for (j=tmp1;j<tmp3;j++)
519
1.88M
         exc[j]=MAC16_16(exc[j],SHL16(gain[2-i],7),exc2[j-pp-pitch]);
520
358k
   }
521
4.89M
   for (i=0;i<nsf;i++)
522
4.77M
   {
523
4.77M
      spx_word32_t tmp = ADD32(ADD32(MULT16_16(gain[0],x[2][i]),MULT16_16(gain[1],x[1][i])),
524
4.77M
                            MULT16_16(gain[2],x[0][i]));
525
4.77M
      new_target[i] = SUB16(new_target[i], EXTRACT16(PSHR32(tmp,6)));
526
4.77M
   }
527
119k
   err = inner_prod(new_target, new_target, nsf);
528
529
119k
   return err;
530
119k
}
531
532
/** Finds the best quantized 3-tap pitch predictor by analysis by synthesis */
533
int pitch_search_3tap(
534
spx_word16_t target[],                 /* Target vector */
535
spx_word16_t *sw,
536
spx_coef_t ak[],                     /* LPCs for this subframe */
537
spx_coef_t awk1[],                   /* Weighted LPCs #1 for this subframe */
538
spx_coef_t awk2[],                   /* Weighted LPCs #2 for this subframe */
539
spx_sig_t exc[],                    /* Excitation */
540
const void *par,
541
int   start,                    /* Smallest pitch value allowed */
542
int   end,                      /* Largest pitch value allowed */
543
spx_word16_t pitch_coef,               /* Voicing (pitch) coefficient */
544
int   p,                        /* Number of LPC coeffs */
545
int   nsf,                      /* Number of samples in subframe */
546
SpeexBits *bits,
547
char *stack,
548
spx_word16_t *exc2,
549
spx_word16_t *r,
550
int complexity,
551
int cdbk_offset,
552
int plc_tuning,
553
spx_word32_t *cumul_gain
554
)
555
28.8k
{
556
28.8k
   int i;
557
28.8k
   int cdbk_index, pitch=0, best_gain_index=0;
558
28.8k
   VARDECL(spx_sig_t *best_exc);
559
28.8k
   VARDECL(spx_word16_t *new_target);
560
28.8k
   VARDECL(spx_word16_t *best_target);
561
28.8k
   int best_pitch=0;
562
28.8k
   spx_word32_t err, best_err=-1;
563
28.8k
   int N;
564
28.8k
   const ltp_params *params;
565
28.8k
   const signed char *gain_cdbk;
566
28.8k
   int   gain_cdbk_size;
567
28.8k
   int scaledown=0;
568
569
28.8k
   VARDECL(int *nbest);
570
571
28.8k
   params = (const ltp_params*) par;
572
28.8k
   gain_cdbk_size = 1<<params->gain_bits;
573
28.8k
   gain_cdbk = params->gain_cdbk + 4*gain_cdbk_size*cdbk_offset;
574
575
28.8k
   N=complexity;
576
28.8k
   if (N>10)
577
0
      N=10;
578
28.8k
   if (N<1)
579
4.66k
      N=1;
580
581
28.8k
   ALLOC(nbest, N, int);
582
28.8k
   params = (const ltp_params*) par;
583
584
28.8k
   if (end<start)
585
1.61k
   {
586
1.61k
      speex_bits_pack(bits, 0, params->pitch_bits);
587
1.61k
      speex_bits_pack(bits, 0, params->gain_bits);
588
1.61k
      SPEEX_MEMSET(exc, 0, nsf);
589
1.61k
      return start;
590
1.61k
   }
591
592
27.2k
#ifdef FIXED_POINT
593
   /* Check if we need to scale everything down in the pitch search to avoid overflows */
594
977k
   for (i=0;i<nsf;i++)
595
954k
   {
596
954k
      if (ABS16(target[i])>16383)
597
4.59k
      {
598
4.59k
         scaledown=1;
599
4.59k
         break;
600
4.59k
      }
601
954k
   }
602
2.91M
   for (i=-end;i<0;i++)
603
2.89M
   {
604
2.89M
      if (ABS16(exc2[i])>16383)
605
4.03k
      {
606
4.03k
         scaledown=1;
607
4.03k
         break;
608
4.03k
      }
609
2.89M
   }
610
27.2k
#endif
611
27.2k
   if (N>end-start+1)
612
3.40k
      N=end-start+1;
613
27.2k
   if (end != start)
614
22.7k
      open_loop_nbest_pitch(sw, start, end, nsf, nbest, NULL, N, stack);
615
4.49k
   else
616
4.49k
      nbest[0] = start;
617
618
27.2k
   ALLOC(best_exc, nsf, spx_sig_t);
619
27.2k
   ALLOC(new_target, nsf, spx_word16_t);
620
27.2k
   ALLOC(best_target, nsf, spx_word16_t);
621
622
146k
   for (i=0;i<N;i++)
623
119k
   {
624
119k
      pitch=nbest[i];
625
119k
      SPEEX_MEMSET(exc, 0, nsf);
626
119k
      err=pitch_gain_search_3tap(target, ak, awk1, awk2, exc, gain_cdbk, gain_cdbk_size, pitch, p, nsf,
627
119k
                                 bits, stack, exc2, r, new_target, &cdbk_index, plc_tuning, *cumul_gain, scaledown);
628
119k
      if (err<best_err || best_err<0)
629
36.5k
      {
630
36.5k
         SPEEX_COPY(best_exc, exc, nsf);
631
36.5k
         SPEEX_COPY(best_target, new_target, nsf);
632
36.5k
         best_err=err;
633
36.5k
         best_pitch=pitch;
634
36.5k
         best_gain_index=cdbk_index;
635
36.5k
      }
636
119k
   }
637
   /*printf ("pitch: %d %d\n", best_pitch, best_gain_index);*/
638
27.2k
   speex_bits_pack(bits, best_pitch-start, params->pitch_bits);
639
27.2k
   speex_bits_pack(bits, best_gain_index, params->gain_bits);
640
27.2k
#ifdef FIXED_POINT
641
27.2k
   *cumul_gain = MULT16_32_Q13(SHL16(params->gain_cdbk[4*best_gain_index+3],8), MAX32(1024,*cumul_gain));
642
#else
643
   *cumul_gain = 0.03125*MAX32(1024,*cumul_gain)*params->gain_cdbk[4*best_gain_index+3];
644
#endif
645
   /*printf ("%f\n", cumul_gain);*/
646
   /*printf ("encode pitch: %d %d\n", best_pitch, best_gain_index);*/
647
27.2k
   SPEEX_COPY(exc, best_exc, nsf);
648
27.2k
   SPEEX_COPY(target, best_target, nsf);
649
27.2k
#ifdef FIXED_POINT
650
   /* Scale target back up if needed */
651
27.2k
   if (scaledown)
652
6.14k
   {
653
251k
      for (i=0;i<nsf;i++)
654
245k
         target[i]=SHL16(target[i],1);
655
6.14k
   }
656
27.2k
#endif
657
27.2k
   return pitch;
658
28.8k
}
659
#endif /* DISABLE_ENCODER */
660
661
#ifndef DISABLE_DECODER
662
void pitch_unquant_3tap(
663
spx_word16_t exc[],             /* Input excitation */
664
spx_word32_t exc_out[],         /* Output excitation */
665
int   start,                    /* Smallest pitch value allowed */
666
int   end,                      /* Largest pitch value allowed */
667
spx_word16_t pitch_coef,        /* Voicing (pitch) coefficient */
668
const void *par,
669
int   nsf,                      /* Number of samples in subframe */
670
int *pitch_val,
671
spx_word16_t *gain_val,
672
SpeexBits *bits,
673
char *stack,
674
int count_lost,
675
int subframe_offset,
676
spx_word16_t last_pitch_gain,
677
int cdbk_offset
678
)
679
0
{
680
0
   int i;
681
0
   int pitch;
682
0
   int gain_index;
683
0
   spx_word16_t gain[3];
684
0
   const signed char *gain_cdbk;
685
0
   int gain_cdbk_size;
686
0
   const ltp_params *params;
687
688
0
   params = (const ltp_params*) par;
689
0
   gain_cdbk_size = 1<<params->gain_bits;
690
0
   gain_cdbk = params->gain_cdbk + 4*gain_cdbk_size*cdbk_offset;
691
692
0
   pitch = speex_bits_unpack_unsigned(bits, params->pitch_bits);
693
0
   pitch += start;
694
0
   gain_index = speex_bits_unpack_unsigned(bits, params->gain_bits);
695
   /*printf ("decode pitch: %d %d\n", pitch, gain_index);*/
696
0
#ifdef FIXED_POINT
697
0
   gain[0] = ADD16(32,(spx_word16_t)gain_cdbk[gain_index*4]);
698
0
   gain[1] = ADD16(32,(spx_word16_t)gain_cdbk[gain_index*4+1]);
699
0
   gain[2] = ADD16(32,(spx_word16_t)gain_cdbk[gain_index*4+2]);
700
#else
701
   gain[0] = 0.015625*gain_cdbk[gain_index*4]+.5;
702
   gain[1] = 0.015625*gain_cdbk[gain_index*4+1]+.5;
703
   gain[2] = 0.015625*gain_cdbk[gain_index*4+2]+.5;
704
#endif
705
706
0
   if (count_lost && pitch > subframe_offset)
707
0
   {
708
0
      spx_word16_t gain_sum;
709
0
      if (1) {
710
0
#ifdef FIXED_POINT
711
0
         spx_word16_t tmp = count_lost < 4 ? last_pitch_gain : SHR16(last_pitch_gain,1);
712
0
         if (tmp>62)
713
0
            tmp=62;
714
#else
715
         spx_word16_t tmp = count_lost < 4 ? last_pitch_gain : 0.5 * last_pitch_gain;
716
         if (tmp>.95)
717
            tmp=.95;
718
#endif
719
0
         gain_sum = gain_3tap_to_1tap(gain);
720
721
0
         if (gain_sum > tmp)
722
0
         {
723
0
            spx_word16_t fact = DIV32_16(SHL32(EXTEND32(tmp),14),gain_sum);
724
0
            for (i=0;i<3;i++)
725
0
               gain[i]=MULT16_16_Q14(fact,gain[i]);
726
0
         }
727
728
0
      }
729
730
0
   }
731
732
0
   *pitch_val = pitch;
733
0
   gain_val[0]=gain[0];
734
0
   gain_val[1]=gain[1];
735
0
   gain_val[2]=gain[2];
736
0
   gain[0] = SHL16(gain[0],7);
737
0
   gain[1] = SHL16(gain[1],7);
738
0
   gain[2] = SHL16(gain[2],7);
739
0
   SPEEX_MEMSET(exc_out, 0, nsf);
740
0
   for (i=0;i<3;i++)
741
0
   {
742
0
      int j;
743
0
      int tmp1, tmp3;
744
0
      int pp=pitch+1-i;
745
0
      tmp1=nsf;
746
0
      if (tmp1>pp)
747
0
         tmp1=pp;
748
0
      for (j=0;j<tmp1;j++)
749
0
         exc_out[j]=MAC16_16(exc_out[j],gain[2-i],exc[j-pp]);
750
0
      tmp3=nsf;
751
0
      if (tmp3>pp+pitch)
752
0
         tmp3=pp+pitch;
753
0
      for (j=tmp1;j<tmp3;j++)
754
0
         exc_out[j]=MAC16_16(exc_out[j],gain[2-i],exc[j-pp-pitch]);
755
0
   }
756
   /*for (i=0;i<nsf;i++)
757
   exc[i]=PSHR32(exc32[i],13);*/
758
0
}
759
#endif /* DISABLE_DECODER */
760
761
#ifndef DISABLE_ENCODER
762
/** Forced pitch delay and gain */
763
int forced_pitch_quant(
764
spx_word16_t target[],                 /* Target vector */
765
spx_word16_t *sw,
766
spx_coef_t ak[],                     /* LPCs for this subframe */
767
spx_coef_t awk1[],                   /* Weighted LPCs #1 for this subframe */
768
spx_coef_t awk2[],                   /* Weighted LPCs #2 for this subframe */
769
spx_sig_t exc[],                    /* Excitation */
770
const void *par,
771
int   start,                    /* Smallest pitch value allowed */
772
int   end,                      /* Largest pitch value allowed */
773
spx_word16_t pitch_coef,               /* Voicing (pitch) coefficient */
774
int   p,                        /* Number of LPC coeffs */
775
int   nsf,                      /* Number of samples in subframe */
776
SpeexBits *bits,
777
char *stack,
778
spx_word16_t *exc2,
779
spx_word16_t *r,
780
int complexity,
781
int cdbk_offset,
782
int plc_tuning,
783
spx_word32_t *cumul_gain
784
)
785
44.8k
{
786
44.8k
   int i;
787
44.8k
   VARDECL(spx_word16_t *res);
788
44.8k
   ALLOC(res, nsf, spx_word16_t);
789
44.8k
#ifdef FIXED_POINT
790
44.8k
   if (pitch_coef>63)
791
464
      pitch_coef=63;
792
#else
793
   if (pitch_coef>.99)
794
      pitch_coef=.99;
795
#endif
796
1.50M
   for (i=0;i<nsf&&i<start;i++)
797
1.46M
   {
798
1.46M
      exc[i]=MULT16_16(SHL16(pitch_coef, 7),exc2[i-start]);
799
1.46M
   }
800
373k
   for (;i<nsf;i++)
801
328k
   {
802
328k
      exc[i]=MULT16_32_Q15(SHL16(pitch_coef, 9),exc[i-start]);
803
328k
   }
804
1.83M
   for (i=0;i<nsf;i++)
805
1.79M
      res[i] = EXTRACT16(PSHR32(exc[i], SIG_SHIFT-1));
806
44.8k
   syn_percep_zero16(res, ak, awk1, awk2, res, nsf, p, stack);
807
1.83M
   for (i=0;i<nsf;i++)
808
1.79M
      target[i]=EXTRACT16(SATURATE(SUB32(EXTEND32(target[i]),EXTEND32(res[i])),32700));
809
44.8k
   return start;
810
44.8k
}
811
#endif /* DISABLE_ENCODER */
812
813
#ifndef DISABLE_DECODER
814
/** Unquantize forced pitch delay and gain */
815
void forced_pitch_unquant(
816
spx_word16_t exc[],             /* Input excitation */
817
spx_word32_t exc_out[],         /* Output excitation */
818
int   start,                    /* Smallest pitch value allowed */
819
int   end,                      /* Largest pitch value allowed */
820
spx_word16_t pitch_coef,        /* Voicing (pitch) coefficient */
821
const void *par,
822
int   nsf,                      /* Number of samples in subframe */
823
int *pitch_val,
824
spx_word16_t *gain_val,
825
SpeexBits *bits,
826
char *stack,
827
int count_lost,
828
int subframe_offset,
829
spx_word16_t last_pitch_gain,
830
int cdbk_offset
831
)
832
0
{
833
0
   int i;
834
0
#ifdef FIXED_POINT
835
0
   if (pitch_coef>63)
836
0
      pitch_coef=63;
837
#else
838
   if (pitch_coef>.99)
839
      pitch_coef=.99;
840
#endif
841
0
   for (i=0;i<nsf;i++)
842
0
   {
843
0
      exc_out[i]=MULT16_16(exc[i-start],SHL16(pitch_coef,7));
844
0
      exc[i] = EXTRACT16(PSHR32(exc_out[i],13));
845
0
   }
846
0
   *pitch_val = start;
847
0
   gain_val[0]=gain_val[2]=0;
848
0
   gain_val[1] = pitch_coef;
849
0
}
850
#endif /* DISABLE_DECODER */