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

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

media/libopus/silk/fixed/burg_modified_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 "SigProc_FIX.h"
michael@0	33	#include "define.h"
michael@0	34	#include "tuning_parameters.h"
michael@0	35	#include "pitch.h"
michael@0	36
michael@0	37	#define MAX_FRAME_SIZE 384 /* subfr_length * nb_subfr = ( 0.005 * 16000 + 16 ) * 4 = 384 */
michael@0	38
michael@0	39	#define QA 25
michael@0	40	#define N_BITS_HEAD_ROOM 2
michael@0	41	#define MIN_RSHIFTS -16
michael@0	42	#define MAX_RSHIFTS (32 - QA)
michael@0	43
michael@0	44	/* Compute reflection coefficients from input signal */
michael@0	45	void silk_burg_modified(
michael@0	46	opus_int32 *res_nrg, /* O Residual energy */
michael@0	47	opus_int *res_nrg_Q, /* O Residual energy Q value */
michael@0	48	opus_int32 A_Q16[], /* O Prediction coefficients (length order) */
michael@0	49	const opus_int16 x[], /* I Input signal, length: nb_subfr * ( D + subfr_length ) */
michael@0	50	const opus_int32 minInvGain_Q30, /* I Inverse of max prediction gain */
michael@0	51	const opus_int subfr_length, /* I Input signal subframe length (incl. D preceding samples) */
michael@0	52	const opus_int nb_subfr, /* I Number of subframes stacked in x */
michael@0	53	const opus_int D, /* I Order */
michael@0	54	int arch /* I Run-time architecture */
michael@0	55	)
michael@0	56	{
michael@0	57	opus_int k, n, s, lz, rshifts, rshifts_extra, reached_max_gain;
michael@0	58	opus_int32 C0, num, nrg, rc_Q31, invGain_Q30, Atmp_QA, Atmp1, tmp1, tmp2, x1, x2;
michael@0	59	const opus_int16 *x_ptr;
michael@0	60	opus_int32 C_first_row[ SILK_MAX_ORDER_LPC ];
michael@0	61	opus_int32 C_last_row[ SILK_MAX_ORDER_LPC ];
michael@0	62	opus_int32 Af_QA[ SILK_MAX_ORDER_LPC ];
michael@0	63	opus_int32 CAf[ SILK_MAX_ORDER_LPC + 1 ];
michael@0	64	opus_int32 CAb[ SILK_MAX_ORDER_LPC + 1 ];
michael@0	65	opus_int32 xcorr[ SILK_MAX_ORDER_LPC ];
michael@0	66
michael@0	67	silk_assert( subfr_length * nb_subfr <= MAX_FRAME_SIZE );
michael@0	68
michael@0	69	/* Compute autocorrelations, added over subframes */
michael@0	70	silk_sum_sqr_shift( &C0, &rshifts, x, nb_subfr * subfr_length );
michael@0	71	if( rshifts > MAX_RSHIFTS ) {
michael@0	72	C0 = silk_LSHIFT32( C0, rshifts - MAX_RSHIFTS );
michael@0	73	silk_assert( C0 > 0 );
michael@0	74	rshifts = MAX_RSHIFTS;
michael@0	75	} else {
michael@0	76	lz = silk_CLZ32( C0 ) - 1;
michael@0	77	rshifts_extra = N_BITS_HEAD_ROOM - lz;
michael@0	78	if( rshifts_extra > 0 ) {
michael@0	79	rshifts_extra = silk_min( rshifts_extra, MAX_RSHIFTS - rshifts );
michael@0	80	C0 = silk_RSHIFT32( C0, rshifts_extra );
michael@0	81	} else {
michael@0	82	rshifts_extra = silk_max( rshifts_extra, MIN_RSHIFTS - rshifts );
michael@0	83	C0 = silk_LSHIFT32( C0, -rshifts_extra );
michael@0	84	}
michael@0	85	rshifts += rshifts_extra;
michael@0	86	}
michael@0	87	CAb[ 0 ] = CAf[ 0 ] = C0 + silk_SMMUL( SILK_FIX_CONST( FIND_LPC_COND_FAC, 32 ), C0 ) + 1; /* Q(-rshifts) */
michael@0	88	silk_memset( C_first_row, 0, SILK_MAX_ORDER_LPC * sizeof( opus_int32 ) );
michael@0	89	if( rshifts > 0 ) {
michael@0	90	for( s = 0; s < nb_subfr; s++ ) {
michael@0	91	x_ptr = x + s * subfr_length;
michael@0	92	for( n = 1; n < D + 1; n++ ) {
michael@0	93	C_first_row[ n - 1 ] += (opus_int32)silk_RSHIFT64(
michael@0	94	silk_inner_prod16_aligned_64( x_ptr, x_ptr + n, subfr_length - n ), rshifts );
michael@0	95	}
michael@0	96	}
michael@0	97	} else {
michael@0	98	for( s = 0; s < nb_subfr; s++ ) {
michael@0	99	int i;
michael@0	100	opus_int32 d;
michael@0	101	x_ptr = x + s * subfr_length;
michael@0	102	celt_pitch_xcorr(x_ptr, x_ptr + 1, xcorr, subfr_length - D, D, arch );
michael@0	103	for( n = 1; n < D + 1; n++ ) {
michael@0	104	for ( i = n + subfr_length - D, d = 0; i < subfr_length; i++ )
michael@0	105	d = MAC16_16( d, x_ptr[ i ], x_ptr[ i - n ] );
michael@0	106	xcorr[ n - 1 ] += d;
michael@0	107	}
michael@0	108	for( n = 1; n < D + 1; n++ ) {
michael@0	109	C_first_row[ n - 1 ] += silk_LSHIFT32( xcorr[ n - 1 ], -rshifts );
michael@0	110	}
michael@0	111	}
michael@0	112	}
michael@0	113	silk_memcpy( C_last_row, C_first_row, SILK_MAX_ORDER_LPC * sizeof( opus_int32 ) );
michael@0	114
michael@0	115	/* Initialize */
michael@0	116	CAb[ 0 ] = CAf[ 0 ] = C0 + silk_SMMUL( SILK_FIX_CONST( FIND_LPC_COND_FAC, 32 ), C0 ) + 1; /* Q(-rshifts) */
michael@0	117
michael@0	118	invGain_Q30 = (opus_int32)1 << 30;
michael@0	119	reached_max_gain = 0;
michael@0	120	for( n = 0; n < D; n++ ) {
michael@0	121	/* Update first row of correlation matrix (without first element) */
michael@0	122	/* Update last row of correlation matrix (without last element, stored in reversed order) */
michael@0	123	/* Update C * Af */
michael@0	124	/* Update C * flipud(Af) (stored in reversed order) */
michael@0	125	if( rshifts > -2 ) {
michael@0	126	for( s = 0; s < nb_subfr; s++ ) {
michael@0	127	x_ptr = x + s * subfr_length;
michael@0	128	x1 = -silk_LSHIFT32( (opus_int32)x_ptr[ n ], 16 - rshifts ); /* Q(16-rshifts) */
michael@0	129	x2 = -silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n - 1 ], 16 - rshifts ); /* Q(16-rshifts) */
michael@0	130	tmp1 = silk_LSHIFT32( (opus_int32)x_ptr[ n ], QA - 16 ); /* Q(QA-16) */
michael@0	131	tmp2 = silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n - 1 ], QA - 16 ); /* Q(QA-16) */
michael@0	132	for( k = 0; k < n; k++ ) {
michael@0	133	C_first_row[ k ] = silk_SMLAWB( C_first_row[ k ], x1, x_ptr[ n - k - 1 ] ); /* Q( -rshifts ) */
michael@0	134	C_last_row[ k ] = silk_SMLAWB( C_last_row[ k ], x2, x_ptr[ subfr_length - n + k ] ); /* Q( -rshifts ) */
michael@0	135	Atmp_QA = Af_QA[ k ];
michael@0	136	tmp1 = silk_SMLAWB( tmp1, Atmp_QA, x_ptr[ n - k - 1 ] ); /* Q(QA-16) */
michael@0	137	tmp2 = silk_SMLAWB( tmp2, Atmp_QA, x_ptr[ subfr_length - n + k ] ); /* Q(QA-16) */
michael@0	138	}
michael@0	139	tmp1 = silk_LSHIFT32( -tmp1, 32 - QA - rshifts ); /* Q(16-rshifts) */
michael@0	140	tmp2 = silk_LSHIFT32( -tmp2, 32 - QA - rshifts ); /* Q(16-rshifts) */
michael@0	141	for( k = 0; k <= n; k++ ) {
michael@0	142	CAf[ k ] = silk_SMLAWB( CAf[ k ], tmp1, x_ptr[ n - k ] ); /* Q( -rshift ) */
michael@0	143	CAb[ k ] = silk_SMLAWB( CAb[ k ], tmp2, x_ptr[ subfr_length - n + k - 1 ] ); /* Q( -rshift ) */
michael@0	144	}
michael@0	145	}
michael@0	146	} else {
michael@0	147	for( s = 0; s < nb_subfr; s++ ) {
michael@0	148	x_ptr = x + s * subfr_length;
michael@0	149	x1 = -silk_LSHIFT32( (opus_int32)x_ptr[ n ], -rshifts ); /* Q( -rshifts ) */
michael@0	150	x2 = -silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n - 1 ], -rshifts ); /* Q( -rshifts ) */
michael@0	151	tmp1 = silk_LSHIFT32( (opus_int32)x_ptr[ n ], 17 ); /* Q17 */
michael@0	152	tmp2 = silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n - 1 ], 17 ); /* Q17 */
michael@0	153	for( k = 0; k < n; k++ ) {
michael@0	154	C_first_row[ k ] = silk_MLA( C_first_row[ k ], x1, x_ptr[ n - k - 1 ] ); /* Q( -rshifts ) */
michael@0	155	C_last_row[ k ] = silk_MLA( C_last_row[ k ], x2, x_ptr[ subfr_length - n + k ] ); /* Q( -rshifts ) */
michael@0	156	Atmp1 = silk_RSHIFT_ROUND( Af_QA[ k ], QA - 17 ); /* Q17 */
michael@0	157	tmp1 = silk_MLA( tmp1, x_ptr[ n - k - 1 ], Atmp1 ); /* Q17 */
michael@0	158	tmp2 = silk_MLA( tmp2, x_ptr[ subfr_length - n + k ], Atmp1 ); /* Q17 */
michael@0	159	}
michael@0	160	tmp1 = -tmp1; /* Q17 */
michael@0	161	tmp2 = -tmp2; /* Q17 */
michael@0	162	for( k = 0; k <= n; k++ ) {
michael@0	163	CAf[ k ] = silk_SMLAWW( CAf[ k ], tmp1,
michael@0	164	silk_LSHIFT32( (opus_int32)x_ptr[ n - k ], -rshifts - 1 ) ); /* Q( -rshift ) */
michael@0	165	CAb[ k ] = silk_SMLAWW( CAb[ k ], tmp2,
michael@0	166	silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n + k - 1 ], -rshifts - 1 ) ); /* Q( -rshift ) */
michael@0	167	}
michael@0	168	}
michael@0	169	}
michael@0	170
michael@0	171	/* Calculate nominator and denominator for the next order reflection (parcor) coefficient */
michael@0	172	tmp1 = C_first_row[ n ]; /* Q( -rshifts ) */
michael@0	173	tmp2 = C_last_row[ n ]; /* Q( -rshifts ) */
michael@0	174	num = 0; /* Q( -rshifts ) */
michael@0	175	nrg = silk_ADD32( CAb[ 0 ], CAf[ 0 ] ); /* Q( 1-rshifts ) */
michael@0	176	for( k = 0; k < n; k++ ) {
michael@0	177	Atmp_QA = Af_QA[ k ];
michael@0	178	lz = silk_CLZ32( silk_abs( Atmp_QA ) ) - 1;
michael@0	179	lz = silk_min( 32 - QA, lz );
michael@0	180	Atmp1 = silk_LSHIFT32( Atmp_QA, lz ); /* Q( QA + lz ) */
michael@0	181
michael@0	182	tmp1 = silk_ADD_LSHIFT32( tmp1, silk_SMMUL( C_last_row[ n - k - 1 ], Atmp1 ), 32 - QA - lz ); /* Q( -rshifts ) */
michael@0	183	tmp2 = silk_ADD_LSHIFT32( tmp2, silk_SMMUL( C_first_row[ n - k - 1 ], Atmp1 ), 32 - QA - lz ); /* Q( -rshifts ) */
michael@0	184	num = silk_ADD_LSHIFT32( num, silk_SMMUL( CAb[ n - k ], Atmp1 ), 32 - QA - lz ); /* Q( -rshifts ) */
michael@0	185	nrg = silk_ADD_LSHIFT32( nrg, silk_SMMUL( silk_ADD32( CAb[ k + 1 ], CAf[ k + 1 ] ),
michael@0	186	Atmp1 ), 32 - QA - lz ); /* Q( 1-rshifts ) */
michael@0	187	}
michael@0	188	CAf[ n + 1 ] = tmp1; /* Q( -rshifts ) */
michael@0	189	CAb[ n + 1 ] = tmp2; /* Q( -rshifts ) */
michael@0	190	num = silk_ADD32( num, tmp2 ); /* Q( -rshifts ) */
michael@0	191	num = silk_LSHIFT32( -num, 1 ); /* Q( 1-rshifts ) */
michael@0	192
michael@0	193	/* Calculate the next order reflection (parcor) coefficient */
michael@0	194	if( silk_abs( num ) < nrg ) {
michael@0	195	rc_Q31 = silk_DIV32_varQ( num, nrg, 31 );
michael@0	196	} else {
michael@0	197	rc_Q31 = ( num > 0 ) ? silk_int32_MAX : silk_int32_MIN;
michael@0	198	}
michael@0	199
michael@0	200	/* Update inverse prediction gain */
michael@0	201	tmp1 = ( (opus_int32)1 << 30 ) - silk_SMMUL( rc_Q31, rc_Q31 );
michael@0	202	tmp1 = silk_LSHIFT( silk_SMMUL( invGain_Q30, tmp1 ), 2 );
michael@0	203	if( tmp1 <= minInvGain_Q30 ) {
michael@0	204	/* Max prediction gain exceeded; set reflection coefficient such that max prediction gain is exactly hit */
michael@0	205	tmp2 = ( (opus_int32)1 << 30 ) - silk_DIV32_varQ( minInvGain_Q30, invGain_Q30, 30 ); /* Q30 */
michael@0	206	rc_Q31 = silk_SQRT_APPROX( tmp2 ); /* Q15 */
michael@0	207	/* Newton-Raphson iteration */
michael@0	208	rc_Q31 = silk_RSHIFT32( rc_Q31 + silk_DIV32( tmp2, rc_Q31 ), 1 ); /* Q15 */
michael@0	209	rc_Q31 = silk_LSHIFT32( rc_Q31, 16 ); /* Q31 */
michael@0	210	if( num < 0 ) {
michael@0	211	/* Ensure adjusted reflection coefficients has the original sign */
michael@0	212	rc_Q31 = -rc_Q31;
michael@0	213	}
michael@0	214	invGain_Q30 = minInvGain_Q30;
michael@0	215	reached_max_gain = 1;
michael@0	216	} else {
michael@0	217	invGain_Q30 = tmp1;
michael@0	218	}
michael@0	219
michael@0	220	/* Update the AR coefficients */
michael@0	221	for( k = 0; k < (n + 1) >> 1; k++ ) {
michael@0	222	tmp1 = Af_QA[ k ]; /* QA */
michael@0	223	tmp2 = Af_QA[ n - k - 1 ]; /* QA */
michael@0	224	Af_QA[ k ] = silk_ADD_LSHIFT32( tmp1, silk_SMMUL( tmp2, rc_Q31 ), 1 ); /* QA */
michael@0	225	Af_QA[ n - k - 1 ] = silk_ADD_LSHIFT32( tmp2, silk_SMMUL( tmp1, rc_Q31 ), 1 ); /* QA */
michael@0	226	}
michael@0	227	Af_QA[ n ] = silk_RSHIFT32( rc_Q31, 31 - QA ); /* QA */
michael@0	228
michael@0	229	if( reached_max_gain ) {
michael@0	230	/* Reached max prediction gain; set remaining coefficients to zero and exit loop */
michael@0	231	for( k = n + 1; k < D; k++ ) {
michael@0	232	Af_QA[ k ] = 0;
michael@0	233	}
michael@0	234	break;
michael@0	235	}
michael@0	236
michael@0	237	/* Update C * Af and C * Ab */
michael@0	238	for( k = 0; k <= n + 1; k++ ) {
michael@0	239	tmp1 = CAf[ k ]; /* Q( -rshifts ) */
michael@0	240	tmp2 = CAb[ n - k + 1 ]; /* Q( -rshifts ) */
michael@0	241	CAf[ k ] = silk_ADD_LSHIFT32( tmp1, silk_SMMUL( tmp2, rc_Q31 ), 1 ); /* Q( -rshifts ) */
michael@0	242	CAb[ n - k + 1 ] = silk_ADD_LSHIFT32( tmp2, silk_SMMUL( tmp1, rc_Q31 ), 1 ); /* Q( -rshifts ) */
michael@0	243	}
michael@0	244	}
michael@0	245
michael@0	246	if( reached_max_gain ) {
michael@0	247	for( k = 0; k < D; k++ ) {
michael@0	248	/* Scale coefficients */
michael@0	249	A_Q16[ k ] = -silk_RSHIFT_ROUND( Af_QA[ k ], QA - 16 );
michael@0	250	}
michael@0	251	/* Subtract energy of preceding samples from C0 */
michael@0	252	if( rshifts > 0 ) {
michael@0	253	for( s = 0; s < nb_subfr; s++ ) {
michael@0	254	x_ptr = x + s * subfr_length;
michael@0	255	C0 -= (opus_int32)silk_RSHIFT64( silk_inner_prod16_aligned_64( x_ptr, x_ptr, D ), rshifts );
michael@0	256	}
michael@0	257	} else {
michael@0	258	for( s = 0; s < nb_subfr; s++ ) {
michael@0	259	x_ptr = x + s * subfr_length;
michael@0	260	C0 -= silk_LSHIFT32( silk_inner_prod_aligned( x_ptr, x_ptr, D ), -rshifts );
michael@0	261	}
michael@0	262	}
michael@0	263	/* Approximate residual energy */
michael@0	264	*res_nrg = silk_LSHIFT( silk_SMMUL( invGain_Q30, C0 ), 2 );
michael@0	265	*res_nrg_Q = -rshifts;
michael@0	266	} else {
michael@0	267	/* Return residual energy */
michael@0	268	nrg = CAf[ 0 ]; /* Q( -rshifts ) */
michael@0	269	tmp1 = (opus_int32)1 << 16; /* Q16 */
michael@0	270	for( k = 0; k < D; k++ ) {
michael@0	271	Atmp1 = silk_RSHIFT_ROUND( Af_QA[ k ], QA - 16 ); /* Q16 */
michael@0	272	nrg = silk_SMLAWW( nrg, CAf[ k + 1 ], Atmp1 ); /* Q( -rshifts ) */
michael@0	273	tmp1 = silk_SMLAWW( tmp1, Atmp1, Atmp1 ); /* Q16 */
michael@0	274	A_Q16[ k ] = -Atmp1;
michael@0	275	}
michael@0	276	*res_nrg = silk_SMLAWW( nrg, silk_SMMUL( SILK_FIX_CONST( FIND_LPC_COND_FAC, 32 ), C0 ), -tmp1 );/* Q( -rshifts ) */
michael@0	277	*res_nrg_Q = -rshifts;
michael@0	278	}
michael@0	279	}