Line data Source code
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 : /* Conversion between prediction filter coefficients and NLSFs */
29 : /* Requires the order to be an even number */
30 : /* A piecewise linear approximation maps LSF <-> cos(LSF) */
31 : /* Therefore the result is not accurate NLSFs, but the two */
32 : /* functions are accurate inverses of each other */
33 :
34 : #ifdef HAVE_CONFIG_H
35 : #include "config.h"
36 : #endif
37 :
38 : #include "SigProc_FIX.h"
39 : #include "tables.h"
40 :
41 : /* Number of binary divisions, when not in low complexity mode */
42 : #define BIN_DIV_STEPS_A2NLSF_FIX 3 /* must be no higher than 16 - log2( LSF_COS_TAB_SZ_FIX ) */
43 : #define MAX_ITERATIONS_A2NLSF_FIX 16
44 :
45 : /* Helper function for A2NLSF(..) */
46 : /* Transforms polynomials from cos(n*f) to cos(f)^n */
47 0 : static OPUS_INLINE void silk_A2NLSF_trans_poly(
48 : opus_int32 *p, /* I/O Polynomial */
49 : const opus_int dd /* I Polynomial order (= filter order / 2 ) */
50 : )
51 : {
52 : opus_int k, n;
53 :
54 0 : for( k = 2; k <= dd; k++ ) {
55 0 : for( n = dd; n > k; n-- ) {
56 0 : p[ n - 2 ] -= p[ n ];
57 : }
58 0 : p[ k - 2 ] -= silk_LSHIFT( p[ k ], 1 );
59 : }
60 0 : }
61 : /* Helper function for A2NLSF(..) */
62 : /* Polynomial evaluation */
63 0 : static OPUS_INLINE opus_int32 silk_A2NLSF_eval_poly( /* return the polynomial evaluation, in Q16 */
64 : opus_int32 *p, /* I Polynomial, Q16 */
65 : const opus_int32 x, /* I Evaluation point, Q12 */
66 : const opus_int dd /* I Order */
67 : )
68 : {
69 : opus_int n;
70 : opus_int32 x_Q16, y32;
71 :
72 0 : y32 = p[ dd ]; /* Q16 */
73 0 : x_Q16 = silk_LSHIFT( x, 4 );
74 :
75 0 : if ( opus_likely( 8 == dd ) )
76 : {
77 0 : y32 = silk_SMLAWW( p[ 7 ], y32, x_Q16 );
78 0 : y32 = silk_SMLAWW( p[ 6 ], y32, x_Q16 );
79 0 : y32 = silk_SMLAWW( p[ 5 ], y32, x_Q16 );
80 0 : y32 = silk_SMLAWW( p[ 4 ], y32, x_Q16 );
81 0 : y32 = silk_SMLAWW( p[ 3 ], y32, x_Q16 );
82 0 : y32 = silk_SMLAWW( p[ 2 ], y32, x_Q16 );
83 0 : y32 = silk_SMLAWW( p[ 1 ], y32, x_Q16 );
84 0 : y32 = silk_SMLAWW( p[ 0 ], y32, x_Q16 );
85 : }
86 : else
87 : {
88 0 : for( n = dd - 1; n >= 0; n-- ) {
89 0 : y32 = silk_SMLAWW( p[ n ], y32, x_Q16 ); /* Q16 */
90 : }
91 : }
92 0 : return y32;
93 : }
94 :
95 0 : static OPUS_INLINE void silk_A2NLSF_init(
96 : const opus_int32 *a_Q16,
97 : opus_int32 *P,
98 : opus_int32 *Q,
99 : const opus_int dd
100 : )
101 : {
102 : opus_int k;
103 :
104 : /* Convert filter coefs to even and odd polynomials */
105 0 : P[dd] = silk_LSHIFT( 1, 16 );
106 0 : Q[dd] = silk_LSHIFT( 1, 16 );
107 0 : for( k = 0; k < dd; k++ ) {
108 0 : P[ k ] = -a_Q16[ dd - k - 1 ] - a_Q16[ dd + k ]; /* Q16 */
109 0 : Q[ k ] = -a_Q16[ dd - k - 1 ] + a_Q16[ dd + k ]; /* Q16 */
110 : }
111 :
112 : /* Divide out zeros as we have that for even filter orders, */
113 : /* z = 1 is always a root in Q, and */
114 : /* z = -1 is always a root in P */
115 0 : for( k = dd; k > 0; k-- ) {
116 0 : P[ k - 1 ] -= P[ k ];
117 0 : Q[ k - 1 ] += Q[ k ];
118 : }
119 :
120 : /* Transform polynomials from cos(n*f) to cos(f)^n */
121 0 : silk_A2NLSF_trans_poly( P, dd );
122 0 : silk_A2NLSF_trans_poly( Q, dd );
123 0 : }
124 :
125 : /* Compute Normalized Line Spectral Frequencies (NLSFs) from whitening filter coefficients */
126 : /* If not all roots are found, the a_Q16 coefficients are bandwidth expanded until convergence. */
127 0 : void silk_A2NLSF(
128 : opus_int16 *NLSF, /* O Normalized Line Spectral Frequencies in Q15 (0..2^15-1) [d] */
129 : opus_int32 *a_Q16, /* I/O Monic whitening filter coefficients in Q16 [d] */
130 : const opus_int d /* I Filter order (must be even) */
131 : )
132 : {
133 : opus_int i, k, m, dd, root_ix, ffrac;
134 : opus_int32 xlo, xhi, xmid;
135 : opus_int32 ylo, yhi, ymid, thr;
136 : opus_int32 nom, den;
137 : opus_int32 P[ SILK_MAX_ORDER_LPC / 2 + 1 ];
138 : opus_int32 Q[ SILK_MAX_ORDER_LPC / 2 + 1 ];
139 : opus_int32 *PQ[ 2 ];
140 : opus_int32 *p;
141 :
142 : /* Store pointers to array */
143 0 : PQ[ 0 ] = P;
144 0 : PQ[ 1 ] = Q;
145 :
146 0 : dd = silk_RSHIFT( d, 1 );
147 :
148 0 : silk_A2NLSF_init( a_Q16, P, Q, dd );
149 :
150 : /* Find roots, alternating between P and Q */
151 0 : p = P; /* Pointer to polynomial */
152 :
153 0 : xlo = silk_LSFCosTab_FIX_Q12[ 0 ]; /* Q12*/
154 0 : ylo = silk_A2NLSF_eval_poly( p, xlo, dd );
155 :
156 0 : if( ylo < 0 ) {
157 : /* Set the first NLSF to zero and move on to the next */
158 0 : NLSF[ 0 ] = 0;
159 0 : p = Q; /* Pointer to polynomial */
160 0 : ylo = silk_A2NLSF_eval_poly( p, xlo, dd );
161 0 : root_ix = 1; /* Index of current root */
162 : } else {
163 0 : root_ix = 0; /* Index of current root */
164 : }
165 0 : k = 1; /* Loop counter */
166 0 : i = 0; /* Counter for bandwidth expansions applied */
167 0 : thr = 0;
168 : while( 1 ) {
169 : /* Evaluate polynomial */
170 0 : xhi = silk_LSFCosTab_FIX_Q12[ k ]; /* Q12 */
171 0 : yhi = silk_A2NLSF_eval_poly( p, xhi, dd );
172 :
173 : /* Detect zero crossing */
174 0 : if( ( ylo <= 0 && yhi >= thr ) || ( ylo >= 0 && yhi <= -thr ) ) {
175 0 : if( yhi == 0 ) {
176 : /* If the root lies exactly at the end of the current */
177 : /* interval, look for the next root in the next interval */
178 0 : thr = 1;
179 : } else {
180 0 : thr = 0;
181 : }
182 : /* Binary division */
183 0 : ffrac = -256;
184 0 : for( m = 0; m < BIN_DIV_STEPS_A2NLSF_FIX; m++ ) {
185 : /* Evaluate polynomial */
186 0 : xmid = silk_RSHIFT_ROUND( xlo + xhi, 1 );
187 0 : ymid = silk_A2NLSF_eval_poly( p, xmid, dd );
188 :
189 : /* Detect zero crossing */
190 0 : if( ( ylo <= 0 && ymid >= 0 ) || ( ylo >= 0 && ymid <= 0 ) ) {
191 : /* Reduce frequency */
192 0 : xhi = xmid;
193 0 : yhi = ymid;
194 : } else {
195 : /* Increase frequency */
196 0 : xlo = xmid;
197 0 : ylo = ymid;
198 0 : ffrac = silk_ADD_RSHIFT( ffrac, 128, m );
199 : }
200 : }
201 :
202 : /* Interpolate */
203 0 : if( silk_abs( ylo ) < 65536 ) {
204 : /* Avoid dividing by zero */
205 0 : den = ylo - yhi;
206 0 : nom = silk_LSHIFT( ylo, 8 - BIN_DIV_STEPS_A2NLSF_FIX ) + silk_RSHIFT( den, 1 );
207 0 : if( den != 0 ) {
208 0 : ffrac += silk_DIV32( nom, den );
209 : }
210 : } else {
211 : /* No risk of dividing by zero because abs(ylo - yhi) >= abs(ylo) >= 65536 */
212 0 : ffrac += silk_DIV32( ylo, silk_RSHIFT( ylo - yhi, 8 - BIN_DIV_STEPS_A2NLSF_FIX ) );
213 : }
214 0 : NLSF[ root_ix ] = (opus_int16)silk_min_32( silk_LSHIFT( (opus_int32)k, 8 ) + ffrac, silk_int16_MAX );
215 :
216 0 : silk_assert( NLSF[ root_ix ] >= 0 );
217 :
218 0 : root_ix++; /* Next root */
219 0 : if( root_ix >= d ) {
220 : /* Found all roots */
221 0 : break;
222 : }
223 : /* Alternate pointer to polynomial */
224 0 : p = PQ[ root_ix & 1 ];
225 :
226 : /* Evaluate polynomial */
227 0 : xlo = silk_LSFCosTab_FIX_Q12[ k - 1 ]; /* Q12*/
228 0 : ylo = silk_LSHIFT( 1 - ( root_ix & 2 ), 12 );
229 : } else {
230 : /* Increment loop counter */
231 0 : k++;
232 0 : xlo = xhi;
233 0 : ylo = yhi;
234 0 : thr = 0;
235 :
236 0 : if( k > LSF_COS_TAB_SZ_FIX ) {
237 0 : i++;
238 0 : if( i > MAX_ITERATIONS_A2NLSF_FIX ) {
239 : /* Set NLSFs to white spectrum and exit */
240 0 : NLSF[ 0 ] = (opus_int16)silk_DIV32_16( 1 << 15, d + 1 );
241 0 : for( k = 1; k < d; k++ ) {
242 0 : NLSF[ k ] = (opus_int16)silk_ADD16( NLSF[ k-1 ], NLSF[ 0 ] );
243 : }
244 0 : return;
245 : }
246 :
247 : /* Error: Apply progressively more bandwidth expansion and run again */
248 0 : silk_bwexpander_32( a_Q16, d, 65536 - silk_LSHIFT( 1, i ) );
249 :
250 0 : silk_A2NLSF_init( a_Q16, P, Q, dd );
251 0 : p = P; /* Pointer to polynomial */
252 0 : xlo = silk_LSFCosTab_FIX_Q12[ 0 ]; /* Q12*/
253 0 : ylo = silk_A2NLSF_eval_poly( p, xlo, dd );
254 0 : if( ylo < 0 ) {
255 : /* Set the first NLSF to zero and move on to the next */
256 0 : NLSF[ 0 ] = 0;
257 0 : p = Q; /* Pointer to polynomial */
258 0 : ylo = silk_A2NLSF_eval_poly( p, xlo, dd );
259 0 : root_ix = 1; /* Index of current root */
260 : } else {
261 0 : root_ix = 0; /* Index of current root */
262 : }
263 0 : k = 1; /* Reset loop counter */
264 : }
265 : }
266 : }
267 : }
|
