Initial Skype commit taken from FreeSwitch, which got it from the IETF draft.
[opus.git] / src / SKP_Silk_LPC_inv_pred_gain.c
1 /***********************************************************************\r
2 Copyright (c) 2006-2010, Skype Limited. All rights reserved. \r
3 Redistribution and use in source and binary forms, with or without \r
4 modification, (subject to the limitations in the disclaimer below) \r
5 are permitted provided that the following conditions are met:\r
6 - Redistributions of source code must retain the above copyright notice,\r
7 this list of conditions and the following disclaimer.\r
8 - Redistributions in binary form must reproduce the above copyright \r
9 notice, this list of conditions and the following disclaimer in the \r
10 documentation and/or other materials provided with the distribution.\r
11 - Neither the name of Skype Limited, nor the names of specific \r
12 contributors, may be used to endorse or promote products derived from \r
13 this software without specific prior written permission.\r
14 NO EXPRESS OR IMPLIED LICENSES TO ANY PARTY'S PATENT RIGHTS ARE GRANTED \r
15 BY THIS LICENSE. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND \r
16 CONTRIBUTORS ''AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,\r
17 BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND \r
18 FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE \r
19 COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, \r
20 INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT\r
21 NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF \r
22 USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON \r
23 ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT \r
24 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE \r
25 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.\r
26 ***********************************************************************/\r
27 \r
28 /*                                                                      *\r
29  * SKP_Silk_LPC_inverse_pred_gain.c                                   *\r
30  *                                                                      *\r
31  * Compute inverse of LPC prediction gain, and                          *\r
32  * test if LPC coefficients are stable (all poles within unit circle)   *\r
33  *                                                                      *\r
34  * Copyright 2008 (c), Skype Limited                                           *\r
35  *                                                                      */\r
36 #include "SKP_Silk_SigProc_FIX.h"\r
37 #define QA          16\r
38 #define A_LIMIT     65520\r
39 \r
40 /* Compute inverse of LPC prediction gain, and                          */\r
41 /* test if LPC coefficients are stable (all poles within unit circle)   */\r
42 SKP_int SKP_Silk_LPC_inverse_pred_gain(       /* O:   Returns 1 if unstable, otherwise 0          */\r
43     SKP_int32           *invGain_Q30,           /* O:   Inverse prediction gain, Q30 energy domain  */\r
44     const SKP_int16     *A_Q12,                 /* I:   Prediction coefficients, Q12 [order]        */\r
45     const SKP_int       order                   /* I:   Prediction order                            */\r
46 )\r
47 {\r
48     SKP_int   k, n, headrm;\r
49     SKP_int32 rc_Q31, rc_mult1_Q30, rc_mult2_Q16;\r
50     SKP_int32 Atmp_QA[ 2 ][ SigProc_MAX_ORDER_LPC ], tmp_QA;\r
51     SKP_int32 *Aold_QA, *Anew_QA;\r
52 \r
53     Anew_QA = Atmp_QA[ order & 1 ];\r
54     /* Increase Q domain of the AR coefficients */\r
55     for( k = 0; k < order; k++ ) {\r
56         Anew_QA[ k ] = SKP_LSHIFT( (SKP_int32)A_Q12[ k ], QA - 12 );\r
57     }\r
58 \r
59     *invGain_Q30 = ( 1 << 30 );\r
60     for( k = order - 1; k > 0; k-- ) {\r
61         /* Check for stability */\r
62         if( ( Anew_QA[ k ] > A_LIMIT ) || ( Anew_QA[ k ] < -A_LIMIT ) ) {\r
63             return 1;\r
64         }\r
65 \r
66         /* Set RC equal to negated AR coef */\r
67         rc_Q31 = -SKP_LSHIFT( Anew_QA[ k ], 31 - QA );\r
68         \r
69         /* rc_mult1_Q30 range: [ 1 : 2^30-1 ] */\r
70         rc_mult1_Q30 = ( SKP_int32_MAX >> 1 ) - SKP_SMMUL( rc_Q31, rc_Q31 );\r
71         SKP_assert( rc_mult1_Q30 > ( 1 << 15 ) );                   /* reduce A_LIMIT if fails */\r
72         SKP_assert( rc_mult1_Q30 < ( 1 << 30 ) );\r
73 \r
74         /* rc_mult2_Q16 range: [ 2^16 : SKP_int32_MAX ] */\r
75         rc_mult2_Q16 = SKP_INVERSE32_varQ( rc_mult1_Q30, 46 );      /* 16 = 46 - 30 */\r
76 \r
77         /* Update inverse gain */\r
78         /* invGain_Q30 range: [ 0 : 2^30 ] */\r
79         *invGain_Q30 = SKP_LSHIFT( SKP_SMMUL( *invGain_Q30, rc_mult1_Q30 ), 2 );\r
80         SKP_assert( *invGain_Q30 >= 0           );\r
81         SKP_assert( *invGain_Q30 <= ( 1 << 30 ) );\r
82 \r
83         /* Swap pointers */\r
84         Aold_QA = Anew_QA;\r
85         Anew_QA = Atmp_QA[ k & 1 ];\r
86         \r
87         /* Update AR coefficient */\r
88         headrm = SKP_Silk_CLZ32( rc_mult2_Q16 ) - 1;\r
89         rc_mult2_Q16 = SKP_LSHIFT( rc_mult2_Q16, headrm );          /* Q: 16 + headrm */\r
90         for( n = 0; n < k; n++ ) {\r
91             tmp_QA = Aold_QA[ n ] - SKP_LSHIFT( SKP_SMMUL( Aold_QA[ k - n - 1 ], rc_Q31 ), 1 );\r
92             Anew_QA[ n ] = SKP_LSHIFT( SKP_SMMUL( tmp_QA, rc_mult2_Q16 ), 16 - headrm );\r
93         }\r
94     }\r
95 \r
96     /* Check for stability */\r
97     if( ( Anew_QA[ 0 ] > A_LIMIT ) || ( Anew_QA[ 0 ] < -A_LIMIT ) ) {\r
98         return 1;\r
99     }\r
100 \r
101     /* Set RC equal to negated AR coef */\r
102     rc_Q31 = -SKP_LSHIFT( Anew_QA[ 0 ], 31 - QA );\r
103 \r
104     /* Range: [ 1 : 2^30 ] */\r
105     rc_mult1_Q30 = ( SKP_int32_MAX >> 1 ) - SKP_SMMUL( rc_Q31, rc_Q31 );\r
106 \r
107     /* Update inverse gain */\r
108     /* Range: [ 0 : 2^30 ] */\r
109     *invGain_Q30 = SKP_LSHIFT( SKP_SMMUL( *invGain_Q30, rc_mult1_Q30 ), 2 );\r
110     SKP_assert( *invGain_Q30 >= 0     );\r
111     SKP_assert( *invGain_Q30 <= 1<<30 );\r
112 \r
113     return 0;\r
114 }\r
115 \r
116 /* For input in Q13 domain */\r
117 SKP_int SKP_Silk_LPC_inverse_pred_gain_Q13(   /* O:   Returns 1 if unstable, otherwise 0          */\r
118     SKP_int32           *invGain_Q30,           /* O:   Inverse prediction gain, Q30 energy domain  */\r
119     const SKP_int16     *A_Q13,                 /* I:   Prediction coefficients, Q13 [order]        */\r
120     const SKP_int       order                   /* I:   Prediction order                            */\r
121 )\r
122 {\r
123     SKP_int   k, n, headrm;\r
124     SKP_int32 rc_Q31, rc_mult1_Q30, rc_mult2_Q16;\r
125     SKP_int32 Atmp_QA[ 2 ][ SigProc_MAX_ORDER_LPC ], tmp_QA;\r
126     SKP_int32 *Aold_QA, *Anew_QA;\r
127 \r
128     Anew_QA = Atmp_QA[ order & 1 ];\r
129     /* Increase Q domain of the AR coefficients */\r
130     for( k = 0; k < order; k++ ) {\r
131         Anew_QA[ k ] = SKP_LSHIFT( (SKP_int32)A_Q13[ k ], QA - 13 );\r
132     }\r
133 \r
134     *invGain_Q30 = ( 1 << 30 );\r
135     for( k = order - 1; k > 0; k-- ) {\r
136         /* Check for stability */\r
137         if( ( Anew_QA[ k ] > A_LIMIT ) || ( Anew_QA[ k ] < -A_LIMIT ) ) {\r
138             return 1;\r
139         }\r
140 \r
141         /* Set RC equal to negated AR coef */\r
142         rc_Q31 = -SKP_LSHIFT( Anew_QA[ k ], 31 - QA );\r
143         \r
144         /* rc_mult1_Q30 range: [ 1 : 2^30-1 ] */\r
145         rc_mult1_Q30 = ( SKP_int32_MAX >> 1 ) - SKP_SMMUL( rc_Q31, rc_Q31 );\r
146         SKP_assert( rc_mult1_Q30 > ( 1 << 15 ) );                   /* reduce A_LIMIT if fails */\r
147         SKP_assert( rc_mult1_Q30 < ( 1 << 30 ) );\r
148 \r
149         /* rc_mult2_Q16 range: [ 2^16 : SKP_int32_MAX ] */\r
150         rc_mult2_Q16 = SKP_INVERSE32_varQ( rc_mult1_Q30, 46 );      /* 16 = 46 - 30 */\r
151 \r
152         /* Update inverse gain */\r
153         /* invGain_Q30 range: [ 0 : 2^30 ] */\r
154         *invGain_Q30 = SKP_LSHIFT( SKP_SMMUL( *invGain_Q30, rc_mult1_Q30 ), 2 );\r
155         SKP_assert( *invGain_Q30 >= 0     );\r
156         SKP_assert( *invGain_Q30 <= 1<<30 );\r
157 \r
158         /* Swap pointers */\r
159         Aold_QA = Anew_QA;\r
160         Anew_QA = Atmp_QA[ k & 1 ];\r
161         \r
162         /* Update AR coefficient */\r
163         headrm = SKP_Silk_CLZ32( rc_mult2_Q16 ) - 1;\r
164         rc_mult2_Q16 = SKP_LSHIFT( rc_mult2_Q16, headrm );          /* Q: 16 + headrm */\r
165         for( n = 0; n < k; n++ ) {\r
166             tmp_QA = Aold_QA[ n ] - SKP_LSHIFT( SKP_SMMUL( Aold_QA[ k - n - 1 ], rc_Q31 ), 1 );\r
167             Anew_QA[ n ] = SKP_LSHIFT( SKP_SMMUL( tmp_QA, rc_mult2_Q16 ), 16 - headrm );\r
168         }\r
169     }\r
170 \r
171     /* Check for stability */\r
172     if( ( Anew_QA[ 0 ] > A_LIMIT ) || ( Anew_QA[ 0 ] < -A_LIMIT ) ) {\r
173         return 1;\r
174     }\r
175 \r
176     /* Set RC equal to negated AR coef */\r
177     rc_Q31 = -SKP_LSHIFT( Anew_QA[ 0 ], 31 - QA );\r
178 \r
179     /* Range: [ 1 : 2^30 ] */\r
180     rc_mult1_Q30 = ( SKP_int32_MAX >> 1 ) - SKP_SMMUL( rc_Q31, rc_Q31 );\r
181 \r
182     /* Update inverse gain */\r
183     /* Range: [ 0 : 2^30 ] */\r
184     *invGain_Q30 = SKP_LSHIFT( SKP_SMMUL( *invGain_Q30, rc_mult1_Q30 ), 2 );\r
185     SKP_assert( *invGain_Q30 >= 0     );\r
186     SKP_assert( *invGain_Q30 <= 1<<30 );\r
187 \r
188     return 0;\r
189 }\r