Coverage Report

Created: 2025-07-23 06:50

/src/speex/libspeex/ltp.c
Line
Count
Source (jump to first uncovered line)
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.61M
{
60
6.61M
   spx_word32_t sum=0;
61
6.61M
   len >>= 2;
62
146M
   while(len--)
63
139M
   {
64
139M
      spx_word32_t part=0;
65
139M
      part = MAC16_16(part,*x++,*y++);
66
139M
      part = MAC16_16(part,*x++,*y++);
67
139M
      part = MAC16_16(part,*x++,*y++);
68
139M
      part = MAC16_16(part,*x++,*y++);
69
      /* HINT: If you had a 40-bit accumulator, you could shift only at the end */
70
139M
      sum = ADD32(sum,SHR32(part,6));
71
139M
   }
72
6.61M
   return sum;
73
6.61M
}
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
43.2k
{
145
43.2k
   int i;
146
5.30M
   for (i=0;i<nb_pitch;i++)
147
5.26M
   {
148
      /* Compute correlation*/
149
5.26M
      corr[nb_pitch-1-i]=inner_prod(_x, _y+i, len);
150
5.26M
   }
151
152
43.2k
}
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.5M
{
159
10.5M
   spx_word32_t sum = 0;
160
10.5M
   sum = ADD32(sum,MULT16_16(MULT16_16_16(g[0],pitch_control),C[0]));
161
10.5M
   sum = ADD32(sum,MULT16_16(MULT16_16_16(g[1],pitch_control),C[1]));
162
10.5M
   sum = ADD32(sum,MULT16_16(MULT16_16_16(g[2],pitch_control),C[2]));
163
10.5M
   sum = SUB32(sum,MULT16_16(MULT16_16_16(g[0],g[1]),C[3]));
164
10.5M
   sum = SUB32(sum,MULT16_16(MULT16_16_16(g[2],g[1]),C[4]));
165
10.5M
   sum = SUB32(sum,MULT16_16(MULT16_16_16(g[2],g[0]),C[5]));
166
10.5M
   sum = SUB32(sum,MULT16_16(MULT16_16_16(g[0],g[0]),C[6]));
167
10.5M
   sum = SUB32(sum,MULT16_16(MULT16_16_16(g[1],g[1]),C[7]));
168
10.5M
   sum = SUB32(sum,MULT16_16(MULT16_16_16(g[2],g[2]),C[8]));
169
10.5M
   return sum;
170
10.5M
}
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
43.2k
{
176
43.2k
   int i,j,k;
177
43.2k
   VARDECL(spx_word32_t *best_score);
178
43.2k
   VARDECL(spx_word32_t *best_ener);
179
43.2k
   spx_word32_t e0;
180
43.2k
   VARDECL(spx_word32_t *corr);
181
43.2k
#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
43.2k
   VARDECL(spx_word16_t *corr16);
185
43.2k
   VARDECL(spx_word16_t *ener16);
186
43.2k
   spx_word32_t *energy;
187
43.2k
   int cshift=0, eshift=0;
188
43.2k
   int scaledown = 0;
189
43.2k
   ALLOC(corr16, end-start+1, spx_word16_t);
190
43.2k
   ALLOC(ener16, end-start+1, spx_word16_t);
191
43.2k
   ALLOC(corr, end-start+1, spx_word32_t);
192
43.2k
   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
43.2k
   ALLOC(best_score, N, spx_word32_t);
205
43.2k
   ALLOC(best_ener, N, spx_word32_t);
206
278k
   for (i=0;i<N;i++)
207
235k
   {
208
235k
        best_score[i]=-1;
209
235k
        best_ener[i]=0;
210
235k
        pitch[i]=start;
211
235k
   }
212
213
43.2k
#ifdef FIXED_POINT
214
8.09M
   for (i=-end;i<len;i++)
215
8.06M
   {
216
8.06M
      if (ABS16(sw[i])>16383)
217
11.7k
      {
218
11.7k
         scaledown=1;
219
11.7k
         break;
220
11.7k
      }
221
8.06M
   }
222
   /* If the weighted input is close to saturation, then we scale it down */
223
43.2k
   if (scaledown)
224
11.7k
   {
225
2.73M
      for (i=-end;i<len;i++)
226
2.72M
      {
227
2.72M
         sw[i]=SHR16(sw[i],1);
228
2.72M
      }
229
11.7k
   }
230
43.2k
#endif
231
43.2k
   energy[0]=inner_prod(sw-start, sw-start, len);
232
43.2k
   e0=inner_prod(sw, sw, len);
233
5.26M
   for (i=start;i<end;i++)
234
5.22M
   {
235
      /* Update energy for next pitch*/
236
5.22M
      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
5.22M
      if (energy[i-start+1] < 0)
238
0
         energy[i-start+1] = 0;
239
5.22M
   }
240
241
43.2k
#ifdef FIXED_POINT
242
43.2k
   eshift = normalize16(energy, ener16, 32766, end-start+1);
243
43.2k
#endif
244
245
   /* In fixed-point, this actually overrites the energy array (aliased to corr) */
246
43.2k
   pitch_xcorr(sw, sw-end, corr, len, end-start+1, stack);
247
248
43.2k
#ifdef FIXED_POINT
249
   /* Normalize to 180 so we can square it and it still fits in 16 bits */
250
43.2k
   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
43.2k
   if (scaledown)
253
11.7k
   {
254
2.73M
      for (i=-end;i<len;i++)
255
2.72M
      {
256
2.72M
         sw[i]=SHL16(sw[i],1);
257
2.72M
      }
258
11.7k
   }
259
43.2k
#endif
260
261
   /* Search for the best pitch prediction gain */
262
5.30M
   for (i=start;i<=end;i++)
263
5.26M
   {
264
5.26M
      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
5.26M
      if (MULT16_16(tmp,best_ener[N-1])>MULT16_16(best_score[N-1],ADD16(1,ener16[i-start])))
267
808k
      {
268
         /* We can safely put it last and then check */
269
808k
         best_score[N-1]=tmp;
270
808k
         best_ener[N-1]=ener16[i-start]+1;
271
808k
         pitch[N-1]=i;
272
         /* Check if it comes in front of others */
273
2.59M
         for (j=0;j<N-1;j++)
274
2.45M
         {
275
2.45M
            if (MULT16_16(tmp,best_ener[j])>MULT16_16(best_score[j],ADD16(1,ener16[i-start])))
276
674k
            {
277
3.09M
               for (k=N-1;k>j;k--)
278
2.41M
               {
279
2.41M
                  best_score[k]=best_score[k-1];
280
2.41M
                  best_ener[k]=best_ener[k-1];
281
2.41M
                  pitch[k]=pitch[k-1];
282
2.41M
               }
283
674k
               best_score[j]=tmp;
284
674k
               best_ener[j]=ener16[i-start]+1;
285
674k
               pitch[j]=i;
286
674k
               break;
287
674k
            }
288
2.45M
         }
289
808k
      }
290
5.26M
   }
291
292
   /* Compute open-loop gain if necessary */
293
43.2k
   if (gain)
294
18.8k
   {
295
131k
      for (j=0;j<N;j++)
296
112k
      {
297
112k
         spx_word16_t g;
298
112k
         i=pitch[j];
299
112k
         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
112k
         if (g<0)
302
46.3k
            g = 0;
303
112k
         gain[j]=g;
304
112k
      }
305
18.8k
   }
306
307
308
43.2k
}
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
126k
{
319
126k
  const signed char *ptr=gain_cdbk;
320
126k
  int                best_cdbk=0;
321
126k
  spx_word32_t       best_sum=-VERY_LARGE32;
322
126k
  spx_word32_t       sum=0;
323
126k
  spx_word16_t       g[3];
324
126k
  spx_word16_t       pitch_control=64;
325
126k
  spx_word16_t       gain_sum;
326
126k
  int                i;
327
328
10.6M
  for (i=0;i<gain_cdbk_size;i++) {
329
330
10.5M
    ptr = gain_cdbk+4*i;
331
10.5M
    g[0]=ADD16((spx_word16_t)ptr[0],32);
332
10.5M
    g[1]=ADD16((spx_word16_t)ptr[1],32);
333
10.5M
    g[2]=ADD16((spx_word16_t)ptr[2],32);
334
10.5M
    gain_sum = (spx_word16_t)ptr[3];
335
336
10.5M
    sum = compute_pitch_error(C16, g, pitch_control);
337
338
10.5M
    if (sum>best_sum && gain_sum<=max_gain) {
339
440k
      best_sum=sum;
340
440k
      best_cdbk=i;
341
440k
    }
342
10.5M
  }
343
344
126k
  return best_cdbk;
345
126k
}
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
126k
{
371
126k
   int i,j;
372
126k
   VARDECL(spx_word16_t *tmp1);
373
126k
   VARDECL(spx_word16_t *e);
374
126k
   spx_word16_t *x[3];
375
126k
   spx_word32_t corr[3];
376
126k
   spx_word32_t A[3][3];
377
126k
   spx_word16_t gain[3];
378
126k
   spx_word32_t err;
379
126k
   spx_word16_t max_gain=128;
380
126k
   int          best_cdbk=0;
381
382
126k
   ALLOC(tmp1, 3*nsf, spx_word16_t);
383
126k
   ALLOC(e, nsf, spx_word16_t);
384
385
126k
   if (cumul_gain > 262144)
386
592
      max_gain = 31;
387
388
126k
   x[0]=tmp1;
389
126k
   x[1]=tmp1+nsf;
390
126k
   x[2]=tmp1+2*nsf;
391
392
5.18M
   for (j=0;j<nsf;j++)
393
5.05M
      new_target[j] = target[j];
394
395
126k
   {
396
126k
      int bound;
397
126k
      VARDECL(spx_mem_t *mm);
398
126k
      int pp=pitch-1;
399
126k
      ALLOC(mm, p, spx_mem_t);
400
126k
      bound = nsf;
401
126k
      if (nsf-pp>0)
402
48.0k
         bound = pp;
403
4.43M
      for (j=0;j<bound;j++)
404
4.31M
         e[j]=exc2[j-pp];
405
126k
      bound = nsf;
406
126k
      if (nsf-pp-pitch>0)
407
16.2k
         bound = pp+pitch;
408
802k
      for (;j<bound;j++)
409
676k
         e[j]=exc2[j-pp-pitch];
410
195k
      for (;j<nsf;j++)
411
68.9k
         e[j]=0;
412
126k
#ifdef FIXED_POINT
413
      /* Scale target and excitation down if needed (avoiding overflow) */
414
126k
      if (scaledown)
415
29.8k
      {
416
1.22M
         for (j=0;j<nsf;j++)
417
1.19M
            e[j] = SHR16(e[j],1);
418
1.22M
         for (j=0;j<nsf;j++)
419
1.19M
            new_target[j] = SHR16(new_target[j],1);
420
29.8k
      }
421
126k
#endif
422
1.39M
      for (j=0;j<p;j++)
423
1.26M
         mm[j] = 0;
424
126k
      iir_mem16(e, ak, e, nsf, p, mm, stack);
425
1.39M
      for (j=0;j<p;j++)
426
1.26M
         mm[j] = 0;
427
126k
      filter10(e, awk1, awk2, e, nsf, mm, stack);
428
5.18M
      for (j=0;j<nsf;j++)
429
5.05M
         x[2][j] = e[j];
430
126k
   }
431
379k
   for (i=1;i>=0;i--)
432
252k
   {
433
252k
      spx_word16_t e0=exc2[-pitch-1+i];
434
252k
#ifdef FIXED_POINT
435
      /* Scale excitation down if needed (avoiding overflow) */
436
252k
      if (scaledown)
437
59.6k
         e0 = SHR16(e0,1);
438
252k
#endif
439
252k
      x[i][0]=MULT16_16_Q14(r[0], e0);
440
10.1M
      for (j=0;j<nsf-1;j++)
441
9.85M
         x[i][j+1]=ADD32(x[i+1][j],MULT16_16_P14(r[j+1], e0));
442
252k
   }
443
444
505k
   for (i=0;i<3;i++)
445
379k
      corr[i]=inner_prod(x[i],new_target,nsf);
446
505k
   for (i=0;i<3;i++)
447
1.13M
      for (j=0;j<=i;j++)
448
758k
         A[i][j]=A[j][i]=inner_prod(x[i],x[j],nsf);
449
450
126k
   {
451
126k
      spx_word32_t C[9];
452
126k
#ifdef FIXED_POINT
453
126k
      spx_word16_t C16[9];
454
#else
455
      spx_word16_t *C16=C;
456
#endif
457
126k
      C[0]=corr[2];
458
126k
      C[1]=corr[1];
459
126k
      C[2]=corr[0];
460
126k
      C[3]=A[1][2];
461
126k
      C[4]=A[0][1];
462
126k
      C[5]=A[0][2];
463
126k
      C[6]=A[2][2];
464
126k
      C[7]=A[1][1];
465
126k
      C[8]=A[0][0];
466
467
      /*plc_tuning *= 2;*/
468
126k
      if (plc_tuning<2)
469
0
         plc_tuning=2;
470
126k
      if (plc_tuning>30)
471
0
         plc_tuning=30;
472
126k
#ifdef FIXED_POINT
473
126k
      C[0] = SHL32(C[0],1);
474
126k
      C[1] = SHL32(C[1],1);
475
126k
      C[2] = SHL32(C[2],1);
476
126k
      C[3] = SHL32(C[3],1);
477
126k
      C[4] = SHL32(C[4],1);
478
126k
      C[5] = SHL32(C[5],1);
479
126k
      C[6] = MAC16_32_Q15(C[6],MULT16_16_16(plc_tuning,655),C[6]);
480
126k
      C[7] = MAC16_32_Q15(C[7],MULT16_16_16(plc_tuning,655),C[7]);
481
126k
      C[8] = MAC16_32_Q15(C[8],MULT16_16_16(plc_tuning,655),C[8]);
482
126k
      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
126k
      best_cdbk = pitch_gain_search_3tap_vq(gain_cdbk, gain_cdbk_size, C16, max_gain);
490
491
126k
#ifdef FIXED_POINT
492
126k
      gain[0] = ADD16(32,(spx_word16_t)gain_cdbk[best_cdbk*4]);
493
126k
      gain[1] = ADD16(32,(spx_word16_t)gain_cdbk[best_cdbk*4+1]);
494
126k
      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
126k
      *cdbk_index=best_cdbk;
502
126k
   }
503
504
126k
   SPEEX_MEMSET(exc, 0, nsf);
505
505k
   for (i=0;i<3;i++)
506
379k
   {
507
379k
      int j;
508
379k
      int tmp1, tmp3;
509
379k
      int pp=pitch+1-i;
510
379k
      tmp1=nsf;
511
379k
      if (tmp1>pp)
512
138k
         tmp1=pp;
513
13.4M
      for (j=0;j<tmp1;j++)
514
13.0M
         exc[j]=MAC16_16(exc[j],SHL16(gain[2-i],7),exc2[j-pp]);
515
379k
      tmp3=nsf;
516
379k
      if (tmp3>pp+pitch)
517
41.4k
         tmp3=pp+pitch;
518
2.31M
      for (j=tmp1;j<tmp3;j++)
519
1.93M
         exc[j]=MAC16_16(exc[j],SHL16(gain[2-i],7),exc2[j-pp-pitch]);
520
379k
   }
521
5.18M
   for (i=0;i<nsf;i++)
522
5.05M
   {
523
5.05M
      spx_word32_t tmp = ADD32(ADD32(MULT16_16(gain[0],x[2][i]),MULT16_16(gain[1],x[1][i])),
524
5.05M
                            MULT16_16(gain[2],x[0][i]));
525
5.05M
      new_target[i] = SUB16(new_target[i], EXTRACT16(PSHR32(tmp,6)));
526
5.05M
   }
527
126k
   err = inner_prod(new_target, new_target, nsf);
528
529
126k
   return err;
530
126k
}
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
30.3k
{
556
30.3k
   int i;
557
30.3k
   int cdbk_index, pitch=0, best_gain_index=0;
558
30.3k
   VARDECL(spx_sig_t *best_exc);
559
30.3k
   VARDECL(spx_word16_t *new_target);
560
30.3k
   VARDECL(spx_word16_t *best_target);
561
30.3k
   int best_pitch=0;
562
30.3k
   spx_word32_t err, best_err=-1;
563
30.3k
   int N;
564
30.3k
   const ltp_params *params;
565
30.3k
   const signed char *gain_cdbk;
566
30.3k
   int   gain_cdbk_size;
567
30.3k
   int scaledown=0;
568
569
30.3k
   VARDECL(int *nbest);
570
571
30.3k
   params = (const ltp_params*) par;
572
30.3k
   gain_cdbk_size = 1<<params->gain_bits;
573
30.3k
   gain_cdbk = params->gain_cdbk + 4*gain_cdbk_size*cdbk_offset;
574
575
30.3k
   N=complexity;
576
30.3k
   if (N>10)
577
0
      N=10;
578
30.3k
   if (N<1)
579
4.84k
      N=1;
580
581
30.3k
   ALLOC(nbest, N, int);
582
30.3k
   params = (const ltp_params*) par;
583
584
30.3k
   if (end<start)
585
1.69k
   {
586
1.69k
      speex_bits_pack(bits, 0, params->pitch_bits);
587
1.69k
      speex_bits_pack(bits, 0, params->gain_bits);
588
1.69k
      SPEEX_MEMSET(exc, 0, nsf);
589
1.69k
      return start;
590
1.69k
   }
591
592
28.6k
#ifdef FIXED_POINT
593
   /* Check if we need to scale everything down in the pitch search to avoid overflows */
594
1.01M
   for (i=0;i<nsf;i++)
595
990k
   {
596
990k
      if (ABS16(target[i])>16383)
597
5.17k
      {
598
5.17k
         scaledown=1;
599
5.17k
         break;
600
5.17k
      }
601
990k
   }
602
3.03M
   for (i=-end;i<0;i++)
603
3.01M
   {
604
3.01M
      if (ABS16(exc2[i])>16383)
605
4.74k
      {
606
4.74k
         scaledown=1;
607
4.74k
         break;
608
4.74k
      }
609
3.01M
   }
610
28.6k
#endif
611
28.6k
   if (N>end-start+1)
612
3.36k
      N=end-start+1;
613
28.6k
   if (end != start)
614
24.4k
      open_loop_nbest_pitch(sw, start, end, nsf, nbest, NULL, N, stack);
615
4.25k
   else
616
4.25k
      nbest[0] = start;
617
618
28.6k
   ALLOC(best_exc, nsf, spx_sig_t);
619
28.6k
   ALLOC(new_target, nsf, spx_word16_t);
620
28.6k
   ALLOC(best_target, nsf, spx_word16_t);
621
622
155k
   for (i=0;i<N;i++)
623
126k
   {
624
126k
      pitch=nbest[i];
625
126k
      SPEEX_MEMSET(exc, 0, nsf);
626
126k
      err=pitch_gain_search_3tap(target, ak, awk1, awk2, exc, gain_cdbk, gain_cdbk_size, pitch, p, nsf,
627
126k
                                 bits, stack, exc2, r, new_target, &cdbk_index, plc_tuning, *cumul_gain, scaledown);
628
126k
      if (err<best_err || best_err<0)
629
38.7k
      {
630
38.7k
         SPEEX_COPY(best_exc, exc, nsf);
631
38.7k
         SPEEX_COPY(best_target, new_target, nsf);
632
38.7k
         best_err=err;
633
38.7k
         best_pitch=pitch;
634
38.7k
         best_gain_index=cdbk_index;
635
38.7k
      }
636
126k
   }
637
   /*printf ("pitch: %d %d\n", best_pitch, best_gain_index);*/
638
28.6k
   speex_bits_pack(bits, best_pitch-start, params->pitch_bits);
639
28.6k
   speex_bits_pack(bits, best_gain_index, params->gain_bits);
640
28.6k
#ifdef FIXED_POINT
641
28.6k
   *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
28.6k
   SPEEX_COPY(exc, best_exc, nsf);
648
28.6k
   SPEEX_COPY(target, best_target, nsf);
649
28.6k
#ifdef FIXED_POINT
650
   /* Scale target back up if needed */
651
28.6k
   if (scaledown)
652
6.83k
   {
653
280k
      for (i=0;i<nsf;i++)
654
273k
         target[i]=SHL16(target[i],1);
655
6.83k
   }
656
28.6k
#endif
657
28.6k
   return pitch;
658
30.3k
}
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
46.7k
{
786
46.7k
   int i;
787
46.7k
   VARDECL(spx_word16_t *res);
788
46.7k
   ALLOC(res, nsf, spx_word16_t);
789
46.7k
#ifdef FIXED_POINT
790
46.7k
   if (pitch_coef>63)
791
436
      pitch_coef=63;
792
#else
793
   if (pitch_coef>.99)
794
      pitch_coef=.99;
795
#endif
796
1.55M
   for (i=0;i<nsf&&i<start;i++)
797
1.51M
   {
798
1.51M
      exc[i]=MULT16_16(SHL16(pitch_coef, 7),exc2[i-start]);
799
1.51M
   }
800
402k
   for (;i<nsf;i++)
801
356k
   {
802
356k
      exc[i]=MULT16_32_Q15(SHL16(pitch_coef, 9),exc[i-start]);
803
356k
   }
804
1.91M
   for (i=0;i<nsf;i++)
805
1.86M
      res[i] = EXTRACT16(PSHR32(exc[i], SIG_SHIFT-1));
806
46.7k
   syn_percep_zero16(res, ak, awk1, awk2, res, nsf, p, stack);
807
1.91M
   for (i=0;i<nsf;i++)
808
1.86M
      target[i]=EXTRACT16(SATURATE(SUB32(EXTEND32(target[i]),EXTEND32(res[i])),32700));
809
46.7k
   return start;
810
46.7k
}
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 */