1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/media/libopus/silk/float/burg_modified_FLP.c Wed Dec 31 06:09:35 2014 +0100
1.3 @@ -0,0 +1,186 @@
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_FLP.h"
1.36 +#include "tuning_parameters.h"
1.37 +#include "define.h"
1.38 +
1.39 +#define MAX_FRAME_SIZE 384 /* subfr_length * nb_subfr = ( 0.005 * 16000 + 16 ) * 4 = 384*/
1.40 +
1.41 +/* Compute reflection coefficients from input signal */
1.42 +silk_float silk_burg_modified_FLP( /* O returns residual energy */
1.43 + silk_float A[], /* O prediction coefficients (length order) */
1.44 + const silk_float x[], /* I input signal, length: nb_subfr*(D+L_sub) */
1.45 + const silk_float minInvGain, /* I minimum inverse prediction gain */
1.46 + const opus_int subfr_length, /* I input signal subframe length (incl. D preceding samples) */
1.47 + const opus_int nb_subfr, /* I number of subframes stacked in x */
1.48 + const opus_int D /* I order */
1.49 +)
1.50 +{
1.51 + opus_int k, n, s, reached_max_gain;
1.52 + double C0, invGain, num, nrg_f, nrg_b, rc, Atmp, tmp1, tmp2;
1.53 + const silk_float *x_ptr;
1.54 + double C_first_row[ SILK_MAX_ORDER_LPC ], C_last_row[ SILK_MAX_ORDER_LPC ];
1.55 + double CAf[ SILK_MAX_ORDER_LPC + 1 ], CAb[ SILK_MAX_ORDER_LPC + 1 ];
1.56 + double Af[ SILK_MAX_ORDER_LPC ];
1.57 +
1.58 + silk_assert( subfr_length * nb_subfr <= MAX_FRAME_SIZE );
1.59 +
1.60 + /* Compute autocorrelations, added over subframes */
1.61 + C0 = silk_energy_FLP( x, nb_subfr * subfr_length );
1.62 + silk_memset( C_first_row, 0, SILK_MAX_ORDER_LPC * sizeof( double ) );
1.63 + for( s = 0; s < nb_subfr; s++ ) {
1.64 + x_ptr = x + s * subfr_length;
1.65 + for( n = 1; n < D + 1; n++ ) {
1.66 + C_first_row[ n - 1 ] += silk_inner_product_FLP( x_ptr, x_ptr + n, subfr_length - n );
1.67 + }
1.68 + }
1.69 + silk_memcpy( C_last_row, C_first_row, SILK_MAX_ORDER_LPC * sizeof( double ) );
1.70 +
1.71 + /* Initialize */
1.72 + CAb[ 0 ] = CAf[ 0 ] = C0 + FIND_LPC_COND_FAC * C0 + 1e-9f;
1.73 + invGain = 1.0f;
1.74 + reached_max_gain = 0;
1.75 + for( n = 0; n < D; n++ ) {
1.76 + /* Update first row of correlation matrix (without first element) */
1.77 + /* Update last row of correlation matrix (without last element, stored in reversed order) */
1.78 + /* Update C * Af */
1.79 + /* Update C * flipud(Af) (stored in reversed order) */
1.80 + for( s = 0; s < nb_subfr; s++ ) {
1.81 + x_ptr = x + s * subfr_length;
1.82 + tmp1 = x_ptr[ n ];
1.83 + tmp2 = x_ptr[ subfr_length - n - 1 ];
1.84 + for( k = 0; k < n; k++ ) {
1.85 + C_first_row[ k ] -= x_ptr[ n ] * x_ptr[ n - k - 1 ];
1.86 + C_last_row[ k ] -= x_ptr[ subfr_length - n - 1 ] * x_ptr[ subfr_length - n + k ];
1.87 + Atmp = Af[ k ];
1.88 + tmp1 += x_ptr[ n - k - 1 ] * Atmp;
1.89 + tmp2 += x_ptr[ subfr_length - n + k ] * Atmp;
1.90 + }
1.91 + for( k = 0; k <= n; k++ ) {
1.92 + CAf[ k ] -= tmp1 * x_ptr[ n - k ];
1.93 + CAb[ k ] -= tmp2 * x_ptr[ subfr_length - n + k - 1 ];
1.94 + }
1.95 + }
1.96 + tmp1 = C_first_row[ n ];
1.97 + tmp2 = C_last_row[ n ];
1.98 + for( k = 0; k < n; k++ ) {
1.99 + Atmp = Af[ k ];
1.100 + tmp1 += C_last_row[ n - k - 1 ] * Atmp;
1.101 + tmp2 += C_first_row[ n - k - 1 ] * Atmp;
1.102 + }
1.103 + CAf[ n + 1 ] = tmp1;
1.104 + CAb[ n + 1 ] = tmp2;
1.105 +
1.106 + /* Calculate nominator and denominator for the next order reflection (parcor) coefficient */
1.107 + num = CAb[ n + 1 ];
1.108 + nrg_b = CAb[ 0 ];
1.109 + nrg_f = CAf[ 0 ];
1.110 + for( k = 0; k < n; k++ ) {
1.111 + Atmp = Af[ k ];
1.112 + num += CAb[ n - k ] * Atmp;
1.113 + nrg_b += CAb[ k + 1 ] * Atmp;
1.114 + nrg_f += CAf[ k + 1 ] * Atmp;
1.115 + }
1.116 + silk_assert( nrg_f > 0.0 );
1.117 + silk_assert( nrg_b > 0.0 );
1.118 +
1.119 + /* Calculate the next order reflection (parcor) coefficient */
1.120 + rc = -2.0 * num / ( nrg_f + nrg_b );
1.121 + silk_assert( rc > -1.0 && rc < 1.0 );
1.122 +
1.123 + /* Update inverse prediction gain */
1.124 + tmp1 = invGain * ( 1.0 - rc * rc );
1.125 + if( tmp1 <= minInvGain ) {
1.126 + /* Max prediction gain exceeded; set reflection coefficient such that max prediction gain is exactly hit */
1.127 + rc = sqrt( 1.0 - minInvGain / invGain );
1.128 + if( num > 0 ) {
1.129 + /* Ensure adjusted reflection coefficients has the original sign */
1.130 + rc = -rc;
1.131 + }
1.132 + invGain = minInvGain;
1.133 + reached_max_gain = 1;
1.134 + } else {
1.135 + invGain = tmp1;
1.136 + }
1.137 +
1.138 + /* Update the AR coefficients */
1.139 + for( k = 0; k < (n + 1) >> 1; k++ ) {
1.140 + tmp1 = Af[ k ];
1.141 + tmp2 = Af[ n - k - 1 ];
1.142 + Af[ k ] = tmp1 + rc * tmp2;
1.143 + Af[ n - k - 1 ] = tmp2 + rc * tmp1;
1.144 + }
1.145 + Af[ n ] = rc;
1.146 +
1.147 + if( reached_max_gain ) {
1.148 + /* Reached max prediction gain; set remaining coefficients to zero and exit loop */
1.149 + for( k = n + 1; k < D; k++ ) {
1.150 + Af[ k ] = 0.0;
1.151 + }
1.152 + break;
1.153 + }
1.154 +
1.155 + /* Update C * Af and C * Ab */
1.156 + for( k = 0; k <= n + 1; k++ ) {
1.157 + tmp1 = CAf[ k ];
1.158 + CAf[ k ] += rc * CAb[ n - k + 1 ];
1.159 + CAb[ n - k + 1 ] += rc * tmp1;
1.160 + }
1.161 + }
1.162 +
1.163 + if( reached_max_gain ) {
1.164 + /* Convert to silk_float */
1.165 + for( k = 0; k < D; k++ ) {
1.166 + A[ k ] = (silk_float)( -Af[ k ] );
1.167 + }
1.168 + /* Subtract energy of preceding samples from C0 */
1.169 + for( s = 0; s < nb_subfr; s++ ) {
1.170 + C0 -= silk_energy_FLP( x + s * subfr_length, D );
1.171 + }
1.172 + /* Approximate residual energy */
1.173 + nrg_f = C0 * invGain;
1.174 + } else {
1.175 + /* Compute residual energy and store coefficients as silk_float */
1.176 + nrg_f = CAf[ 0 ];
1.177 + tmp1 = 1.0;
1.178 + for( k = 0; k < D; k++ ) {
1.179 + Atmp = Af[ k ];
1.180 + nrg_f += CAf[ k + 1 ] * Atmp;
1.181 + tmp1 += Atmp * Atmp;
1.182 + A[ k ] = (silk_float)(-Atmp);
1.183 + }
1.184 + nrg_f -= FIND_LPC_COND_FAC * C0 * tmp1;
1.185 + }
1.186 +
1.187 + /* Return residual energy */
1.188 + return (silk_float)nrg_f;
1.189 +}
