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 <math.h>
33 #include <string.h>
34 #include <assert.h>
35
36 #include "arch.h"
37 #include "burg.h"
38
39 #define MAX_FRAME_SIZE 384 /* subfr_length * nb_subfr = ( 0.005 * 16000 + 16 ) * 4 = 384*/
40 #define SILK_MAX_ORDER_LPC 16
41 #define FIND_LPC_COND_FAC 1e-5f
42
43 /* sum of squares of a silk_float array, with result as double */
silk_energy_FLP(const float * data,int dataSize)44 static double silk_energy_FLP(
45 const float *data,
46 int dataSize
47 )
48 {
49 int i;
50 double result;
51
52 /* 4x unrolled loop */
53 result = 0.0;
54 for( i = 0; i < dataSize - 3; i += 4 ) {
55 result += data[ i + 0 ] * (double)data[ i + 0 ] +
56 data[ i + 1 ] * (double)data[ i + 1 ] +
57 data[ i + 2 ] * (double)data[ i + 2 ] +
58 data[ i + 3 ] * (double)data[ i + 3 ];
59 }
60
61 /* add any remaining products */
62 for( ; i < dataSize; i++ ) {
63 result += data[ i ] * (double)data[ i ];
64 }
65
66 assert( result >= 0.0 );
67 return result;
68 }
69
70 /* inner product of two silk_float arrays, with result as double */
silk_inner_product_FLP(const float * data1,const float * data2,int dataSize)71 static double silk_inner_product_FLP(
72 const float *data1,
73 const float *data2,
74 int dataSize
75 )
76 {
77 int i;
78 double result;
79
80 /* 4x unrolled loop */
81 result = 0.0;
82 for( i = 0; i < dataSize - 3; i += 4 ) {
83 result += data1[ i + 0 ] * (double)data2[ i + 0 ] +
84 data1[ i + 1 ] * (double)data2[ i + 1 ] +
85 data1[ i + 2 ] * (double)data2[ i + 2 ] +
86 data1[ i + 3 ] * (double)data2[ i + 3 ];
87 }
88
89 /* add any remaining products */
90 for( ; i < dataSize; i++ ) {
91 result += data1[ i ] * (double)data2[ i ];
92 }
93
94 return result;
95 }
96
97
98 /* Compute reflection coefficients from input signal */
silk_burg_analysis(float A[],const float x[],const float minInvGain,const int subfr_length,const int nb_subfr,const int D)99 float silk_burg_analysis( /* O returns residual energy */
100 float A[], /* O prediction coefficients (length order) */
101 const float x[], /* I input signal, length: nb_subfr*(D+L_sub) */
102 const float minInvGain, /* I minimum inverse prediction gain */
103 const int subfr_length, /* I input signal subframe length (incl. D preceding samples) */
104 const int nb_subfr, /* I number of subframes stacked in x */
105 const int D /* I order */
106 )
107 {
108 int k, n, s, reached_max_gain;
109 double C0, invGain, num, nrg_f, nrg_b, rc, Atmp, tmp1, tmp2;
110 const float *x_ptr;
111 double C_first_row[ SILK_MAX_ORDER_LPC ], C_last_row[ SILK_MAX_ORDER_LPC ];
112 double CAf[ SILK_MAX_ORDER_LPC + 1 ], CAb[ SILK_MAX_ORDER_LPC + 1 ];
113 double Af[ SILK_MAX_ORDER_LPC ];
114
115 assert( subfr_length * nb_subfr <= MAX_FRAME_SIZE );
116
117 /* Compute autocorrelations, added over subframes */
118 C0 = silk_energy_FLP( x, nb_subfr * subfr_length );
119 memset( C_first_row, 0, SILK_MAX_ORDER_LPC * sizeof( double ) );
120 for( s = 0; s < nb_subfr; s++ ) {
121 x_ptr = x + s * subfr_length;
122 for( n = 1; n < D + 1; n++ ) {
123 C_first_row[ n - 1 ] += silk_inner_product_FLP( x_ptr, x_ptr + n, subfr_length - n );
124 }
125 }
126 memcpy( C_last_row, C_first_row, SILK_MAX_ORDER_LPC * sizeof( double ) );
127
128 /* Initialize */
129 CAb[ 0 ] = CAf[ 0 ] = C0 + FIND_LPC_COND_FAC * C0 + 1e-9f;
130 invGain = 1.0f;
131 reached_max_gain = 0;
132 for( n = 0; n < D; n++ ) {
133 /* Update first row of correlation matrix (without first element) */
134 /* Update last row of correlation matrix (without last element, stored in reversed order) */
135 /* Update C * Af */
136 /* Update C * flipud(Af) (stored in reversed order) */
137 for( s = 0; s < nb_subfr; s++ ) {
138 x_ptr = x + s * subfr_length;
139 tmp1 = x_ptr[ n ];
140 tmp2 = x_ptr[ subfr_length - n - 1 ];
141 for( k = 0; k < n; k++ ) {
142 C_first_row[ k ] -= x_ptr[ n ] * x_ptr[ n - k - 1 ];
143 C_last_row[ k ] -= x_ptr[ subfr_length - n - 1 ] * x_ptr[ subfr_length - n + k ];
144 Atmp = Af[ k ];
145 tmp1 += x_ptr[ n - k - 1 ] * Atmp;
146 tmp2 += x_ptr[ subfr_length - n + k ] * Atmp;
147 }
148 for( k = 0; k <= n; k++ ) {
149 CAf[ k ] -= tmp1 * x_ptr[ n - k ];
150 CAb[ k ] -= tmp2 * x_ptr[ subfr_length - n + k - 1 ];
151 }
152 }
153 tmp1 = C_first_row[ n ];
154 tmp2 = C_last_row[ n ];
155 for( k = 0; k < n; k++ ) {
156 Atmp = Af[ k ];
157 tmp1 += C_last_row[ n - k - 1 ] * Atmp;
158 tmp2 += C_first_row[ n - k - 1 ] * Atmp;
159 }
160 CAf[ n + 1 ] = tmp1;
161 CAb[ n + 1 ] = tmp2;
162
163 /* Calculate nominator and denominator for the next order reflection (parcor) coefficient */
164 num = CAb[ n + 1 ];
165 nrg_b = CAb[ 0 ];
166 nrg_f = CAf[ 0 ];
167 for( k = 0; k < n; k++ ) {
168 Atmp = Af[ k ];
169 num += CAb[ n - k ] * Atmp;
170 nrg_b += CAb[ k + 1 ] * Atmp;
171 nrg_f += CAf[ k + 1 ] * Atmp;
172 }
173 assert( nrg_f > 0.0 );
174 assert( nrg_b > 0.0 );
175
176 /* Calculate the next order reflection (parcor) coefficient */
177 rc = -2.0 * num / ( nrg_f + nrg_b );
178 assert( rc > -1.0 && rc < 1.0 );
179
180 /* Update inverse prediction gain */
181 tmp1 = invGain * ( 1.0 - rc * rc );
182 if( tmp1 <= minInvGain ) {
183 /* Max prediction gain exceeded; set reflection coefficient such that max prediction gain is exactly hit */
184 rc = sqrt( 1.0 - minInvGain / invGain );
185 if( num > 0 ) {
186 /* Ensure adjusted reflection coefficients has the original sign */
187 rc = -rc;
188 }
189 invGain = minInvGain;
190 reached_max_gain = 1;
191 } else {
192 invGain = tmp1;
193 }
194
195 /* Update the AR coefficients */
196 for( k = 0; k < (n + 1) >> 1; k++ ) {
197 tmp1 = Af[ k ];
198 tmp2 = Af[ n - k - 1 ];
199 Af[ k ] = tmp1 + rc * tmp2;
200 Af[ n - k - 1 ] = tmp2 + rc * tmp1;
201 }
202 Af[ n ] = rc;
203
204 if( reached_max_gain ) {
205 /* Reached max prediction gain; set remaining coefficients to zero and exit loop */
206 for( k = n + 1; k < D; k++ ) {
207 Af[ k ] = 0.0;
208 }
209 break;
210 }
211
212 /* Update C * Af and C * Ab */
213 for( k = 0; k <= n + 1; k++ ) {
214 tmp1 = CAf[ k ];
215 CAf[ k ] += rc * CAb[ n - k + 1 ];
216 CAb[ n - k + 1 ] += rc * tmp1;
217 }
218 }
219
220 if( reached_max_gain ) {
221 /* Convert to float */
222 for( k = 0; k < D; k++ ) {
223 A[ k ] = (float)( -Af[ k ] );
224 }
225 /* Subtract energy of preceding samples from C0 */
226 for( s = 0; s < nb_subfr; s++ ) {
227 C0 -= silk_energy_FLP( x + s * subfr_length, D );
228 }
229 /* Approximate residual energy */
230 nrg_f = C0 * invGain;
231 } else {
232 /* Compute residual energy and store coefficients as float */
233 nrg_f = CAf[ 0 ];
234 tmp1 = 1.0;
235 for( k = 0; k < D; k++ ) {
236 Atmp = Af[ k ];
237 nrg_f += CAf[ k + 1 ] * Atmp;
238 tmp1 += Atmp * Atmp;
239 A[ k ] = (float)(-Atmp);
240 }
241 nrg_f -= FIND_LPC_COND_FAC * C0 * tmp1;
242 }
243
244 /* Return residual energy */
245 return MAX32(0, (float)nrg_f);
246 }
247