1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/media/libopus/silk/A2NLSF.c Wed Dec 31 06:09:35 2014 +0100
1.3 @@ -0,0 +1,252 @@
1.4 +/***********************************************************************
1.5 +Copyright (c) 2006-2011, Skype Limited. All rights reserved.
1.6 +Redistribution and use in source and binary forms, with or without
1.7 +modification, are permitted provided that the following conditions
1.8 +are met:
1.9 +- Redistributions of source code must retain the above copyright notice,
1.10 +this list of conditions and the following disclaimer.
1.11 +- Redistributions in binary form must reproduce the above copyright
1.12 +notice, this list of conditions and the following disclaimer in the
1.13 +documentation and/or other materials provided with the distribution.
1.14 +- Neither the name of Internet Society, IETF or IETF Trust, nor the
1.15 +names of specific contributors, may be used to endorse or promote
1.16 +products derived from this software without specific prior written
1.17 +permission.
1.18 +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
1.19 +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
1.20 +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
1.21 +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
1.22 +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
1.23 +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
1.24 +SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
1.25 +INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
1.26 +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
1.27 +ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
1.28 +POSSIBILITY OF SUCH DAMAGE.
1.29 +***********************************************************************/
1.30 +
1.31 +/* Conversion between prediction filter coefficients and NLSFs */
1.32 +/* Requires the order to be an even number */
1.33 +/* A piecewise linear approximation maps LSF <-> cos(LSF) */
1.34 +/* Therefore the result is not accurate NLSFs, but the two */
1.35 +/* functions are accurate inverses of each other */
1.36 +
1.37 +#ifdef HAVE_CONFIG_H
1.38 +#include "config.h"
1.39 +#endif
1.40 +
1.41 +#include "SigProc_FIX.h"
1.42 +#include "tables.h"
1.43 +
1.44 +/* Number of binary divisions, when not in low complexity mode */
1.45 +#define BIN_DIV_STEPS_A2NLSF_FIX 3 /* must be no higher than 16 - log2( LSF_COS_TAB_SZ_FIX ) */
1.46 +#define MAX_ITERATIONS_A2NLSF_FIX 30
1.47 +
1.48 +/* Helper function for A2NLSF(..) */
1.49 +/* Transforms polynomials from cos(n*f) to cos(f)^n */
1.50 +static OPUS_INLINE void silk_A2NLSF_trans_poly(
1.51 + opus_int32 *p, /* I/O Polynomial */
1.52 + const opus_int dd /* I Polynomial order (= filter order / 2 ) */
1.53 +)
1.54 +{
1.55 + opus_int k, n;
1.56 +
1.57 + for( k = 2; k <= dd; k++ ) {
1.58 + for( n = dd; n > k; n-- ) {
1.59 + p[ n - 2 ] -= p[ n ];
1.60 + }
1.61 + p[ k - 2 ] -= silk_LSHIFT( p[ k ], 1 );
1.62 + }
1.63 +}
1.64 +/* Helper function for A2NLSF(..) */
1.65 +/* Polynomial evaluation */
1.66 +static OPUS_INLINE opus_int32 silk_A2NLSF_eval_poly( /* return the polynomial evaluation, in Q16 */
1.67 + opus_int32 *p, /* I Polynomial, Q16 */
1.68 + const opus_int32 x, /* I Evaluation point, Q12 */
1.69 + const opus_int dd /* I Order */
1.70 +)
1.71 +{
1.72 + opus_int n;
1.73 + opus_int32 x_Q16, y32;
1.74 +
1.75 + y32 = p[ dd ]; /* Q16 */
1.76 + x_Q16 = silk_LSHIFT( x, 4 );
1.77 + for( n = dd - 1; n >= 0; n-- ) {
1.78 + y32 = silk_SMLAWW( p[ n ], y32, x_Q16 ); /* Q16 */
1.79 + }
1.80 + return y32;
1.81 +}
1.82 +
1.83 +static OPUS_INLINE void silk_A2NLSF_init(
1.84 + const opus_int32 *a_Q16,
1.85 + opus_int32 *P,
1.86 + opus_int32 *Q,
1.87 + const opus_int dd
1.88 +)
1.89 +{
1.90 + opus_int k;
1.91 +
1.92 + /* Convert filter coefs to even and odd polynomials */
1.93 + P[dd] = silk_LSHIFT( 1, 16 );
1.94 + Q[dd] = silk_LSHIFT( 1, 16 );
1.95 + for( k = 0; k < dd; k++ ) {
1.96 + P[ k ] = -a_Q16[ dd - k - 1 ] - a_Q16[ dd + k ]; /* Q16 */
1.97 + Q[ k ] = -a_Q16[ dd - k - 1 ] + a_Q16[ dd + k ]; /* Q16 */
1.98 + }
1.99 +
1.100 + /* Divide out zeros as we have that for even filter orders, */
1.101 + /* z = 1 is always a root in Q, and */
1.102 + /* z = -1 is always a root in P */
1.103 + for( k = dd; k > 0; k-- ) {
1.104 + P[ k - 1 ] -= P[ k ];
1.105 + Q[ k - 1 ] += Q[ k ];
1.106 + }
1.107 +
1.108 + /* Transform polynomials from cos(n*f) to cos(f)^n */
1.109 + silk_A2NLSF_trans_poly( P, dd );
1.110 + silk_A2NLSF_trans_poly( Q, dd );
1.111 +}
1.112 +
1.113 +/* Compute Normalized Line Spectral Frequencies (NLSFs) from whitening filter coefficients */
1.114 +/* If not all roots are found, the a_Q16 coefficients are bandwidth expanded until convergence. */
1.115 +void silk_A2NLSF(
1.116 + opus_int16 *NLSF, /* O Normalized Line Spectral Frequencies in Q15 (0..2^15-1) [d] */
1.117 + opus_int32 *a_Q16, /* I/O Monic whitening filter coefficients in Q16 [d] */
1.118 + const opus_int d /* I Filter order (must be even) */
1.119 +)
1.120 +{
1.121 + opus_int i, k, m, dd, root_ix, ffrac;
1.122 + opus_int32 xlo, xhi, xmid;
1.123 + opus_int32 ylo, yhi, ymid, thr;
1.124 + opus_int32 nom, den;
1.125 + opus_int32 P[ SILK_MAX_ORDER_LPC / 2 + 1 ];
1.126 + opus_int32 Q[ SILK_MAX_ORDER_LPC / 2 + 1 ];
1.127 + opus_int32 *PQ[ 2 ];
1.128 + opus_int32 *p;
1.129 +
1.130 + /* Store pointers to array */
1.131 + PQ[ 0 ] = P;
1.132 + PQ[ 1 ] = Q;
1.133 +
1.134 + dd = silk_RSHIFT( d, 1 );
1.135 +
1.136 + silk_A2NLSF_init( a_Q16, P, Q, dd );
1.137 +
1.138 + /* Find roots, alternating between P and Q */
1.139 + p = P; /* Pointer to polynomial */
1.140 +
1.141 + xlo = silk_LSFCosTab_FIX_Q12[ 0 ]; /* Q12*/
1.142 + ylo = silk_A2NLSF_eval_poly( p, xlo, dd );
1.143 +
1.144 + if( ylo < 0 ) {
1.145 + /* Set the first NLSF to zero and move on to the next */
1.146 + NLSF[ 0 ] = 0;
1.147 + p = Q; /* Pointer to polynomial */
1.148 + ylo = silk_A2NLSF_eval_poly( p, xlo, dd );
1.149 + root_ix = 1; /* Index of current root */
1.150 + } else {
1.151 + root_ix = 0; /* Index of current root */
1.152 + }
1.153 + k = 1; /* Loop counter */
1.154 + i = 0; /* Counter for bandwidth expansions applied */
1.155 + thr = 0;
1.156 + while( 1 ) {
1.157 + /* Evaluate polynomial */
1.158 + xhi = silk_LSFCosTab_FIX_Q12[ k ]; /* Q12 */
1.159 + yhi = silk_A2NLSF_eval_poly( p, xhi, dd );
1.160 +
1.161 + /* Detect zero crossing */
1.162 + if( ( ylo <= 0 && yhi >= thr ) || ( ylo >= 0 && yhi <= -thr ) ) {
1.163 + if( yhi == 0 ) {
1.164 + /* If the root lies exactly at the end of the current */
1.165 + /* interval, look for the next root in the next interval */
1.166 + thr = 1;
1.167 + } else {
1.168 + thr = 0;
1.169 + }
1.170 + /* Binary division */
1.171 + ffrac = -256;
1.172 + for( m = 0; m < BIN_DIV_STEPS_A2NLSF_FIX; m++ ) {
1.173 + /* Evaluate polynomial */
1.174 + xmid = silk_RSHIFT_ROUND( xlo + xhi, 1 );
1.175 + ymid = silk_A2NLSF_eval_poly( p, xmid, dd );
1.176 +
1.177 + /* Detect zero crossing */
1.178 + if( ( ylo <= 0 && ymid >= 0 ) || ( ylo >= 0 && ymid <= 0 ) ) {
1.179 + /* Reduce frequency */
1.180 + xhi = xmid;
1.181 + yhi = ymid;
1.182 + } else {
1.183 + /* Increase frequency */
1.184 + xlo = xmid;
1.185 + ylo = ymid;
1.186 + ffrac = silk_ADD_RSHIFT( ffrac, 128, m );
1.187 + }
1.188 + }
1.189 +
1.190 + /* Interpolate */
1.191 + if( silk_abs( ylo ) < 65536 ) {
1.192 + /* Avoid dividing by zero */
1.193 + den = ylo - yhi;
1.194 + nom = silk_LSHIFT( ylo, 8 - BIN_DIV_STEPS_A2NLSF_FIX ) + silk_RSHIFT( den, 1 );
1.195 + if( den != 0 ) {
1.196 + ffrac += silk_DIV32( nom, den );
1.197 + }
1.198 + } else {
1.199 + /* No risk of dividing by zero because abs(ylo - yhi) >= abs(ylo) >= 65536 */
1.200 + ffrac += silk_DIV32( ylo, silk_RSHIFT( ylo - yhi, 8 - BIN_DIV_STEPS_A2NLSF_FIX ) );
1.201 + }
1.202 + NLSF[ root_ix ] = (opus_int16)silk_min_32( silk_LSHIFT( (opus_int32)k, 8 ) + ffrac, silk_int16_MAX );
1.203 +
1.204 + silk_assert( NLSF[ root_ix ] >= 0 );
1.205 +
1.206 + root_ix++; /* Next root */
1.207 + if( root_ix >= d ) {
1.208 + /* Found all roots */
1.209 + break;
1.210 + }
1.211 + /* Alternate pointer to polynomial */
1.212 + p = PQ[ root_ix & 1 ];
1.213 +
1.214 + /* Evaluate polynomial */
1.215 + xlo = silk_LSFCosTab_FIX_Q12[ k - 1 ]; /* Q12*/
1.216 + ylo = silk_LSHIFT( 1 - ( root_ix & 2 ), 12 );
1.217 + } else {
1.218 + /* Increment loop counter */
1.219 + k++;
1.220 + xlo = xhi;
1.221 + ylo = yhi;
1.222 + thr = 0;
1.223 +
1.224 + if( k > LSF_COS_TAB_SZ_FIX ) {
1.225 + i++;
1.226 + if( i > MAX_ITERATIONS_A2NLSF_FIX ) {
1.227 + /* Set NLSFs to white spectrum and exit */
1.228 + NLSF[ 0 ] = (opus_int16)silk_DIV32_16( 1 << 15, d + 1 );
1.229 + for( k = 1; k < d; k++ ) {
1.230 + NLSF[ k ] = (opus_int16)silk_SMULBB( k + 1, NLSF[ 0 ] );
1.231 + }
1.232 + return;
1.233 + }
1.234 +
1.235 + /* Error: Apply progressively more bandwidth expansion and run again */
1.236 + silk_bwexpander_32( a_Q16, d, 65536 - silk_SMULBB( 10 + i, i ) ); /* 10_Q16 = 0.00015*/
1.237 +
1.238 + silk_A2NLSF_init( a_Q16, P, Q, dd );
1.239 + p = P; /* Pointer to polynomial */
1.240 + xlo = silk_LSFCosTab_FIX_Q12[ 0 ]; /* Q12*/
1.241 + ylo = silk_A2NLSF_eval_poly( p, xlo, dd );
1.242 + if( ylo < 0 ) {
1.243 + /* Set the first NLSF to zero and move on to the next */
1.244 + NLSF[ 0 ] = 0;
1.245 + p = Q; /* Pointer to polynomial */
1.246 + ylo = silk_A2NLSF_eval_poly( p, xlo, dd );
1.247 + root_ix = 1; /* Index of current root */
1.248 + } else {
1.249 + root_ix = 0; /* Index of current root */
1.250 + }
1.251 + k = 1; /* Reset loop counter */
1.252 + }
1.253 + }
1.254 + }
1.255 +}
