The Tor Browser / file comparison
comparison: media/libopus/silk/float/burg_modified_FLP.c
media/libopus/silk/float/burg_modified_FLP.c
- changeset 0
- 6474c204b198
equal
deleted
inserted
replaced
| -1:000000000000 | 0:43abeb1a4849 |
|---|---|
| 1 /*********************************************************************** | |
| 2 Copyright (c) 2006-2011, Skype Limited. All rights reserved. | |
| 3 Redistribution and use in source and binary forms, with or without | |
| 4 modification, are permitted provided that the following conditions | |
| 5 are met: | |
| 6 - Redistributions of source code must retain the above copyright notice, | |
| 7 this list of conditions and the following disclaimer. | |
| 8 - Redistributions in binary form must reproduce the above copyright | |
| 9 notice, this list of conditions and the following disclaimer in the | |
| 10 documentation and/or other materials provided with the distribution. | |
| 11 - Neither the name of Internet Society, IETF or IETF Trust, nor the | |
| 12 names of specific contributors, may be used to endorse or promote | |
| 13 products derived from this software without specific prior written | |
| 14 permission. | |
| 15 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" | |
| 16 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE | |
| 17 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE | |
| 18 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE | |
| 19 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR | |
| 20 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF | |
| 21 SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS | |
| 22 INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN | |
| 23 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) | |
| 24 ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE | |
| 25 POSSIBILITY OF SUCH DAMAGE. | |
| 26 ***********************************************************************/ | |
| 27 | |
| 28 #ifdef HAVE_CONFIG_H | |
| 29 #include "config.h" | |
| 30 #endif | |
| 31 | |
| 32 #include "SigProc_FLP.h" | |
| 33 #include "tuning_parameters.h" | |
| 34 #include "define.h" | |
| 35 | |
| 36 #define MAX_FRAME_SIZE 384 /* subfr_length * nb_subfr = ( 0.005 * 16000 + 16 ) * 4 = 384*/ | |
| 37 | |
| 38 /* Compute reflection coefficients from input signal */ | |
| 39 silk_float silk_burg_modified_FLP( /* O returns residual energy */ | |
| 40 silk_float A[], /* O prediction coefficients (length order) */ | |
| 41 const silk_float x[], /* I input signal, length: nb_subfr*(D+L_sub) */ | |
| 42 const silk_float minInvGain, /* I minimum inverse prediction gain */ | |
| 43 const opus_int subfr_length, /* I input signal subframe length (incl. D preceding samples) */ | |
| 44 const opus_int nb_subfr, /* I number of subframes stacked in x */ | |
| 45 const opus_int D /* I order */ | |
| 46 ) | |
| 47 { | |
| 48 opus_int k, n, s, reached_max_gain; | |
| 49 double C0, invGain, num, nrg_f, nrg_b, rc, Atmp, tmp1, tmp2; | |
| 50 const silk_float *x_ptr; | |
| 51 double C_first_row[ SILK_MAX_ORDER_LPC ], C_last_row[ SILK_MAX_ORDER_LPC ]; | |
| 52 double CAf[ SILK_MAX_ORDER_LPC + 1 ], CAb[ SILK_MAX_ORDER_LPC + 1 ]; | |
| 53 double Af[ SILK_MAX_ORDER_LPC ]; | |
| 54 | |
| 55 silk_assert( subfr_length * nb_subfr <= MAX_FRAME_SIZE ); | |
| 56 | |
| 57 /* Compute autocorrelations, added over subframes */ | |
| 58 C0 = silk_energy_FLP( x, nb_subfr * subfr_length ); | |
| 59 silk_memset( C_first_row, 0, SILK_MAX_ORDER_LPC * sizeof( double ) ); | |
| 60 for( s = 0; s < nb_subfr; s++ ) { | |
| 61 x_ptr = x + s * subfr_length; | |
| 62 for( n = 1; n < D + 1; n++ ) { | |
| 63 C_first_row[ n - 1 ] += silk_inner_product_FLP( x_ptr, x_ptr + n, subfr_length - n ); | |
| 64 } | |
| 65 } | |
| 66 silk_memcpy( C_last_row, C_first_row, SILK_MAX_ORDER_LPC * sizeof( double ) ); | |
| 67 | |
| 68 /* Initialize */ | |
| 69 CAb[ 0 ] = CAf[ 0 ] = C0 + FIND_LPC_COND_FAC * C0 + 1e-9f; | |
| 70 invGain = 1.0f; | |
| 71 reached_max_gain = 0; | |
| 72 for( n = 0; n < D; n++ ) { | |
| 73 /* Update first row of correlation matrix (without first element) */ | |
| 74 /* Update last row of correlation matrix (without last element, stored in reversed order) */ | |
| 75 /* Update C * Af */ | |
| 76 /* Update C * flipud(Af) (stored in reversed order) */ | |
| 77 for( s = 0; s < nb_subfr; s++ ) { | |
| 78 x_ptr = x + s * subfr_length; | |
| 79 tmp1 = x_ptr[ n ]; | |
| 80 tmp2 = x_ptr[ subfr_length - n - 1 ]; | |
| 81 for( k = 0; k < n; k++ ) { | |
| 82 C_first_row[ k ] -= x_ptr[ n ] * x_ptr[ n - k - 1 ]; | |
| 83 C_last_row[ k ] -= x_ptr[ subfr_length - n - 1 ] * x_ptr[ subfr_length - n + k ]; | |
| 84 Atmp = Af[ k ]; | |
| 85 tmp1 += x_ptr[ n - k - 1 ] * Atmp; | |
| 86 tmp2 += x_ptr[ subfr_length - n + k ] * Atmp; | |
| 87 } | |
| 88 for( k = 0; k <= n; k++ ) { | |
| 89 CAf[ k ] -= tmp1 * x_ptr[ n - k ]; | |
| 90 CAb[ k ] -= tmp2 * x_ptr[ subfr_length - n + k - 1 ]; | |
| 91 } | |
| 92 } | |
| 93 tmp1 = C_first_row[ n ]; | |
| 94 tmp2 = C_last_row[ n ]; | |
| 95 for( k = 0; k < n; k++ ) { | |
| 96 Atmp = Af[ k ]; | |
| 97 tmp1 += C_last_row[ n - k - 1 ] * Atmp; | |
| 98 tmp2 += C_first_row[ n - k - 1 ] * Atmp; | |
| 99 } | |
| 100 CAf[ n + 1 ] = tmp1; | |
| 101 CAb[ n + 1 ] = tmp2; | |
| 102 | |
| 103 /* Calculate nominator and denominator for the next order reflection (parcor) coefficient */ | |
| 104 num = CAb[ n + 1 ]; | |
| 105 nrg_b = CAb[ 0 ]; | |
| 106 nrg_f = CAf[ 0 ]; | |
| 107 for( k = 0; k < n; k++ ) { | |
| 108 Atmp = Af[ k ]; | |
| 109 num += CAb[ n - k ] * Atmp; | |
| 110 nrg_b += CAb[ k + 1 ] * Atmp; | |
| 111 nrg_f += CAf[ k + 1 ] * Atmp; | |
| 112 } | |
| 113 silk_assert( nrg_f > 0.0 ); | |
| 114 silk_assert( nrg_b > 0.0 ); | |
| 115 | |
| 116 /* Calculate the next order reflection (parcor) coefficient */ | |
| 117 rc = -2.0 * num / ( nrg_f + nrg_b ); | |
| 118 silk_assert( rc > -1.0 && rc < 1.0 ); | |
| 119 | |
| 120 /* Update inverse prediction gain */ | |
| 121 tmp1 = invGain * ( 1.0 - rc * rc ); | |
| 122 if( tmp1 <= minInvGain ) { | |
| 123 /* Max prediction gain exceeded; set reflection coefficient such that max prediction gain is exactly hit */ | |
| 124 rc = sqrt( 1.0 - minInvGain / invGain ); | |
| 125 if( num > 0 ) { | |
| 126 /* Ensure adjusted reflection coefficients has the original sign */ | |
| 127 rc = -rc; | |
| 128 } | |
| 129 invGain = minInvGain; | |
| 130 reached_max_gain = 1; | |
| 131 } else { | |
| 132 invGain = tmp1; | |
| 133 } | |
| 134 | |
| 135 /* Update the AR coefficients */ | |
| 136 for( k = 0; k < (n + 1) >> 1; k++ ) { | |
| 137 tmp1 = Af[ k ]; | |
| 138 tmp2 = Af[ n - k - 1 ]; | |
| 139 Af[ k ] = tmp1 + rc * tmp2; | |
| 140 Af[ n - k - 1 ] = tmp2 + rc * tmp1; | |
| 141 } | |
| 142 Af[ n ] = rc; | |
| 143 | |
| 144 if( reached_max_gain ) { | |
| 145 /* Reached max prediction gain; set remaining coefficients to zero and exit loop */ | |
| 146 for( k = n + 1; k < D; k++ ) { | |
| 147 Af[ k ] = 0.0; | |
| 148 } | |
| 149 break; | |
| 150 } | |
| 151 | |
| 152 /* Update C * Af and C * Ab */ | |
| 153 for( k = 0; k <= n + 1; k++ ) { | |
| 154 tmp1 = CAf[ k ]; | |
| 155 CAf[ k ] += rc * CAb[ n - k + 1 ]; | |
| 156 CAb[ n - k + 1 ] += rc * tmp1; | |
| 157 } | |
| 158 } | |
| 159 | |
| 160 if( reached_max_gain ) { | |
| 161 /* Convert to silk_float */ | |
| 162 for( k = 0; k < D; k++ ) { | |
| 163 A[ k ] = (silk_float)( -Af[ k ] ); | |
| 164 } | |
| 165 /* Subtract energy of preceding samples from C0 */ | |
| 166 for( s = 0; s < nb_subfr; s++ ) { | |
| 167 C0 -= silk_energy_FLP( x + s * subfr_length, D ); | |
| 168 } | |
| 169 /* Approximate residual energy */ | |
| 170 nrg_f = C0 * invGain; | |
| 171 } else { | |
| 172 /* Compute residual energy and store coefficients as silk_float */ | |
| 173 nrg_f = CAf[ 0 ]; | |
| 174 tmp1 = 1.0; | |
| 175 for( k = 0; k < D; k++ ) { | |
| 176 Atmp = Af[ k ]; | |
| 177 nrg_f += CAf[ k + 1 ] * Atmp; | |
| 178 tmp1 += Atmp * Atmp; | |
| 179 A[ k ] = (silk_float)(-Atmp); | |
| 180 } | |
| 181 nrg_f -= FIND_LPC_COND_FAC * C0 * tmp1; | |
| 182 } | |
| 183 | |
| 184 /* Return residual energy */ | |
| 185 return (silk_float)nrg_f; | |
| 186 } |
