Mercurial > hg > sv-dependency-builds

annotate src/opus-1.3/silk/fixed/burg_modified_FIX.c @ 79:91c729825bca pa_catalina

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Update build for AUDIO_COMPONENT_FIX
author Chris Cannam
date Wed, 30 Oct 2019 12:40:34 +0000
parents 7aeed7906520
children
rev   line source
Chris@69 1 /***********************************************************************
Chris@69 2 Copyright (c) 2006-2011, Skype Limited. All rights reserved.
Chris@69 3 Redistribution and use in source and binary forms, with or without
Chris@69 4 modification, are permitted provided that the following conditions
Chris@69 5 are met:
Chris@69 6 - Redistributions of source code must retain the above copyright notice,
Chris@69 7 this list of conditions and the following disclaimer.
Chris@69 8 - Redistributions in binary form must reproduce the above copyright
Chris@69 9 notice, this list of conditions and the following disclaimer in the
Chris@69 10 documentation and/or other materials provided with the distribution.
Chris@69 11 - Neither the name of Internet Society, IETF or IETF Trust, nor the
Chris@69 12 names of specific contributors, may be used to endorse or promote
Chris@69 13 products derived from this software without specific prior written
Chris@69 14 permission.
Chris@69 15 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
Chris@69 16 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
Chris@69 17 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
Chris@69 18 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
Chris@69 19 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
Chris@69 20 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
Chris@69 21 SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
Chris@69 22 INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
Chris@69 23 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
Chris@69 24 ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
Chris@69 25 POSSIBILITY OF SUCH DAMAGE.
Chris@69 26 ***********************************************************************/
Chris@69 27
Chris@69 28 #ifdef HAVE_CONFIG_H
Chris@69 29 #include "config.h"
Chris@69 30 #endif
Chris@69 31
Chris@69 32 #include "SigProc_FIX.h"
Chris@69 33 #include "define.h"
Chris@69 34 #include "tuning_parameters.h"
Chris@69 35 #include "pitch.h"
Chris@69 36
Chris@69 37 #define MAX_FRAME_SIZE 384 /* subfr_length * nb_subfr = ( 0.005 * 16000 + 16 ) * 4 = 384 */
Chris@69 38
Chris@69 39 #define QA 25
Chris@69 40 #define N_BITS_HEAD_ROOM 3
Chris@69 41 #define MIN_RSHIFTS -16
Chris@69 42 #define MAX_RSHIFTS (32 - QA)
Chris@69 43
Chris@69 44 /* Compute reflection coefficients from input signal */
Chris@69 45 void silk_burg_modified_c(
Chris@69 46 opus_int32 *res_nrg, /* O Residual energy */
Chris@69 47 opus_int *res_nrg_Q, /* O Residual energy Q value */
Chris@69 48 opus_int32 A_Q16[], /* O Prediction coefficients (length order) */
Chris@69 49 const opus_int16 x[], /* I Input signal, length: nb_subfr * ( D + subfr_length ) */
Chris@69 50 const opus_int32 minInvGain_Q30, /* I Inverse of max prediction gain */
Chris@69 51 const opus_int subfr_length, /* I Input signal subframe length (incl. D preceding samples) */
Chris@69 52 const opus_int nb_subfr, /* I Number of subframes stacked in x */
Chris@69 53 const opus_int D, /* I Order */
Chris@69 54 int arch /* I Run-time architecture */
Chris@69 55 )
Chris@69 56 {
Chris@69 57 opus_int k, n, s, lz, rshifts, reached_max_gain;
Chris@69 58 opus_int32 C0, num, nrg, rc_Q31, invGain_Q30, Atmp_QA, Atmp1, tmp1, tmp2, x1, x2;
Chris@69 59 const opus_int16 *x_ptr;
Chris@69 60 opus_int32 C_first_row[ SILK_MAX_ORDER_LPC ];
Chris@69 61 opus_int32 C_last_row[ SILK_MAX_ORDER_LPC ];
Chris@69 62 opus_int32 Af_QA[ SILK_MAX_ORDER_LPC ];
Chris@69 63 opus_int32 CAf[ SILK_MAX_ORDER_LPC + 1 ];
Chris@69 64 opus_int32 CAb[ SILK_MAX_ORDER_LPC + 1 ];
Chris@69 65 opus_int32 xcorr[ SILK_MAX_ORDER_LPC ];
Chris@69 66 opus_int64 C0_64;
Chris@69 67
Chris@69 68 celt_assert( subfr_length * nb_subfr <= MAX_FRAME_SIZE );
Chris@69 69
Chris@69 70 /* Compute autocorrelations, added over subframes */
Chris@69 71 C0_64 = silk_inner_prod16_aligned_64( x, x, subfr_length*nb_subfr, arch );
Chris@69 72 lz = silk_CLZ64(C0_64);
Chris@69 73 rshifts = 32 + 1 + N_BITS_HEAD_ROOM - lz;
Chris@69 74 if (rshifts > MAX_RSHIFTS) rshifts = MAX_RSHIFTS;
Chris@69 75 if (rshifts < MIN_RSHIFTS) rshifts = MIN_RSHIFTS;
Chris@69 76
Chris@69 77 if (rshifts > 0) {
Chris@69 78 C0 = (opus_int32)silk_RSHIFT64(C0_64, rshifts );
Chris@69 79 } else {
Chris@69 80 C0 = silk_LSHIFT32((opus_int32)C0_64, -rshifts );
Chris@69 81 }
Chris@69 82
Chris@69 83 CAb[ 0 ] = CAf[ 0 ] = C0 + silk_SMMUL( SILK_FIX_CONST( FIND_LPC_COND_FAC, 32 ), C0 ) + 1; /* Q(-rshifts) */
Chris@69 84 silk_memset( C_first_row, 0, SILK_MAX_ORDER_LPC * sizeof( opus_int32 ) );
Chris@69 85 if( rshifts > 0 ) {
Chris@69 86 for( s = 0; s < nb_subfr; s++ ) {
Chris@69 87 x_ptr = x + s * subfr_length;
Chris@69 88 for( n = 1; n < D + 1; n++ ) {
Chris@69 89 C_first_row[ n - 1 ] += (opus_int32)silk_RSHIFT64(
Chris@69 90 silk_inner_prod16_aligned_64( x_ptr, x_ptr + n, subfr_length - n, arch ), rshifts );
Chris@69 91 }
Chris@69 92 }
Chris@69 93 } else {
Chris@69 94 for( s = 0; s < nb_subfr; s++ ) {
Chris@69 95 int i;
Chris@69 96 opus_int32 d;
Chris@69 97 x_ptr = x + s * subfr_length;
Chris@69 98 celt_pitch_xcorr(x_ptr, x_ptr + 1, xcorr, subfr_length - D, D, arch );
Chris@69 99 for( n = 1; n < D + 1; n++ ) {
Chris@69 100 for ( i = n + subfr_length - D, d = 0; i < subfr_length; i++ )
Chris@69 101 d = MAC16_16( d, x_ptr[ i ], x_ptr[ i - n ] );
Chris@69 102 xcorr[ n - 1 ] += d;
Chris@69 103 }
Chris@69 104 for( n = 1; n < D + 1; n++ ) {
Chris@69 105 C_first_row[ n - 1 ] += silk_LSHIFT32( xcorr[ n - 1 ], -rshifts );
Chris@69 106 }
Chris@69 107 }
Chris@69 108 }
Chris@69 109 silk_memcpy( C_last_row, C_first_row, SILK_MAX_ORDER_LPC * sizeof( opus_int32 ) );
Chris@69 110
Chris@69 111 /* Initialize */
Chris@69 112 CAb[ 0 ] = CAf[ 0 ] = C0 + silk_SMMUL( SILK_FIX_CONST( FIND_LPC_COND_FAC, 32 ), C0 ) + 1; /* Q(-rshifts) */
Chris@69 113
Chris@69 114 invGain_Q30 = (opus_int32)1 << 30;
Chris@69 115 reached_max_gain = 0;
Chris@69 116 for( n = 0; n < D; n++ ) {
Chris@69 117 /* Update first row of correlation matrix (without first element) */
Chris@69 118 /* Update last row of correlation matrix (without last element, stored in reversed order) */
Chris@69 119 /* Update C * Af */
Chris@69 120 /* Update C * flipud(Af) (stored in reversed order) */
Chris@69 121 if( rshifts > -2 ) {
Chris@69 122 for( s = 0; s < nb_subfr; s++ ) {
Chris@69 123 x_ptr = x + s * subfr_length;
Chris@69 124 x1 = -silk_LSHIFT32( (opus_int32)x_ptr[ n ], 16 - rshifts ); /* Q(16-rshifts) */
Chris@69 125 x2 = -silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n - 1 ], 16 - rshifts ); /* Q(16-rshifts) */
Chris@69 126 tmp1 = silk_LSHIFT32( (opus_int32)x_ptr[ n ], QA - 16 ); /* Q(QA-16) */
Chris@69 127 tmp2 = silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n - 1 ], QA - 16 ); /* Q(QA-16) */
Chris@69 128 for( k = 0; k < n; k++ ) {
Chris@69 129 C_first_row[ k ] = silk_SMLAWB( C_first_row[ k ], x1, x_ptr[ n - k - 1 ] ); /* Q( -rshifts ) */
Chris@69 130 C_last_row[ k ] = silk_SMLAWB( C_last_row[ k ], x2, x_ptr[ subfr_length - n + k ] ); /* Q( -rshifts ) */
Chris@69 131 Atmp_QA = Af_QA[ k ];
Chris@69 132 tmp1 = silk_SMLAWB( tmp1, Atmp_QA, x_ptr[ n - k - 1 ] ); /* Q(QA-16) */
Chris@69 133 tmp2 = silk_SMLAWB( tmp2, Atmp_QA, x_ptr[ subfr_length - n + k ] ); /* Q(QA-16) */
Chris@69 134 }
Chris@69 135 tmp1 = silk_LSHIFT32( -tmp1, 32 - QA - rshifts ); /* Q(16-rshifts) */
Chris@69 136 tmp2 = silk_LSHIFT32( -tmp2, 32 - QA - rshifts ); /* Q(16-rshifts) */
Chris@69 137 for( k = 0; k <= n; k++ ) {
Chris@69 138 CAf[ k ] = silk_SMLAWB( CAf[ k ], tmp1, x_ptr[ n - k ] ); /* Q( -rshift ) */
Chris@69 139 CAb[ k ] = silk_SMLAWB( CAb[ k ], tmp2, x_ptr[ subfr_length - n + k - 1 ] ); /* Q( -rshift ) */
Chris@69 140 }
Chris@69 141 }
Chris@69 142 } else {
Chris@69 143 for( s = 0; s < nb_subfr; s++ ) {
Chris@69 144 x_ptr = x + s * subfr_length;
Chris@69 145 x1 = -silk_LSHIFT32( (opus_int32)x_ptr[ n ], -rshifts ); /* Q( -rshifts ) */
Chris@69 146 x2 = -silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n - 1 ], -rshifts ); /* Q( -rshifts ) */
Chris@69 147 tmp1 = silk_LSHIFT32( (opus_int32)x_ptr[ n ], 17 ); /* Q17 */
Chris@69 148 tmp2 = silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n - 1 ], 17 ); /* Q17 */
Chris@69 149 for( k = 0; k < n; k++ ) {
Chris@69 150 C_first_row[ k ] = silk_MLA( C_first_row[ k ], x1, x_ptr[ n - k - 1 ] ); /* Q( -rshifts ) */
Chris@69 151 C_last_row[ k ] = silk_MLA( C_last_row[ k ], x2, x_ptr[ subfr_length - n + k ] ); /* Q( -rshifts ) */
Chris@69 152 Atmp1 = silk_RSHIFT_ROUND( Af_QA[ k ], QA - 17 ); /* Q17 */
Chris@69 153 /* We sometimes have get overflows in the multiplications (even beyond +/- 2^32),
Chris@69 154 but they cancel each other and the real result seems to always fit in a 32-bit
Chris@69 155 signed integer. This was determined experimentally, not theoretically (unfortunately). */
Chris@69 156 tmp1 = silk_MLA_ovflw( tmp1, x_ptr[ n - k - 1 ], Atmp1 ); /* Q17 */
Chris@69 157 tmp2 = silk_MLA_ovflw( tmp2, x_ptr[ subfr_length - n + k ], Atmp1 ); /* Q17 */
Chris@69 158 }
Chris@69 159 tmp1 = -tmp1; /* Q17 */
Chris@69 160 tmp2 = -tmp2; /* Q17 */
Chris@69 161 for( k = 0; k <= n; k++ ) {
Chris@69 162 CAf[ k ] = silk_SMLAWW( CAf[ k ], tmp1,
Chris@69 163 silk_LSHIFT32( (opus_int32)x_ptr[ n - k ], -rshifts - 1 ) ); /* Q( -rshift ) */
Chris@69 164 CAb[ k ] = silk_SMLAWW( CAb[ k ], tmp2,
Chris@69 165 silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n + k - 1 ], -rshifts - 1 ) ); /* Q( -rshift ) */
Chris@69 166 }
Chris@69 167 }
Chris@69 168 }
Chris@69 169
Chris@69 170 /* Calculate nominator and denominator for the next order reflection (parcor) coefficient */
Chris@69 171 tmp1 = C_first_row[ n ]; /* Q( -rshifts ) */
Chris@69 172 tmp2 = C_last_row[ n ]; /* Q( -rshifts ) */
Chris@69 173 num = 0; /* Q( -rshifts ) */
Chris@69 174 nrg = silk_ADD32( CAb[ 0 ], CAf[ 0 ] ); /* Q( 1-rshifts ) */
Chris@69 175 for( k = 0; k < n; k++ ) {
Chris@69 176 Atmp_QA = Af_QA[ k ];
Chris@69 177 lz = silk_CLZ32( silk_abs( Atmp_QA ) ) - 1;
Chris@69 178 lz = silk_min( 32 - QA, lz );
Chris@69 179 Atmp1 = silk_LSHIFT32( Atmp_QA, lz ); /* Q( QA + lz ) */
Chris@69 180
Chris@69 181 tmp1 = silk_ADD_LSHIFT32( tmp1, silk_SMMUL( C_last_row[ n - k - 1 ], Atmp1 ), 32 - QA - lz ); /* Q( -rshifts ) */
Chris@69 182 tmp2 = silk_ADD_LSHIFT32( tmp2, silk_SMMUL( C_first_row[ n - k - 1 ], Atmp1 ), 32 - QA - lz ); /* Q( -rshifts ) */
Chris@69 183 num = silk_ADD_LSHIFT32( num, silk_SMMUL( CAb[ n - k ], Atmp1 ), 32 - QA - lz ); /* Q( -rshifts ) */
Chris@69 184 nrg = silk_ADD_LSHIFT32( nrg, silk_SMMUL( silk_ADD32( CAb[ k + 1 ], CAf[ k + 1 ] ),
Chris@69 185 Atmp1 ), 32 - QA - lz ); /* Q( 1-rshifts ) */
Chris@69 186 }
Chris@69 187 CAf[ n + 1 ] = tmp1; /* Q( -rshifts ) */
Chris@69 188 CAb[ n + 1 ] = tmp2; /* Q( -rshifts ) */
Chris@69 189 num = silk_ADD32( num, tmp2 ); /* Q( -rshifts ) */
Chris@69 190 num = silk_LSHIFT32( -num, 1 ); /* Q( 1-rshifts ) */
Chris@69 191
Chris@69 192 /* Calculate the next order reflection (parcor) coefficient */
Chris@69 193 if( silk_abs( num ) < nrg ) {
Chris@69 194 rc_Q31 = silk_DIV32_varQ( num, nrg, 31 );
Chris@69 195 } else {
Chris@69 196 rc_Q31 = ( num > 0 ) ? silk_int32_MAX : silk_int32_MIN;
Chris@69 197 }
Chris@69 198
Chris@69 199 /* Update inverse prediction gain */
Chris@69 200 tmp1 = ( (opus_int32)1 << 30 ) - silk_SMMUL( rc_Q31, rc_Q31 );
Chris@69 201 tmp1 = silk_LSHIFT( silk_SMMUL( invGain_Q30, tmp1 ), 2 );
Chris@69 202 if( tmp1 <= minInvGain_Q30 ) {
Chris@69 203 /* Max prediction gain exceeded; set reflection coefficient such that max prediction gain is exactly hit */
Chris@69 204 tmp2 = ( (opus_int32)1 << 30 ) - silk_DIV32_varQ( minInvGain_Q30, invGain_Q30, 30 ); /* Q30 */
Chris@69 205 rc_Q31 = silk_SQRT_APPROX( tmp2 ); /* Q15 */
Chris@69 206 if( rc_Q31 > 0 ) {
Chris@69 207 /* Newton-Raphson iteration */
Chris@69 208 rc_Q31 = silk_RSHIFT32( rc_Q31 + silk_DIV32( tmp2, rc_Q31 ), 1 ); /* Q15 */
Chris@69 209 rc_Q31 = silk_LSHIFT32( rc_Q31, 16 ); /* Q31 */
Chris@69 210 if( num < 0 ) {
Chris@69 211 /* Ensure adjusted reflection coefficients has the original sign */
Chris@69 212 rc_Q31 = -rc_Q31;
Chris@69 213 }
Chris@69 214 }
Chris@69 215 invGain_Q30 = minInvGain_Q30;
Chris@69 216 reached_max_gain = 1;
Chris@69 217 } else {
Chris@69 218 invGain_Q30 = tmp1;
Chris@69 219 }
Chris@69 220
Chris@69 221 /* Update the AR coefficients */
Chris@69 222 for( k = 0; k < (n + 1) >> 1; k++ ) {
Chris@69 223 tmp1 = Af_QA[ k ]; /* QA */
Chris@69 224 tmp2 = Af_QA[ n - k - 1 ]; /* QA */
Chris@69 225 Af_QA[ k ] = silk_ADD_LSHIFT32( tmp1, silk_SMMUL( tmp2, rc_Q31 ), 1 ); /* QA */
Chris@69 226 Af_QA[ n - k - 1 ] = silk_ADD_LSHIFT32( tmp2, silk_SMMUL( tmp1, rc_Q31 ), 1 ); /* QA */
Chris@69 227 }
Chris@69 228 Af_QA[ n ] = silk_RSHIFT32( rc_Q31, 31 - QA ); /* QA */
Chris@69 229
Chris@69 230 if( reached_max_gain ) {
Chris@69 231 /* Reached max prediction gain; set remaining coefficients to zero and exit loop */
Chris@69 232 for( k = n + 1; k < D; k++ ) {
Chris@69 233 Af_QA[ k ] = 0;
Chris@69 234 }
Chris@69 235 break;
Chris@69 236 }
Chris@69 237
Chris@69 238 /* Update C * Af and C * Ab */
Chris@69 239 for( k = 0; k <= n + 1; k++ ) {
Chris@69 240 tmp1 = CAf[ k ]; /* Q( -rshifts ) */
Chris@69 241 tmp2 = CAb[ n - k + 1 ]; /* Q( -rshifts ) */
Chris@69 242 CAf[ k ] = silk_ADD_LSHIFT32( tmp1, silk_SMMUL( tmp2, rc_Q31 ), 1 ); /* Q( -rshifts ) */
Chris@69 243 CAb[ n - k + 1 ] = silk_ADD_LSHIFT32( tmp2, silk_SMMUL( tmp1, rc_Q31 ), 1 ); /* Q( -rshifts ) */
Chris@69 244 }
Chris@69 245 }
Chris@69 246
Chris@69 247 if( reached_max_gain ) {
Chris@69 248 for( k = 0; k < D; k++ ) {
Chris@69 249 /* Scale coefficients */
Chris@69 250 A_Q16[ k ] = -silk_RSHIFT_ROUND( Af_QA[ k ], QA - 16 );
Chris@69 251 }
Chris@69 252 /* Subtract energy of preceding samples from C0 */
Chris@69 253 if( rshifts > 0 ) {
Chris@69 254 for( s = 0; s < nb_subfr; s++ ) {
Chris@69 255 x_ptr = x + s * subfr_length;
Chris@69 256 C0 -= (opus_int32)silk_RSHIFT64( silk_inner_prod16_aligned_64( x_ptr, x_ptr, D, arch ), rshifts );
Chris@69 257 }
Chris@69 258 } else {
Chris@69 259 for( s = 0; s < nb_subfr; s++ ) {
Chris@69 260 x_ptr = x + s * subfr_length;
Chris@69 261 C0 -= silk_LSHIFT32( silk_inner_prod_aligned( x_ptr, x_ptr, D, arch), -rshifts);
Chris@69 262 }
Chris@69 263 }
Chris@69 264 /* Approximate residual energy */
Chris@69 265 *res_nrg = silk_LSHIFT( silk_SMMUL( invGain_Q30, C0 ), 2 );
Chris@69 266 *res_nrg_Q = -rshifts;
Chris@69 267 } else {
Chris@69 268 /* Return residual energy */
Chris@69 269 nrg = CAf[ 0 ]; /* Q( -rshifts ) */
Chris@69 270 tmp1 = (opus_int32)1 << 16; /* Q16 */
Chris@69 271 for( k = 0; k < D; k++ ) {
Chris@69 272 Atmp1 = silk_RSHIFT_ROUND( Af_QA[ k ], QA - 16 ); /* Q16 */
Chris@69 273 nrg = silk_SMLAWW( nrg, CAf[ k + 1 ], Atmp1 ); /* Q( -rshifts ) */
Chris@69 274 tmp1 = silk_SMLAWW( tmp1, Atmp1, Atmp1 ); /* Q16 */
Chris@69 275 A_Q16[ k ] = -Atmp1;
Chris@69 276 }
Chris@69 277 *res_nrg = silk_SMLAWW( nrg, silk_SMMUL( SILK_FIX_CONST( FIND_LPC_COND_FAC, 32 ), C0 ), -tmp1 );/* Q( -rshifts ) */
Chris@69 278 *res_nrg_Q = -rshifts;
Chris@69 279 }
Chris@69 280 }