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