Chris@69
|
1 /***********************************************************************
|
Chris@69
|
2 Copyright (c) 2006-2011, Skype Limited. All rights reserved.
|
Chris@69
|
3 Redistribution and use in source and binary forms, with or without
|
Chris@69
|
4 modification, are permitted provided that the following conditions
|
Chris@69
|
5 are met:
|
Chris@69
|
6 - Redistributions of source code must retain the above copyright notice,
|
Chris@69
|
7 this list of conditions and the following disclaimer.
|
Chris@69
|
8 - Redistributions in binary form must reproduce the above copyright
|
Chris@69
|
9 notice, this list of conditions and the following disclaimer in the
|
Chris@69
|
10 documentation and/or other materials provided with the distribution.
|
Chris@69
|
11 - Neither the name of Internet Society, IETF or IETF Trust, nor the
|
Chris@69
|
12 names of specific contributors, may be used to endorse or promote
|
Chris@69
|
13 products derived from this software without specific prior written
|
Chris@69
|
14 permission.
|
Chris@69
|
15 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
|
Chris@69
|
16 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
|
Chris@69
|
17 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
|
Chris@69
|
18 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
|
Chris@69
|
19 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
|
Chris@69
|
20 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
|
Chris@69
|
21 SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
|
Chris@69
|
22 INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
|
Chris@69
|
23 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
|
Chris@69
|
24 ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
|
Chris@69
|
25 POSSIBILITY OF SUCH DAMAGE.
|
Chris@69
|
26 ***********************************************************************/
|
Chris@69
|
27
|
Chris@69
|
28 #ifdef HAVE_CONFIG_H
|
Chris@69
|
29 #include "config.h"
|
Chris@69
|
30 #endif
|
Chris@69
|
31
|
Chris@69
|
32 /**********************************************************************
|
Chris@69
|
33 * Correlation Matrix Computations for LS estimate.
|
Chris@69
|
34 **********************************************************************/
|
Chris@69
|
35
|
Chris@69
|
36 #include "main_FIX.h"
|
Chris@69
|
37
|
Chris@69
|
38 /* Calculates correlation vector X'*t */
|
Chris@69
|
39 void silk_corrVector_FIX(
|
Chris@69
|
40 const opus_int16 *x, /* I x vector [L + order - 1] used to form data matrix X */
|
Chris@69
|
41 const opus_int16 *t, /* I Target vector [L] */
|
Chris@69
|
42 const opus_int L, /* I Length of vectors */
|
Chris@69
|
43 const opus_int order, /* I Max lag for correlation */
|
Chris@69
|
44 opus_int32 *Xt, /* O Pointer to X'*t correlation vector [order] */
|
Chris@69
|
45 const opus_int rshifts, /* I Right shifts of correlations */
|
Chris@69
|
46 int arch /* I Run-time architecture */
|
Chris@69
|
47 )
|
Chris@69
|
48 {
|
Chris@69
|
49 opus_int lag, i;
|
Chris@69
|
50 const opus_int16 *ptr1, *ptr2;
|
Chris@69
|
51 opus_int32 inner_prod;
|
Chris@69
|
52
|
Chris@69
|
53 ptr1 = &x[ order - 1 ]; /* Points to first sample of column 0 of X: X[:,0] */
|
Chris@69
|
54 ptr2 = t;
|
Chris@69
|
55 /* Calculate X'*t */
|
Chris@69
|
56 if( rshifts > 0 ) {
|
Chris@69
|
57 /* Right shifting used */
|
Chris@69
|
58 for( lag = 0; lag < order; lag++ ) {
|
Chris@69
|
59 inner_prod = 0;
|
Chris@69
|
60 for( i = 0; i < L; i++ ) {
|
Chris@69
|
61 inner_prod = silk_ADD_RSHIFT32( inner_prod, silk_SMULBB( ptr1[ i ], ptr2[i] ), rshifts );
|
Chris@69
|
62 }
|
Chris@69
|
63 Xt[ lag ] = inner_prod; /* X[:,lag]'*t */
|
Chris@69
|
64 ptr1--; /* Go to next column of X */
|
Chris@69
|
65 }
|
Chris@69
|
66 } else {
|
Chris@69
|
67 silk_assert( rshifts == 0 );
|
Chris@69
|
68 for( lag = 0; lag < order; lag++ ) {
|
Chris@69
|
69 Xt[ lag ] = silk_inner_prod_aligned( ptr1, ptr2, L, arch ); /* X[:,lag]'*t */
|
Chris@69
|
70 ptr1--; /* Go to next column of X */
|
Chris@69
|
71 }
|
Chris@69
|
72 }
|
Chris@69
|
73 }
|
Chris@69
|
74
|
Chris@69
|
75 /* Calculates correlation matrix X'*X */
|
Chris@69
|
76 void silk_corrMatrix_FIX(
|
Chris@69
|
77 const opus_int16 *x, /* I x vector [L + order - 1] used to form data matrix X */
|
Chris@69
|
78 const opus_int L, /* I Length of vectors */
|
Chris@69
|
79 const opus_int order, /* I Max lag for correlation */
|
Chris@69
|
80 opus_int32 *XX, /* O Pointer to X'*X correlation matrix [ order x order ] */
|
Chris@69
|
81 opus_int32 *nrg, /* O Energy of x vector */
|
Chris@69
|
82 opus_int *rshifts, /* O Right shifts of correlations and energy */
|
Chris@69
|
83 int arch /* I Run-time architecture */
|
Chris@69
|
84 )
|
Chris@69
|
85 {
|
Chris@69
|
86 opus_int i, j, lag;
|
Chris@69
|
87 opus_int32 energy;
|
Chris@69
|
88 const opus_int16 *ptr1, *ptr2;
|
Chris@69
|
89
|
Chris@69
|
90 /* Calculate energy to find shift used to fit in 32 bits */
|
Chris@69
|
91 silk_sum_sqr_shift( nrg, rshifts, x, L + order - 1 );
|
Chris@69
|
92 energy = *nrg;
|
Chris@69
|
93
|
Chris@69
|
94 /* Calculate energy of first column (0) of X: X[:,0]'*X[:,0] */
|
Chris@69
|
95 /* Remove contribution of first order - 1 samples */
|
Chris@69
|
96 for( i = 0; i < order - 1; i++ ) {
|
Chris@69
|
97 energy -= silk_RSHIFT32( silk_SMULBB( x[ i ], x[ i ] ), *rshifts );
|
Chris@69
|
98 }
|
Chris@69
|
99
|
Chris@69
|
100 /* Calculate energy of remaining columns of X: X[:,j]'*X[:,j] */
|
Chris@69
|
101 /* Fill out the diagonal of the correlation matrix */
|
Chris@69
|
102 matrix_ptr( XX, 0, 0, order ) = energy;
|
Chris@69
|
103 silk_assert( energy >= 0 );
|
Chris@69
|
104 ptr1 = &x[ order - 1 ]; /* First sample of column 0 of X */
|
Chris@69
|
105 for( j = 1; j < order; j++ ) {
|
Chris@69
|
106 energy = silk_SUB32( energy, silk_RSHIFT32( silk_SMULBB( ptr1[ L - j ], ptr1[ L - j ] ), *rshifts ) );
|
Chris@69
|
107 energy = silk_ADD32( energy, silk_RSHIFT32( silk_SMULBB( ptr1[ -j ], ptr1[ -j ] ), *rshifts ) );
|
Chris@69
|
108 matrix_ptr( XX, j, j, order ) = energy;
|
Chris@69
|
109 silk_assert( energy >= 0 );
|
Chris@69
|
110 }
|
Chris@69
|
111
|
Chris@69
|
112 ptr2 = &x[ order - 2 ]; /* First sample of column 1 of X */
|
Chris@69
|
113 /* Calculate the remaining elements of the correlation matrix */
|
Chris@69
|
114 if( *rshifts > 0 ) {
|
Chris@69
|
115 /* Right shifting used */
|
Chris@69
|
116 for( lag = 1; lag < order; lag++ ) {
|
Chris@69
|
117 /* Inner product of column 0 and column lag: X[:,0]'*X[:,lag] */
|
Chris@69
|
118 energy = 0;
|
Chris@69
|
119 for( i = 0; i < L; i++ ) {
|
Chris@69
|
120 energy += silk_RSHIFT32( silk_SMULBB( ptr1[ i ], ptr2[i] ), *rshifts );
|
Chris@69
|
121 }
|
Chris@69
|
122 /* Calculate remaining off diagonal: X[:,j]'*X[:,j + lag] */
|
Chris@69
|
123 matrix_ptr( XX, lag, 0, order ) = energy;
|
Chris@69
|
124 matrix_ptr( XX, 0, lag, order ) = energy;
|
Chris@69
|
125 for( j = 1; j < ( order - lag ); j++ ) {
|
Chris@69
|
126 energy = silk_SUB32( energy, silk_RSHIFT32( silk_SMULBB( ptr1[ L - j ], ptr2[ L - j ] ), *rshifts ) );
|
Chris@69
|
127 energy = silk_ADD32( energy, silk_RSHIFT32( silk_SMULBB( ptr1[ -j ], ptr2[ -j ] ), *rshifts ) );
|
Chris@69
|
128 matrix_ptr( XX, lag + j, j, order ) = energy;
|
Chris@69
|
129 matrix_ptr( XX, j, lag + j, order ) = energy;
|
Chris@69
|
130 }
|
Chris@69
|
131 ptr2--; /* Update pointer to first sample of next column (lag) in X */
|
Chris@69
|
132 }
|
Chris@69
|
133 } else {
|
Chris@69
|
134 for( lag = 1; lag < order; lag++ ) {
|
Chris@69
|
135 /* Inner product of column 0 and column lag: X[:,0]'*X[:,lag] */
|
Chris@69
|
136 energy = silk_inner_prod_aligned( ptr1, ptr2, L, arch );
|
Chris@69
|
137 matrix_ptr( XX, lag, 0, order ) = energy;
|
Chris@69
|
138 matrix_ptr( XX, 0, lag, order ) = energy;
|
Chris@69
|
139 /* Calculate remaining off diagonal: X[:,j]'*X[:,j + lag] */
|
Chris@69
|
140 for( j = 1; j < ( order - lag ); j++ ) {
|
Chris@69
|
141 energy = silk_SUB32( energy, silk_SMULBB( ptr1[ L - j ], ptr2[ L - j ] ) );
|
Chris@69
|
142 energy = silk_SMLABB( energy, ptr1[ -j ], ptr2[ -j ] );
|
Chris@69
|
143 matrix_ptr( XX, lag + j, j, order ) = energy;
|
Chris@69
|
144 matrix_ptr( XX, j, lag + j, order ) = energy;
|
Chris@69
|
145 }
|
Chris@69
|
146 ptr2--;/* Update pointer to first sample of next column (lag) in X */
|
Chris@69
|
147 }
|
Chris@69
|
148 }
|
Chris@69
|
149 }
|
Chris@69
|
150
|