1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/media/libopus/silk/fixed/burg_modified_FIX.c Wed Dec 31 06:09:35 2014 +0100
1.3 @@ -0,0 +1,279 @@
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 +#ifdef HAVE_CONFIG_H
1.32 +#include "config.h"
1.33 +#endif
1.34 +
1.35 +#include "SigProc_FIX.h"
1.36 +#include "define.h"
1.37 +#include "tuning_parameters.h"
1.38 +#include "pitch.h"
1.39 +
1.40 +#define MAX_FRAME_SIZE 384 /* subfr_length * nb_subfr = ( 0.005 * 16000 + 16 ) * 4 = 384 */
1.41 +
1.42 +#define QA 25
1.43 +#define N_BITS_HEAD_ROOM 2
1.44 +#define MIN_RSHIFTS -16
1.45 +#define MAX_RSHIFTS (32 - QA)
1.46 +
1.47 +/* Compute reflection coefficients from input signal */
1.48 +void silk_burg_modified(
1.49 + opus_int32 *res_nrg, /* O Residual energy */
1.50 + opus_int *res_nrg_Q, /* O Residual energy Q value */
1.51 + opus_int32 A_Q16[], /* O Prediction coefficients (length order) */
1.52 + const opus_int16 x[], /* I Input signal, length: nb_subfr * ( D + subfr_length ) */
1.53 + const opus_int32 minInvGain_Q30, /* I Inverse of max prediction gain */
1.54 + const opus_int subfr_length, /* I Input signal subframe length (incl. D preceding samples) */
1.55 + const opus_int nb_subfr, /* I Number of subframes stacked in x */
1.56 + const opus_int D, /* I Order */
1.57 + int arch /* I Run-time architecture */
1.58 +)
1.59 +{
1.60 + opus_int k, n, s, lz, rshifts, rshifts_extra, reached_max_gain;
1.61 + opus_int32 C0, num, nrg, rc_Q31, invGain_Q30, Atmp_QA, Atmp1, tmp1, tmp2, x1, x2;
1.62 + const opus_int16 *x_ptr;
1.63 + opus_int32 C_first_row[ SILK_MAX_ORDER_LPC ];
1.64 + opus_int32 C_last_row[ SILK_MAX_ORDER_LPC ];
1.65 + opus_int32 Af_QA[ SILK_MAX_ORDER_LPC ];
1.66 + opus_int32 CAf[ SILK_MAX_ORDER_LPC + 1 ];
1.67 + opus_int32 CAb[ SILK_MAX_ORDER_LPC + 1 ];
1.68 + opus_int32 xcorr[ SILK_MAX_ORDER_LPC ];
1.69 +
1.70 + silk_assert( subfr_length * nb_subfr <= MAX_FRAME_SIZE );
1.71 +
1.72 + /* Compute autocorrelations, added over subframes */
1.73 + silk_sum_sqr_shift( &C0, &rshifts, x, nb_subfr * subfr_length );
1.74 + if( rshifts > MAX_RSHIFTS ) {
1.75 + C0 = silk_LSHIFT32( C0, rshifts - MAX_RSHIFTS );
1.76 + silk_assert( C0 > 0 );
1.77 + rshifts = MAX_RSHIFTS;
1.78 + } else {
1.79 + lz = silk_CLZ32( C0 ) - 1;
1.80 + rshifts_extra = N_BITS_HEAD_ROOM - lz;
1.81 + if( rshifts_extra > 0 ) {
1.82 + rshifts_extra = silk_min( rshifts_extra, MAX_RSHIFTS - rshifts );
1.83 + C0 = silk_RSHIFT32( C0, rshifts_extra );
1.84 + } else {
1.85 + rshifts_extra = silk_max( rshifts_extra, MIN_RSHIFTS - rshifts );
1.86 + C0 = silk_LSHIFT32( C0, -rshifts_extra );
1.87 + }
1.88 + rshifts += rshifts_extra;
1.89 + }
1.90 + CAb[ 0 ] = CAf[ 0 ] = C0 + silk_SMMUL( SILK_FIX_CONST( FIND_LPC_COND_FAC, 32 ), C0 ) + 1; /* Q(-rshifts) */
1.91 + silk_memset( C_first_row, 0, SILK_MAX_ORDER_LPC * sizeof( opus_int32 ) );
1.92 + if( rshifts > 0 ) {
1.93 + for( s = 0; s < nb_subfr; s++ ) {
1.94 + x_ptr = x + s * subfr_length;
1.95 + for( n = 1; n < D + 1; n++ ) {
1.96 + C_first_row[ n - 1 ] += (opus_int32)silk_RSHIFT64(
1.97 + silk_inner_prod16_aligned_64( x_ptr, x_ptr + n, subfr_length - n ), rshifts );
1.98 + }
1.99 + }
1.100 + } else {
1.101 + for( s = 0; s < nb_subfr; s++ ) {
1.102 + int i;
1.103 + opus_int32 d;
1.104 + x_ptr = x + s * subfr_length;
1.105 + celt_pitch_xcorr(x_ptr, x_ptr + 1, xcorr, subfr_length - D, D, arch );
1.106 + for( n = 1; n < D + 1; n++ ) {
1.107 + for ( i = n + subfr_length - D, d = 0; i < subfr_length; i++ )
1.108 + d = MAC16_16( d, x_ptr[ i ], x_ptr[ i - n ] );
1.109 + xcorr[ n - 1 ] += d;
1.110 + }
1.111 + for( n = 1; n < D + 1; n++ ) {
1.112 + C_first_row[ n - 1 ] += silk_LSHIFT32( xcorr[ n - 1 ], -rshifts );
1.113 + }
1.114 + }
1.115 + }
1.116 + silk_memcpy( C_last_row, C_first_row, SILK_MAX_ORDER_LPC * sizeof( opus_int32 ) );
1.117 +
1.118 + /* Initialize */
1.119 + CAb[ 0 ] = CAf[ 0 ] = C0 + silk_SMMUL( SILK_FIX_CONST( FIND_LPC_COND_FAC, 32 ), C0 ) + 1; /* Q(-rshifts) */
1.120 +
1.121 + invGain_Q30 = (opus_int32)1 << 30;
1.122 + reached_max_gain = 0;
1.123 + for( n = 0; n < D; n++ ) {
1.124 + /* Update first row of correlation matrix (without first element) */
1.125 + /* Update last row of correlation matrix (without last element, stored in reversed order) */
1.126 + /* Update C * Af */
1.127 + /* Update C * flipud(Af) (stored in reversed order) */
1.128 + if( rshifts > -2 ) {
1.129 + for( s = 0; s < nb_subfr; s++ ) {
1.130 + x_ptr = x + s * subfr_length;
1.131 + x1 = -silk_LSHIFT32( (opus_int32)x_ptr[ n ], 16 - rshifts ); /* Q(16-rshifts) */
1.132 + x2 = -silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n - 1 ], 16 - rshifts ); /* Q(16-rshifts) */
1.133 + tmp1 = silk_LSHIFT32( (opus_int32)x_ptr[ n ], QA - 16 ); /* Q(QA-16) */
1.134 + tmp2 = silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n - 1 ], QA - 16 ); /* Q(QA-16) */
1.135 + for( k = 0; k < n; k++ ) {
1.136 + C_first_row[ k ] = silk_SMLAWB( C_first_row[ k ], x1, x_ptr[ n - k - 1 ] ); /* Q( -rshifts ) */
1.137 + C_last_row[ k ] = silk_SMLAWB( C_last_row[ k ], x2, x_ptr[ subfr_length - n + k ] ); /* Q( -rshifts ) */
1.138 + Atmp_QA = Af_QA[ k ];
1.139 + tmp1 = silk_SMLAWB( tmp1, Atmp_QA, x_ptr[ n - k - 1 ] ); /* Q(QA-16) */
1.140 + tmp2 = silk_SMLAWB( tmp2, Atmp_QA, x_ptr[ subfr_length - n + k ] ); /* Q(QA-16) */
1.141 + }
1.142 + tmp1 = silk_LSHIFT32( -tmp1, 32 - QA - rshifts ); /* Q(16-rshifts) */
1.143 + tmp2 = silk_LSHIFT32( -tmp2, 32 - QA - rshifts ); /* Q(16-rshifts) */
1.144 + for( k = 0; k <= n; k++ ) {
1.145 + CAf[ k ] = silk_SMLAWB( CAf[ k ], tmp1, x_ptr[ n - k ] ); /* Q( -rshift ) */
1.146 + CAb[ k ] = silk_SMLAWB( CAb[ k ], tmp2, x_ptr[ subfr_length - n + k - 1 ] ); /* Q( -rshift ) */
1.147 + }
1.148 + }
1.149 + } else {
1.150 + for( s = 0; s < nb_subfr; s++ ) {
1.151 + x_ptr = x + s * subfr_length;
1.152 + x1 = -silk_LSHIFT32( (opus_int32)x_ptr[ n ], -rshifts ); /* Q( -rshifts ) */
1.153 + x2 = -silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n - 1 ], -rshifts ); /* Q( -rshifts ) */
1.154 + tmp1 = silk_LSHIFT32( (opus_int32)x_ptr[ n ], 17 ); /* Q17 */
1.155 + tmp2 = silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n - 1 ], 17 ); /* Q17 */
1.156 + for( k = 0; k < n; k++ ) {
1.157 + C_first_row[ k ] = silk_MLA( C_first_row[ k ], x1, x_ptr[ n - k - 1 ] ); /* Q( -rshifts ) */
1.158 + C_last_row[ k ] = silk_MLA( C_last_row[ k ], x2, x_ptr[ subfr_length - n + k ] ); /* Q( -rshifts ) */
1.159 + Atmp1 = silk_RSHIFT_ROUND( Af_QA[ k ], QA - 17 ); /* Q17 */
1.160 + tmp1 = silk_MLA( tmp1, x_ptr[ n - k - 1 ], Atmp1 ); /* Q17 */
1.161 + tmp2 = silk_MLA( tmp2, x_ptr[ subfr_length - n + k ], Atmp1 ); /* Q17 */
1.162 + }
1.163 + tmp1 = -tmp1; /* Q17 */
1.164 + tmp2 = -tmp2; /* Q17 */
1.165 + for( k = 0; k <= n; k++ ) {
1.166 + CAf[ k ] = silk_SMLAWW( CAf[ k ], tmp1,
1.167 + silk_LSHIFT32( (opus_int32)x_ptr[ n - k ], -rshifts - 1 ) ); /* Q( -rshift ) */
1.168 + CAb[ k ] = silk_SMLAWW( CAb[ k ], tmp2,
1.169 + silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n + k - 1 ], -rshifts - 1 ) ); /* Q( -rshift ) */
1.170 + }
1.171 + }
1.172 + }
1.173 +
1.174 + /* Calculate nominator and denominator for the next order reflection (parcor) coefficient */
1.175 + tmp1 = C_first_row[ n ]; /* Q( -rshifts ) */
1.176 + tmp2 = C_last_row[ n ]; /* Q( -rshifts ) */
1.177 + num = 0; /* Q( -rshifts ) */
1.178 + nrg = silk_ADD32( CAb[ 0 ], CAf[ 0 ] ); /* Q( 1-rshifts ) */
1.179 + for( k = 0; k < n; k++ ) {
1.180 + Atmp_QA = Af_QA[ k ];
1.181 + lz = silk_CLZ32( silk_abs( Atmp_QA ) ) - 1;
1.182 + lz = silk_min( 32 - QA, lz );
1.183 + Atmp1 = silk_LSHIFT32( Atmp_QA, lz ); /* Q( QA + lz ) */
1.184 +
1.185 + tmp1 = silk_ADD_LSHIFT32( tmp1, silk_SMMUL( C_last_row[ n - k - 1 ], Atmp1 ), 32 - QA - lz ); /* Q( -rshifts ) */
1.186 + tmp2 = silk_ADD_LSHIFT32( tmp2, silk_SMMUL( C_first_row[ n - k - 1 ], Atmp1 ), 32 - QA - lz ); /* Q( -rshifts ) */
1.187 + num = silk_ADD_LSHIFT32( num, silk_SMMUL( CAb[ n - k ], Atmp1 ), 32 - QA - lz ); /* Q( -rshifts ) */
1.188 + nrg = silk_ADD_LSHIFT32( nrg, silk_SMMUL( silk_ADD32( CAb[ k + 1 ], CAf[ k + 1 ] ),
1.189 + Atmp1 ), 32 - QA - lz ); /* Q( 1-rshifts ) */
1.190 + }
1.191 + CAf[ n + 1 ] = tmp1; /* Q( -rshifts ) */
1.192 + CAb[ n + 1 ] = tmp2; /* Q( -rshifts ) */
1.193 + num = silk_ADD32( num, tmp2 ); /* Q( -rshifts ) */
1.194 + num = silk_LSHIFT32( -num, 1 ); /* Q( 1-rshifts ) */
1.195 +
1.196 + /* Calculate the next order reflection (parcor) coefficient */
1.197 + if( silk_abs( num ) < nrg ) {
1.198 + rc_Q31 = silk_DIV32_varQ( num, nrg, 31 );
1.199 + } else {
1.200 + rc_Q31 = ( num > 0 ) ? silk_int32_MAX : silk_int32_MIN;
1.201 + }
1.202 +
1.203 + /* Update inverse prediction gain */
1.204 + tmp1 = ( (opus_int32)1 << 30 ) - silk_SMMUL( rc_Q31, rc_Q31 );
1.205 + tmp1 = silk_LSHIFT( silk_SMMUL( invGain_Q30, tmp1 ), 2 );
1.206 + if( tmp1 <= minInvGain_Q30 ) {
1.207 + /* Max prediction gain exceeded; set reflection coefficient such that max prediction gain is exactly hit */
1.208 + tmp2 = ( (opus_int32)1 << 30 ) - silk_DIV32_varQ( minInvGain_Q30, invGain_Q30, 30 ); /* Q30 */
1.209 + rc_Q31 = silk_SQRT_APPROX( tmp2 ); /* Q15 */
1.210 + /* Newton-Raphson iteration */
1.211 + rc_Q31 = silk_RSHIFT32( rc_Q31 + silk_DIV32( tmp2, rc_Q31 ), 1 ); /* Q15 */
1.212 + rc_Q31 = silk_LSHIFT32( rc_Q31, 16 ); /* Q31 */
1.213 + if( num < 0 ) {
1.214 + /* Ensure adjusted reflection coefficients has the original sign */
1.215 + rc_Q31 = -rc_Q31;
1.216 + }
1.217 + invGain_Q30 = minInvGain_Q30;
1.218 + reached_max_gain = 1;
1.219 + } else {
1.220 + invGain_Q30 = tmp1;
1.221 + }
1.222 +
1.223 + /* Update the AR coefficients */
1.224 + for( k = 0; k < (n + 1) >> 1; k++ ) {
1.225 + tmp1 = Af_QA[ k ]; /* QA */
1.226 + tmp2 = Af_QA[ n - k - 1 ]; /* QA */
1.227 + Af_QA[ k ] = silk_ADD_LSHIFT32( tmp1, silk_SMMUL( tmp2, rc_Q31 ), 1 ); /* QA */
1.228 + Af_QA[ n - k - 1 ] = silk_ADD_LSHIFT32( tmp2, silk_SMMUL( tmp1, rc_Q31 ), 1 ); /* QA */
1.229 + }
1.230 + Af_QA[ n ] = silk_RSHIFT32( rc_Q31, 31 - QA ); /* QA */
1.231 +
1.232 + if( reached_max_gain ) {
1.233 + /* Reached max prediction gain; set remaining coefficients to zero and exit loop */
1.234 + for( k = n + 1; k < D; k++ ) {
1.235 + Af_QA[ k ] = 0;
1.236 + }
1.237 + break;
1.238 + }
1.239 +
1.240 + /* Update C * Af and C * Ab */
1.241 + for( k = 0; k <= n + 1; k++ ) {
1.242 + tmp1 = CAf[ k ]; /* Q( -rshifts ) */
1.243 + tmp2 = CAb[ n - k + 1 ]; /* Q( -rshifts ) */
1.244 + CAf[ k ] = silk_ADD_LSHIFT32( tmp1, silk_SMMUL( tmp2, rc_Q31 ), 1 ); /* Q( -rshifts ) */
1.245 + CAb[ n - k + 1 ] = silk_ADD_LSHIFT32( tmp2, silk_SMMUL( tmp1, rc_Q31 ), 1 ); /* Q( -rshifts ) */
1.246 + }
1.247 + }
1.248 +
1.249 + if( reached_max_gain ) {
1.250 + for( k = 0; k < D; k++ ) {
1.251 + /* Scale coefficients */
1.252 + A_Q16[ k ] = -silk_RSHIFT_ROUND( Af_QA[ k ], QA - 16 );
1.253 + }
1.254 + /* Subtract energy of preceding samples from C0 */
1.255 + if( rshifts > 0 ) {
1.256 + for( s = 0; s < nb_subfr; s++ ) {
1.257 + x_ptr = x + s * subfr_length;
1.258 + C0 -= (opus_int32)silk_RSHIFT64( silk_inner_prod16_aligned_64( x_ptr, x_ptr, D ), rshifts );
1.259 + }
1.260 + } else {
1.261 + for( s = 0; s < nb_subfr; s++ ) {
1.262 + x_ptr = x + s * subfr_length;
1.263 + C0 -= silk_LSHIFT32( silk_inner_prod_aligned( x_ptr, x_ptr, D ), -rshifts );
1.264 + }
1.265 + }
1.266 + /* Approximate residual energy */
1.267 + *res_nrg = silk_LSHIFT( silk_SMMUL( invGain_Q30, C0 ), 2 );
1.268 + *res_nrg_Q = -rshifts;
1.269 + } else {
1.270 + /* Return residual energy */
1.271 + nrg = CAf[ 0 ]; /* Q( -rshifts ) */
1.272 + tmp1 = (opus_int32)1 << 16; /* Q16 */
1.273 + for( k = 0; k < D; k++ ) {
1.274 + Atmp1 = silk_RSHIFT_ROUND( Af_QA[ k ], QA - 16 ); /* Q16 */
1.275 + nrg = silk_SMLAWW( nrg, CAf[ k + 1 ], Atmp1 ); /* Q( -rshifts ) */
1.276 + tmp1 = silk_SMLAWW( tmp1, Atmp1, Atmp1 ); /* Q16 */
1.277 + A_Q16[ k ] = -Atmp1;
1.278 + }
1.279 + *res_nrg = silk_SMLAWW( nrg, silk_SMMUL( SILK_FIX_CONST( FIND_LPC_COND_FAC, 32 ), C0 ), -tmp1 );/* Q( -rshifts ) */
1.280 + *res_nrg_Q = -rshifts;
1.281 + }
1.282 +}
