3rdparty/iaxclient/lib/libspeex/lsp.c

   1 /*---------------------------------------------------------------------------*\
   2 Original copyright
   3         FILE........: AKSLSPD.C
   4         TYPE........: Turbo C
   5         COMPANY.....: Voicetronix
   6         AUTHOR......: David Rowe
   7         DATE CREATED: 24/2/93
   8
   9 Heavily modified by Jean-Marc Valin (fixed-point, optimizations,
  10                                      additional functions, ...)
  11
  12    This file contains functions for converting Linear Prediction
  13    Coefficients (LPC) to Line Spectral Pair (LSP) and back. Note that the
  14    LSP coefficients are not in radians format but in the x domain of the
  15    unit circle.
  16
  17    Speex License:
  18
  19    Redistribution and use in source and binary forms, with or without
  20    modification, are permitted provided that the following conditions
  21    are met:
  22
  23    - Redistributions of source code must retain the above copyright
  24    notice, this list of conditions and the following disclaimer.
  25
  26    - Redistributions in binary form must reproduce the above copyright
  27    notice, this list of conditions and the following disclaimer in the
  28    documentation and/or other materials provided with the distribution.
  29
  30    - Neither the name of the Xiph.org Foundation nor the names of its
  31    contributors may be used to endorse or promote products derived from
  32    this software without specific prior written permission.
  33
  34    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  35    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  36    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  37    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
  38    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
  39    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
  40    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
  41    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
  42    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
  43    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  44    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  45 */
  46
  47 #ifdef HAVE_CONFIG_H
  48 #include "config.h"
  49 #endif
  50 #ifdef _MSC_VER
  51 #include "winpoop.h"
  52 #endif
  53 #include <math.h>
  54 #include "lsp.h"
  55 #include "stack_alloc.h"
  56 #include "math_approx.h"
  57
  58 #ifndef M_PI
  59 #define M_PI           3.14159265358979323846  /* pi */
  60 #endif
  61
  62 #ifndef NULL
  63 #define NULL 0
  64 #endif
  65
  66 #ifdef FIXED_POINT
  67
  68 #define C1 8192
  69 #define C2 -4096
  70 #define C3 340
  71 #define C4 -10
  72
  73 static spx_word16_t spx_cos(spx_word16_t x)
  74 {
  75    spx_word16_t x2;
  76
  77    if (x<12868)
  78    {
  79       x2 = MULT16_16_P13(x,x);
  80       return ADD32(C1, MULT16_16_P13(x2, ADD32(C2, MULT16_16_P13(x2, ADD32(C3, MULT16_16_P13(C4, x2))))));
  81    } else {
  82       x = SUB16(25736,x);
  83       x2 = MULT16_16_P13(x,x);
  84       return SUB32(-C1, MULT16_16_P13(x2, ADD32(C2, MULT16_16_P13(x2, ADD32(C3, MULT16_16_P13(C4, x2))))));
  85       /*return SUB32(-C1, MULT16_16_Q13(x2, ADD32(C2, MULT16_16_Q13(C3, x2))));*/
  86    }
  87 }
  88
  89
  90 #define FREQ_SCALE 16384
  91
  92 /*#define ANGLE2X(a) (32768*cos(((a)/8192.)))*/
  93 #define ANGLE2X(a) (SHL16(spx_cos(a),2))
  94
  95 /*#define X2ANGLE(x) (acos(.00006103515625*(x))*LSP_SCALING)*/
  96 #define X2ANGLE(x) (spx_acos(x))
  97
  98 #else
  99
 100 /*#define C1 0.99940307
 101 #define C2 -0.49558072
 102 #define C3 0.03679168*/
 103
 104 #define C1 0.9999932946f
 105 #define C2 -0.4999124376f
 106 #define C3 0.0414877472f
 107 #define C4 -0.0012712095f
 108
 109
 110 #define SPX_PI_2 1.5707963268
 111 static inline spx_word16_t spx_cos(spx_word16_t x)
 112 {
 113    if (x<SPX_PI_2)
 114    {
 115       x *= x;
 116       return C1 + x*(C2+x*(C3+C4*x));
 117    } else {
 118       x = M_PI-x;
 119       x *= x;
 120       return NEG16(C1 + x*(C2+x*(C3+C4*x)));
 121    }
 122 }
 123 #define FREQ_SCALE 1.
 124 #define ANGLE2X(a) (spx_cos(a))
 125 #define X2ANGLE(x) (acos(x))
 126
 127 #endif
 128
 129
 130 /*---------------------------------------------------------------------------*\
 131
 132         FUNCTION....: cheb_poly_eva()
 133
 134         AUTHOR......: David Rowe
 135         DATE CREATED: 24/2/93
 136
 137     This function evaluates a series of Chebyshev polynomials
 138
 139 \*---------------------------------------------------------------------------*/
 140
 141 #ifdef FIXED_POINT
 142
 143 static inline spx_word32_t cheb_poly_eva(spx_word32_t *coef,spx_word16_t x,int m,char *stack)
 144 /*  float coef[]        coefficients of the polynomial to be evaluated  */
 145 /*  float x             the point where polynomial is to be evaluated   */
 146 /*  int m               order of the polynomial                         */
 147 {
 148     int i;
 149     VARDECL(spx_word16_t *T);
 150     spx_word32_t sum;
 151     int m2=m>>1;
 152     VARDECL(spx_word16_t *coefn);
 153
 154     /*Prevents overflows*/
 155     if (x>16383)
 156        x = 16383;
 157     if (x<-16383)
 158        x = -16383;
 159
 160     /* Allocate memory for Chebyshev series formulation */
 161     ALLOC(T, m2+1, spx_word16_t);
 162     ALLOC(coefn, m2+1, spx_word16_t);
 163
 164     for (i=0;i<m2+1;i++)
 165     {
 166        coefn[i] = coef[i];
 167        /*printf ("%f ", coef[i]);*/
 168     }
 169     /*printf ("\n");*/
 170
 171     /* Initialise values */
 172     T[0]=16384;
 173     T[1]=x;
 174
 175     /* Evaluate Chebyshev series formulation using iterative approach  */
 176     /* Evaluate polynomial and return value also free memory space */
 177     sum = ADD32(coefn[m2], MULT16_16_P14(coefn[m2-1],x));
 178     /*x *= 2;*/
 179     for(i=2;i<=m2;i++)
 180     {
 181        T[i] = SUB16(MULT16_16_Q13(x,T[i-1]), T[i-2]);
 182        sum = ADD32(sum, MULT16_16_P14(coefn[m2-i],T[i]));
 183        /*printf ("%f ", sum);*/
 184     }
 185
 186     /*printf ("\n");*/
 187     return sum;
 188 }
 189 #else
 190 static float cheb_poly_eva(spx_word32_t *coef,float x,int m,char *stack)
 191 /*  float coef[]        coefficients of the polynomial to be evaluated  */
 192 /*  float x             the point where polynomial is to be evaluated   */
 193 /*  int m               order of the polynomial                         */
 194 {
 195     int i;
 196     VARDECL(float *T);
 197     float sum;
 198     int m2=m>>1;
 199
 200     /* Allocate memory for Chebyshev series formulation */
 201     ALLOC(T, m2+1, float);
 202
 203     /* Initialise values */
 204     T[0]=1;
 205     T[1]=x;
 206
 207     /* Evaluate Chebyshev series formulation using iterative approach  */
 208     /* Evaluate polynomial and return value also free memory space */
 209     sum = coef[m2] + coef[m2-1]*x;
 210     x *= 2;
 211     for(i=2;i<=m2;i++)
 212     {
 213        T[i] = x*T[i-1] - T[i-2];
 214        sum += coef[m2-i] * T[i];
 215     }
 216
 217     return sum;
 218 }
 219 #endif
 220
 221 /*---------------------------------------------------------------------------*\
 222
 223         FUNCTION....: lpc_to_lsp()
 224
 225         AUTHOR......: David Rowe
 226         DATE CREATED: 24/2/93
 227
 228     This function converts LPC coefficients to LSP
 229     coefficients.
 230
 231 \*---------------------------------------------------------------------------*/
 232
 233 #ifdef FIXED_POINT
 234 #define SIGN_CHANGE(a,b) (((a)&0x70000000)^((b)&0x70000000)||(b==0))
 235 #else
 236 #define SIGN_CHANGE(a,b) (((a)*(b))<0.0)
 237 #endif
 238
 239
 240 int lpc_to_lsp (spx_coef_t *a,int lpcrdr,spx_lsp_t *freq,int nb,spx_word16_t delta, char *stack)
 241 /*  float *a                    lpc coefficients                        */
 242 /*  int lpcrdr                  order of LPC coefficients (10)          */
 243 /*  float *freq                 LSP frequencies in the x domain         */
 244 /*  int nb                      number of sub-intervals (4)             */
 245 /*  float delta                 grid spacing interval (0.02)            */
 246
 247
 248 {
 249     spx_word16_t temp_xr,xl,xr,xm=0;
 250     spx_word32_t psuml,psumr,psumm,temp_psumr/*,temp_qsumr*/;
 251     int i,j,m,flag,k;
 252     VARDECL(spx_word32_t *Q);                   /* ptrs for memory allocation           */
 253     VARDECL(spx_word32_t *P);
 254     spx_word32_t *px;                   /* ptrs of respective P'(z) & Q'(z)     */
 255     spx_word32_t *qx;
 256     spx_word32_t *p;
 257     spx_word32_t *q;
 258     spx_word32_t *pt;                   /* ptr used for cheb_poly_eval()
 259                                 whether P' or Q'                        */
 260     int roots=0;                /* DR 8/2/94: number of roots found     */
 261     flag = 1;                   /*  program is searching for a root when,
 262                                 1 else has found one                    */
 263     m = lpcrdr/2;               /* order of P'(z) & Q'(z) polynomials   */
 264
 265     /* Allocate memory space for polynomials */
 266     ALLOC(Q, (m+1), spx_word32_t);
 267     ALLOC(P, (m+1), spx_word32_t);
 268
 269     /* determine P'(z)'s and Q'(z)'s coefficients where
 270       P'(z) = P(z)/(1 + z^(-1)) and Q'(z) = Q(z)/(1-z^(-1)) */
 271
 272     px = P;                      /* initialise ptrs                     */
 273     qx = Q;
 274     p = px;
 275     q = qx;
 276
 277 #ifdef FIXED_POINT
 278     *px++ = LPC_SCALING;
 279     *qx++ = LPC_SCALING;
 280     for(i=1;i<=m;i++){
 281        *px++ = SUB32(ADD32(EXTEND32(a[i]),EXTEND32(a[lpcrdr+1-i])), *p++);
 282        *qx++ = ADD32(SUB32(EXTEND32(a[i]),EXTEND32(a[lpcrdr+1-i])), *q++);
 283     }
 284     px = P;
 285     qx = Q;
 286     for(i=0;i<m;i++)
 287     {
 288        /*if (fabs(*px)>=32768)
 289           speex_warning_int("px", *px);
 290        if (fabs(*qx)>=32768)
 291        speex_warning_int("qx", *qx);*/
 292        *px = PSHR32(*px,2);
 293        *qx = PSHR32(*qx,2);
 294        px++;
 295        qx++;
 296     }
 297     /* The reason for this lies in the way cheb_poly_eva() is implemented for fixed-point */
 298     P[m] = PSHR32(P[m],3);
 299     Q[m] = PSHR32(Q[m],3);
 300 #else
 301     *px++ = LPC_SCALING;
 302     *qx++ = LPC_SCALING;
 303     for(i=1;i<=m;i++){
 304        *px++ = (a[i]+a[lpcrdr+1-i]) - *p++;
 305        *qx++ = (a[i]-a[lpcrdr+1-i]) + *q++;
 306     }
 307     px = P;
 308     qx = Q;
 309     for(i=0;i<m;i++){
 310        *px = 2**px;
 311        *qx = 2**qx;
 312        px++;
 313        qx++;
 314     }
 315 #endif
 316
 317     px = P;                     /* re-initialise ptrs                   */
 318     qx = Q;
 319
 320     /* Search for a zero in P'(z) polynomial first and then alternate to Q'(z).
 321     Keep alternating between the two polynomials as each zero is found  */
 322
 323     xr = 0;                     /* initialise xr to zero                */
 324     xl = FREQ_SCALE;                    /* start at point xl = 1                */
 325
 326
 327     for(j=0;j<lpcrdr;j++){
 328         if(j&1)                 /* determines whether P' or Q' is eval. */
 329             pt = qx;
 330         else
 331             pt = px;
 332
 333         psuml = cheb_poly_eva(pt,xl,lpcrdr,stack);      /* evals poly. at xl    */
 334         flag = 1;
 335         while(flag && (xr >= -FREQ_SCALE)){
 336            spx_word16_t dd;
 337            /* Modified by JMV to provide smaller steps around x=+-1 */
 338 #ifdef FIXED_POINT
 339            dd = MULT16_16_Q15(delta,SUB16(FREQ_SCALE, MULT16_16_Q14(MULT16_16_Q14(xl,xl),14000)));
 340            if (psuml<512 && psuml>-512)
 341               dd = PSHR16(dd,1);
 342 #else
 343            dd=delta*(1-.9*xl*xl);
 344            if (fabs(psuml)<.2)
 345               dd *= .5;
 346 #endif
 347            xr = SUB16(xl, dd);                          /* interval spacing     */
 348             psumr = cheb_poly_eva(pt,xr,lpcrdr,stack);/* poly(xl-delta_x)       */
 349             temp_psumr = psumr;
 350             temp_xr = xr;
 351
 352     /* if no sign change increment xr and re-evaluate poly(xr). Repeat til
 353     sign change.
 354     if a sign change has occurred the interval is bisected and then
 355     checked again for a sign change which determines in which
 356     interval the zero lies in.
 357     If there is no sign change between poly(xm) and poly(xl) set interval
 358     between xm and xr else set interval between xl and xr and repeat till
 359     root is located within the specified limits                         */
 360
 361             if(SIGN_CHANGE(psumr,psuml))
 362             {
 363                 roots++;
 364
 365                 psumm=psuml;
 366                 for(k=0;k<=nb;k++){
 367 #ifdef FIXED_POINT
 368                     xm = ADD16(PSHR16(xl,1),PSHR16(xr,1));              /* bisect the interval  */
 369 #else
 370                     xm = .5*(xl+xr);            /* bisect the interval  */
 371 #endif
 372                     psumm=cheb_poly_eva(pt,xm,lpcrdr,stack);
 373                     /*if(psumm*psuml>0.)*/
 374                     if(!SIGN_CHANGE(psumm,psuml))
 375                     {
 376                         psuml=psumm;
 377                         xl=xm;
 378                     } else {
 379                         psumr=psumm;
 380                         xr=xm;
 381                     }
 382                 }
 383
 384                /* once zero is found, reset initial interval to xr      */
 385                freq[j] = X2ANGLE(xm);
 386                xl = xm;
 387                flag = 0;                /* reset flag for next search   */
 388             }
 389             else{
 390                 psuml=temp_psumr;
 391                 xl=temp_xr;
 392             }
 393         }
 394     }
 395     return(roots);
 396 }
 397
 398
 399 /*---------------------------------------------------------------------------*\
 400
 401         FUNCTION....: lsp_to_lpc()
 402
 403         AUTHOR......: David Rowe
 404         DATE CREATED: 24/2/93
 405
 406     lsp_to_lpc: This function converts LSP coefficients to LPC
 407     coefficients.
 408
 409 \*---------------------------------------------------------------------------*/
 410
 411 #ifdef FIXED_POINT
 412
 413 void lsp_to_lpc(spx_lsp_t *freq,spx_coef_t *ak,int lpcrdr, char *stack)
 414 /*  float *freq         array of LSP frequencies in the x domain        */
 415 /*  float *ak           array of LPC coefficients                       */
 416 /*  int lpcrdr          order of LPC coefficients                       */
 417
 418
 419 {
 420     int i,j;
 421     spx_word32_t xout1,xout2,xin1,xin2;
 422     VARDECL(spx_word32_t *Wp);
 423     spx_word32_t *pw,*n1,*n2,*n3,*n4=NULL;
 424     VARDECL(spx_word16_t *freqn);
 425     int m = lpcrdr>>1;
 426
 427     ALLOC(freqn, lpcrdr, spx_word16_t);
 428     for (i=0;i<lpcrdr;i++)
 429        freqn[i] = ANGLE2X(freq[i]);
 430
 431     ALLOC(Wp, 4*m+2, spx_word32_t);
 432     pw = Wp;
 433
 434
 435     /* initialise contents of array */
 436
 437     for(i=0;i<=4*m+1;i++){              /* set contents of buffer to 0 */
 438         *pw++ = 0;
 439     }
 440
 441     /* Set pointers up */
 442
 443     pw = Wp;
 444     xin1 = 1048576;
 445     xin2 = 1048576;
 446
 447     /* reconstruct P(z) and Q(z) by  cascading second order
 448       polynomials in form 1 - 2xz(-1) +z(-2), where x is the
 449       LSP coefficient */
 450
 451     for(j=0;j<=lpcrdr;j++){
 452        spx_word16_t *fr=freqn;
 453         for(i=0;i<m;i++){
 454             n1 = pw+(i<<2);
 455             n2 = n1 + 1;
 456             n3 = n2 + 1;
 457             n4 = n3 + 1;
 458             xout1 = ADD32(SUB32(xin1, MULT16_32_Q14(*fr,*n1)), *n2);
 459             fr++;
 460             xout2 = ADD32(SUB32(xin2, MULT16_32_Q14(*fr,*n3)), *n4);
 461             fr++;
 462             *n2 = *n1;
 463             *n4 = *n3;
 464             *n1 = xin1;
 465             *n3 = xin2;
 466             xin1 = xout1;
 467             xin2 = xout2;
 468         }
 469         xout1 = xin1 + *(n4+1);
 470         xout2 = xin2 - *(n4+2);
 471         /* FIXME: perhaps apply bandwidth expansion in case of overflow? */
 472         /*FIXME: Is it OK to have a long constant? */
 473         if (xout1 + xout2>SHL(32766,8))
 474            ak[j] = 32767;
 475         else if (xout1 + xout2 < -SHL(32766,8))
 476            ak[j] = -32767;
 477         else
 478            ak[j] = EXTRACT16(PSHR32(ADD32(xout1,xout2),8));
 479         *(n4+1) = xin1;
 480         *(n4+2) = xin2;
 481
 482         xin1 = 0;
 483         xin2 = 0;
 484     }
 485 }
 486 #else
 487
 488 void lsp_to_lpc(spx_lsp_t *freq,spx_coef_t *ak,int lpcrdr, char *stack)
 489 /*  float *freq         array of LSP frequencies in the x domain        */
 490 /*  float *ak           array of LPC coefficients                       */
 491 /*  int lpcrdr          order of LPC coefficients                       */
 492
 493
 494 {
 495     int i,j;
 496     float xout1,xout2,xin1,xin2;
 497     VARDECL(float *Wp);
 498     float *pw,*n1,*n2,*n3,*n4=NULL;
 499     VARDECL(float *x_freq);
 500     int m = lpcrdr>>1;
 501
 502     ALLOC(Wp, 4*m+2, float);
 503     pw = Wp;
 504
 505     /* initialise contents of array */
 506
 507     for(i=0;i<=4*m+1;i++){              /* set contents of buffer to 0 */
 508         *pw++ = 0.0;
 509     }
 510
 511     /* Set pointers up */
 512
 513     pw = Wp;
 514     xin1 = 1.0;
 515     xin2 = 1.0;
 516
 517     ALLOC(x_freq, lpcrdr, float);
 518     for (i=0;i<lpcrdr;i++)
 519        x_freq[i] = ANGLE2X(freq[i]);
 520
 521     /* reconstruct P(z) and Q(z) by  cascading second order
 522       polynomials in form 1 - 2xz(-1) +z(-2), where x is the
 523       LSP coefficient */
 524
 525     for(j=0;j<=lpcrdr;j++){
 526        int i2=0;
 527         for(i=0;i<m;i++,i2+=2){
 528             n1 = pw+(i*4);
 529             n2 = n1 + 1;
 530             n3 = n2 + 1;
 531             n4 = n3 + 1;
 532             xout1 = xin1 - 2.f*x_freq[i2] * *n1 + *n2;
 533             xout2 = xin2 - 2.f*x_freq[i2+1] * *n3 + *n4;
 534             *n2 = *n1;
 535             *n4 = *n3;
 536             *n1 = xin1;
 537             *n3 = xin2;
 538             xin1 = xout1;
 539             xin2 = xout2;
 540         }
 541         xout1 = xin1 + *(n4+1);
 542         xout2 = xin2 - *(n4+2);
 543         ak[j] = (xout1 + xout2)*0.5f;
 544         *(n4+1) = xin1;
 545         *(n4+2) = xin2;
 546
 547         xin1 = 0.0;
 548         xin2 = 0.0;
 549     }
 550
 551 }
 552 #endif
 553
 554
 555 #ifdef FIXED_POINT
 556
 557 /*Makes sure the LSPs are stable*/
 558 void lsp_enforce_margin(spx_lsp_t *lsp, int len, spx_word16_t margin)
 559 {
 560    int i;
 561    spx_word16_t m = margin;
 562    spx_word16_t m2 = 25736-margin;
 563
 564    if (lsp[0]<m)
 565       lsp[0]=m;
 566    if (lsp[len-1]>m2)
 567       lsp[len-1]=m2;
 568    for (i=1;i<len-1;i++)
 569    {
 570       if (lsp[i]<lsp[i-1]+m)
 571          lsp[i]=lsp[i-1]+m;
 572
 573       if (lsp[i]>lsp[i+1]-m)
 574          lsp[i]= SHR16(lsp[i],1) + SHR16(lsp[i+1]-m,1);
 575    }
 576 }
 577
 578
 579 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)
 580 {
 581    int i;
 582    spx_word16_t tmp = DIV32_16(SHL32(1 + subframe,14),nb_subframes);
 583    spx_word16_t tmp2 = 16384-tmp;
 584    for (i=0;i<len;i++)
 585    {
 586       interp_lsp[i] = MULT16_16_P14(tmp2,old_lsp[i]) + MULT16_16_P14(tmp,new_lsp[i]);
 587    }
 588 }
 589
 590 #else
 591
 592 /*Makes sure the LSPs are stable*/
 593 void lsp_enforce_margin(spx_lsp_t *lsp, int len, spx_word16_t margin)
 594 {
 595    int i;
 596    if (lsp[0]<LSP_SCALING*margin)
 597       lsp[0]=LSP_SCALING*margin;
 598    if (lsp[len-1]>LSP_SCALING*(M_PI-margin))
 599       lsp[len-1]=LSP_SCALING*(M_PI-margin);
 600    for (i=1;i<len-1;i++)
 601    {
 602       if (lsp[i]<lsp[i-1]+LSP_SCALING*margin)
 603          lsp[i]=lsp[i-1]+LSP_SCALING*margin;
 604
 605       if (lsp[i]>lsp[i+1]-LSP_SCALING*margin)
 606          lsp[i]= .5f* (lsp[i] + lsp[i+1]-LSP_SCALING*margin);
 607    }
 608 }
 609
 610
 611 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)
 612 {
 613    int i;
 614    float tmp = (1.0f + subframe)/nb_subframes;
 615    for (i=0;i<len;i++)
 616    {
 617       interp_lsp[i] = (1-tmp)*old_lsp[i] + tmp*new_lsp[i];
 618    }
 619 }
 620
 621 #endif