Line | Count | Source (jump to first uncovered line) |
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 | 0 | #define BIN_DIV_STEPS_A2NLSF_FIX 3 /* must be no higher than 16 - log2( LSF_COS_TAB_SZ_FIX ) */ |
43 | 0 | #define MAX_ITERATIONS_A2NLSF_FIX 16 |
44 | | |
45 | | /* Helper function for A2NLSF(..) */ |
46 | | /* Transforms polynomials from cos(n*f) to cos(f)^n */ |
47 | | 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 | 0 | { |
52 | 0 | 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 | 0 | } |
58 | 0 | p[ k - 2 ] -= silk_LSHIFT( p[ k ], 1 ); |
59 | 0 | } |
60 | 0 | } |
61 | | /* Helper function for A2NLSF(..) */ |
62 | | /* Polynomial evaluation */ |
63 | | 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 | 0 | { |
69 | 0 | opus_int n; |
70 | 0 | 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 | 0 | { |
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 | 0 | } |
86 | 0 | else |
87 | 0 | { |
88 | 0 | for( n = dd - 1; n >= 0; n-- ) { |
89 | 0 | y32 = silk_SMLAWW( p[ n ], y32, x_Q16 ); /* Q16 */ |
90 | 0 | } |
91 | 0 | } |
92 | 0 | return y32; |
93 | 0 | } |
94 | | |
95 | | 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 | 0 | { |
102 | 0 | 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 | 0 | } |
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 | 0 | } |
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 | | 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 | 0 | { |
133 | 0 | opus_int i, k, m, dd, root_ix, ffrac; |
134 | 0 | opus_int32 xlo, xhi, xmid; |
135 | 0 | opus_int32 ylo, yhi, ymid, thr; |
136 | 0 | opus_int32 nom, den; |
137 | 0 | opus_int32 P[ SILK_MAX_ORDER_LPC / 2 + 1 ]; |
138 | 0 | opus_int32 Q[ SILK_MAX_ORDER_LPC / 2 + 1 ]; |
139 | 0 | opus_int32 *PQ[ 2 ]; |
140 | 0 | 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 | 0 | } else { |
163 | 0 | root_ix = 0; /* Index of current root */ |
164 | 0 | } |
165 | 0 | k = 1; /* Loop counter */ |
166 | 0 | i = 0; /* Counter for bandwidth expansions applied */ |
167 | 0 | thr = 0; |
168 | 0 | 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 | 0 | } else { |
180 | 0 | thr = 0; |
181 | 0 | } |
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 | 0 | } else { |
195 | | /* Increase frequency */ |
196 | 0 | xlo = xmid; |
197 | 0 | ylo = ymid; |
198 | 0 | ffrac = silk_ADD_RSHIFT( ffrac, 128, m ); |
199 | 0 | } |
200 | 0 | } |
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 | 0 | } |
210 | 0 | } 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 | 0 | } |
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 | 0 | } |
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 | 0 | } 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 | 0 | } |
244 | 0 | return; |
245 | 0 | } |
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 | 0 | } else { |
261 | 0 | root_ix = 0; /* Index of current root */ |
262 | 0 | } |
263 | 0 | k = 1; /* Reset loop counter */ |
264 | 0 | } |
265 | 0 | } |
266 | 0 | } |
267 | 0 | } |