Ignore:
Timestamp:
Jul 26, 2006 5:04:54 PM (18 years ago)
Author:
bennylp
Message:
  • Bring speex codec up to date with their SVN trunk
  • Speex codec should work in FIXED_POINT mode when PJ_HAS_FLOATING_POINT is set to zero.
  • ulaw2linear will return zero if zero is given (this would make the VAD works better, and it also fixed click noise when call is established/hangup).
File:
1 edited

Legend:

Unmodified
Added
Removed

  • pjproject/trunk/pjmedia/src/pjmedia-codec/speex/lsp.c

    r278 r628  
    11/*---------------------------------------------------------------------------*\ 
    22Original copyright 
    3         FILE........: AKSLSPD.C 
    4         TYPE........: Turbo C 
    5         COMPANY.....: Voicetronix 
     3        FILE........: lsp.c 
    64        AUTHOR......: David Rowe 
    75        DATE CREATED: 24/2/93 
    … …  
    4543*/ 
    4644 
     45/*---------------------------------------------------------------------------*\ 
     46 
     47  Introduction to Line Spectrum Pairs (LSPs) 
     48  ------------------------------------------ 
     49 
     50  LSPs are used to encode the LPC filter coefficients {ak} for 
     51  transmission over the channel.  LSPs have several properties (like 
     52  less sensitivity to quantisation noise) that make them superior to 
     53  direct quantisation of {ak}. 
     54 
     55  A(z) is a polynomial of order lpcrdr with {ak} as the coefficients. 
     56 
     57  A(z) is transformed to P(z) and Q(z) (using a substitution and some 
     58  algebra), to obtain something like: 
     59 
     60    A(z) = 0.5[P(z)(z+z^-1) + Q(z)(z-z^-1)]  (1) 
     61 
     62  As you can imagine A(z) has complex zeros all over the z-plane. P(z) 
     63  and Q(z) have the very neat property of only having zeros _on_ the 
     64  unit circle.  So to find them we take a test point z=exp(jw) and 
     65  evaluate P (exp(jw)) and Q(exp(jw)) using a grid of points between 0 
     66  and pi. 
     67 
     68  The zeros (roots) of P(z) also happen to alternate, which is why we 
     69  swap coefficients as we find roots.  So the process of finding the 
     70  LSP frequencies is basically finding the roots of 5th order 
     71  polynomials. 
     72 
     73  The root so P(z) and Q(z) occur in symmetrical pairs at +/-w, hence 
     74  the name Line Spectrum Pairs (LSPs). 
     75 
     76  To convert back to ak we just evaluate (1), "clocking" an impulse 
     77  thru it lpcrdr times gives us the impulse response of A(z) which is 
     78  {ak}. 
     79 
     80\*---------------------------------------------------------------------------*/ 
     81 
    4782#ifdef HAVE_CONFIG_H 
    4883#include "config.h" 
    … …  
    6499#ifdef FIXED_POINT 
    65100 
    66  
    67  
    68101#define FREQ_SCALE 16384 
    69102 
    … …  
    73106/*#define X2ANGLE(x) (acos(.00006103515625*(x))*LSP_SCALING)*/ 
    74107#define X2ANGLE(x) (spx_acos(x)) 
     108 
     109#ifdef BFIN_ASM 
     110#include "lsp_bfin.h" 
     111#endif 
    75112 
    76113#else 
    … …  
    89126/*---------------------------------------------------------------------------*\ 
    90127 
    91         FUNCTION....: cheb_poly_eva() 
    92  
    93         AUTHOR......: David Rowe 
    94         DATE CREATED: 24/2/93 
    95  
    96     This function evaluates a series of Chebyshev polynomials 
     128   FUNCTION....: cheb_poly_eva() 
     129 
     130   AUTHOR......: David Rowe 
     131   DATE CREATED: 24/2/93 
     132 
     133   This function evaluates a series of Chebyshev polynomials 
    97134 
    98135\*---------------------------------------------------------------------------*/ 
    … …  
    100137#ifdef FIXED_POINT 
    101138 
    102 static inline spx_word32_t cheb_poly_eva(spx_word32_t *coef,spx_word16_t x,int m,char *stack) 
    103 /*  float coef[]        coefficients of the polynomial to be evaluated  */ 
    104 /*  float x             the point where polynomial is to be evaluated   */ 
    105 /*  int m               order of the polynomial                         */ 
     139#ifndef OVERRIDE_CHEB_POLY_EVA 
     140static inline spx_word32_t cheb_poly_eva( 
     141  spx_word16_t *coef, /* P or Q coefs in Q13 format               */ 
     142  spx_word16_t     x, /* cos of freq (-1.0 to 1.0) in Q14 format  */ 
     143  int              m, /* LPC order/2                              */ 
     144  char         *stack 
     145) 
    106146{ 
    107147    int i; 
    108     VARDECL(spx_word16_t *T); 
     148    spx_word16_t b0, b1; 
    109149    spx_word32_t sum; 
    110     int m2=m>>1; 
    111     VARDECL(spx_word16_t *coefn); 
    112150 
    113151    /*Prevents overflows*/ 
    … …  
    117155       x = -16383; 
    118156 
    119     /* Allocate memory for Chebyshev series formulation */ 
    120     ALLOC(T, m2+1, spx_word16_t); 
    121     ALLOC(coefn, m2+1, spx_word16_t); 
    122  
    123     for (i=0;i<m2+1;i++) 
     157    /* Initialise values */ 
     158    b1=16384; 
     159    b0=x; 
     160 
     161    /* Evaluate Chebyshev series formulation usin g iterative approach  */ 
     162    sum = ADD32(EXTEND32(coef[m]), EXTEND32(MULT16_16_P14(coef[m-1],x))); 
     163    for(i=2;i<=m;i++) 
    124164    { 
    125        coefn[i] = coef[i]; 
    126        /*printf ("%f ", coef[i]);*/ 
    127     } 
    128     /*printf ("\n");*/ 
    129  
    130     /* Initialise values */ 
    131     T[0]=16384; 
    132     T[1]=x; 
    133  
    134     /* Evaluate Chebyshev series formulation using iterative approach  */ 
    135     /* Evaluate polynomial and return value also free memory space */ 
    136     sum = ADD32(EXTEND32(coefn[m2]), EXTEND32(MULT16_16_P14(coefn[m2-1],x))); 
    137     /*x *= 2;*/ 
    138     for(i=2;i<=m2;i++) 
    139     { 
    140        T[i] = SUB16(MULT16_16_Q13(x,T[i-1]), T[i-2]); 
    141        sum = ADD32(sum, EXTEND32(MULT16_16_P14(coefn[m2-i],T[i]))); 
    142        /*printf ("%f ", sum);*/ 
    143     } 
    144      
    145     /*printf ("\n");*/ 
    146     return sum; 
    147 } 
    148 #else 
    149 static float cheb_poly_eva(spx_word32_t *coef,float x,int m,char *stack) 
    150 /*  float coef[]        coefficients of the polynomial to be evaluated  */ 
    151 /*  float x             the point where polynomial is to be evaluated   */ 
    152 /*  int m               order of the polynomial                         */ 
    153 { 
    154     int i; 
    155     VARDECL(float *T); 
    156     float sum; 
    157     int m2=m>>1; 
    158  
    159     /* Allocate memory for Chebyshev series formulation */ 
    160     ALLOC(T, m2+1, float); 
    161  
    162     /* Initialise values */ 
    163     T[0]=1; 
    164     T[1]=x; 
    165  
    166     /* Evaluate Chebyshev series formulation using iterative approach  */ 
    167     /* Evaluate polynomial and return value also free memory space */ 
    168     sum = coef[m2] + coef[m2-1]*x; 
    169     x *= 2; 
    170     for(i=2;i<=m2;i++) 
    171     { 
    172        T[i] = x*T[i-1] - T[i-2]; 
    173        sum += coef[m2-i] * T[i]; 
     165       spx_word16_t tmp=b0; 
     166       b0 = SUB16(MULT16_16_Q13(x,b0), b1); 
     167       b1 = tmp; 
     168       sum = ADD32(sum, EXTEND32(MULT16_16_P14(coef[m-i],b0))); 
    174169    } 
    175170     
    … …  
    178173#endif 
    179174 
     175#else 
     176 
     177static float cheb_poly_eva(spx_word32_t *coef, spx_word16_t x, int m, char *stack) 
     178{ 
     179   int k; 
     180   float b0, b1, tmp; 
     181 
     182   /* Initial conditions */ 
     183   b0=0; /* b_(m+1) */ 
     184   b1=0; /* b_(m+2) */ 
     185 
     186   x*=2; 
     187 
     188   /* Calculate the b_(k) */ 
     189   for(k=m;k>0;k--) 
     190   { 
     191      tmp=b0;                           /* tmp holds the previous value of b0 */ 
     192      b0=x*b0-b1+coef[m-k];    /* b0 holds its new value based on b0 and b1 */ 
     193      b1=tmp;                           /* b1 holds the previous value of b0 */ 
     194   } 
     195 
     196   return(-b1+.5*x*b0+coef[m]); 
     197} 
     198#endif 
     199 
    180200/*---------------------------------------------------------------------------*\ 
    181201 
    182         FUNCTION....: lpc_to_lsp() 
    183  
    184         AUTHOR......: David Rowe 
    185         DATE CREATED: 24/2/93 
     202    FUNCTION....: lpc_to_lsp() 
     203 
     204    AUTHOR......: David Rowe 
     205    DATE CREATED: 24/2/93 
    186206 
    187207    This function converts LPC coefficients to LSP 
    … …  
    211231    VARDECL(spx_word32_t *Q);                   /* ptrs for memory allocation           */ 
    212232    VARDECL(spx_word32_t *P); 
     233    VARDECL(spx_word16_t *Q16);         /* ptrs for memory allocation           */ 
     234    VARDECL(spx_word16_t *P16); 
    213235    spx_word32_t *px;                   /* ptrs of respective P'(z) & Q'(z)     */ 
    214236    spx_word32_t *qx; 
    215237    spx_word32_t *p; 
    216238    spx_word32_t *q; 
    217     spx_word32_t *pt;                   /* ptr used for cheb_poly_eval() 
     239    spx_word16_t *pt;                   /* ptr used for cheb_poly_eval() 
    218240                                whether P' or Q'                        */ 
    219241    int roots=0;                /* DR 8/2/94: number of roots found     */ 
    … …  
    277299    qx = Q; 
    278300 
     301    /* now that we have computed P and Q convert to 16 bits to 
     302       speed up cheb_poly_eval */ 
     303 
     304    ALLOC(P16, m+1, spx_word16_t); 
     305    ALLOC(Q16, m+1, spx_word16_t); 
     306 
     307    for (i=0;i<m+1;i++) 
     308    { 
     309       P16[i] = P[i]; 
     310       Q16[i] = Q[i]; 
     311    } 
     312 
    279313    /* Search for a zero in P'(z) polynomial first and then alternate to Q'(z). 
    280314    Keep alternating between the two polynomials as each zero is found  */ 
    … …  
    283317    xl = FREQ_SCALE;                    /* start at point xl = 1                */ 
    284318 
    285  
    286319    for(j=0;j<lpcrdr;j++){ 
    287320        if(j&1)                 /* determines whether P' or Q' is eval. */ 
    288             pt = qx; 
     321            pt = Q16; 
    289322        else 
    290             pt = px; 
    291  
    292         psuml = cheb_poly_eva(pt,xl,lpcrdr,stack);      /* evals poly. at xl    */ 
     323            pt = P16; 
     324 
     325        psuml = cheb_poly_eva(pt,xl,m,stack);   /* evals poly. at xl    */ 
    293326        flag = 1; 
    294327        while(flag && (xr >= -FREQ_SCALE)){ 
    … …  
    305338#endif 
    306339           xr = SUB16(xl, dd);                          /* interval spacing     */ 
    307             psumr = cheb_poly_eva(pt,xr,lpcrdr,stack);/* poly(xl-delta_x)       */ 
     340            psumr = cheb_poly_eva(pt,xr,m,stack);/* poly(xl-delta_x)    */ 
    308341            temp_psumr = psumr; 
    309342            temp_xr = xr; 
    … …  
    329362                    xm = .5*(xl+xr);            /* bisect the interval  */ 
    330363#endif 
    331                     psumm=cheb_poly_eva(pt,xm,lpcrdr,stack); 
     364                    psumm=cheb_poly_eva(pt,xm,m,stack); 
    332365                    /*if(psumm*psuml>0.)*/ 
    333366                    if(!SIGN_CHANGE(psumm,psuml)) 
    … …  
    355388} 
    356389 
    357  
    358390/*---------------------------------------------------------------------------*\ 
    359391 
    … …  
    363395        DATE CREATED: 24/2/93 
    364396 
    365     lsp_to_lpc: This function converts LSP coefficients to LPC 
    366     coefficients. 
     397        Converts LSP coefficients to LPC coefficients. 
    367398 
    368399\*---------------------------------------------------------------------------*/ 
    … …  
    374405/*  float *ak           array of LPC coefficients                       */ 
    375406/*  int lpcrdr          order of LPC coefficients                       */ 
    376  
    377  
    378407{ 
    379408    int i,j; 
    380     spx_word32_t xout1,xout2,xin1,xin2; 
    381     VARDECL(spx_word32_t *Wp); 
    382     spx_word32_t *pw,*n1,*n2,*n3,*n4=NULL; 
     409    spx_word32_t xout1,xout2,xin; 
     410    spx_word32_t mult, a; 
    383411    VARDECL(spx_word16_t *freqn); 
     412    VARDECL(spx_word32_t **xp); 
     413    VARDECL(spx_word32_t *xpmem); 
     414    VARDECL(spx_word32_t **xq); 
     415    VARDECL(spx_word32_t *xqmem); 
    384416    int m = lpcrdr>>1; 
     417 
     418    /*  
    385419     
     420       Reconstruct P(z) and Q(z) by cascading second order polynomials 
     421       in form 1 - 2cos(w)z(-1) + z(-2), where w is the LSP frequency. 
     422       In the time domain this is: 
     423 
     424       y(n) = x(n) - 2cos(w)x(n-1) + x(n-2) 
     425     
     426       This is what the ALLOCS below are trying to do: 
     427 
     428         int xp[m+1][lpcrdr+1+2]; // P matrix in QIMP 
     429         int xq[m+1][lpcrdr+1+2]; // Q matrix in QIMP 
     430 
     431       These matrices store the output of each stage on each row.  The 
     432       final (m-th) row has the output of the final (m-th) cascaded 
     433       2nd order filter.  The first row is the impulse input to the 
     434       system (not written as it is known). 
     435 
     436       The version below takes advantage of the fact that a lot of the 
     437       outputs are zero or known, for example if we put an inpulse 
     438       into the first section the "clock" it 10 times only the first 3 
     439       outputs samples are non-zero (it's an FIR filter). 
     440    */ 
     441 
     442    ALLOC(xp, (m+1), spx_word32_t*); 
     443    ALLOC(xpmem, (m+1)*(lpcrdr+1+2), spx_word32_t); 
     444 
     445    ALLOC(xq, (m+1), spx_word32_t*); 
     446    ALLOC(xqmem, (m+1)*(lpcrdr+1+2), spx_word32_t); 
     447     
     448    for(i=0; i<=m; i++) { 
     449      xp[i] = xpmem + i*(lpcrdr+1+2); 
     450      xq[i] = xqmem + i*(lpcrdr+1+2); 
     451    } 
     452 
     453    /* work out 2cos terms in Q14 */ 
     454 
    386455    ALLOC(freqn, lpcrdr, spx_word16_t); 
    387     for (i=0;i<lpcrdr;i++) 
     456    for (i=0;i<lpcrdr;i++)  
    388457       freqn[i] = ANGLE2X(freq[i]); 
    389458 
    390     ALLOC(Wp, 4*m+2, spx_word32_t); 
    391     pw = Wp; 
    392  
    393  
    394     /* initialise contents of array */ 
    395  
    396     for(i=0;i<=4*m+1;i++){              /* set contents of buffer to 0 */ 
    397         *pw++ = 0; 
    398     } 
    399  
    400     /* Set pointers up */ 
    401  
    402     pw = Wp; 
    403     xin1 = 1048576; 
    404     xin2 = 1048576; 
    405  
    406     /* reconstruct P(z) and Q(z) by  cascading second order 
    407       polynomials in form 1 - 2xz(-1) +z(-2), where x is the 
    408       LSP coefficient */ 
    409  
    410     for(j=0;j<=lpcrdr;j++){ 
    411        spx_word16_t *fr=freqn; 
    412         for(i=0;i<m;i++){ 
    413             n1 = pw+(i<<2); 
    414             n2 = n1 + 1; 
    415             n3 = n2 + 1; 
    416             n4 = n3 + 1; 
    417             xout1 = ADD32(SUB32(xin1, MULT16_32_Q14(*fr,*n1)), *n2); 
    418             fr++; 
    419             xout2 = ADD32(SUB32(xin2, MULT16_32_Q14(*fr,*n3)), *n4); 
    420             fr++; 
    421             *n2 = *n1; 
    422             *n4 = *n3; 
    423             *n1 = xin1; 
    424             *n3 = xin2; 
    425             xin1 = xout1; 
    426             xin2 = xout2; 
    427         } 
    428         xout1 = xin1 + *(n4+1); 
    429         xout2 = xin2 - *(n4+2); 
    430         /* FIXME: perhaps apply bandwidth expansion in case of overflow? */ 
    431         if (j>0) 
    432         { 
    433         if (xout1 + xout2>SHL32(EXTEND32(32766),8)) 
    434            ak[j-1] = 32767; 
    435         else if (xout1 + xout2 < -SHL32(EXTEND32(32766),8)) 
    436            ak[j-1] = -32767; 
    437         else 
    438            ak[j-1] = EXTRACT16(PSHR32(ADD32(xout1,xout2),8)); 
    439         } else {/*speex_warning_int("ak[0] = ", EXTRACT16(PSHR32(ADD32(xout1,xout2),8)));*/} 
    440         *(n4+1) = xin1; 
    441         *(n4+2) = xin2; 
    442  
    443         xin1 = 0; 
    444         xin2 = 0; 
    445     } 
    446 } 
     459    #define QIMP  21   /* scaling for impulse */ 
     460 
     461    xin = SHL32(EXTEND32(1), (QIMP-1)); /* 0.5 in QIMP format */ 
     462    
     463    /* first col and last non-zero values of each row are trivial */ 
     464     
     465    for(i=0;i<=m;i++) { 
     466     xp[i][1] = 0; 
     467     xp[i][2] = xin; 
     468     xp[i][2+2*i] = xin; 
     469     xq[i][1] = 0; 
     470     xq[i][2] = xin; 
     471     xq[i][2+2*i] = xin; 
     472    } 
     473 
     474    /* 2nd row (first output row) is trivial */ 
     475 
     476    xp[1][3] = -MULT16_32_Q14(freqn[0],xp[0][2]); 
     477    xq[1][3] = -MULT16_32_Q14(freqn[1],xq[0][2]); 
     478 
     479    xout1 = xout2 = 0; 
     480 
     481    /* now generate remaining rows */ 
     482 
     483    for(i=1;i<m;i++) { 
     484 
     485      for(j=1;j<2*(i+1)-1;j++) { 
     486        mult = MULT16_32_Q14(freqn[2*i],xp[i][j+1]); 
     487        xp[i+1][j+2] = ADD32(SUB32(xp[i][j+2], mult), xp[i][j]); 
     488        mult = MULT16_32_Q14(freqn[2*i+1],xq[i][j+1]); 
     489        xq[i+1][j+2] = ADD32(SUB32(xq[i][j+2], mult), xq[i][j]); 
     490      } 
     491 
     492      /* for last col xp[i][j+2] = xq[i][j+2] = 0 */ 
     493 
     494      mult = MULT16_32_Q14(freqn[2*i],xp[i][j+1]); 
     495      xp[i+1][j+2] = SUB32(xp[i][j], mult); 
     496      mult = MULT16_32_Q14(freqn[2*i+1],xq[i][j+1]); 
     497      xq[i+1][j+2] = SUB32(xq[i][j], mult); 
     498    } 
     499 
     500    /* process last row to extra a{k} */ 
     501 
     502    for(j=1;j<=lpcrdr;j++) { 
     503      int shift = QIMP-13; 
     504 
     505      /* final filter sections */ 
     506      a = PSHR32(xp[m][j+2] + xout1 + xq[m][j+2] - xout2, shift);  
     507      xout1 = xp[m][j+2]; 
     508      xout2 = xq[m][j+2]; 
     509       
     510      /* hard limit ak's to +/- 32767 */ 
     511 
     512      if (a < -32767) a = 32767; 
     513      if (a > 32767) a = 32767; 
     514      ak[j-1] = (short)a; 
     515      
     516    } 
     517 
     518} 
     519 
    447520#else 
    448521 
Note: See TracChangeset for help on using the changeset viewer.