← Back to branch summary

~ubuntu-branches/ubuntu/hardy/speex/hardy-security

« back to all changes in this revision

Viewing changes to libspeex/lsp.c

  • Committer: Bazaar Package Importer
  • Author(s): Mark Purcell
  • Date: 2005-12-07 23:22:21 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20051207232221-nme7vf9m182p7dpe
Tags: 1.1.11.1-1
New upstream release

Show diffs side-by-side

added added

removed removed

Lines of Context:
libspeex/lsp.c
6
6
        AUTHOR......: David Rowe
7
7
        DATE CREATED: 24/2/93
8
8
 
9
 
Modified by Jean-Marc Valin
 
9
Heavily modified by Jean-Marc Valin (fixed-point, optimizations, 
 
10
                                     additional functions, ...)
10
11
 
11
12
   This file contains functions for converting Linear Prediction
12
13
   Coefficients (LPC) to Line Spectral Pair (LSP) and back. Note that the
76
77
      x2 = MULT16_16_P13(x,x);
77
78
      return ADD32(C1, MULT16_16_P13(x2, ADD32(C2, MULT16_16_P13(x2, ADD32(C3, MULT16_16_P13(C4, x2))))));
78
79
   } else {
79
 
      x = 25736-x;
 
80
      x = SUB16(25736,x);
80
81
      x2 = MULT16_16_P13(x,x);
81
82
      return SUB32(-C1, MULT16_16_P13(x2, ADD32(C2, MULT16_16_P13(x2, ADD32(C3, MULT16_16_P13(C4, x2))))));
82
83
      /*return SUB32(-C1, MULT16_16_Q13(x2, ADD32(C2, MULT16_16_Q13(C3, x2))));*/
87
88
#define FREQ_SCALE 16384
88
89
 
89
90
/*#define ANGLE2X(a) (32768*cos(((a)/8192.)))*/
90
 
#define ANGLE2X(a) (SHL(spx_cos(a),2))
 
91
#define ANGLE2X(a) (SHL16(spx_cos(a),2))
91
92
 
92
93
/*#define X2ANGLE(x) (acos(.00006103515625*(x))*LSP_SCALING)*/
93
94
#define X2ANGLE(x) (spx_acos(x))
98
99
#define C2 -0.49558072
99
100
#define C3 0.03679168*/
100
101
 
101
 
#define C1 0.9999932946
102
 
#define C2 -0.4999124376
103
 
#define C3 0.0414877472
104
 
#define C4 -0.0012712095
 
102
#define C1 0.9999932946f
 
103
#define C2 -0.4999124376f
 
104
#define C3 0.0414877472f
 
105
#define C4 -0.0012712095f
105
106
 
106
107
 
107
108
#define SPX_PI_2 1.5707963268
114
115
   } else {
115
116
      x = M_PI-x;
116
117
      x *= x;
117
 
      return -(C1 + x*(C2+x*(C3+C4*x)));
 
118
      return NEG16(C1 + x*(C2+x*(C3+C4*x)));
118
119
   }
119
120
}
120
121
#define FREQ_SCALE 1.
137
138
 
138
139
#ifdef FIXED_POINT
139
140
 
140
 
static spx_word32_t cheb_poly_eva(spx_word32_t *coef,spx_word16_t x,int m,char *stack)
 
