Adding a unit test for LPC_inv_pred_gain()

It checks that no clearly unstable filter passes the LPC_inv_pred_gain()
test. Also, this will make it possible to check assembly for correctness.

Modified from an original patch from Linfeng Zhang <linfengz@google.com>.

.gitignore

@@ -78,6 +78,7 @@ silk/fixed/x64
 silk/float/Debug
 silk/float/Release
 silk/float/x64
+silk/tests/test_unit_LPC_inv_pred_gain
 src/Debug
 src/Release
 src/x64

Makefile.am

@@ -84,9 +84,36 @@ pkginclude_HEADERS = include/opus.h include/opus_multistream.h include/opus_type
 noinst_HEADERS = $(OPUS_HEAD) $(SILK_HEAD) $(CELT_HEAD)
 
 if EXTRA_PROGRAMS
-noinst_PROGRAMS = opus_demo repacketizer_demo opus_compare tests/test_opus_api tests/test_opus_encode tests/test_opus_decode tests/test_opus_padding celt/tests/test_unit_cwrs32 celt/tests/test_unit_dft celt/tests/test_unit_entropy celt/tests/test_unit_laplace celt/tests/test_unit_mathops celt/tests/test_unit_mdct celt/tests/test_unit_rotation celt/tests/test_unit_types
-
-TESTS = celt/tests/test_unit_types celt/tests/test_unit_mathops celt/tests/test_unit_entropy celt/tests/test_unit_laplace celt/tests/test_unit_dft celt/tests/test_unit_mdct celt/tests/test_unit_rotation celt/tests/test_unit_cwrs32 tests/test_opus_api tests/test_opus_decode tests/test_opus_encode tests/test_opus_padding
+noinst_PROGRAMS = celt/tests/test_unit_cwrs32 \
+                  celt/tests/test_unit_dft \
+                  celt/tests/test_unit_entropy \
+                  celt/tests/test_unit_laplace \
+                  celt/tests/test_unit_mathops \
+                  celt/tests/test_unit_mdct \
+                  celt/tests/test_unit_rotation \
+                  celt/tests/test_unit_types \
+                  opus_compare \
+                  opus_demo \
+                  repacketizer_demo \
+                  silk/tests/test_unit_LPC_inv_pred_gain \
+                  tests/test_opus_api \
+                  tests/test_opus_decode \
+                  tests/test_opus_encode \
+                  tests/test_opus_padding
+
+TESTS = celt/tests/test_unit_cwrs32 \
+        celt/tests/test_unit_dft \
+        celt/tests/test_unit_entropy \
+        celt/tests/test_unit_laplace \
+        celt/tests/test_unit_mathops \
+        celt/tests/test_unit_mdct \
+        celt/tests/test_unit_rotation \
+        celt/tests/test_unit_types \
+        silk/tests/test_unit_LPC_inv_pred_gain \
+        tests/test_opus_api \
+        tests/test_opus_decode \
+        tests/test_opus_encode \
+        tests/test_opus_padding
 
 opus_demo_SOURCES = src/opus_demo.c
 
@@ -111,6 +138,11 @@ tests_test_opus_decode_LDADD = libopus.la $(NE10_LIBS) $(LIBM)
 tests_test_opus_padding_SOURCES = tests/test_opus_padding.c tests/test_opus_common.h
 tests_test_opus_padding_LDADD = libopus.la $(NE10_LIBS) $(LIBM)
 
+silk_tests_test_unit_LPC_inv_pred_gain_SOURCES = silk/tests/test_unit_LPC_inv_pred_gain.c
+silk_tests_test_unit_LPC_inv_pred_gain_LDADD = libopus.la $(NE10_LIBS) $(LIBM)
+# this target requires hidden symbols
+silk_tests_test_unit_LPC_inv_pred_gain_LDFLAGS = -static
+
 celt_tests_test_unit_cwrs32_SOURCES = celt/tests/test_unit_cwrs32.c
 celt_tests_test_unit_cwrs32_LDADD = $(LIBM)
 
