|
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 #ifdef HAVE_CONFIG_H |
|
29 #include "config.h" |
|
30 #endif |
|
31 |
|
32 /********************************************************************** |
|
33 * Correlation Matrix Computations for LS estimate. |
|
34 **********************************************************************/ |
|
35 |
|
36 #include "main_FIX.h" |
|
37 |
|
38 /* Calculates correlation vector X'*t */ |
|
39 void silk_corrVector_FIX( |
|
40 const opus_int16 *x, /* I x vector [L + order - 1] used to form data matrix X */ |
|
41 const opus_int16 *t, /* I Target vector [L] */ |
|
42 const opus_int L, /* I Length of vectors */ |
|
43 const opus_int order, /* I Max lag for correlation */ |
|
44 opus_int32 *Xt, /* O Pointer to X'*t correlation vector [order] */ |
|
45 const opus_int rshifts /* I Right shifts of correlations */ |
|
46 ) |
|
47 { |
|
48 opus_int lag, i; |
|
49 const opus_int16 *ptr1, *ptr2; |
|
50 opus_int32 inner_prod; |
|
51 |
|
52 ptr1 = &x[ order - 1 ]; /* Points to first sample of column 0 of X: X[:,0] */ |
|
53 ptr2 = t; |
|
54 /* Calculate X'*t */ |
|
55 if( rshifts > 0 ) { |
|
56 /* Right shifting used */ |
|
57 for( lag = 0; lag < order; lag++ ) { |
|
58 inner_prod = 0; |
|
59 for( i = 0; i < L; i++ ) { |
|
60 inner_prod += silk_RSHIFT32( silk_SMULBB( ptr1[ i ], ptr2[i] ), rshifts ); |
|
61 } |
|
62 Xt[ lag ] = inner_prod; /* X[:,lag]'*t */ |
|
63 ptr1--; /* Go to next column of X */ |
|
64 } |
|
65 } else { |
|
66 silk_assert( rshifts == 0 ); |
|
67 for( lag = 0; lag < order; lag++ ) { |
|
68 Xt[ lag ] = silk_inner_prod_aligned( ptr1, ptr2, L ); /* X[:,lag]'*t */ |
|
69 ptr1--; /* Go to next column of X */ |
|
70 } |
|
71 } |
|
72 } |
|
73 |
|
74 /* Calculates correlation matrix X'*X */ |
|
75 void silk_corrMatrix_FIX( |
|
76 const opus_int16 *x, /* I x vector [L + order - 1] used to form data matrix X */ |
|
77 const opus_int L, /* I Length of vectors */ |
|
78 const opus_int order, /* I Max lag for correlation */ |
|
79 const opus_int head_room, /* I Desired headroom */ |
|
80 opus_int32 *XX, /* O Pointer to X'*X correlation matrix [ order x order ] */ |
|
81 opus_int *rshifts /* I/O Right shifts of correlations */ |
|
82 ) |
|
83 { |
|
84 opus_int i, j, lag, rshifts_local, head_room_rshifts; |
|
85 opus_int32 energy; |
|
86 const opus_int16 *ptr1, *ptr2; |
|
87 |
|
88 /* Calculate energy to find shift used to fit in 32 bits */ |
|
89 silk_sum_sqr_shift( &energy, &rshifts_local, x, L + order - 1 ); |
|
90 /* Add shifts to get the desired head room */ |
|
91 head_room_rshifts = silk_max( head_room - silk_CLZ32( energy ), 0 ); |
|
92 |
|
93 energy = silk_RSHIFT32( energy, head_room_rshifts ); |
|
94 rshifts_local += head_room_rshifts; |
|
95 |
|
96 /* Calculate energy of first column (0) of X: X[:,0]'*X[:,0] */ |
|
97 /* Remove contribution of first order - 1 samples */ |
|
98 for( i = 0; i < order - 1; i++ ) { |
|
99 energy -= silk_RSHIFT32( silk_SMULBB( x[ i ], x[ i ] ), rshifts_local ); |
|
100 } |
|
101 if( rshifts_local < *rshifts ) { |
|
102 /* Adjust energy */ |
|
103 energy = silk_RSHIFT32( energy, *rshifts - rshifts_local ); |
|
104 rshifts_local = *rshifts; |
|
105 } |
|
106 |
|
107 /* Calculate energy of remaining columns of X: X[:,j]'*X[:,j] */ |
|
108 /* Fill out the diagonal of the correlation matrix */ |
|
109 matrix_ptr( XX, 0, 0, order ) = energy; |
|
110 ptr1 = &x[ order - 1 ]; /* First sample of column 0 of X */ |
|
111 for( j = 1; j < order; j++ ) { |
|
112 energy = silk_SUB32( energy, silk_RSHIFT32( silk_SMULBB( ptr1[ L - j ], ptr1[ L - j ] ), rshifts_local ) ); |
|
113 energy = silk_ADD32( energy, silk_RSHIFT32( silk_SMULBB( ptr1[ -j ], ptr1[ -j ] ), rshifts_local ) ); |
|
114 matrix_ptr( XX, j, j, order ) = energy; |
|
115 } |
|
116 |
|
117 ptr2 = &x[ order - 2 ]; /* First sample of column 1 of X */ |
|
118 /* Calculate the remaining elements of the correlation matrix */ |
|
119 if( rshifts_local > 0 ) { |
|
120 /* Right shifting used */ |
|
121 for( lag = 1; lag < order; lag++ ) { |
|
122 /* Inner product of column 0 and column lag: X[:,0]'*X[:,lag] */ |
|
123 energy = 0; |
|
124 for( i = 0; i < L; i++ ) { |
|
125 energy += silk_RSHIFT32( silk_SMULBB( ptr1[ i ], ptr2[i] ), rshifts_local ); |
|
126 } |
|
127 /* Calculate remaining off diagonal: X[:,j]'*X[:,j + lag] */ |
|
128 matrix_ptr( XX, lag, 0, order ) = energy; |
|
129 matrix_ptr( XX, 0, lag, order ) = energy; |
|
130 for( j = 1; j < ( order - lag ); j++ ) { |
|
131 energy = silk_SUB32( energy, silk_RSHIFT32( silk_SMULBB( ptr1[ L - j ], ptr2[ L - j ] ), rshifts_local ) ); |
|
132 energy = silk_ADD32( energy, silk_RSHIFT32( silk_SMULBB( ptr1[ -j ], ptr2[ -j ] ), rshifts_local ) ); |
|
133 matrix_ptr( XX, lag + j, j, order ) = energy; |
|
134 matrix_ptr( XX, j, lag + j, order ) = energy; |
|
135 } |
|
136 ptr2--; /* Update pointer to first sample of next column (lag) in X */ |
|
137 } |
|
138 } else { |
|
139 for( lag = 1; lag < order; lag++ ) { |
|
140 /* Inner product of column 0 and column lag: X[:,0]'*X[:,lag] */ |
|
141 energy = silk_inner_prod_aligned( ptr1, ptr2, L ); |
|
142 matrix_ptr( XX, lag, 0, order ) = energy; |
|
143 matrix_ptr( XX, 0, lag, order ) = energy; |
|
144 /* Calculate remaining off diagonal: X[:,j]'*X[:,j + lag] */ |
|
145 for( j = 1; j < ( order - lag ); j++ ) { |
|
146 energy = silk_SUB32( energy, silk_SMULBB( ptr1[ L - j ], ptr2[ L - j ] ) ); |
|
147 energy = silk_SMLABB( energy, ptr1[ -j ], ptr2[ -j ] ); |
|
148 matrix_ptr( XX, lag + j, j, order ) = energy; |
|
149 matrix_ptr( XX, j, lag + j, order ) = energy; |
|
150 } |
|
151 ptr2--;/* Update pointer to first sample of next column (lag) in X */ |
|
152 } |
|
153 } |
|
154 *rshifts = rshifts_local; |
|
155 } |
|
156 |