Chris@69: /***********************************************************************
Chris@69: Copyright (c) 2017 Google Inc., Jean-Marc Valin
Chris@69: Redistribution and use in source and binary forms, with or without
Chris@69: modification, are permitted provided that the following conditions
Chris@69: are met:
Chris@69: - Redistributions of source code must retain the above copyright notice,
Chris@69: this list of conditions and the following disclaimer.
Chris@69: - Redistributions in binary form must reproduce the above copyright
Chris@69: notice, this list of conditions and the following disclaimer in the
Chris@69: documentation and/or other materials provided with the distribution.
Chris@69: - Neither the name of Internet Society, IETF or IETF Trust, nor the
Chris@69: names of specific contributors, may be used to endorse or promote
Chris@69: products derived from this software without specific prior written
Chris@69: permission.
Chris@69: THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
Chris@69: AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
Chris@69: IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
Chris@69: ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
Chris@69: LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
Chris@69: CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
Chris@69: SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
Chris@69: INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
Chris@69: CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
Chris@69: ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
Chris@69: POSSIBILITY OF SUCH DAMAGE.
Chris@69: ***********************************************************************/
Chris@69: 
Chris@69: #ifdef HAVE_CONFIG_H
Chris@69: #include "config.h"
Chris@69: #endif
Chris@69: 
Chris@69: #include <stdio.h>
Chris@69: #include <stdlib.h>
Chris@69: #include "celt/stack_alloc.h"
Chris@69: #include "cpu_support.h"
Chris@69: #include "SigProc_FIX.h"
Chris@69: 
Chris@69: /* Computes the impulse response of the filter so we
Chris@69:    can catch filters that are definitely unstable. Some
Chris@69:    unstable filters may be classified as stable, but not
Chris@69:    the other way around. */
Chris@69: int check_stability(opus_int16 *A_Q12, int order) {
Chris@69:     int i;
Chris@69:     int j;
Chris@69:     int sum_a, sum_abs_a;
Chris@69:     sum_a = sum_abs_a = 0;
Chris@69:     for( j = 0; j < order; j++ ) {
Chris@69:         sum_a += A_Q12[ j ];
Chris@69:         sum_abs_a += silk_abs( A_Q12[ j ] );
Chris@69:     }
Chris@69:     /* Check DC stability. */
Chris@69:     if( sum_a >= 4096 ) {
Chris@69:         return 0;
Chris@69:     }
Chris@69:     /* If the sum of absolute values is less than 1, the filter
Chris@69:        has to be stable. */
Chris@69:     if( sum_abs_a < 4096 ) {
Chris@69:         return 1;
Chris@69:     }
Chris@69:     double y[SILK_MAX_ORDER_LPC] = {0};
Chris@69:     y[0] = 1;
Chris@69:     for( i = 0; i < 10000; i++ ) {
Chris@69:         double sum = 0;
Chris@69:         for( j = 0; j < order; j++ ) {
Chris@69:             sum += y[ j ]*A_Q12[ j ];
Chris@69:         }
Chris@69:         for( j = order - 1; j > 0; j-- ) {
Chris@69:             y[ j ] = y[ j - 1 ];
Chris@69:         }
Chris@69:         y[ 0 ] = sum*(1./4096);
Chris@69:         /* If impulse response reaches +/- 10000, the filter
Chris@69:            is definitely unstable. */
Chris@69:         if( !(y[ 0 ] < 10000 && y[ 0 ] > -10000) ) {
Chris@69:             return 0;
Chris@69:         }
Chris@69:         /* Test every 8 sample for low amplitude. */
Chris@69:         if( ( i & 0x7 ) == 0 ) {
Chris@69:             double amp = 0;
Chris@69:             for( j = 0; j < order; j++ ) {
Chris@69:                 amp += fabs(y[j]);
Chris@69:             }
Chris@69:             if( amp < 0.00001 ) {
Chris@69:                 return 1;
Chris@69:             }
Chris@69:         }
Chris@69:     }
Chris@69:     return 1;
Chris@69: }
Chris@69: 
Chris@69: int main(void) {
Chris@69:     const int arch = opus_select_arch();
Chris@69:     /* Set to 10000 so all branches in C function are triggered */
Chris@69:     const int loop_num = 10000;
Chris@69:     int count = 0;
Chris@69:     ALLOC_STACK;
Chris@69: 
Chris@69:     /* FIXME: Make the seed random (with option to set it explicitly)
Chris@69:        so we get wider coverage. */
Chris@69:     srand(0);
Chris@69: 
Chris@69:     printf("Testing silk_LPC_inverse_pred_gain() optimization ...\n");
Chris@69:     for( count = 0; count < loop_num; count++ ) {
Chris@69:         unsigned int i;
Chris@69:         opus_int     order;
Chris@69:         unsigned int shift;
Chris@69:         opus_int16   A_Q12[ SILK_MAX_ORDER_LPC ];
Chris@69:         opus_int32 gain;
Chris@69: 
Chris@69:         for( order = 2; order <= SILK_MAX_ORDER_LPC; order += 2 ) { /* order must be even. */
Chris@69:             for( shift = 0; shift < 16; shift++ ) { /* Different dynamic range. */
Chris@69:                 for( i = 0; i < SILK_MAX_ORDER_LPC; i++ ) {
Chris@69:                     A_Q12[i] = ((opus_int16)rand()) >> shift;
Chris@69:                 }
Chris@69:                 gain = silk_LPC_inverse_pred_gain(A_Q12, order, arch);
Chris@69:                 /* Look for filters that silk_LPC_inverse_pred_gain() thinks are
Chris@69:                    stable but definitely aren't. */
Chris@69:                 if( gain != 0 && !check_stability(A_Q12, order) ) {
Chris@69:                     fprintf(stderr, "**Loop %4d failed!**\n", count);
Chris@69:                     return 1;
Chris@69:                 }
Chris@69:             }
Chris@69:         }
Chris@69:         if( !(count % 500) ) {
Chris@69:             printf("Loop %4d passed\n", count);
Chris@69:         }
Chris@69:     }
Chris@69:     printf("silk_LPC_inverse_pred_gain() optimization passed\n");
Chris@69:     return 0;
Chris@69: }