diff options
Diffstat (limited to 'lib/rbcodec/codecs/libopus/silk/float/burg_modified_FLP.c')
-rw-r--r-- | lib/rbcodec/codecs/libopus/silk/float/burg_modified_FLP.c | 186 |
1 files changed, 186 insertions, 0 deletions
diff --git a/lib/rbcodec/codecs/libopus/silk/float/burg_modified_FLP.c b/lib/rbcodec/codecs/libopus/silk/float/burg_modified_FLP.c new file mode 100644 index 0000000000..756b76a35b --- /dev/null +++ b/lib/rbcodec/codecs/libopus/silk/float/burg_modified_FLP.c | |||
@@ -0,0 +1,186 @@ | |||
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 | celt_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 | } | ||