| 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 #include "main_FLP.h" | 
| Chris@69 | 33 | 
| Chris@69 | 34 #define MAX_ITERATIONS_RESIDUAL_NRG         10 | 
| Chris@69 | 35 #define REGULARIZATION_FACTOR               1e-8f | 
| Chris@69 | 36 | 
| Chris@69 | 37 /* Residual energy: nrg = wxx - 2 * wXx * c + c' * wXX * c */ | 
| Chris@69 | 38 silk_float silk_residual_energy_covar_FLP(                              /* O    Weighted residual energy                    */ | 
| Chris@69 | 39     const silk_float                *c,                                 /* I    Filter coefficients                         */ | 
| Chris@69 | 40     silk_float                      *wXX,                               /* I/O  Weighted correlation matrix, reg. out       */ | 
| Chris@69 | 41     const silk_float                *wXx,                               /* I    Weighted correlation vector                 */ | 
| Chris@69 | 42     const silk_float                wxx,                                /* I    Weighted correlation value                  */ | 
| Chris@69 | 43     const opus_int                  D                                   /* I    Dimension                                   */ | 
| Chris@69 | 44 ) | 
| Chris@69 | 45 { | 
| Chris@69 | 46     opus_int   i, j, k; | 
| Chris@69 | 47     silk_float tmp, nrg = 0.0f, regularization; | 
| Chris@69 | 48 | 
| Chris@69 | 49     /* Safety checks */ | 
| Chris@69 | 50     celt_assert( D >= 0 ); | 
| Chris@69 | 51 | 
| Chris@69 | 52     regularization = REGULARIZATION_FACTOR * ( wXX[ 0 ] + wXX[ D * D - 1 ] ); | 
| Chris@69 | 53     for( k = 0; k < MAX_ITERATIONS_RESIDUAL_NRG; k++ ) { | 
| Chris@69 | 54         nrg = wxx; | 
| Chris@69 | 55 | 
| Chris@69 | 56         tmp = 0.0f; | 
| Chris@69 | 57         for( i = 0; i < D; i++ ) { | 
| Chris@69 | 58             tmp += wXx[ i ] * c[ i ]; | 
| Chris@69 | 59         } | 
| Chris@69 | 60         nrg -= 2.0f * tmp; | 
| Chris@69 | 61 | 
| Chris@69 | 62         /* compute c' * wXX * c, assuming wXX is symmetric */ | 
| Chris@69 | 63         for( i = 0; i < D; i++ ) { | 
| Chris@69 | 64             tmp = 0.0f; | 
| Chris@69 | 65             for( j = i + 1; j < D; j++ ) { | 
| Chris@69 | 66                 tmp += matrix_c_ptr( wXX, i, j, D ) * c[ j ]; | 
| Chris@69 | 67             } | 
| Chris@69 | 68             nrg += c[ i ] * ( 2.0f * tmp + matrix_c_ptr( wXX, i, i, D ) * c[ i ] ); | 
| Chris@69 | 69         } | 
| Chris@69 | 70         if( nrg > 0 ) { | 
| Chris@69 | 71             break; | 
| Chris@69 | 72         } else { | 
| Chris@69 | 73             /* Add white noise */ | 
| Chris@69 | 74             for( i = 0; i < D; i++ ) { | 
| Chris@69 | 75                 matrix_c_ptr( wXX, i, i, D ) +=  regularization; | 
| Chris@69 | 76             } | 
| Chris@69 | 77             /* Increase noise for next run */ | 
| Chris@69 | 78             regularization *= 2.0f; | 
| Chris@69 | 79         } | 
| Chris@69 | 80     } | 
| Chris@69 | 81     if( k == MAX_ITERATIONS_RESIDUAL_NRG ) { | 
| Chris@69 | 82         silk_assert( nrg == 0 ); | 
| Chris@69 | 83         nrg = 1.0f; | 
| Chris@69 | 84     } | 
| Chris@69 | 85 | 
| Chris@69 | 86     return nrg; | 
| Chris@69 | 87 } | 
| Chris@69 | 88 | 
| Chris@69 | 89 /* Calculates residual energies of input subframes where all subframes have LPC_order   */ | 
| Chris@69 | 90 /* of preceding samples                                                                 */ | 
| Chris@69 | 91 void silk_residual_energy_FLP( | 
| Chris@69 | 92     silk_float                      nrgs[ MAX_NB_SUBFR ],               /* O    Residual energy per subframe                */ | 
| Chris@69 | 93     const silk_float                x[],                                /* I    Input signal                                */ | 
| Chris@69 | 94     silk_float                      a[ 2 ][ MAX_LPC_ORDER ],            /* I    AR coefs for each frame half                */ | 
| Chris@69 | 95     const silk_float                gains[],                            /* I    Quantization gains                          */ | 
| Chris@69 | 96     const opus_int                  subfr_length,                       /* I    Subframe length                             */ | 
| Chris@69 | 97     const opus_int                  nb_subfr,                           /* I    number of subframes                         */ | 
| Chris@69 | 98     const opus_int                  LPC_order                           /* I    LPC order                                   */ | 
| Chris@69 | 99 ) | 
| Chris@69 | 100 { | 
| Chris@69 | 101     opus_int     shift; | 
| Chris@69 | 102     silk_float   *LPC_res_ptr, LPC_res[ ( MAX_FRAME_LENGTH + MAX_NB_SUBFR * MAX_LPC_ORDER ) / 2 ]; | 
| Chris@69 | 103 | 
| Chris@69 | 104     LPC_res_ptr = LPC_res + LPC_order; | 
| Chris@69 | 105     shift = LPC_order + subfr_length; | 
| Chris@69 | 106 | 
| Chris@69 | 107     /* Filter input to create the LPC residual for each frame half, and measure subframe energies */ | 
| Chris@69 | 108     silk_LPC_analysis_filter_FLP( LPC_res, a[ 0 ], x + 0 * shift, 2 * shift, LPC_order ); | 
| Chris@69 | 109     nrgs[ 0 ] = ( silk_float )( gains[ 0 ] * gains[ 0 ] * silk_energy_FLP( LPC_res_ptr + 0 * shift, subfr_length ) ); | 
| Chris@69 | 110     nrgs[ 1 ] = ( silk_float )( gains[ 1 ] * gains[ 1 ] * silk_energy_FLP( LPC_res_ptr + 1 * shift, subfr_length ) ); | 
| Chris@69 | 111 | 
| Chris@69 | 112     if( nb_subfr == MAX_NB_SUBFR ) { | 
| Chris@69 | 113         silk_LPC_analysis_filter_FLP( LPC_res, a[ 1 ], x + 2 * shift, 2 * shift, LPC_order ); | 
| Chris@69 | 114         nrgs[ 2 ] = ( silk_float )( gains[ 2 ] * gains[ 2 ] * silk_energy_FLP( LPC_res_ptr + 0 * shift, subfr_length ) ); | 
| Chris@69 | 115         nrgs[ 3 ] = ( silk_float )( gains[ 3 ] * gains[ 3 ] * silk_energy_FLP( LPC_res_ptr + 1 * shift, subfr_length ) ); | 
| Chris@69 | 116     } | 
| Chris@69 | 117 } |