@@ -267,7 +299,8 @@ $(CELT_SOURCES_ARM_ASM:%.s=%-gnu.S): $(top_srcdir)/celt/arm/arm2gnu.pl
 OPT_UNIT_TEST_OBJ = $(celt_tests_test_unit_mathops_SOURCES:.c=.o) \
                     $(celt_tests_test_unit_rotation_SOURCES:.c=.o) \
                     $(celt_tests_test_unit_mdct_SOURCES:.c=.o) \
-                    $(celt_tests_test_unit_dft_SOURCES:.c=.o)
+                    $(celt_tests_test_unit_dft_SOURCES:.c=.o) \
+                    $(silk_tests_test_unit_LPC_inv_pred_gain_SOURCES:.c=.o)
 
 if HAVE_SSE
 SSE_OBJ = $(CELT_SOURCES_SSE:.c=.lo)

silk/tests/test_unit_LPC_inv_pred_gain.c

+/***********************************************************************
+Copyright (c) 2017 Google Inc., Jean-Marc Valin
+Redistribution and use in source and binary forms, with or without
+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 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
+#include "config.h"
+#endif
+
+#include <stdio.h>
+#include <stdlib.h>
+#include "celt/stack_alloc.h"
+#include "cpu_support.h"
+#include "SigProc_FIX.h"
+
+/* Computes the impulse response of the filter so we
+   can catch filters that are definitely unstable. Some
+   unstable filters may be classified as stable, but not
+   the other way around. */
+int check_stability(opus_int16 *A_Q12, int order) {
+    int i;
+    int j;
+    int sum_a, sum_abs_a;
+    sum_a = sum_abs_a = 0;
+    for( j = 0; j < order; j++ ) {
+        sum_a += A_Q12[ j ];
+        sum_abs_a += silk_abs( A_Q12[ j ] );
+    }
+    /* Check DC stability. */
+    if( sum_a >= 4096 ) {
+        return 0;
+    }
+    /* If the sum of absolute values is less than 1, the filter
+       has to be stable. */
+    if( sum_abs_a < 4096 ) {
+        return 1;
+    }
+    double y[SILK_MAX_ORDER_LPC] = {0};
+    y[0] = 1;
+    for( i = 0; i < 10000; i++ ) {
+        double sum = 0;
+        for( j = 0; j < order; j++ ) {
+            sum += y[ j ]*A_Q12[ j ];
+        }
+        for( j = order - 1; j > 0; j-- ) {
+            y[ j ] = y[ j - 1 ];
+        }
+        y[ 0 ] = sum*(1./4096);
+        /* If impulse response reaches +/- 10000, the filter
+           is definitely unstable. */
+        if( !(y[ 0 ] < 10000 && y[ 0 ] > -10000) ) {
+            return 0;
+        }
+    }
+    return 1;
+}
+
+int main(void) {
+    /* Set to 10000 so all branches in C function are triggered */
+    const int loop_num = 10000;
+    int count = 0;
+    ALLOC_STACK;
+
+    /* FIXME: Make the seed random (with option to set it explicitly)
+       so we get wider coverage. */
+    srand(0);
+
+    printf("Testing silk_LPC_inverse_pred_gain() optimization ...\n");
+    for( count = 0; count < loop_num; count++ ) {
+        unsigned int i;
+        opus_int     order;
+        unsigned int shift;
+        opus_int16   A_Q12[ SILK_MAX_ORDER_LPC ];
+        opus_int32 gain;
+
+        for( order = 2; order <= SILK_MAX_ORDER_LPC; order += 2 ) { /* order must be even. */
+            for( shift = 0; shift < 16; shift++ ) { /* Different dynamic range. */
+                for( i = 0; i < SILK_MAX_ORDER_LPC; i++ ) {
+                    A_Q12[i] = ((opus_int16)rand()) >> shift;
+                }
+                gain = silk_LPC_inverse_pred_gain(A_Q12, order);
+                /* Look for filters that silk_LPC_inverse_pred_gain() thinks are
+                   stable but definitely aren't. */
+                if( gain != 0 && !check_stability(A_Q12, order) ) {
+                    fprintf(stderr, "**Loop %4d failed!**\n", count);
+                    return 1;
+                }
+            }
+        }
+        if( !(count % 500) ) {
+            printf("Loop %4d passed\n", count);
+        }
+    }
+    printf("silk_LPC_inverse_pred_gain() optimization passed\n");
+    return 0;
+}