License update using the IETF Trust flavour of the BSD on the Silk code
[opus.git] / silk / float / burg_modified_FLP.c
index afda450..0add3a6 100644 (file)
@@ -1,28 +1,28 @@
 /***********************************************************************
 Copyright (c) 2006-2011, Skype Limited. All rights reserved.
 Redistribution and use in source and binary forms, with or without
-modification, (subject to the limitations in the disclaimer below)
-are permitted provided that the following conditions are met:
+modification, are permitted provided that the following conditions
+are met:
 - Redistributions of source code must retain the above copyright notice,
 this list of conditions and the following disclaimer.
 - Redistributions in binary form must reproduce the above copyright
 notice, this list of conditions and the following disclaimer in the
 documentation and/or other materials provided with the distribution.
-- Neither the name of Skype Limited, nor the names of specific
-contributors, may be used to endorse or promote products derived from
-this software without specific prior written permission.
-NO EXPRESS OR IMPLIED LICENSES TO ANY PARTY'S PATENT RIGHTS ARE GRANTED
-BY THIS LICENSE. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
-CONTRIBUTORS ''AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
-BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
-FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
-COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
-INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
-NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
-USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
-ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
-(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
-OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+- Neither the name of Internet Society, IETF or IETF Trust, nor the 
+names of specific contributors, may be used to endorse or promote
+products derived from this software without specific prior written
+permission.
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS”
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGE.
 ***********************************************************************/
 
 #ifdef HAVE_CONFIG_H
@@ -30,29 +30,29 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #endif
 
 #include "SigProc_FLP.h"
+#include "tuning_parameters.h"
+#include "define.h"
 
 #define MAX_FRAME_SIZE              384 /* subfr_length * nb_subfr = ( 0.005 * 16000 + 16 ) * 4 = 384*/
-#define MAX_NB_SUBFR                4
 
 /* Compute reflection coefficients from input signal */
