• file
• revisions
• annotate
• diff
• comparison
• raw

media/libopus/silk/fixed/solve_LS_FIX.c@b8a032363ba2 (annotated)

media/libopus/silk/fixed/solve_LS_FIX.c

Thu, 22 Jan 2015 13:21:57 +0100

author

Michael Schloh von Bennewitz <michael@schloh.com>

date

Thu, 22 Jan 2015 13:21:57 +0100

branch

TOR_BUG_9701

changeset 15

b8a032363ba2

permissions

-rw-r--r--

Incorporate requested changes from Mozilla in review:
https://bugzilla.mozilla.org/show_bug.cgi?id=1123480#c6

michael@0	1	/***********************************************************************
michael@0	2	Copyright (c) 2006-2011, Skype Limited. All rights reserved.
michael@0	3	Redistribution and use in source and binary forms, with or without
michael@0	4	modification, are permitted provided that the following conditions
michael@0	5	are met:
michael@0	6	- Redistributions of source code must retain the above copyright notice,
michael@0	7	this list of conditions and the following disclaimer.
michael@0	8	- Redistributions in binary form must reproduce the above copyright
michael@0	9	notice, this list of conditions and the following disclaimer in the
michael@0	10	documentation and/or other materials provided with the distribution.
michael@0	11	- Neither the name of Internet Society, IETF or IETF Trust, nor the
michael@0	12	names of specific contributors, may be used to endorse or promote
michael@0	13	products derived from this software without specific prior written
michael@0	14	permission.
michael@0	15	THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
michael@0	16	AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
michael@0	17	IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
michael@0	18	ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
michael@0	19	LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
michael@0	20	CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
michael@0	21	SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
michael@0	22	INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
michael@0	23	CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
michael@0	24	ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
michael@0	25	POSSIBILITY OF SUCH DAMAGE.
michael@0	26	***********************************************************************/
michael@0	27
michael@0	28	#ifdef HAVE_CONFIG_H
michael@0	29	#include "config.h"
michael@0	30	#endif
michael@0	31
michael@0	32	#include "main_FIX.h"
michael@0	33	#include "stack_alloc.h"
michael@0	34	#include "tuning_parameters.h"
michael@0	35
michael@0	36	/*****************************/
michael@0	37	/* Internal function headers */
michael@0	38	/*****************************/
michael@0	39
michael@0	40	typedef struct {
michael@0	41	opus_int32 Q36_part;
michael@0	42	opus_int32 Q48_part;
michael@0	43	} inv_D_t;
michael@0	44
michael@0	45	/* Factorize square matrix A into LDL form */
michael@0	46	static OPUS_INLINE void silk_LDL_factorize_FIX(
michael@0	47	opus_int32 *A, /* I/O Pointer to Symetric Square Matrix */
michael@0	48	opus_int M, /* I Size of Matrix */
michael@0	49	opus_int32 *L_Q16, /* I/O Pointer to Square Upper triangular Matrix */
michael@0	50	inv_D_t *inv_D /* I/O Pointer to vector holding inverted diagonal elements of D */
michael@0	51	);
michael@0	52
michael@0	53	/* Solve Lx = b, when L is lower triangular and has ones on the diagonal */
michael@0	54	static OPUS_INLINE void silk_LS_SolveFirst_FIX(
michael@0	55	const opus_int32 *L_Q16, /* I Pointer to Lower Triangular Matrix */
michael@0	56	opus_int M, /* I Dim of Matrix equation */
michael@0	57	const opus_int32 *b, /* I b Vector */
michael@0	58	opus_int32 *x_Q16 /* O x Vector */
michael@0	59	);
michael@0	60
michael@0	61	/* Solve L^t*x = b, where L is lower triangular with ones on the diagonal */
michael@0	62	static OPUS_INLINE void silk_LS_SolveLast_FIX(
michael@0	63	const opus_int32 *L_Q16, /* I Pointer to Lower Triangular Matrix */
michael@0	64	const opus_int M, /* I Dim of Matrix equation */
michael@0	65	const opus_int32 *b, /* I b Vector */
michael@0	66	opus_int32 *x_Q16 /* O x Vector */
michael@0	67	);
michael@0	68
michael@0	69	static OPUS_INLINE void silk_LS_divide_Q16_FIX(
michael@0	70	opus_int32 T[], /* I/O Numenator vector */
michael@0	71	inv_D_t *inv_D, /* I 1 / D vector */
michael@0	72	opus_int M /* I dimension */
michael@0	73	);
michael@0	74
michael@0	75	/* Solves Ax = b, assuming A is symmetric */
michael@0	76	void silk_solve_LDL_FIX(
michael@0	77	opus_int32 *A, /* I Pointer to symetric square matrix A */
michael@0	78	opus_int M, /* I Size of matrix */
michael@0	79	const opus_int32 *b, /* I Pointer to b vector */
michael@0	80	opus_int32 *x_Q16 /* O Pointer to x solution vector */
michael@0	81	)
michael@0	82	{
michael@0	83	VARDECL( opus_int32, L_Q16 );
michael@0	84	opus_int32 Y[ MAX_MATRIX_SIZE ];
michael@0	85	inv_D_t inv_D[ MAX_MATRIX_SIZE ];
michael@0	86	SAVE_STACK;
michael@0	87
michael@0	88	silk_assert( M <= MAX_MATRIX_SIZE );
michael@0	89	ALLOC( L_Q16, M * M, opus_int32 );
michael@0	90
michael@0	91	/***************************************************
michael@0	92	Factorize A by LDL such that A = L*D*L',
michael@0	93	where L is lower triangular with ones on diagonal
michael@0	94	****************************************************/
michael@0	95	silk_LDL_factorize_FIX( A, M, L_Q16, inv_D );
michael@0	96
michael@0	97	/****************************************************
michael@0	98	* substitute D*L'*x = Y. ie:
michael@0	99	L*D*L'*x = b => L*Y = b <=> Y = inv(L)*b
michael@0	100	******************************************************/
michael@0	101	silk_LS_SolveFirst_FIX( L_Q16, M, b, Y );
michael@0	102
michael@0	103	/****************************************************
michael@0	104	D*L'*x = Y <=> L'*x = inv(D)*Y, because D is
michael@0	105	diagonal just multiply with 1/d_i
michael@0	106	****************************************************/
michael@0	107	silk_LS_divide_Q16_FIX( Y, inv_D, M );
michael@0	108
michael@0	109	/****************************************************
michael@0	110	x = inv(L') * inv(D) * Y
michael@0	111	*****************************************************/
michael@0	112	silk_LS_SolveLast_FIX( L_Q16, M, Y, x_Q16 );
michael@0	113	RESTORE_STACK;
michael@0	114	}
michael@0	115
michael@0	116	static OPUS_INLINE void silk_LDL_factorize_FIX(
michael@0	117	opus_int32 *A, /* I/O Pointer to Symetric Square Matrix */
michael@0	118	opus_int M, /* I Size of Matrix */
michael@0	119	opus_int32 *L_Q16, /* I/O Pointer to Square Upper triangular Matrix */
michael@0	120	inv_D_t *inv_D /* I/O Pointer to vector holding inverted diagonal elements of D */
michael@0	121	)
michael@0	122	{
michael@0	123	opus_int i, j, k, status, loop_count;
michael@0	124	const opus_int32 *ptr1, *ptr2;
michael@0	125	opus_int32 diag_min_value, tmp_32, err;
michael@0	126	opus_int32 v_Q0[ MAX_MATRIX_SIZE ], D_Q0[ MAX_MATRIX_SIZE ];
michael@0	127	opus_int32 one_div_diag_Q36, one_div_diag_Q40, one_div_diag_Q48;
michael@0	128
michael@0	129	silk_assert( M <= MAX_MATRIX_SIZE );
michael@0	130
michael@0	131	status = 1;
michael@0	132	diag_min_value = silk_max_32( silk_SMMUL( silk_ADD_SAT32( A[ 0 ], A[ silk_SMULBB( M, M ) - 1 ] ), SILK_FIX_CONST( FIND_LTP_COND_FAC, 31 ) ), 1 << 9 );
michael@0	133	for( loop_count = 0; loop_count < M && status == 1; loop_count++ ) {
michael@0	134	status = 0;
michael@0	135	for( j = 0; j < M; j++ ) {
michael@0	136	ptr1 = matrix_adr( L_Q16, j, 0, M );
michael@0	137	tmp_32 = 0;
michael@0	138	for( i = 0; i < j; i++ ) {
michael@0	139	v_Q0[ i ] = silk_SMULWW( D_Q0[ i ], ptr1[ i ] ); /* Q0 */
michael@0	140	tmp_32 = silk_SMLAWW( tmp_32, v_Q0[ i ], ptr1[ i ] ); /* Q0 */
michael@0	141	}
michael@0	142	tmp_32 = silk_SUB32( matrix_ptr( A, j, j, M ), tmp_32 );
michael@0	143
michael@0	144	if( tmp_32 < diag_min_value ) {
michael@0	145	tmp_32 = silk_SUB32( silk_SMULBB( loop_count + 1, diag_min_value ), tmp_32 );
michael@0	146	/* Matrix not positive semi-definite, or ill conditioned */
michael@0	147	for( i = 0; i < M; i++ ) {
michael@0	148	matrix_ptr( A, i, i, M ) = silk_ADD32( matrix_ptr( A, i, i, M ), tmp_32 );
michael@0	149	}
michael@0	150	status = 1;
michael@0	151	break;
michael@0	152	}
michael@0	153	D_Q0[ j ] = tmp_32; /* always < max(Correlation) */
michael@0	154
michael@0	155	/* two-step division */
michael@0	156	one_div_diag_Q36 = silk_INVERSE32_varQ( tmp_32, 36 ); /* Q36 */
michael@0	157	one_div_diag_Q40 = silk_LSHIFT( one_div_diag_Q36, 4 ); /* Q40 */
michael@0	158	err = silk_SUB32( (opus_int32)1 << 24, silk_SMULWW( tmp_32, one_div_diag_Q40 ) ); /* Q24 */
michael@0	159	one_div_diag_Q48 = silk_SMULWW( err, one_div_diag_Q40 ); /* Q48 */
michael@0	160
michael@0	161	/* Save 1/Ds */
michael@0	162	inv_D[ j ].Q36_part = one_div_diag_Q36;
michael@0	163	inv_D[ j ].Q48_part = one_div_diag_Q48;
michael@0	164
michael@0	165	matrix_ptr( L_Q16, j, j, M ) = 65536; /* 1.0 in Q16 */
michael@0	166	ptr1 = matrix_adr( A, j, 0, M );
michael@0	167	ptr2 = matrix_adr( L_Q16, j + 1, 0, M );
michael@0	168	for( i = j + 1; i < M; i++ ) {
michael@0	169	tmp_32 = 0;
michael@0	170	for( k = 0; k < j; k++ ) {
michael@0	171	tmp_32 = silk_SMLAWW( tmp_32, v_Q0[ k ], ptr2[ k ] ); /* Q0 */
michael@0	172	}
michael@0	173	tmp_32 = silk_SUB32( ptr1[ i ], tmp_32 ); /* always < max(Correlation) */
michael@0	174
michael@0	175	/* tmp_32 / D_Q0[j] : Divide to Q16 */
michael@0	176	matrix_ptr( L_Q16, i, j, M ) = silk_ADD32( silk_SMMUL( tmp_32, one_div_diag_Q48 ),
michael@0	177	silk_RSHIFT( silk_SMULWW( tmp_32, one_div_diag_Q36 ), 4 ) );
michael@0	178
michael@0	179	/* go to next column */
michael@0	180	ptr2 += M;
michael@0	181	}
michael@0	182	}
michael@0	183	}
michael@0	184
michael@0	185	silk_assert( status == 0 );
michael@0	186	}
michael@0	187
michael@0	188	static OPUS_INLINE void silk_LS_divide_Q16_FIX(
michael@0	189	opus_int32 T[], /* I/O Numenator vector */
michael@0	190	inv_D_t *inv_D, /* I 1 / D vector */
michael@0	191	opus_int M /* I dimension */
michael@0	192	)
michael@0	193	{
michael@0	194	opus_int i;
michael@0	195	opus_int32 tmp_32;
michael@0	196	opus_int32 one_div_diag_Q36, one_div_diag_Q48;
michael@0	197
michael@0	198	for( i = 0; i < M; i++ ) {
michael@0	199	one_div_diag_Q36 = inv_D[ i ].Q36_part;
michael@0	200	one_div_diag_Q48 = inv_D[ i ].Q48_part;
michael@0	201
michael@0	202	tmp_32 = T[ i ];
michael@0	203	T[ i ] = silk_ADD32( silk_SMMUL( tmp_32, one_div_diag_Q48 ), silk_RSHIFT( silk_SMULWW( tmp_32, one_div_diag_Q36 ), 4 ) );
michael@0	204	}
michael@0	205	}
michael@0	206
michael@0	207	/* Solve Lx = b, when L is lower triangular and has ones on the diagonal */
michael@0	208	static OPUS_INLINE void silk_LS_SolveFirst_FIX(
michael@0	209	const opus_int32 *L_Q16, /* I Pointer to Lower Triangular Matrix */
michael@0	210	opus_int M, /* I Dim of Matrix equation */
michael@0	211	const opus_int32 *b, /* I b Vector */
michael@0	212	opus_int32 *x_Q16 /* O x Vector */
michael@0	213	)
michael@0	214	{
michael@0	215	opus_int i, j;
michael@0	216	const opus_int32 *ptr32;
michael@0	217	opus_int32 tmp_32;
michael@0	218
michael@0	219	for( i = 0; i < M; i++ ) {
michael@0	220	ptr32 = matrix_adr( L_Q16, i, 0, M );
michael@0	221	tmp_32 = 0;
michael@0	222	for( j = 0; j < i; j++ ) {
michael@0	223	tmp_32 = silk_SMLAWW( tmp_32, ptr32[ j ], x_Q16[ j ] );
michael@0	224	}
michael@0	225	x_Q16[ i ] = silk_SUB32( b[ i ], tmp_32 );
michael@0	226	}
michael@0	227	}
michael@0	228
michael@0	229	/* Solve L^t*x = b, where L is lower triangular with ones on the diagonal */
michael@0	230	static OPUS_INLINE void silk_LS_SolveLast_FIX(
michael@0	231	const opus_int32 *L_Q16, /* I Pointer to Lower Triangular Matrix */
michael@0	232	const opus_int M, /* I Dim of Matrix equation */
michael@0	233	const opus_int32 *b, /* I b Vector */
michael@0	234	opus_int32 *x_Q16 /* O x Vector */
michael@0	235	)
michael@0	236	{
michael@0	237	opus_int i, j;
michael@0	238	const opus_int32 *ptr32;
michael@0	239	opus_int32 tmp_32;
michael@0	240
michael@0	241	for( i = M - 1; i >= 0; i-- ) {
michael@0	242	ptr32 = matrix_adr( L_Q16, 0, i, M );
michael@0	243	tmp_32 = 0;
michael@0	244	for( j = M - 1; j > i; j-- ) {
michael@0	245	tmp_32 = silk_SMLAWW( tmp_32, ptr32[ silk_SMULBB( j, M ) ], x_Q16[ j ] );
michael@0	246	}
michael@0	247	x_Q16[ i ] = silk_SUB32( b[ i ], tmp_32 );
michael@0	248	}
michael@0	249	}