OLD | NEW |
1 /*********************************************************************** | 1 /*********************************************************************** |
2 Copyright (c) 2006-2011, Skype Limited. All rights reserved. | 2 Copyright (c) 2006-2011, Skype Limited. All rights reserved. |
3 Redistribution and use in source and binary forms, with or without | 3 Redistribution and use in source and binary forms, with or without |
4 modification, are permitted provided that the following conditions | 4 modification, are permitted provided that the following conditions |
5 are met: | 5 are met: |
6 - Redistributions of source code must retain the above copyright notice, | 6 - Redistributions of source code must retain the above copyright notice, |
7 this list of conditions and the following disclaimer. | 7 this list of conditions and the following disclaimer. |
8 - Redistributions in binary form must reproduce the above copyright | 8 - Redistributions in binary form must reproduce the above copyright |
9 notice, this list of conditions and the following disclaimer in the | 9 notice, this list of conditions and the following disclaimer in the |
10 documentation and/or other materials provided with the distribution. | 10 documentation and/or other materials provided with the distribution. |
11 - Neither the name of Internet Society, IETF or IETF Trust, nor the | 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 | 12 names of specific contributors, may be used to endorse or promote |
13 products derived from this software without specific prior written | 13 products derived from this software without specific prior written |
14 permission. | 14 permission. |
15 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” | 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 | 16 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE |
17 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE | 17 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE |
18 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE | 18 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE |
19 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR | 19 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR |
20 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF | 20 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF |
21 SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS | 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 | 22 INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN |
23 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) | 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 | 24 ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE |
25 POSSIBILITY OF SUCH DAMAGE. | 25 POSSIBILITY OF SUCH DAMAGE. |
26 ***********************************************************************/ | 26 ***********************************************************************/ |
27 | 27 |
28 #ifdef HAVE_CONFIG_H | 28 #ifdef HAVE_CONFIG_H |
29 #include "config.h" | 29 #include "config.h" |
30 #endif | 30 #endif |
31 | 31 |
32 #include "main_FIX.h" | 32 #include "main_FIX.h" |
| 33 #include "stack_alloc.h" |
33 #include "tuning_parameters.h" | 34 #include "tuning_parameters.h" |
34 | 35 |
35 /*****************************/ | 36 /*****************************/ |
36 /* Internal function headers */ | 37 /* Internal function headers */ |
37 /*****************************/ | 38 /*****************************/ |
38 | 39 |
39 typedef struct { | 40 typedef struct { |
40 opus_int32 Q36_part; | 41 opus_int32 Q36_part; |
41 opus_int32 Q48_part; | 42 opus_int32 Q48_part; |
42 } inv_D_t; | 43 } inv_D_t; |
(...skipping 29 matching lines...) Expand all Loading... |
72 ); | 73 ); |
73 | 74 |
74 /* Solves Ax = b, assuming A is symmetric */ | 75 /* Solves Ax = b, assuming A is symmetric */ |
75 void silk_solve_LDL_FIX( | 76 void silk_solve_LDL_FIX( |
76 opus_int32 *A, /* I
Pointer to symetric square matrix A
*/ | 77 opus_int32 *A, /* I
Pointer to symetric square matrix A
*/ |
77 opus_int M, /* I
Size of matrix
*/ | 78 opus_int M, /* I
Size of matrix
*/ |
78 const opus_int32 *b, /* I
Pointer to b vector
*/ | 79 const opus_int32 *b, /* I
Pointer to b vector
*/ |
79 opus_int32 *x_Q16 /* O
Pointer to x solution vector
*/ | 80 opus_int32 *x_Q16 /* O
Pointer to x solution vector
*/ |
80 ) | 81 ) |
81 { | 82 { |
82 opus_int32 L_Q16[ MAX_MATRIX_SIZE * MAX_MATRIX_SIZE ]; | 83 VARDECL( opus_int32, L_Q16 ); |
83 opus_int32 Y[ MAX_MATRIX_SIZE ]; | 84 opus_int32 Y[ MAX_MATRIX_SIZE ]; |
84 inv_D_t inv_D[ MAX_MATRIX_SIZE ]; | 85 inv_D_t inv_D[ MAX_MATRIX_SIZE ]; |
| 86 SAVE_STACK; |
85 | 87 |
86 silk_assert( M <= MAX_MATRIX_SIZE ); | 88 silk_assert( M <= MAX_MATRIX_SIZE ); |
| 89 ALLOC( L_Q16, M * M, opus_int32 ); |
87 | 90 |
88 /*************************************************** | 91 /*************************************************** |
89 Factorize A by LDL such that A = L*D*L', | 92 Factorize A by LDL such that A = L*D*L', |
90 where L is lower triangular with ones on diagonal | 93 where L is lower triangular with ones on diagonal |
91 ****************************************************/ | 94 ****************************************************/ |
92 silk_LDL_factorize_FIX( A, M, L_Q16, inv_D ); | 95 silk_LDL_factorize_FIX( A, M, L_Q16, inv_D ); |
93 | 96 |
94 /**************************************************** | 97 /**************************************************** |
95 * substitute D*L'*x = Y. ie: | 98 * substitute D*L'*x = Y. ie: |
96 L*D*L'*x = b => L*Y = b <=> Y = inv(L)*b | 99 L*D*L'*x = b => L*Y = b <=> Y = inv(L)*b |
97 ******************************************************/ | 100 ******************************************************/ |
98 silk_LS_SolveFirst_FIX( L_Q16, M, b, Y ); | 101 silk_LS_SolveFirst_FIX( L_Q16, M, b, Y ); |
99 | 102 |
100 /**************************************************** | 103 /**************************************************** |
101 D*L'*x = Y <=> L'*x = inv(D)*Y, because D is | 104 D*L'*x = Y <=> L'*x = inv(D)*Y, because D is |
102 diagonal just multiply with 1/d_i | 105 diagonal just multiply with 1/d_i |
103 ****************************************************/ | 106 ****************************************************/ |
104 silk_LS_divide_Q16_FIX( Y, inv_D, M ); | 107 silk_LS_divide_Q16_FIX( Y, inv_D, M ); |
105 | 108 |
106 /**************************************************** | 109 /**************************************************** |
107 x = inv(L') * inv(D) * Y | 110 x = inv(L') * inv(D) * Y |
108 *****************************************************/ | 111 *****************************************************/ |
109 silk_LS_SolveLast_FIX( L_Q16, M, Y, x_Q16 ); | 112 silk_LS_SolveLast_FIX( L_Q16, M, Y, x_Q16 ); |
| 113 RESTORE_STACK; |
110 } | 114 } |
111 | 115 |
112 static inline void silk_LDL_factorize_FIX( | 116 static inline void silk_LDL_factorize_FIX( |
113 opus_int32 *A, /* I/O Pointer to Symetric Square Matrix
*/ | 117 opus_int32 *A, /* I/O Pointer to Symetric Square Matrix
*/ |
114 opus_int M, /* I Size of Matrix
*/ | 118 opus_int M, /* I Size of Matrix
*/ |
115 opus_int32 *L_Q16, /* I/O Pointer to Square Upper triangular Ma
trix */ | 119 opus_int32 *L_Q16, /* I/O Pointer to Square Upper triangular Ma
trix */ |
116 inv_D_t *inv_D /* I/O Pointer to vector holding inverted di
agonal elements of D */ | 120 inv_D_t *inv_D /* I/O Pointer to vector holding inverted di
agonal elements of D */ |
117 ) | 121 ) |
118 { | 122 { |
119 opus_int i, j, k, status, loop_count; | 123 opus_int i, j, k, status, loop_count; |
(...skipping 116 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
236 | 240 |
237 for( i = M - 1; i >= 0; i-- ) { | 241 for( i = M - 1; i >= 0; i-- ) { |
238 ptr32 = matrix_adr( L_Q16, 0, i, M ); | 242 ptr32 = matrix_adr( L_Q16, 0, i, M ); |
239 tmp_32 = 0; | 243 tmp_32 = 0; |
240 for( j = M - 1; j > i; j-- ) { | 244 for( j = M - 1; j > i; j-- ) { |
241 tmp_32 = silk_SMLAWW( tmp_32, ptr32[ silk_SMULBB( j, M ) ], x_Q16[ j
] ); | 245 tmp_32 = silk_SMLAWW( tmp_32, ptr32[ silk_SMULBB( j, M ) ], x_Q16[ j
] ); |
242 } | 246 } |
243 x_Q16[ i ] = silk_SUB32( b[ i ], tmp_32 ); | 247 x_Q16[ i ] = silk_SUB32( b[ i ], tmp_32 ); |
244 } | 248 } |
245 } | 249 } |
OLD | NEW |