media/libopus/silk/A2NLSF.c

Wed, 31 Dec 2014 06:09:35 +0100

author
Michael Schloh von Bennewitz <michael@schloh.com>
date
Wed, 31 Dec 2014 06:09:35 +0100
changeset 0
6474c204b198
permissions
-rw-r--r--

Cloned upstream origin tor-browser at tor-browser-31.3.0esr-4.5-1-build1
revision ID fc1c9ff7c1b2defdbc039f12214767608f46423f for hacking purpose.

michael@0 1 /***********************************************************************
michael@0 2 Copyright (c) 2006-2011, Skype Limited. All rights reserved.
michael@0 3 Redistribution and use in source and binary forms, with or without
michael@0 4 modification, are permitted provided that the following conditions
michael@0 5 are met:
michael@0 6 - Redistributions of source code must retain the above copyright notice,
michael@0 7 this list of conditions and the following disclaimer.
michael@0 8 - Redistributions in binary form must reproduce the above copyright
michael@0 9 notice, this list of conditions and the following disclaimer in the
michael@0 10 documentation and/or other materials provided with the distribution.
michael@0 11 - Neither the name of Internet Society, IETF or IETF Trust, nor the
michael@0 12 names of specific contributors, may be used to endorse or promote
michael@0 13 products derived from this software without specific prior written
michael@0 14 permission.
michael@0 15 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
michael@0 16 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
michael@0 17 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
michael@0 18 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
michael@0 19 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
michael@0 20 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
michael@0 21 SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
michael@0 22 INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
michael@0 23 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
michael@0 24 ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
michael@0 25 POSSIBILITY OF SUCH DAMAGE.
michael@0 26 ***********************************************************************/
michael@0 27
michael@0 28 /* Conversion between prediction filter coefficients and NLSFs */
michael@0 29 /* Requires the order to be an even number */
michael@0 30 /* A piecewise linear approximation maps LSF <-> cos(LSF) */
michael@0 31 /* Therefore the result is not accurate NLSFs, but the two */
michael@0 32 /* functions are accurate inverses of each other */
michael@0 33
michael@0 34 #ifdef HAVE_CONFIG_H
michael@0 35 #include "config.h"
michael@0 36 #endif
michael@0 37
michael@0 38 #include "SigProc_FIX.h"
michael@0 39 #include "tables.h"
michael@0 40
michael@0 41 /* Number of binary divisions, when not in low complexity mode */
michael@0 42 #define BIN_DIV_STEPS_A2NLSF_FIX 3 /* must be no higher than 16 - log2( LSF_COS_TAB_SZ_FIX ) */
michael@0 43 #define MAX_ITERATIONS_A2NLSF_FIX 30
michael@0 44
michael@0 45 /* Helper function for A2NLSF(..) */
michael@0 46 /* Transforms polynomials from cos(n*f) to cos(f)^n */
michael@0 47 static OPUS_INLINE void silk_A2NLSF_trans_poly(
michael@0 48 opus_int32 *p, /* I/O Polynomial */
michael@0 49 const opus_int dd /* I Polynomial order (= filter order / 2 ) */
michael@0 50 )
michael@0 51 {
michael@0 52 opus_int k, n;
michael@0 53
michael@0 54 for( k = 2; k <= dd; k++ ) {
michael@0 55 for( n = dd; n > k; n-- ) {
michael@0 56 p[ n - 2 ] -= p[ n ];
michael@0 57 }
michael@0 58 p[ k - 2 ] -= silk_LSHIFT( p[ k ], 1 );
michael@0 59 }
michael@0 60 }
michael@0 61 /* Helper function for A2NLSF(..) */
michael@0 62 /* Polynomial evaluation */
michael@0 63 static OPUS_INLINE opus_int32 silk_A2NLSF_eval_poly( /* return the polynomial evaluation, in Q16 */
michael@0 64 opus_int32 *p, /* I Polynomial, Q16 */
michael@0 65 const opus_int32 x, /* I Evaluation point, Q12 */
michael@0 66 const opus_int dd /* I Order */
michael@0 67 )
michael@0 68 {
michael@0 69 opus_int n;
michael@0 70 opus_int32 x_Q16, y32;
michael@0 71
michael@0 72 y32 = p[ dd ]; /* Q16 */
michael@0 73 x_Q16 = silk_LSHIFT( x, 4 );
michael@0 74 for( n = dd - 1; n >= 0; n-- ) {
michael@0 75 y32 = silk_SMLAWW( p[ n ], y32, x_Q16 ); /* Q16 */
michael@0 76 }
michael@0 77 return y32;
michael@0 78 }
michael@0 79
michael@0 80 static OPUS_INLINE void silk_A2NLSF_init(
michael@0 81 const opus_int32 *a_Q16,
michael@0 82 opus_int32 *P,
michael@0 83 opus_int32 *Q,
michael@0 84 const opus_int dd
michael@0 85 )
michael@0 86 {
michael@0 87 opus_int k;
michael@0 88
michael@0 89 /* Convert filter coefs to even and odd polynomials */
michael@0 90 P[dd] = silk_LSHIFT( 1, 16 );
michael@0 91 Q[dd] = silk_LSHIFT( 1, 16 );
michael@0 92 for( k = 0; k < dd; k++ ) {
michael@0 93 P[ k ] = -a_Q16[ dd - k - 1 ] - a_Q16[ dd + k ]; /* Q16 */
michael@0 94 Q[ k ] = -a_Q16[ dd - k - 1 ] + a_Q16[ dd + k ]; /* Q16 */
michael@0 95 }
michael@0 96
michael@0 97 /* Divide out zeros as we have that for even filter orders, */
michael@0 98 /* z = 1 is always a root in Q, and */
michael@0 99 /* z = -1 is always a root in P */
michael@0 100 for( k = dd; k > 0; k-- ) {
michael@0 101 P[ k - 1 ] -= P[ k ];
michael@0 102 Q[ k - 1 ] += Q[ k ];
michael@0 103 }
michael@0 104
michael@0 105 /* Transform polynomials from cos(n*f) to cos(f)^n */
michael@0 106 silk_A2NLSF_trans_poly( P, dd );
michael@0 107 silk_A2NLSF_trans_poly( Q, dd );
michael@0 108 }
michael@0 109
michael@0 110 /* Compute Normalized Line Spectral Frequencies (NLSFs) from whitening filter coefficients */
michael@0 111 /* If not all roots are found, the a_Q16 coefficients are bandwidth expanded until convergence. */
michael@0 112 void silk_A2NLSF(
michael@0 113 opus_int16 *NLSF, /* O Normalized Line Spectral Frequencies in Q15 (0..2^15-1) [d] */
michael@0 114 opus_int32 *a_Q16, /* I/O Monic whitening filter coefficients in Q16 [d] */
michael@0 115 const opus_int d /* I Filter order (must be even) */
michael@0 116 )
michael@0 117 {
michael@0 118 opus_int i, k, m, dd, root_ix, ffrac;
michael@0 119 opus_int32 xlo, xhi, xmid;
michael@0 120 opus_int32 ylo, yhi, ymid, thr;
michael@0 121 opus_int32 nom, den;
michael@0 122 opus_int32 P[ SILK_MAX_ORDER_LPC / 2 + 1 ];
michael@0 123 opus_int32 Q[ SILK_MAX_ORDER_LPC / 2 + 1 ];
michael@0 124 opus_int32 *PQ[ 2 ];
michael@0 125 opus_int32 *p;
michael@0 126
michael@0 127 /* Store pointers to array */
michael@0 128 PQ[ 0 ] = P;
michael@0 129 PQ[ 1 ] = Q;
michael@0 130
michael@0 131 dd = silk_RSHIFT( d, 1 );
michael@0 132
michael@0 133 silk_A2NLSF_init( a_Q16, P, Q, dd );
michael@0 134
michael@0 135 /* Find roots, alternating between P and Q */
michael@0 136 p = P; /* Pointer to polynomial */
michael@0 137
michael@0 138 xlo = silk_LSFCosTab_FIX_Q12[ 0 ]; /* Q12*/
michael@0 139 ylo = silk_A2NLSF_eval_poly( p, xlo, dd );
michael@0 140
michael@0 141 if( ylo < 0 ) {
michael@0 142 /* Set the first NLSF to zero and move on to the next */
michael@0 143 NLSF[ 0 ] = 0;
michael@0 144 p = Q; /* Pointer to polynomial */
michael@0 145 ylo = silk_A2NLSF_eval_poly( p, xlo, dd );
michael@0 146 root_ix = 1; /* Index of current root */
michael@0 147 } else {
michael@0 148 root_ix = 0; /* Index of current root */
michael@0 149 }
michael@0 150 k = 1; /* Loop counter */
michael@0 151 i = 0; /* Counter for bandwidth expansions applied */
michael@0 152 thr = 0;
michael@0 153 while( 1 ) {
michael@0 154 /* Evaluate polynomial */
michael@0 155 xhi = silk_LSFCosTab_FIX_Q12[ k ]; /* Q12 */
michael@0 156 yhi = silk_A2NLSF_eval_poly( p, xhi, dd );
michael@0 157
michael@0 158 /* Detect zero crossing */
michael@0 159 if( ( ylo <= 0 && yhi >= thr ) || ( ylo >= 0 && yhi <= -thr ) ) {
michael@0 160 if( yhi == 0 ) {
michael@0 161 /* If the root lies exactly at the end of the current */
michael@0 162 /* interval, look for the next root in the next interval */
michael@0 163 thr = 1;
michael@0 164 } else {
michael@0 165 thr = 0;
michael@0 166 }
michael@0 167 /* Binary division */
michael@0 168 ffrac = -256;
michael@0 169 for( m = 0; m < BIN_DIV_STEPS_A2NLSF_FIX; m++ ) {
michael@0 170 /* Evaluate polynomial */
michael@0 171 xmid = silk_RSHIFT_ROUND( xlo + xhi, 1 );
michael@0 172 ymid = silk_A2NLSF_eval_poly( p, xmid, dd );
michael@0 173
michael@0 174 /* Detect zero crossing */
michael@0 175 if( ( ylo <= 0 && ymid >= 0 ) || ( ylo >= 0 && ymid <= 0 ) ) {
michael@0 176 /* Reduce frequency */
michael@0 177 xhi = xmid;
michael@0 178 yhi = ymid;
michael@0 179 } else {
michael@0 180 /* Increase frequency */
michael@0 181 xlo = xmid;
michael@0 182 ylo = ymid;
michael@0 183 ffrac = silk_ADD_RSHIFT( ffrac, 128, m );
michael@0 184 }
michael@0 185 }
michael@0 186
michael@0 187 /* Interpolate */
michael@0 188 if( silk_abs( ylo ) < 65536 ) {
michael@0 189 /* Avoid dividing by zero */
michael@0 190 den = ylo - yhi;
michael@0 191 nom = silk_LSHIFT( ylo, 8 - BIN_DIV_STEPS_A2NLSF_FIX ) + silk_RSHIFT( den, 1 );
michael@0 192 if( den != 0 ) {
michael@0 193 ffrac += silk_DIV32( nom, den );
michael@0 194 }
michael@0 195 } else {
michael@0 196 /* No risk of dividing by zero because abs(ylo - yhi) >= abs(ylo) >= 65536 */
michael@0 197 ffrac += silk_DIV32( ylo, silk_RSHIFT( ylo - yhi, 8 - BIN_DIV_STEPS_A2NLSF_FIX ) );
michael@0 198 }
michael@0 199 NLSF[ root_ix ] = (opus_int16)silk_min_32( silk_LSHIFT( (opus_int32)k, 8 ) + ffrac, silk_int16_MAX );
michael@0 200
michael@0 201 silk_assert( NLSF[ root_ix ] >= 0 );
michael@0 202
michael@0 203 root_ix++; /* Next root */
michael@0 204 if( root_ix >= d ) {
michael@0 205 /* Found all roots */
michael@0 206 break;
michael@0 207 }
michael@0 208 /* Alternate pointer to polynomial */
michael@0 209 p = PQ[ root_ix & 1 ];
michael@0 210
michael@0 211 /* Evaluate polynomial */
michael@0 212 xlo = silk_LSFCosTab_FIX_Q12[ k - 1 ]; /* Q12*/
michael@0 213 ylo = silk_LSHIFT( 1 - ( root_ix & 2 ), 12 );
michael@0 214 } else {
michael@0 215 /* Increment loop counter */
michael@0 216 k++;
michael@0 217 xlo = xhi;
michael@0 218 ylo = yhi;
michael@0 219 thr = 0;
michael@0 220
michael@0 221 if( k > LSF_COS_TAB_SZ_FIX ) {
michael@0 222 i++;
michael@0 223 if( i > MAX_ITERATIONS_A2NLSF_FIX ) {
michael@0 224 /* Set NLSFs to white spectrum and exit */
michael@0 225 NLSF[ 0 ] = (opus_int16)silk_DIV32_16( 1 << 15, d + 1 );
michael@0 226 for( k = 1; k < d; k++ ) {
michael@0 227 NLSF[ k ] = (opus_int16)silk_SMULBB( k + 1, NLSF[ 0 ] );
michael@0 228 }
michael@0 229 return;
michael@0 230 }
michael@0 231
michael@0 232 /* Error: Apply progressively more bandwidth expansion and run again */
michael@0 233 silk_bwexpander_32( a_Q16, d, 65536 - silk_SMULBB( 10 + i, i ) ); /* 10_Q16 = 0.00015*/
michael@0 234
michael@0 235 silk_A2NLSF_init( a_Q16, P, Q, dd );
michael@0 236 p = P; /* Pointer to polynomial */
michael@0 237 xlo = silk_LSFCosTab_FIX_Q12[ 0 ]; /* Q12*/
michael@0 238 ylo = silk_A2NLSF_eval_poly( p, xlo, dd );
michael@0 239 if( ylo < 0 ) {
michael@0 240 /* Set the first NLSF to zero and move on to the next */
michael@0 241 NLSF[ 0 ] = 0;
michael@0 242 p = Q; /* Pointer to polynomial */
michael@0 243 ylo = silk_A2NLSF_eval_poly( p, xlo, dd );
michael@0 244 root_ix = 1; /* Index of current root */
michael@0 245 } else {
michael@0 246 root_ix = 0; /* Index of current root */
michael@0 247 }
michael@0 248 k = 1; /* Reset loop counter */
michael@0 249 }
michael@0 250 }
michael@0 251 }
michael@0 252 }

mercurial