1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/media/libopus/silk/float/find_LTP_FLP.c Wed Dec 31 06:09:35 2014 +0100
1.3 @@ -0,0 +1,132 @@
1.4 +/***********************************************************************
1.5 +Copyright (c) 2006-2011, Skype Limited. All rights reserved.
1.6 +Redistribution and use in source and binary forms, with or without
1.7 +modification, are permitted provided that the following conditions
1.8 +are met:
1.9 +- Redistributions of source code must retain the above copyright notice,
1.10 +this list of conditions and the following disclaimer.
1.11 +- Redistributions in binary form must reproduce the above copyright
1.12 +notice, this list of conditions and the following disclaimer in the
1.13 +documentation and/or other materials provided with the distribution.
1.14 +- Neither the name of Internet Society, IETF or IETF Trust, nor the
1.15 +names of specific contributors, may be used to endorse or promote
1.16 +products derived from this software without specific prior written
1.17 +permission.
1.18 +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
1.19 +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
1.20 +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
1.21 +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
1.22 +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
1.23 +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
1.24 +SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
1.25 +INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
1.26 +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
1.27 +ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
1.28 +POSSIBILITY OF SUCH DAMAGE.
1.29 +***********************************************************************/
1.30 +
1.31 +#ifdef HAVE_CONFIG_H
1.32 +#include "config.h"
1.33 +#endif
1.34 +
1.35 +#include "main_FLP.h"
1.36 +#include "tuning_parameters.h"
1.37 +
1.38 +void silk_find_LTP_FLP(
1.39 + silk_float b[ MAX_NB_SUBFR * LTP_ORDER ], /* O LTP coefs */
1.40 + silk_float WLTP[ MAX_NB_SUBFR * LTP_ORDER * LTP_ORDER ], /* O Weight for LTP quantization */
1.41 + silk_float *LTPredCodGain, /* O LTP coding gain */
1.42 + const silk_float r_lpc[], /* I LPC residual */
1.43 + const opus_int lag[ MAX_NB_SUBFR ], /* I LTP lags */
1.44 + const silk_float Wght[ MAX_NB_SUBFR ], /* I Weights */
1.45 + const opus_int subfr_length, /* I Subframe length */
1.46 + const opus_int nb_subfr, /* I number of subframes */
1.47 + const opus_int mem_offset /* I Number of samples in LTP memory */
1.48 +)
1.49 +{
1.50 + opus_int i, k;
1.51 + silk_float *b_ptr, temp, *WLTP_ptr;
1.52 + silk_float LPC_res_nrg, LPC_LTP_res_nrg;
1.53 + silk_float d[ MAX_NB_SUBFR ], m, g, delta_b[ LTP_ORDER ];
1.54 + silk_float w[ MAX_NB_SUBFR ], nrg[ MAX_NB_SUBFR ], regu;
1.55 + silk_float Rr[ LTP_ORDER ], rr[ MAX_NB_SUBFR ];
1.56 + const silk_float *r_ptr, *lag_ptr;
1.57 +
1.58 + b_ptr = b;
1.59 + WLTP_ptr = WLTP;
1.60 + r_ptr = &r_lpc[ mem_offset ];
1.61 + for( k = 0; k < nb_subfr; k++ ) {
1.62 + lag_ptr = r_ptr - ( lag[ k ] + LTP_ORDER / 2 );
1.63 +
1.64 + silk_corrMatrix_FLP( lag_ptr, subfr_length, LTP_ORDER, WLTP_ptr );
1.65 + silk_corrVector_FLP( lag_ptr, r_ptr, subfr_length, LTP_ORDER, Rr );
1.66 +
1.67 + rr[ k ] = ( silk_float )silk_energy_FLP( r_ptr, subfr_length );
1.68 + regu = 1.0f + rr[ k ] +
1.69 + matrix_ptr( WLTP_ptr, 0, 0, LTP_ORDER ) +
1.70 + matrix_ptr( WLTP_ptr, LTP_ORDER-1, LTP_ORDER-1, LTP_ORDER );
1.71 + regu *= LTP_DAMPING / 3;
1.72 + silk_regularize_correlations_FLP( WLTP_ptr, &rr[ k ], regu, LTP_ORDER );
1.73 + silk_solve_LDL_FLP( WLTP_ptr, LTP_ORDER, Rr, b_ptr );
1.74 +
1.75 + /* Calculate residual energy */
1.76 + nrg[ k ] = silk_residual_energy_covar_FLP( b_ptr, WLTP_ptr, Rr, rr[ k ], LTP_ORDER );
1.77 +
1.78 + temp = Wght[ k ] / ( nrg[ k ] * Wght[ k ] + 0.01f * subfr_length );
1.79 + silk_scale_vector_FLP( WLTP_ptr, temp, LTP_ORDER * LTP_ORDER );
1.80 + w[ k ] = matrix_ptr( WLTP_ptr, LTP_ORDER / 2, LTP_ORDER / 2, LTP_ORDER );
1.81 +
1.82 + r_ptr += subfr_length;
1.83 + b_ptr += LTP_ORDER;
1.84 + WLTP_ptr += LTP_ORDER * LTP_ORDER;
1.85 + }
1.86 +
1.87 + /* Compute LTP coding gain */
1.88 + if( LTPredCodGain != NULL ) {
1.89 + LPC_LTP_res_nrg = 1e-6f;
1.90 + LPC_res_nrg = 0.0f;
1.91 + for( k = 0; k < nb_subfr; k++ ) {
1.92 + LPC_res_nrg += rr[ k ] * Wght[ k ];
1.93 + LPC_LTP_res_nrg += nrg[ k ] * Wght[ k ];
1.94 + }
1.95 +
1.96 + silk_assert( LPC_LTP_res_nrg > 0 );
1.97 + *LTPredCodGain = 3.0f * silk_log2( LPC_res_nrg / LPC_LTP_res_nrg );
1.98 + }
1.99 +
1.100 + /* Smoothing */
1.101 + /* d = sum( B, 1 ); */
1.102 + b_ptr = b;
1.103 + for( k = 0; k < nb_subfr; k++ ) {
1.104 + d[ k ] = 0;
1.105 + for( i = 0; i < LTP_ORDER; i++ ) {
1.106 + d[ k ] += b_ptr[ i ];
1.107 + }
1.108 + b_ptr += LTP_ORDER;
1.109 + }
1.110 + /* m = ( w * d' ) / ( sum( w ) + 1e-3 ); */
1.111 + temp = 1e-3f;
1.112 + for( k = 0; k < nb_subfr; k++ ) {
1.113 + temp += w[ k ];
1.114 + }
1.115 + m = 0;
1.116 + for( k = 0; k < nb_subfr; k++ ) {
1.117 + m += d[ k ] * w[ k ];
1.118 + }
1.119 + m = m / temp;
1.120 +
1.121 + b_ptr = b;
1.122 + for( k = 0; k < nb_subfr; k++ ) {
1.123 + g = LTP_SMOOTHING / ( LTP_SMOOTHING + w[ k ] ) * ( m - d[ k ] );
1.124 + temp = 0;
1.125 + for( i = 0; i < LTP_ORDER; i++ ) {
1.126 + delta_b[ i ] = silk_max_float( b_ptr[ i ], 0.1f );
1.127 + temp += delta_b[ i ];
1.128 + }
1.129 + temp = g / temp;
1.130 + for( i = 0; i < LTP_ORDER; i++ ) {
1.131 + b_ptr[ i ] = b_ptr[ i ] + delta_b[ i ] * temp;
1.132 + }
1.133 + b_ptr += LTP_ORDER;
1.134 + }
1.135 +}