141
static inline spx_word32_t cheb_poly_eva(spx_word32_t *coef,spx_word16_t x,int m,char *stack)
141
142
/*  float coef[]        coefficients of the polynomial to be evaluated  */
142
143
/*  float x             the point where polynomial is to be evaluated   */
143
144
/*  int m               order of the polynomial                         */
144
145
{
145
146
    int i;
146
 
    spx_word16_t *T;
 
147
    VARDECL(spx_word16_t *T);
147
148
    spx_word32_t sum;
148
149
    int m2=m>>1;
149
 
    spx_word16_t *coefn;
 
150
    VARDECL(spx_word16_t *coefn);
150
151
 
151
152
    /*Prevents overflows*/
152
153
    if (x>16383)
155
156
       x = -16383;
156
157
 
157
158
    /* Allocate memory for Chebyshev series formulation */
158
 
    T=PUSH(stack, m2+1, spx_word16_t);
159
 
    coefn=PUSH(stack, m2+1, spx_word16_t);
 
159
    ALLOC(T, m2+1, spx_word16_t);
 
160
    ALLOC(coefn, m2+1, spx_word16_t);
160
161
 
161
162
    for (i=0;i<m2+1;i++)
162
163
    {
171
172
 
172
173
    /* Evaluate Chebyshev series formulation using iterative approach  */
173
174
    /* Evaluate polynomial and return value also free memory space */
174
 
    sum = ADD32(coefn[m2], MULT16_16_P14(coefn[m2-1],x));
 
175
    sum = ADD32(EXTEND32(coefn[m2]), EXTEND32(MULT16_16_P14(coefn[m2-1],x)));
175
176
    /*x *= 2;*/
176
177
    for(i=2;i<=m2;i++)
177
178
    {
178
 
       T[i] = MULT16_16_Q13(x,T[i-1]) - T[i-2];
179
 
       sum = ADD32(sum, MULT16_16_P14(coefn[m2-i],T[i]));
 
179
       T[i] = SUB16(MULT16_16_Q13(x,T[i-1]), T[i-2]);
 
180
       sum = ADD32(sum, EXTEND32(MULT16_16_P14(coefn[m2-i],T[i])));
180
181
       /*printf ("%f ", sum);*/
181
182
    }
182
183
    
190
191
/*  int m               order of the polynomial                         */
191
192
{
192
193
    int i;
193
 
    float *T,sum;
 
194
    VARDECL(float *T);
 
195
    float sum;
194
196
    int m2=m>>1;
195
197
 
196
198
    /* Allocate memory for Chebyshev series formulation */
197
 
    T=PUSH(stack, m2+1, float);
 
199
    ALLOC(T, m2+1, float);
198
200
 
199
201
    /* Initialise values */
200
202
    T[0]=1;
245
247
    spx_word16_t temp_xr,xl,xr,xm=0;
246
248
    spx_word32_t psuml,psumr,psumm,temp_psumr/*,temp_qsumr*/;
247
249
    int i,j,m,flag,k;
248
 
    spx_word32_t *Q;                    /* ptrs for memory allocation           */
249
 
    spx_word32_t *P;
 
250
    VARDECL(spx_word32_t *Q);                   /* ptrs for memory allocation           */
 
251
    VARDECL(spx_word32_t *P);
250
252
    spx_word32_t *px;                   /* ptrs of respective P'(z) & Q'(z)     */
251
253
    spx_word32_t *qx;
252
254
    spx_word32_t *p;
259
261
    m = lpcrdr/2;               /* order of P'(z) & Q'(z) polynomials   */
260
262
 
261
263
    /* Allocate memory space for polynomials */
262
 
    Q = PUSH(stack, (m+1), spx_word32_t);
263
 
    P = PUSH(stack, (m+1), spx_word32_t);
 
264
    ALLOC(Q, (m+1), spx_word32_t);
 
265
    ALLOC(P, (m+1), spx_word32_t);
264
266
 
265
267
    /* determine P'(z)'s and Q'(z)'s coefficients where
266
268
      P'(z) = P(z)/(1 + z^(-1)) and Q'(z) = Q(z)/(1-z^(-1)) */
273
275
#ifdef FIXED_POINT
274
276
    *px++ = LPC_SCALING;
275
277
    *qx++ = LPC_SCALING;
276
 
    for(i=1;i<=m;i++){
277
 
        *px++ = (a[i]+a[lpcrdr+1-i]) - *p++;
278
 
        *qx++ = (a[i]-a[lpcrdr+1-i]) + *q++;
 
278
    for(i=0;i<m;i++){
 
279
       *px++ = SUB32(ADD32(EXTEND32(a[i]),EXTEND32(a[lpcrdr-i-1])), *p++);
 
280
       *qx++ = ADD32(SUB32(EXTEND32(a[i]),EXTEND32(a[lpcrdr-i-1])), *q++);
279
281
    }
280
282
    px = P;
281
283
    qx = Q;
285
287
          speex_warning_int("px", *px);
286
288
       if (fabs(*qx)>=32768)
287
289
       speex_warning_int("qx", *qx);*/
288
 
       *px = (2+*px)>>2;
289
 
       *qx = (2+*qx)>>2;
 
290
       *px = PSHR32(*px,2);
 
291
       *qx = PSHR32(*qx,2);
290
292
       px++;
291
293
       qx++;
292
294
    }
293
 
    P[m] = PSHR(P[m],3);
294
 
    Q[m] = PSHR(Q[m],3);
 
295
    /* The reason for this lies in the way cheb_poly_eva() is implemented for fixed-point */
 
296
    P[m] = PSHR32(P[m],3);
 
297
    Q[m] = PSHR32(Q[m],3);
295
298
#else
296
299
    *px++ = LPC_SCALING;
297
300
    *qx++ = LPC_SCALING;
298
 
    for(i=1;i<=m;i++){
299
 
        *px++ = (a[i]+a[lpcrdr+1-i]) - *p++;
300
 
        *qx++ = (a[i]-a[lpcrdr+1-i]) + *q++;
 
301
    for(i=0;i<m;i++){
 
302
       *px++ = (a[i]+a[lpcrdr-1-i]) - *p++;
 
303
       *qx++ = (a[i]-a[lpcrdr-1-i]) + *q++;
301
304
    }
302
305
    px = P;
303
306
    qx = Q;
304
307
    for(i=0;i<m;i++){
305
 
        *px = 2**px;
306
 
        *qx = 2**qx;
307
 
         px++;
308
 
         qx++;
 
308
       *px = 2**px;
 
309
       *qx = 2**qx;
 
310
       px++;
 
311
       qx++;
309
312
    }
310
313
#endif
311
314
 
331
334
           spx_word16_t dd;
332
335
           /* Modified by JMV to provide smaller steps around x=+-1 */
333
336
#ifdef FIXED_POINT
334
 
           dd = MULT16_16_Q15(delta,(FREQ_SCALE - MULT16_16_Q14(MULT16_16_Q14(xl,xl),14000)));
 
337
           dd = MULT16_16_Q15(delta,SUB16(FREQ_SCALE, MULT16_16_Q14(MULT16_16_Q14(xl,xl),14000)));
335
338
           if (psuml<512 && psuml>-512)
336
 
              dd = PSHR(dd,1);
 
339
              dd = PSHR16(dd,1);
337
340
#else
338
341
           dd=delta*(1-.9*xl*xl);
339
342
           if (fabs(psuml)<.2)
340
343
              dd *= .5;
341
344
#endif
342
 
           xr = xl - dd;                                /* interval spacing     */
 
345
           xr = SUB16(xl, dd);                          /* interval spacing     */
343
346
            psumr = cheb_poly_eva(pt,xr,lpcrdr,stack);/* poly(xl-delta_x)       */
344
347
            temp_psumr = psumr;
345
348
            temp_xr = xr;
360
363
                psumm=psuml;
361
364
                for(k=0;k<=nb;k++){
362
365
#ifdef FIXED_POINT
363
 
                    xm = ADD16(PSHR(xl,1),PSHR(xr,1));          /* bisect the interval  */
 
366
                    xm = ADD16(PSHR16(xl,1),PSHR16(xr,1));              /* bisect the interval  */
364
367
#else
365
368
                    xm = .5*(xl+xr);            /* bisect the interval  */
366
369
#endif
414
417
{
415
418
    int i,j;
416
419
    spx_word32_t xout1,xout2,xin1,xin2;
417
 
    spx_word32_t *Wp;
 
420
    VARDECL(spx_word32_t *Wp);
418
421
    spx_word32_t *pw,*n1,*n2,*n3,*n4=NULL;
419
 
    spx_word16_t *freqn;
 
422
    VARDECL(spx_word16_t *freqn);
420
423
    int m = lpcrdr>>1;
421
424
    
422
 
    freqn = PUSH(stack, lpcrdr, spx_word16_t);
 
425
    ALLOC(freqn, lpcrdr, spx_word16_t);
423
426
    for (i=0;i<lpcrdr;i++)
424
427
       freqn[i] = ANGLE2X(freq[i]);
425
428
 
426
 
    Wp = PUSH(stack, 4*m+2, spx_word32_t);
 
429
    ALLOC(Wp, 4*m+2, spx_word32_t);
427
430
    pw = Wp;
428
431
 
429
432
 
430
433
    /* initialise contents of array */
431
434
 
432
435
    for(i=0;i<=4*m+1;i++){              /* set contents of buffer to 0 */
433
 
        *pw++ = 0.0;
 
436
        *pw++ = 0;
434
437
    }
435
438
 
436
439
    /* Set pointers up */
444
447
      LSP coefficient */
445
448
 
446
449
    for(j=0;j<=lpcrdr;j++){
447
 
       int i2=0;
448
 
        for(i=0;i<m;i++,i2+=2){
449
 
            n1 = pw+(i*4);
 
450
       spx_word16_t *fr=freqn;
 
451
        for(i=0;i<m;i++){
 
452
            n1 = pw+(i<<2);
450
453
            n2 = n1 + 1;
451
454
            n3 = n2 + 1;
452
455
            n4 = n3 + 1;
453
 
            xout1 = ADD32(SUB32(xin1, MULT16_32_Q14(freqn[i2],*n1)), *n2);
454
 
            xout2 = ADD32(SUB32(xin2, MULT16_32_Q14(freqn[i2+1],*n3)), *n4);
 
456
            xout1 = ADD32(SUB32(xin1, MULT16_32_Q14(*fr,*n1)), *n2);
 
457
            fr++;
 
458
            xout2 = ADD32(SUB32(xin2, MULT16_32_Q14(*fr,*n3)), *n4);
 
459
            fr++;
455
460
            *n2 = *n1;
456
461
            *n4 = *n3;
457
462
            *n1 = xin1;
462
467
        xout1 = xin1 + *(n4+1);
463
468
        xout2 = xin2 - *(n4+2);
464
469
        /* FIXME: perhaps apply bandwidth expansion in case of overflow? */
465
 
        if (xout1 + xout2>256*32766)
466
 
           ak[j] = 32767;
467
 
        else if (xout1 + xout2 < -256*32767)
468
 
           ak[j] = -32768;
 
470
        if (j>0)
 
471
        {
 
472
        if (xout1 + xout2>SHL32(EXTEND32(32766),8))
 
473
           ak[j-1] = 32767;
 
474
        else if (xout1 + xout2 < -SHL32(EXTEND32(32766),8))
 
475
           ak[j-1] = -32767;
469
476
        else
470
 
           ak[j] = PSHR(ADD32(xout1,xout2),8);
 
477
           ak[j-1] = EXTRACT16(PSHR32(ADD32(xout1,xout2),8));
 
478
        } else {/*speex_warning_int("ak[0] = ", EXTRACT16(PSHR32(ADD32(xout1,xout2),8)));*/}
471
479
        *(n4+1) = xin1;
472
480
        *(n4+2) = xin2;
473
481
 
474
 
        xin1 = 0.0;
475
 
        xin2 = 0.0;
 
482
        xin1 = 0;
 
483
        xin2 = 0;
476
484
    }
477
485
}
478
486
#else
486
494
{
487
495
    int i,j;
488
496
    float xout1,xout2,xin1,xin2;
489
 
    float *Wp;
 
497
    VARDECL(float *Wp);
490
498
    float *pw,*n1,*n2,*n3,*n4=NULL;
491
 
    float *x_freq;
 
499
    VARDECL(float *x_freq);
492
500
    int m = lpcrdr>>1;
493
501
 
494
 
    Wp = PUSH(stack, 4*m+2, float);
 
502
    ALLOC(Wp, 4*m+2, float);
495
503
    pw = Wp;
496
504
 
497
505
    /* initialise contents of array */
506
514
    xin1 = 1.0;
507
515
    xin2 = 1.0;
508
516
 
509
 
    x_freq=PUSH(stack, lpcrdr, float);
 
517
    ALLOC(x_freq, lpcrdr, float);
510
518
    for (i=0;i<lpcrdr;i++)
511
519
       x_freq[i] = ANGLE2X(freq[i]);
512
520
 
532
540
        }
533
541
        xout1 = xin1 + *(n4+1);
534
542
        xout2 = xin2 - *(n4+2);
535
 
        ak[j] = (xout1 + xout2)*0.5f;
 
543
        if (j>0)
 
544
           ak[j-1] = (xout1 + xout2)*0.5f;
536
545
        *(n4+1) = xin1;
537
546
        *(n4+2) = xin2;
538
547
 
563
572
         lsp[i]=lsp[i-1]+m;
564
573
 
565
574
      if (lsp[i]>lsp[i+1]-m)
566
 
         lsp[i]= SHR(lsp[i],1) + SHR(lsp[i+1]-m,1);
 
575
         lsp[i]= SHR16(lsp[i],1) + SHR16(lsp[i+1]-m,1);
567
576
   }
568
577
}
569
578
 
571
580
void lsp_interpolate(spx_lsp_t *old_lsp, spx_lsp_t *new_lsp, spx_lsp_t *interp_lsp, int len, int subframe, int nb_subframes)
572
581
{
573
582
   int i;
574
 
   spx_word16_t tmp = DIV32_16(SHL(1 + subframe,14),nb_subframes);
 
583
   spx_word16_t tmp = DIV32_16(SHL32(EXTEND32(1 + subframe),14),nb_subframes);
575
584
   spx_word16_t tmp2 = 16384-tmp;
576
585
   for (i=0;i<len;i++)
577
586
   {