-silk_float silk_burg_modified_FLP(   /* O    returns residual energy                                         */
-    silk_float       A[],                /* O    prediction coefficients (length order)                          */
-    const silk_float x[],                /* I    input signal, length: nb_subfr*(D+L_sub)                        */
-    const opus_int   subfr_length,       /* I    input signal subframe length (including D preceeding samples)   */
-    const opus_int   nb_subfr,           /* I    number of subframes stacked in x                                */
-    const silk_float WhiteNoiseFrac,     /* I    fraction added to zero-lag autocorrelation                      */
-    const opus_int   D                   /* I    order                                                           */
+silk_float silk_burg_modified_FLP(          /* O    returns residual energy                                     */
+    silk_float          A[],                /* O    prediction coefficients (length order)                      */
+    const silk_float    x[],                /* I    input signal, length: nb_subfr*(D+L_sub)                    */
+    const silk_float    minInvGain,         /* I    minimum inverse prediction gain                             */
+    const opus_int      subfr_length,       /* I    input signal subframe length (incl. D preceeding samples)   */
+    const opus_int      nb_subfr,           /* I    number of subframes stacked in x                            */
+    const opus_int      D                   /* I    order                                                       */
 )
 {
-    opus_int         k, n, s;
-    double          C0, num, nrg_f, nrg_b, rc, Atmp, tmp1, tmp2;
+    opus_int         k, n, s, reached_max_gain;
+    double           C0, invGain, num, nrg_f, nrg_b, rc, Atmp, tmp1, tmp2;
     const silk_float *x_ptr;
-    double          C_first_row[ SILK_MAX_ORDER_LPC ], C_last_row[ SILK_MAX_ORDER_LPC ];
-    double          CAf[ SILK_MAX_ORDER_LPC + 1 ], CAb[ SILK_MAX_ORDER_LPC + 1 ];
-    double          Af[ SILK_MAX_ORDER_LPC ];
+    double           C_first_row[ SILK_MAX_ORDER_LPC ], C_last_row[ SILK_MAX_ORDER_LPC ];
+    double           CAf[ SILK_MAX_ORDER_LPC + 1 ], CAb[ SILK_MAX_ORDER_LPC + 1 ];
+    double           Af[ SILK_MAX_ORDER_LPC ];
 
     silk_assert( subfr_length * nb_subfr <= MAX_FRAME_SIZE );
-    silk_assert( nb_subfr <= MAX_NB_SUBFR );
 
     /* Compute autocorrelations, added over subframes */
     C0 = silk_energy_FLP( x, nb_subfr * subfr_length );
@@ -66,8 +66,9 @@ silk_float silk_burg_modified_FLP(   /* O    returns residual energy
     silk_memcpy( C_last_row, C_first_row, SILK_MAX_ORDER_LPC * sizeof( double ) );
 
     /* Initialize */
-    CAb[ 0 ] = CAf[ 0 ] = C0 + WhiteNoiseFrac * C0 + 1e-9f;
-
+    CAb[ 0 ] = CAf[ 0 ] = C0 + FIND_LPC_COND_FAC * C0 + 1e-9f;
+    invGain = 1.0f;
+    reached_max_gain = 0;
     for( n = 0; n < D; n++ ) {
         /* Update first row of correlation matrix (without first element) */
         /* Update last row of correlation matrix (without last element, stored in reversed order) */
@@ -93,7 +94,7 @@ silk_float silk_burg_modified_FLP(   /* O    returns residual energy
         tmp2 = C_last_row[ n ];
         for( k = 0; k < n; k++ ) {
             Atmp = Af[ k ];
-            tmp1 += C_last_row[ n - k - 1 ]  * Atmp;
+            tmp1 += C_last_row[  n - k - 1 ] * Atmp;
             tmp2 += C_first_row[ n - k - 1 ] * Atmp;
         }
         CAf[ n + 1 ] = tmp1;
@@ -116,6 +117,21 @@ silk_float silk_burg_modified_FLP(   /* O    returns residual energy
         rc = -2.0 * num / ( nrg_f + nrg_b );
         silk_assert( rc > -1.0 && rc < 1.0 );
 
+        /* Update inverse prediction gain */
+        tmp1 = invGain * ( 1.0 - rc * rc );
+        if( tmp1 <= minInvGain ) {
+            /* Max prediction gain exceeded; set reflection coefficient such that max prediction gain is exactly hit */
+            rc = sqrt( 1.0 - minInvGain / invGain );
+            if( num > 0 ) {
+                /* Ensure adjusted reflection coefficients has the original sign */
+                rc = -rc;
+            }
+            invGain = minInvGain;
+            reached_max_gain = 1;
+        } else {
+            invGain = tmp1;
+        }
+
         /* Update the AR coefficients */
         for( k = 0; k < (n + 1) >> 1; k++ ) {
             tmp1 = Af[ k ];
@@ -125,6 +141,14 @@ silk_float silk_burg_modified_FLP(   /* O    returns residual energy
         }
         Af[ n ] = rc;
 
+        if( reached_max_gain ) {
+            /* Reached max prediction gain; set remaining coefficients to zero and exit loop */
+            for( k = n + 1; k < D; k++ ) {
+                Af[ k ] = 0.0;
+            }
+            break;
+        }
+
         /* Update C * Af and C * Ab */
         for( k = 0; k <= n + 1; k++ ) {
             tmp1 = CAf[ k ];
@@ -133,16 +157,30 @@ silk_float silk_burg_modified_FLP(   /* O    returns residual energy
         }
     }
 
-    /* Return residual energy */
-    nrg_f = CAf[ 0 ];
-    tmp1 = 1.0;
-    for( k = 0; k < D; k++ ) {
-        Atmp = Af[ k ];
-        nrg_f += CAf[ k + 1 ] * Atmp;
-        tmp1  += Atmp * Atmp;
-        A[ k ] = (silk_float)(-Atmp);
+    if( reached_max_gain ) {
+        /* Convert to silk_float */
+        for( k = 0; k < D; k++ ) {
+            A[ k ] = (silk_float)( -Af[ k ] );
+        }
+        /* Subtract energy of preceeding samples from C0 */
+        for( s = 0; s < nb_subfr; s++ ) {
+            C0 -= silk_energy_FLP( x + s * subfr_length, D );
+        }
+        /* Approximate residual energy */
+        nrg_f = C0 * invGain;
+    } else {
+        /* Compute residual energy and store coefficients as silk_float */
+        nrg_f = CAf[ 0 ];
+        tmp1 = 1.0;
+        for( k = 0; k < D; k++ ) {
+            Atmp = Af[ k ];
+            nrg_f += CAf[ k + 1 ] * Atmp;
+            tmp1  += Atmp * Atmp;
+            A[ k ] = (silk_float)(-Atmp);
+        }
+        nrg_f -= FIND_LPC_COND_FAC * C0 * tmp1;
     }
-    nrg_f -= WhiteNoiseFrac * C0 * tmp1;
 
+    /* Return residual energy */
     return (silk_float)nrg_f;
 }