Coverage Report

Created: 2025-07-18 06:42

/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
{
60
   spx_word32_t sum=0;
61
   len >>= 2;
62
   while(len--)
63
   {
64
      spx_word32_t part=0;
65
      part = MAC16_16(part,*x++,*y++);
66
      part = MAC16_16(part,*x++,*y++);
67
      part = MAC16_16(part,*x++,*y++);
68
      part = MAC16_16(part,*x++,*y++);
69
      /* HINT: If you had a 40-bit accumulator, you could shift only at the end */
70
      sum = ADD32(sum,SHR32(part,6));
71
   }
72
   return sum;
73
}
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
{
145
   int i;
146
   for (i=0;i<nb_pitch;i++)
147
   {
148
      /* Compute correlation*/
149
      corr[nb_pitch-1-i]=inner_prod(_x, _y+i, len);
150
   }
151
152
}
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
4.32M
{
159
4.32M
   spx_word32_t sum = 0;
160
4.32M
   sum = ADD32(sum,MULT16_16(MULT16_16_16(g[0],pitch_control),C[0]));
161
4.32M
   sum = ADD32(sum,MULT16_16(MULT16_16_16(g[1],pitch_control),C[1]));
162
4.32M
   sum = ADD32(sum,MULT16_16(MULT16_16_16(g[2],pitch_control),C[2]));
163
4.32M
   sum = SUB32(sum,MULT16_16(MULT16_16_16(g[0],g[1]),C[3]));
164
4.32M
   sum = SUB32(sum,MULT16_16(MULT16_16_16(g[2],g[1]),C[4]));
165
4.32M
   sum = SUB32(sum,MULT16_16(MULT16_16_16(g[2],g[0]),C[5]));
166
4.32M
   sum = SUB32(sum,MULT16_16(MULT16_16_16(g[0],g[0]),C[6]));
167
4.32M
   sum = SUB32(sum,MULT16_16(MULT16_16_16(g[1],g[1]),C[7]));
168
4.32M
   sum = SUB32(sum,MULT16_16(MULT16_16_16(g[2],g[2]),C[8]));
169
4.32M
   return sum;
170
4.32M
}
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
16.2k
{
176
16.2k
   int i,j,k;
177
16.2k
   VARDECL(spx_word32_t *best_score);
178
16.2k
   VARDECL(spx_word32_t *best_ener);
179
16.2k
   spx_word32_t e0;
180
16.2k
   VARDECL(spx_word32_t *corr);
181
#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
   VARDECL(spx_word16_t *corr16);
185
   VARDECL(spx_word16_t *ener16);
186
   spx_word32_t *energy;
187
   int cshift=0, eshift=0;
188
   int scaledown = 0;
189
   ALLOC(corr16, end-start+1, spx_word16_t);
190
   ALLOC(ener16, end-start+1, spx_word16_t);
191
   ALLOC(corr, end-start+1, spx_word32_t);
192
   energy = corr;
193
#else
194
   /* In floating-point, we need to float arrays and no normalized copies */
195
16.2k
   VARDECL(spx_word32_t *energy);
196
16.2k
   spx_word16_t *corr16;
197
16.2k
   spx_word16_t *ener16;
198
16.2k
   ALLOC(energy, end-start+2, spx_word32_t);
199
16.2k
   ALLOC(corr, end-start+1, spx_word32_t);
200
16.2k
   corr16 = corr;
201
16.2k
   ener16 = energy;
202
16.2k
#endif
203
204
16.2k
   ALLOC(best_score, N, spx_word32_t);
205
16.2k
   ALLOC(best_ener, N, spx_word32_t);
206
112k
   for (i=0;i<N;i++)
207
95.9k
   {
208
95.9k
        best_score[i]=-1;
209
95.9k
        best_ener[i]=0;
210
95.9k
        pitch[i]=start;
211
95.9k
   }
212
213
#ifdef FIXED_POINT
214
   for (i=-end;i<len;i++)
215
   {
216
      if (ABS16(sw[i])>16383)
217
      {
218
         scaledown=1;
219
         break;
220
      }
221
   }
222
   /* If the weighted input is close to saturation, then we scale it down */
223
   if (scaledown)
224
   {
225
      for (i=-end;i<len;i++)
226
      {
227
         sw[i]=SHR16(sw[i],1);
228
      }
229
   }
230
#endif
231
16.2k
   energy[0]=inner_prod(sw-start, sw-start, len);
232
16.2k
   e0=inner_prod(sw, sw, len);
233
1.97M
   for (i=start;i<end;i++)
234
1.95M
   {
235
      /* Update energy for next pitch*/
236
1.95M
      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
1.95M
      if (energy[i-start+1] < 0)
238
6.92k
         energy[i-start+1] = 0;
239
1.95M
   }
240
241
#ifdef FIXED_POINT
242
   eshift = normalize16(energy, ener16, 32766, end-start+1);
243
#endif
244
245
   /* In fixed-point, this actually overrites the energy array (aliased to corr) */
246
16.2k
   pitch_xcorr(sw, sw-end, corr, len, end-start+1, stack);
247
248
#ifdef FIXED_POINT
249
   /* Normalize to 180 so we can square it and it still fits in 16 bits */
250
   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
   if (scaledown)
253
   {
254
      for (i=-end;i<len;i++)
255
      {
256
         sw[i]=SHL16(sw[i],1);
257
      }
258
   }
259
#endif
260
261
   /* Search for the best pitch prediction gain */
262
1.98M
   for (i=start;i<=end;i++)
263
1.97M
   {
264
1.97M
      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
1.97M
      if (MULT16_16(tmp,best_ener[N-1])>MULT16_16(best_score[N-1],ADD16(1,ener16[i-start])))
267
354k
      {
268
         /* We can safely put it last and then check */
269
354k
         best_score[N-1]=tmp;
270
354k
         best_ener[N-1]=ener16[i-start]+1;
271
354k
         pitch[N-1]=i;
272
         /* Check if it comes in front of others */
273
962k
         for (j=0;j<N-1;j++)
274
917k
         {
275
917k
            if (MULT16_16(tmp,best_ener[j])>MULT16_16(best_score[j],ADD16(1,ener16[i-start])))
276
309k
            {
277
1.64M
               for (k=N-1;k>j;k--)
278
1.33M
               {
279
1.33M
                  best_score[k]=best_score[k-1];
280
1.33M
                  best_ener[k]=best_ener[k-1];
281
1.33M
                  pitch[k]=pitch[k-1];
282
1.33M
               }
283
309k
               best_score[j]=tmp;
284
309k
               best_ener[j]=ener16[i-start]+1;
285
309k
               pitch[j]=i;
286
309k
               break;
287
309k
            }
288
917k
         }
289
354k
      }
290
1.97M
   }
291
292
   /* Compute open-loop gain if necessary */
293
16.2k
   if (gain)
294
7.98k
   {
295
55.8k
      for (j=0;j<N;j++)
296
47.9k
      {
297
47.9k
         spx_word16_t g;
298
47.9k
         i=pitch[j];
299
47.9k
         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
47.9k
         if (g<0)
302
13.9k
            g = 0;
303
47.9k
         gain[j]=g;
304
47.9k
      }
305
7.98k
   }
306
307
308
16.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
49.5k
{
319
49.5k
  const signed char *ptr=gain_cdbk;
320
49.5k
  int                best_cdbk=0;
321
49.5k
  spx_word32_t       best_sum=-VERY_LARGE32;
322
49.5k
  spx_word32_t       sum=0;
323
49.5k
  spx_word16_t       g[3];
324
49.5k
  spx_word16_t       pitch_control=64;
325
49.5k
  spx_word16_t       gain_sum;
326
49.5k
  int                i;
327
328
4.37M
  for (i=0;i<gain_cdbk_size;i++) {
329
330
4.32M
    ptr = gain_cdbk+4*i;
331
4.32M
    g[0]=ADD16((spx_word16_t)ptr[0],32);
332
4.32M
    g[1]=ADD16((spx_word16_t)ptr[1],32);
333
4.32M
    g[2]=ADD16((spx_word16_t)ptr[2],32);
334
4.32M
    gain_sum = (spx_word16_t)ptr[3];
335
336
4.32M
    sum = compute_pitch_error(C16, g, pitch_control);
337
338
4.32M
    if (sum>best_sum && gain_sum<=max_gain) {
339
113k
      best_sum=sum;
340
113k
      best_cdbk=i;
341
113k
    }
342
4.32M
  }
343
344
49.5k
  return best_cdbk;
345
49.5k
}
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
49.5k
{
371
49.5k
   int i,j;
372
49.5k
   VARDECL(spx_word16_t *tmp1);
373
49.5k
   VARDECL(spx_word16_t *e);
374
49.5k
   spx_word16_t *x[3];
375
49.5k
   spx_word32_t corr[3];
376
49.5k
   spx_word32_t A[3][3];
377
49.5k
   spx_word16_t gain[3];
378
49.5k
   spx_word32_t err;
379
49.5k
   spx_word16_t max_gain=128;
380
49.5k
   int          best_cdbk=0;
381
382
49.5k
   ALLOC(tmp1, 3*nsf, spx_word16_t);
383
49.5k
   ALLOC(e, nsf, spx_word16_t);
384
385
49.5k
   if (cumul_gain > 262144)
386
298
      max_gain = 31;
387
388
49.5k
   x[0]=tmp1;
389
49.5k
   x[1]=tmp1+nsf;
390
49.5k
   x[2]=tmp1+2*nsf;
391
392
2.03M
   for (j=0;j<nsf;j++)
393
1.98M
      new_target[j] = target[j];
394
395
49.5k
   {
396
49.5k
      int bound;
397
49.5k
      VARDECL(spx_mem_t *mm);
398
49.5k
      int pp=pitch-1;
399
49.5k
      ALLOC(mm, p, spx_mem_t);
400
49.5k
      bound = nsf;
401
49.5k
      if (nsf-pp>0)
402
24.2k
         bound = pp;
403
1.63M
      for (j=0;j<bound;j++)
404
1.58M
         e[j]=exc2[j-pp];
405
49.5k
      bound = nsf;
406
49.5k
      if (nsf-pp-pitch>0)
407
8.50k
         bound = pp+pitch;
408
405k
      for (;j<bound;j++)
409
355k
         e[j]=exc2[j-pp-pitch];
410
88.9k
      for (;j<nsf;j++)
411
39.4k
         e[j]=0;
412
#ifdef FIXED_POINT
413
      /* Scale target and excitation down if needed (avoiding overflow) */
414
      if (scaledown)
415
      {
416
         for (j=0;j<nsf;j++)
417
            e[j] = SHR16(e[j],1);
418
         for (j=0;j<nsf;j++)
419
            new_target[j] = SHR16(new_target[j],1);
420
      }
421
#endif
422
544k
      for (j=0;j<p;j++)
423
495k
         mm[j] = 0;
424
49.5k
      iir_mem16(e, ak, e, nsf, p, mm, stack);
425
544k
      for (j=0;j<p;j++)
426
495k
         mm[j] = 0;
427
49.5k
      filter10(e, awk1, awk2, e, nsf, mm, stack);
428
2.03M
      for (j=0;j<nsf;j++)
429
1.98M
         x[2][j] = e[j];
430
49.5k
   }
431
148k
   for (i=1;i>=0;i--)
432
99.0k
   {
433
99.0k
      spx_word16_t e0=exc2[-pitch-1+i];
434
#ifdef FIXED_POINT
435
      /* Scale excitation down if needed (avoiding overflow) */
436
      if (scaledown)
437
         e0 = SHR16(e0,1);
438
#endif
439
99.0k
      x[i][0]=MULT16_16_Q14(r[0], e0);
440
3.96M
      for (j=0;j<nsf-1;j++)
441
3.86M
         x[i][j+1]=ADD32(x[i+1][j],MULT16_16_P14(r[j+1], e0));
442
99.0k
   }
443
444
198k
   for (i=0;i<3;i++)
445
148k
      corr[i]=inner_prod(x[i],new_target,nsf);
446
198k
   for (i=0;i<3;i++)
447
445k
      for (j=0;j<=i;j++)
448
297k
         A[i][j]=A[j][i]=inner_prod(x[i],x[j],nsf);
449
450
49.5k
   {
451
49.5k
      spx_word32_t C[9];
452
#ifdef FIXED_POINT
453
      spx_word16_t C16[9];
454
#else
455
49.5k
      spx_word16_t *C16=C;
456
49.5k
#endif
457
49.5k
      C[0]=corr[2];
458
49.5k
      C[1]=corr[1];
459
49.5k
      C[2]=corr[0];
460
49.5k
      C[3]=A[1][2];
461
49.5k
      C[4]=A[0][1];
462
49.5k
      C[5]=A[0][2];
463
49.5k
      C[6]=A[2][2];
464
49.5k
      C[7]=A[1][1];
465
49.5k
      C[8]=A[0][0];
466
467
      /*plc_tuning *= 2;*/
468
49.5k
      if (plc_tuning<2)
469
0
         plc_tuning=2;
470
49.5k
      if (plc_tuning>30)
471
0
         plc_tuning=30;
472
#ifdef FIXED_POINT
473
      C[0] = SHL32(C[0],1);
474
      C[1] = SHL32(C[1],1);
475
      C[2] = SHL32(C[2],1);
476
      C[3] = SHL32(C[3],1);
477
      C[4] = SHL32(C[4],1);
478
      C[5] = SHL32(C[5],1);
479
      C[6] = MAC16_32_Q15(C[6],MULT16_16_16(plc_tuning,655),C[6]);
480
      C[7] = MAC16_32_Q15(C[7],MULT16_16_16(plc_tuning,655),C[7]);
481
      C[8] = MAC16_32_Q15(C[8],MULT16_16_16(plc_tuning,655),C[8]);
482
      normalize16(C, C16, 32767, 9);
483
#else
484
49.5k
      C[6]*=.5*(1+.02*plc_tuning);
485
49.5k
      C[7]*=.5*(1+.02*plc_tuning);
486
49.5k
      C[8]*=.5*(1+.02*plc_tuning);
487
49.5k
#endif
488
489
49.5k
      best_cdbk = pitch_gain_search_3tap_vq(gain_cdbk, gain_cdbk_size, C16, max_gain);
490
491
#ifdef FIXED_POINT
492
      gain[0] = ADD16(32,(spx_word16_t)gain_cdbk[best_cdbk*4]);
493
      gain[1] = ADD16(32,(spx_word16_t)gain_cdbk[best_cdbk*4+1]);
494
      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
49.5k
      gain[0] = 0.015625*gain_cdbk[best_cdbk*4]  + .5;
498
49.5k
      gain[1] = 0.015625*gain_cdbk[best_cdbk*4+1]+ .5;
499
49.5k
      gain[2] = 0.015625*gain_cdbk[best_cdbk*4+2]+ .5;
500
49.5k
#endif
501
49.5k
      *cdbk_index=best_cdbk;
502
49.5k
   }
503
504
49.5k
   SPEEX_MEMSET(exc, 0, nsf);
505
198k
   for (i=0;i<3;i++)
506
148k
   {
507
148k
      int j;
508
148k
      int tmp1, tmp3;
509
148k
      int pp=pitch+1-i;
510
148k
      tmp1=nsf;
511
148k
      if (tmp1>pp)
512
71.5k
         tmp1=pp;
513
4.97M
      for (j=0;j<tmp1;j++)
514
4.83M
         exc[j]=MAC16_16(exc[j],SHL16(gain[2-i],7),exc2[j-pp]);
515
148k
      tmp3=nsf;
516
148k
      if (tmp3>pp+pitch)
517
22.6k
         tmp3=pp+pitch;
518
1.16M
      for (j=tmp1;j<tmp3;j++)
519
1.01M
         exc[j]=MAC16_16(exc[j],SHL16(gain[2-i],7),exc2[j-pp-pitch]);
520
148k
   }
521
2.03M
   for (i=0;i<nsf;i++)
522
1.98M
   {
523
1.98M
      spx_word32_t tmp = ADD32(ADD32(MULT16_16(gain[0],x[2][i]),MULT16_16(gain[1],x[1][i])),
524
1.98M
                            MULT16_16(gain[2],x[0][i]));
525
1.98M
      new_target[i] = SUB16(new_target[i], EXTRACT16(PSHR32(tmp,6)));
526
1.98M
   }
527
49.5k
   err = inner_prod(new_target, new_target, nsf);
528
529
49.5k
   return err;
530
49.5k
}
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
10.4k
{
556
10.4k
   int i;
557
10.4k
   int cdbk_index, pitch=0, best_gain_index=0;
558
10.4k
   VARDECL(spx_sig_t *best_exc);
559
10.4k
   VARDECL(spx_word16_t *new_target);
560
10.4k
   VARDECL(spx_word16_t *best_target);
561
10.4k
   int best_pitch=0;
562
10.4k
   spx_word32_t err, best_err=-1;
563
10.4k
   int N;
564
10.4k
   const ltp_params *params;
565
10.4k
   const signed char *gain_cdbk;
566
10.4k
   int   gain_cdbk_size;
567
10.4k
   int scaledown=0;
568
569
10.4k
   VARDECL(int *nbest);
570
571
10.4k
   params = (const ltp_params*) par;
572
10.4k
   gain_cdbk_size = 1<<params->gain_bits;
573
10.4k
   gain_cdbk = params->gain_cdbk + 4*gain_cdbk_size*cdbk_offset;
574
575
10.4k
   N=complexity;
576
10.4k
   if (N>10)
577
0
      N=10;
578
10.4k
   if (N<1)
579
1.29k
      N=1;
580
581
10.4k
   ALLOC(nbest, N, int);
582
10.4k
   params = (const ltp_params*) par;
583
584
10.4k
   if (end<start)
585
743
   {
586
743
      speex_bits_pack(bits, 0, params->pitch_bits);
587
743
      speex_bits_pack(bits, 0, params->gain_bits);
588
743
      SPEEX_MEMSET(exc, 0, nsf);
589
743
      return start;
590
743
   }
591
592
#ifdef FIXED_POINT
593
   /* Check if we need to scale everything down in the pitch search to avoid overflows */
594
   for (i=0;i<nsf;i++)
595
   {
596
      if (ABS16(target[i])>16383)
597
      {
598
         scaledown=1;
599
         break;
600
      }
601
   }
602
   for (i=-end;i<0;i++)
603
   {
604
      if (ABS16(exc2[i])>16383)
605
      {
606
         scaledown=1;
607
         break;
608
      }
609
   }
610
#endif
611
9.73k
   if (N>end-start+1)
612
1.00k
      N=end-start+1;
613
9.73k
   if (end != start)
614
8.28k
      open_loop_nbest_pitch(sw, start, end, nsf, nbest, NULL, N, stack);
615
1.44k
   else
616
1.44k
      nbest[0] = start;
617
618
9.73k
   ALLOC(best_exc, nsf, spx_sig_t);
619
9.73k
   ALLOC(new_target, nsf, spx_word16_t);
620
9.73k
   ALLOC(best_target, nsf, spx_word16_t);
621
622
59.2k
   for (i=0;i<N;i++)
623
49.5k
   {
624
49.5k
      pitch=nbest[i];
625
49.5k
      SPEEX_MEMSET(exc, 0, nsf);
626
49.5k
      err=pitch_gain_search_3tap(target, ak, awk1, awk2, exc, gain_cdbk, gain_cdbk_size, pitch, p, nsf,
627
49.5k
                                 bits, stack, exc2, r, new_target, &cdbk_index, plc_tuning, *cumul_gain, scaledown);
628
49.5k
      if (err<best_err || best_err<0)
629
12.2k
      {
630
12.2k
         SPEEX_COPY(best_exc, exc, nsf);
631
12.2k
         SPEEX_COPY(best_target, new_target, nsf);
632
12.2k
         best_err=err;
633
12.2k
         best_pitch=pitch;
634
12.2k
         best_gain_index=cdbk_index;
635
12.2k
      }
636
49.5k
   }
637
   /*printf ("pitch: %d %d\n", best_pitch, best_gain_index);*/
638
9.73k
   speex_bits_pack(bits, best_pitch-start, params->pitch_bits);
639
9.73k
   speex_bits_pack(bits, best_gain_index, params->gain_bits);
640
#ifdef FIXED_POINT
641
   *cumul_gain = MULT16_32_Q13(SHL16(params->gain_cdbk[4*best_gain_index+3],8), MAX32(1024,*cumul_gain));
642
#else
643
9.73k
   *cumul_gain = 0.03125*MAX32(1024,*cumul_gain)*params->gain_cdbk[4*best_gain_index+3];
644
9.73k
#endif
645
   /*printf ("%f\n", cumul_gain);*/
646
   /*printf ("encode pitch: %d %d\n", best_pitch, best_gain_index);*/
647
9.73k
   SPEEX_COPY(exc, best_exc, nsf);
648
9.73k
   SPEEX_COPY(target, best_target, nsf);
649
#ifdef FIXED_POINT
650
   /* Scale target back up if needed */
651
   if (scaledown)
652
   {
653
      for (i=0;i<nsf;i++)
654
         target[i]=SHL16(target[i],1);
655
   }
656
#endif
657
9.73k
   return pitch;
658
10.4k
}
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
#ifdef FIXED_POINT
697
   gain[0] = ADD16(32,(spx_word16_t)gain_cdbk[gain_index*4]);
698
   gain[1] = ADD16(32,(spx_word16_t)gain_cdbk[gain_index*4+1]);
699
   gain[2] = ADD16(32,(spx_word16_t)gain_cdbk[gain_index*4+2]);
700
#else
701
0
   gain[0] = 0.015625*gain_cdbk[gain_index*4]+.5;
702
0
   gain[1] = 0.015625*gain_cdbk[gain_index*4+1]+.5;
703
0
   gain[2] = 0.015625*gain_cdbk[gain_index*4+2]+.5;
704
0
#endif
705
706
0
   if (count_lost && pitch > subframe_offset)
707
0
   {
708
0
      spx_word16_t gain_sum;
709
0
      if (1) {
710
#ifdef FIXED_POINT
711
         spx_word16_t tmp = count_lost < 4 ? last_pitch_gain : SHR16(last_pitch_gain,1);
712
         if (tmp>62)
713
            tmp=62;
714
#else
715
0
         spx_word16_t tmp = count_lost < 4 ? last_pitch_gain : 0.5 * last_pitch_gain;
716
0
         if (tmp>.95)
717
0
            tmp=.95;
718
0
#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
19.1k
{
786
19.1k
   int i;
787
19.1k
   VARDECL(spx_word16_t *res);
788
19.1k
   ALLOC(res, nsf, spx_word16_t);
789
#ifdef FIXED_POINT
790
   if (pitch_coef>63)
791
      pitch_coef=63;
792
#else
793
19.1k
   if (pitch_coef>.99)
794
156
      pitch_coef=.99;
795
19.1k
#endif
796
525k
   for (i=0;i<nsf&&i<start;i++)
797
506k
   {
798
506k
      exc[i]=MULT16_16(SHL16(pitch_coef, 7),exc2[i-start]);
799
506k
   }
800
276k
   for (;i<nsf;i++)
801
257k
   {
802
257k
      exc[i]=MULT16_32_Q15(SHL16(pitch_coef, 9),exc[i-start]);
803
257k
   }
804
783k
   for (i=0;i<nsf;i++)
805
764k
      res[i] = EXTRACT16(PSHR32(exc[i], SIG_SHIFT-1));
806
19.1k
   syn_percep_zero16(res, ak, awk1, awk2, res, nsf, p, stack);
807
783k
   for (i=0;i<nsf;i++)
808
764k
      target[i]=EXTRACT16(SATURATE(SUB32(EXTEND32(target[i]),EXTEND32(res[i])),32700));
809
19.1k
   return start;
810
19.1k
}
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
#ifdef FIXED_POINT
835
   if (pitch_coef>63)
836
      pitch_coef=63;
837
#else
838
0
   if (pitch_coef>.99)
839
0
      pitch_coef=.99;
840
0
#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 */