Fixing silk fixed point
[opus.git] / silk / silk_burg_modified.c
1 /***********************************************************************\r
2 Copyright (c) 2006-2011, 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 #include "silk_SigProc_FIX.h"\r
29 \r
30 #define MAX_FRAME_SIZE              384 // subfr_length * nb_subfr = ( 0.005 * 16000 + 16 ) * 4 = 384\r
31 #define MAX_NB_SUBFR                4\r
32 \r
33 #define QA                          25\r
34 #define N_BITS_HEAD_ROOM            2\r
35 #define MIN_RSHIFTS                 -16\r
36 #define MAX_RSHIFTS                 (32 - QA)\r
37 \r
38 /* Compute reflection coefficients from input signal */\r
39 void silk_burg_modified(\r
40     SKP_int32       *res_nrg,           /* O    residual energy                                                 */\r
41     SKP_int         *res_nrg_Q,         /* O    residual energy Q value                                         */\r
42     SKP_int32       A_Q16[],            /* O    prediction coefficients (length order)                          */\r
43     const SKP_int16 x[],                /* I    input signal, length: nb_subfr * ( D + subfr_length )           */\r
44     const SKP_int   subfr_length,       /* I    input signal subframe length (including D preceeding samples)   */\r
45     const SKP_int   nb_subfr,           /* I    number of subframes stacked in x                                */\r
46     const SKP_int32 WhiteNoiseFrac_Q32, /* I    fraction added to zero-lag autocorrelation                      */\r
47     const SKP_int   D                   /* I    order                                                           */\r
48 )\r
49 {\r
50     SKP_int         k, n, s, lz, rshifts, rshifts_extra;\r
51     SKP_int32       C0, num, nrg, rc_Q31, Atmp_QA, Atmp1, tmp1, tmp2, x1, x2;\r
52     const SKP_int16 *x_ptr;\r
53 \r
54     SKP_int32       C_first_row[ SILK_MAX_ORDER_LPC ];\r
55     SKP_int32       C_last_row[  SILK_MAX_ORDER_LPC ];\r
56     SKP_int32       Af_QA[       SILK_MAX_ORDER_LPC ];\r
57 \r
58     SKP_int32       CAf[ SILK_MAX_ORDER_LPC + 1 ];\r
59     SKP_int32       CAb[ SILK_MAX_ORDER_LPC + 1 ];\r
60 \r
61     SKP_assert( subfr_length * nb_subfr <= MAX_FRAME_SIZE );\r
62     SKP_assert( nb_subfr <= MAX_NB_SUBFR );\r
63 \r
64 \r
65     /* Compute autocorrelations, added over subframes */\r
66     silk_sum_sqr_shift( &C0, &rshifts, x, nb_subfr * subfr_length );\r
67     if( rshifts > MAX_RSHIFTS ) {\r
68         C0 = SKP_LSHIFT32( C0, rshifts - MAX_RSHIFTS );\r
69         SKP_assert( C0 > 0 );\r
70         rshifts = MAX_RSHIFTS;\r
71     } else {\r
72         lz = silk_CLZ32( C0 ) - 1;\r
73         rshifts_extra = N_BITS_HEAD_ROOM - lz;\r
74         if( rshifts_extra > 0 ) {\r
75             rshifts_extra = SKP_min( rshifts_extra, MAX_RSHIFTS - rshifts );\r
76             C0 = SKP_RSHIFT32( C0, rshifts_extra );\r
77         } else {\r
78             rshifts_extra = SKP_max( rshifts_extra, MIN_RSHIFTS - rshifts );\r
79             C0 = SKP_LSHIFT32( C0, -rshifts_extra );\r
80         }\r
81         rshifts += rshifts_extra;\r
82     }\r
83     SKP_memset( C_first_row, 0, SILK_MAX_ORDER_LPC * sizeof( SKP_int32 ) );\r
84     if( rshifts > 0 ) {\r
85         for( s = 0; s < nb_subfr; s++ ) {\r
86             x_ptr = x + s * subfr_length;\r
87             for( n = 1; n < D + 1; n++ ) {\r
88                 C_first_row[ n - 1 ] += (SKP_int32)SKP_RSHIFT64( \r
89                     silk_inner_prod16_aligned_64( x_ptr, x_ptr + n, subfr_length - n ), rshifts );\r
90             }\r
91         }\r
92     } else {\r
93         for( s = 0; s < nb_subfr; s++ ) {\r
94             x_ptr = x + s * subfr_length;\r
95             for( n = 1; n < D + 1; n++ ) {\r
96                 C_first_row[ n - 1 ] += SKP_LSHIFT32( \r
97                     silk_inner_prod_aligned( x_ptr, x_ptr + n, subfr_length - n ), -rshifts );\r
98             }\r
99         }\r
100     }\r
101     SKP_memcpy( C_last_row, C_first_row, SILK_MAX_ORDER_LPC * sizeof( SKP_int32 ) );\r
102     \r
103     /* Initialize */\r
104     CAb[ 0 ] = CAf[ 0 ] = C0 + SKP_SMMUL( WhiteNoiseFrac_Q32, C0 ) + 1;         // Q(-rshifts)\r
105 \r
106     for( n = 0; n < D; n++ ) {\r
107         /* Update first row of correlation matrix (without first element) */\r
108         /* Update last row of correlation matrix (without last element, stored in reversed order) */\r
109         /* Update C * Af */\r
110         /* Update C * flipud(Af) (stored in reversed order) */\r
111         if( rshifts > -2 ) {\r
112             for( s = 0; s < nb_subfr; s++ ) {\r
113                 x_ptr = x + s * subfr_length;\r
114                 x1  = -SKP_LSHIFT32( (SKP_int32)x_ptr[ n ],                    16 - rshifts );      // Q(16-rshifts)\r
115                 x2  = -SKP_LSHIFT32( (SKP_int32)x_ptr[ subfr_length - n - 1 ], 16 - rshifts );      // Q(16-rshifts)\r
116                 tmp1 = SKP_LSHIFT32( (SKP_int32)x_ptr[ n ],                    QA - 16 );           // Q(QA-16)\r
117                 tmp2 = SKP_LSHIFT32( (SKP_int32)x_ptr[ subfr_length - n - 1 ], QA - 16 );           // Q(QA-16)\r
118                 for( k = 0; k < n; k++ ) {\r
119                     C_first_row[ k ] = SKP_SMLAWB( C_first_row[ k ], x1, x_ptr[ n - k - 1 ]            ); // Q( -rshifts )\r
120                     C_last_row[ k ]  = SKP_SMLAWB( C_last_row[ k ],  x2, x_ptr[ subfr_length - n + k ] ); // Q( -rshifts )\r
121                     Atmp_QA = Af_QA[ k ];\r
122                     tmp1 = SKP_SMLAWB( tmp1, Atmp_QA, x_ptr[ n - k - 1 ]            );              // Q(QA-16)\r
123                     tmp2 = SKP_SMLAWB( tmp2, Atmp_QA, x_ptr[ subfr_length - n + k ] );              // Q(QA-16)\r
124                 }\r
125                 tmp1 = SKP_LSHIFT32( -tmp1, 32 - QA - rshifts );                                    // Q(16-rshifts)\r
126                 tmp2 = SKP_LSHIFT32( -tmp2, 32 - QA - rshifts );                                    // Q(16-rshifts)\r
127                 for( k = 0; k <= n; k++ ) {\r
128                     CAf[ k ] = SKP_SMLAWB( CAf[ k ], tmp1, x_ptr[ n - k ]                    );     // Q( -rshift )\r
129                     CAb[ k ] = SKP_SMLAWB( CAb[ k ], tmp2, x_ptr[ subfr_length - n + k - 1 ] );     // Q( -rshift )\r
130                 }\r
131             }\r
132         } else {\r
133             for( s = 0; s < nb_subfr; s++ ) {\r
134                 x_ptr = x + s * subfr_length;\r
135                 x1  = -SKP_LSHIFT32( (SKP_int32)x_ptr[ n ],                    -rshifts );          // Q( -rshifts )\r
136                 x2  = -SKP_LSHIFT32( (SKP_int32)x_ptr[ subfr_length - n - 1 ], -rshifts );          // Q( -rshifts )\r
137                 tmp1 = SKP_LSHIFT32( (SKP_int32)x_ptr[ n ],                    17 );                // Q17\r
138                 tmp2 = SKP_LSHIFT32( (SKP_int32)x_ptr[ subfr_length - n - 1 ], 17 );                // Q17\r
139                 for( k = 0; k < n; k++ ) {\r
140                     C_first_row[ k ] = SKP_MLA( C_first_row[ k ], x1, x_ptr[ n - k - 1 ]            ); // Q( -rshifts )\r
141                     C_last_row[ k ]  = SKP_MLA( C_last_row[ k ],  x2, x_ptr[ subfr_length - n + k ] ); // Q( -rshifts )\r
142                     Atmp1 = SKP_RSHIFT_ROUND( Af_QA[ k ], QA - 17 );                                // Q17\r
143                     tmp1 = SKP_MLA( tmp1, x_ptr[ n - k - 1 ],            Atmp1 );                   // Q17\r
144                     tmp2 = SKP_MLA( tmp2, x_ptr[ subfr_length - n + k ], Atmp1 );                   // Q17\r
145                 }\r
146                 tmp1 = -tmp1;                                                                       // Q17\r
147                 tmp2 = -tmp2;                                                                       // Q17\r
148                 for( k = 0; k <= n; k++ ) {\r
149                     CAf[ k ] = SKP_SMLAWW( CAf[ k ], tmp1, \r
150                         SKP_LSHIFT32( (SKP_int32)x_ptr[ n - k ], -rshifts - 1 ) );                  // Q( -rshift )\r
151                     CAb[ k ] = SKP_SMLAWW( CAb[ k ], tmp2, \r
152                         SKP_LSHIFT32( (SKP_int32)x_ptr[ subfr_length - n + k - 1 ], -rshifts - 1 ) );// Q( -rshift )\r
153                 }\r
154             }\r
155         }\r
156 \r
157         /* Calculate nominator and denominator for the next order reflection (parcor) coefficient */\r
158         tmp1 = C_first_row[ n ];                                                            // Q( -rshifts )\r
159         tmp2 = C_last_row[ n ];                                                             // Q( -rshifts )\r
160         num  = 0;                                                                           // Q( -rshifts )\r
161         nrg  = SKP_ADD32( CAb[ 0 ], CAf[ 0 ] );                                             // Q( 1-rshifts )\r
162         for( k = 0; k < n; k++ ) {\r
163             Atmp_QA = Af_QA[ k ];\r
164             lz = silk_CLZ32( SKP_abs( Atmp_QA ) ) - 1;\r
165             lz = SKP_min( 32 - QA, lz );\r
166             Atmp1 = SKP_LSHIFT32( Atmp_QA, lz );                                            // Q( QA + lz )\r
167 \r
168             tmp1 = SKP_ADD_LSHIFT32( tmp1, SKP_SMMUL( C_last_row[  n - k - 1 ], Atmp1 ), 32 - QA - lz );    // Q( -rshifts )\r
169             tmp2 = SKP_ADD_LSHIFT32( tmp2, SKP_SMMUL( C_first_row[ n - k - 1 ], Atmp1 ), 32 - QA - lz );    // Q( -rshifts )\r
170             num  = SKP_ADD_LSHIFT32( num,  SKP_SMMUL( CAb[ n - k ],             Atmp1 ), 32 - QA - lz );    // Q( -rshifts )\r
171             nrg  = SKP_ADD_LSHIFT32( nrg,  SKP_SMMUL( SKP_ADD32( CAb[ k + 1 ], CAf[ k + 1 ] ), \r
172                                                                                 Atmp1 ), 32 - QA - lz );    // Q( 1-rshifts )\r
173         }\r
174         CAf[ n + 1 ] = tmp1;                                                                // Q( -rshifts )\r
175         CAb[ n + 1 ] = tmp2;                                                                // Q( -rshifts )\r
176         num = SKP_ADD32( num, tmp2 );                                                       // Q( -rshifts )\r
177         num = SKP_LSHIFT32( -num, 1 );                                                      // Q( 1-rshifts )\r
178 \r
179         /* Calculate the next order reflection (parcor) coefficient */\r
180         if( SKP_abs( num ) < nrg ) {\r
181             rc_Q31 = silk_DIV32_varQ( num, nrg, 31 );\r
182         } else {\r
183             /* Negative energy or ratio too high; set remaining coefficients to zero and exit loop */\r
184             SKP_memset( &Af_QA[ n ], 0, ( D - n ) * sizeof( SKP_int32 ) );\r
185             SKP_assert( 0 );\r
186             break;\r
187         }\r
188 \r
189         /* Update the AR coefficients */\r
190         for( k = 0; k < (n + 1) >> 1; k++ ) {\r
191             tmp1 = Af_QA[ k ];                                                              // QA\r
192             tmp2 = Af_QA[ n - k - 1 ];                                                      // QA\r
193             Af_QA[ k ]         = SKP_ADD_LSHIFT32( tmp1, SKP_SMMUL( tmp2, rc_Q31 ), 1 );    // QA\r
194             Af_QA[ n - k - 1 ] = SKP_ADD_LSHIFT32( tmp2, SKP_SMMUL( tmp1, rc_Q31 ), 1 );    // QA\r
195         }\r
196         Af_QA[ n ] = SKP_RSHIFT32( rc_Q31, 31 - QA );                                       // QA\r
197 \r
198         /* Update C * Af and C * Ab */\r
199         for( k = 0; k <= n + 1; k++ ) {\r
200             tmp1 = CAf[ k ];                                                                // Q( -rshifts )\r
201             tmp2 = CAb[ n - k + 1 ];                                                        // Q( -rshifts )\r
202             CAf[ k ]         = SKP_ADD_LSHIFT32( tmp1, SKP_SMMUL( tmp2, rc_Q31 ), 1 );      // Q( -rshifts )\r
203             CAb[ n - k + 1 ] = SKP_ADD_LSHIFT32( tmp2, SKP_SMMUL( tmp1, rc_Q31 ), 1 );      // Q( -rshifts )\r
204         }\r
205     }\r
206 \r
207     /* Return residual energy */\r
208     nrg  = CAf[ 0 ];                                                                        // Q( -rshifts )\r
209     tmp1 = 1 << 16;                                                                         // Q16\r
210     for( k = 0; k < D; k++ ) {\r
211         Atmp1 = SKP_RSHIFT_ROUND( Af_QA[ k ], QA - 16 );                                    // Q16\r
212         nrg  = SKP_SMLAWW( nrg, CAf[ k + 1 ], Atmp1 );                                      // Q( -rshifts )\r
213         tmp1 = SKP_SMLAWW( tmp1, Atmp1, Atmp1 );                                            // Q16\r
214         A_Q16[ k ] = -Atmp1;\r
215     }\r
216     *res_nrg = SKP_SMLAWW( nrg, SKP_SMMUL( WhiteNoiseFrac_Q32, C0 ), -tmp1 );               // Q( -rshifts )\r
217     *res_nrg_Q = -rshifts;\r
218 }\r