Alexandre Lision | 744f742 | 2013-09-25 11:39:37 -0400 | [diff] [blame] | 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 | } |