michael@0: /***********************************************************************
michael@0: Copyright (c) 2006-2011, Skype Limited. All rights reserved.
michael@0: Redistribution and use in source and binary forms, with or without
michael@0: modification, are permitted provided that the following conditions
michael@0: are met:
michael@0: - Redistributions of source code must retain the above copyright notice,
michael@0: this list of conditions and the following disclaimer.
michael@0: - Redistributions in binary form must reproduce the above copyright
michael@0: notice, this list of conditions and the following disclaimer in the
michael@0: documentation and/or other materials provided with the distribution.
michael@0: - Neither the name of Internet Society, IETF or IETF Trust, nor the
michael@0: names of specific contributors, may be used to endorse or promote
michael@0: products derived from this software without specific prior written
michael@0: permission.
michael@0: THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
michael@0: AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
michael@0: IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
michael@0: ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
michael@0: LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
michael@0: CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
michael@0: SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
michael@0: INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
michael@0: CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
michael@0: ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
michael@0: POSSIBILITY OF SUCH DAMAGE.
michael@0: ***********************************************************************/
michael@0: 
michael@0: /* Conversion between prediction filter coefficients and NLSFs  */
michael@0: /* Requires the order to be an even number                      */
michael@0: /* A piecewise linear approximation maps LSF <-> cos(LSF)       */
michael@0: /* Therefore the result is not accurate NLSFs, but the two      */
michael@0: /* functions are accurate inverses of each other                */
michael@0: 
michael@0: #ifdef HAVE_CONFIG_H
michael@0: #include "config.h"
michael@0: #endif
michael@0: 
michael@0: #include "SigProc_FIX.h"
michael@0: #include "tables.h"
michael@0: 
michael@0: /* Number of binary divisions, when not in low complexity mode */
michael@0: #define BIN_DIV_STEPS_A2NLSF_FIX      3 /* must be no higher than 16 - log2( LSF_COS_TAB_SZ_FIX ) */
michael@0: #define MAX_ITERATIONS_A2NLSF_FIX    30
michael@0: 
michael@0: /* Helper function for A2NLSF(..)                    */
michael@0: /* Transforms polynomials from cos(n*f) to cos(f)^n  */
michael@0: static OPUS_INLINE void silk_A2NLSF_trans_poly(
michael@0:     opus_int32          *p,                     /* I/O    Polynomial                                */
michael@0:     const opus_int      dd                      /* I      Polynomial order (= filter order / 2 )    */
michael@0: )
michael@0: {
michael@0:     opus_int k, n;
michael@0: 
michael@0:     for( k = 2; k <= dd; k++ ) {
michael@0:         for( n = dd; n > k; n-- ) {
michael@0:             p[ n - 2 ] -= p[ n ];
michael@0:         }
michael@0:         p[ k - 2 ] -= silk_LSHIFT( p[ k ], 1 );
michael@0:     }
michael@0: }
michael@0: /* Helper function for A2NLSF(..) */
michael@0: /* Polynomial evaluation          */
michael@0: static OPUS_INLINE opus_int32 silk_A2NLSF_eval_poly( /* return the polynomial evaluation, in Q16     */
michael@0:     opus_int32          *p,                     /* I    Polynomial, Q16                         */
michael@0:     const opus_int32    x,                      /* I    Evaluation point, Q12                   */
michael@0:     const opus_int      dd                      /* I    Order                                   */
michael@0: )
michael@0: {
michael@0:     opus_int   n;
michael@0:     opus_int32 x_Q16, y32;
michael@0: 
michael@0:     y32 = p[ dd ];                                  /* Q16 */
michael@0:     x_Q16 = silk_LSHIFT( x, 4 );
michael@0:     for( n = dd - 1; n >= 0; n-- ) {
michael@0:         y32 = silk_SMLAWW( p[ n ], y32, x_Q16 );    /* Q16 */
michael@0:     }
michael@0:     return y32;
michael@0: }
michael@0: 
michael@0: static OPUS_INLINE void silk_A2NLSF_init(
michael@0:      const opus_int32    *a_Q16,
michael@0:      opus_int32          *P,
michael@0:      opus_int32          *Q,
michael@0:      const opus_int      dd
michael@0: )
michael@0: {
michael@0:     opus_int k;
michael@0: 
michael@0:     /* Convert filter coefs to even and odd polynomials */
michael@0:     P[dd] = silk_LSHIFT( 1, 16 );
michael@0:     Q[dd] = silk_LSHIFT( 1, 16 );
michael@0:     for( k = 0; k < dd; k++ ) {
michael@0:         P[ k ] = -a_Q16[ dd - k - 1 ] - a_Q16[ dd + k ];    /* Q16 */
michael@0:         Q[ k ] = -a_Q16[ dd - k - 1 ] + a_Q16[ dd + k ];    /* Q16 */
michael@0:     }
michael@0: 
michael@0:     /* Divide out zeros as we have that for even filter orders, */
michael@0:     /* z =  1 is always a root in Q, and                        */
michael@0:     /* z = -1 is always a root in P                             */
michael@0:     for( k = dd; k > 0; k-- ) {
michael@0:         P[ k - 1 ] -= P[ k ];
michael@0:         Q[ k - 1 ] += Q[ k ];
michael@0:     }
michael@0: 
michael@0:     /* Transform polynomials from cos(n*f) to cos(f)^n */
michael@0:     silk_A2NLSF_trans_poly( P, dd );
michael@0:     silk_A2NLSF_trans_poly( Q, dd );
michael@0: }
michael@0: 
michael@0: /* Compute Normalized Line Spectral Frequencies (NLSFs) from whitening filter coefficients      */
michael@0: /* If not all roots are found, the a_Q16 coefficients are bandwidth expanded until convergence. */
michael@0: void silk_A2NLSF(
michael@0:     opus_int16                  *NLSF,              /* O    Normalized Line Spectral Frequencies in Q15 (0..2^15-1) [d] */
michael@0:     opus_int32                  *a_Q16,             /* I/O  Monic whitening filter coefficients in Q16 [d]              */
michael@0:     const opus_int              d                   /* I    Filter order (must be even)                                 */
michael@0: )
michael@0: {
michael@0:     opus_int      i, k, m, dd, root_ix, ffrac;
michael@0:     opus_int32 xlo, xhi, xmid;
michael@0:     opus_int32 ylo, yhi, ymid, thr;
michael@0:     opus_int32 nom, den;
michael@0:     opus_int32 P[ SILK_MAX_ORDER_LPC / 2 + 1 ];
michael@0:     opus_int32 Q[ SILK_MAX_ORDER_LPC / 2 + 1 ];
michael@0:     opus_int32 *PQ[ 2 ];
michael@0:     opus_int32 *p;
michael@0: 
michael@0:     /* Store pointers to array */
michael@0:     PQ[ 0 ] = P;
michael@0:     PQ[ 1 ] = Q;
michael@0: 
michael@0:     dd = silk_RSHIFT( d, 1 );
michael@0: 
michael@0:     silk_A2NLSF_init( a_Q16, P, Q, dd );
michael@0: 
michael@0:     /* Find roots, alternating between P and Q */
michael@0:     p = P;                          /* Pointer to polynomial */
michael@0: 
michael@0:     xlo = silk_LSFCosTab_FIX_Q12[ 0 ]; /* Q12*/
michael@0:     ylo = silk_A2NLSF_eval_poly( p, xlo, dd );
michael@0: 
michael@0:     if( ylo < 0 ) {
michael@0:         /* Set the first NLSF to zero and move on to the next */
michael@0:         NLSF[ 0 ] = 0;
michael@0:         p = Q;                      /* Pointer to polynomial */
michael@0:         ylo = silk_A2NLSF_eval_poly( p, xlo, dd );
michael@0:         root_ix = 1;                /* Index of current root */
michael@0:     } else {
michael@0:         root_ix = 0;                /* Index of current root */
michael@0:     }
michael@0:     k = 1;                          /* Loop counter */
michael@0:     i = 0;                          /* Counter for bandwidth expansions applied */
michael@0:     thr = 0;
michael@0:     while( 1 ) {
michael@0:         /* Evaluate polynomial */
michael@0:         xhi = silk_LSFCosTab_FIX_Q12[ k ]; /* Q12 */
michael@0:         yhi = silk_A2NLSF_eval_poly( p, xhi, dd );
michael@0: 
michael@0:         /* Detect zero crossing */
michael@0:         if( ( ylo <= 0 && yhi >= thr ) || ( ylo >= 0 && yhi <= -thr ) ) {
michael@0:             if( yhi == 0 ) {
michael@0:                 /* If the root lies exactly at the end of the current       */
michael@0:                 /* interval, look for the next root in the next interval    */
michael@0:                 thr = 1;
michael@0:             } else {
michael@0:                 thr = 0;
michael@0:             }
michael@0:             /* Binary division */
michael@0:             ffrac = -256;
michael@0:             for( m = 0; m < BIN_DIV_STEPS_A2NLSF_FIX; m++ ) {
michael@0:                 /* Evaluate polynomial */
michael@0:                 xmid = silk_RSHIFT_ROUND( xlo + xhi, 1 );
michael@0:                 ymid = silk_A2NLSF_eval_poly( p, xmid, dd );
michael@0: 
michael@0:                 /* Detect zero crossing */
michael@0:                 if( ( ylo <= 0 && ymid >= 0 ) || ( ylo >= 0 && ymid <= 0 ) ) {
michael@0:                     /* Reduce frequency */
michael@0:                     xhi = xmid;
michael@0:                     yhi = ymid;
michael@0:                 } else {
michael@0:                     /* Increase frequency */
michael@0:                     xlo = xmid;
michael@0:                     ylo = ymid;
michael@0:                     ffrac = silk_ADD_RSHIFT( ffrac, 128, m );
michael@0:                 }
michael@0:             }
michael@0: 
michael@0:             /* Interpolate */
michael@0:             if( silk_abs( ylo ) < 65536 ) {
michael@0:                 /* Avoid dividing by zero */
michael@0:                 den = ylo - yhi;
michael@0:                 nom = silk_LSHIFT( ylo, 8 - BIN_DIV_STEPS_A2NLSF_FIX ) + silk_RSHIFT( den, 1 );
michael@0:                 if( den != 0 ) {
michael@0:                     ffrac += silk_DIV32( nom, den );
michael@0:                 }
michael@0:             } else {
michael@0:                 /* No risk of dividing by zero because abs(ylo - yhi) >= abs(ylo) >= 65536 */
michael@0:                 ffrac += silk_DIV32( ylo, silk_RSHIFT( ylo - yhi, 8 - BIN_DIV_STEPS_A2NLSF_FIX ) );
michael@0:             }
michael@0:             NLSF[ root_ix ] = (opus_int16)silk_min_32( silk_LSHIFT( (opus_int32)k, 8 ) + ffrac, silk_int16_MAX );
michael@0: 
michael@0:             silk_assert( NLSF[ root_ix ] >= 0 );
michael@0: 
michael@0:             root_ix++;        /* Next root */
michael@0:             if( root_ix >= d ) {
michael@0:                 /* Found all roots */
michael@0:                 break;
michael@0:             }
michael@0:             /* Alternate pointer to polynomial */
michael@0:             p = PQ[ root_ix & 1 ];
michael@0: 
michael@0:             /* Evaluate polynomial */
michael@0:             xlo = silk_LSFCosTab_FIX_Q12[ k - 1 ]; /* Q12*/
michael@0:             ylo = silk_LSHIFT( 1 - ( root_ix & 2 ), 12 );
michael@0:         } else {
michael@0:             /* Increment loop counter */
michael@0:             k++;
michael@0:             xlo = xhi;
michael@0:             ylo = yhi;
michael@0:             thr = 0;
michael@0: 
michael@0:             if( k > LSF_COS_TAB_SZ_FIX ) {
michael@0:                 i++;
michael@0:                 if( i > MAX_ITERATIONS_A2NLSF_FIX ) {
michael@0:                     /* Set NLSFs to white spectrum and exit */
michael@0:                     NLSF[ 0 ] = (opus_int16)silk_DIV32_16( 1 << 15, d + 1 );
michael@0:                     for( k = 1; k < d; k++ ) {
michael@0:                         NLSF[ k ] = (opus_int16)silk_SMULBB( k + 1, NLSF[ 0 ] );
michael@0:                     }
michael@0:                     return;
michael@0:                 }
michael@0: 
michael@0:                 /* Error: Apply progressively more bandwidth expansion and run again */
michael@0:                 silk_bwexpander_32( a_Q16, d, 65536 - silk_SMULBB( 10 + i, i ) ); /* 10_Q16 = 0.00015*/
michael@0: 
michael@0:                 silk_A2NLSF_init( a_Q16, P, Q, dd );
michael@0:                 p = P;                            /* Pointer to polynomial */
michael@0:                 xlo = silk_LSFCosTab_FIX_Q12[ 0 ]; /* Q12*/
michael@0:                 ylo = silk_A2NLSF_eval_poly( p, xlo, dd );
michael@0:                 if( ylo < 0 ) {
michael@0:                     /* Set the first NLSF to zero and move on to the next */
michael@0:                     NLSF[ 0 ] = 0;
michael@0:                     p = Q;                        /* Pointer to polynomial */
michael@0:                     ylo = silk_A2NLSF_eval_poly( p, xlo, dd );
michael@0:                     root_ix = 1;                  /* Index of current root */
michael@0:                 } else {
michael@0:                     root_ix = 0;                  /* Index of current root */
michael@0:                 }
michael@0:                 k = 1;                            /* Reset loop counter */
michael@0:             }
michael@0:         }
michael@0:     }
michael@0: }