/src/opus/silk/fixed/burg_modified_FIX.c

Source
/***********************************************************************
Copyright (c) 2006-2011, Skype Limited. All rights reserved.
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 "SigProc_FIX.h"
#include "define.h"
#include "tuning_parameters.h"
#include "pitch.h"

#define MAX_FRAME_SIZE              384             /* subfr_length * nb_subfr = ( 0.005 * 16000 + 16 ) * 4 = 384 */

#define QA                          25
#define N_BITS_HEAD_ROOM            3
#define MIN_RSHIFTS                 -16
#define MAX_RSHIFTS                 (32 - QA)

/* Compute reflection coefficients from input signal */
void silk_burg_modified_c(
    opus_int32                  *res_nrg,           /* O    Residual energy                                             */
    opus_int                    *res_nrg_Q,         /* O    Residual energy Q value                                     */
    opus_int32                  A_Q16[],            /* O    Prediction coefficients (length order)                      */
    const opus_int16            x[],                /* I    Input signal, length: nb_subfr * ( D + subfr_length )       */
    const opus_int32            minInvGain_Q30,     /* I    Inverse of max prediction gain                              */
    const opus_int              subfr_length,       /* I    Input signal subframe length (incl. D preceding samples)    */
    const opus_int              nb_subfr,           /* I    Number of subframes stacked in x                            */
    const opus_int              D,                  /* I    Order                                                       */
    int                         arch                /* I    Run-time architecture                                       */
)
{
    opus_int         k, n, s, lz, rshifts, reached_max_gain;
    opus_int32       C0, num, nrg, rc_Q31, invGain_Q30, Atmp_QA, Atmp1, tmp1, tmp2, x1, x2;
    const opus_int16 *x_ptr;
    opus_int32       C_first_row[ SILK_MAX_ORDER_LPC ];
    opus_int32       C_last_row[  SILK_MAX_ORDER_LPC ];
    opus_int32       Af_QA[       SILK_MAX_ORDER_LPC ];
    opus_int32       CAf[ SILK_MAX_ORDER_LPC + 1 ];
    opus_int32       CAb[ SILK_MAX_ORDER_LPC + 1 ];
    opus_int32       xcorr[ SILK_MAX_ORDER_LPC ];
    opus_int64       C0_64;

    celt_assert( subfr_length * nb_subfr <= MAX_FRAME_SIZE );

    /* Compute autocorrelations, added over subframes */
    C0_64 = silk_inner_prod16( x, x, subfr_length*nb_subfr, arch );
    lz = silk_CLZ64(C0_64);
    rshifts = 32 + 1 + N_BITS_HEAD_ROOM - lz;
    if (rshifts > MAX_RSHIFTS) rshifts = MAX_RSHIFTS;
    if (rshifts < MIN_RSHIFTS) rshifts = MIN_RSHIFTS;

    if (rshifts > 0) {
        C0 = (opus_int32)silk_RSHIFT64(C0_64, rshifts );
    } else {
        C0 = silk_LSHIFT32((opus_int32)C0_64, -rshifts );
    }

    CAb[ 0 ] = CAf[ 0 ] = C0 + silk_SMMUL( SILK_FIX_CONST( FIND_LPC_COND_FAC, 32 ), C0 ) + 1;                                /* Q(-rshifts) */
    silk_memset( C_first_row, 0, SILK_MAX_ORDER_LPC * sizeof( opus_int32 ) );
    if( rshifts > 0 ) {
        for( s = 0; s < nb_subfr; s++ ) {
            x_ptr = x + s * subfr_length;
            for( n = 1; n < D + 1; n++ ) {
                C_first_row[ n - 1 ] += (opus_int32)silk_RSHIFT64(
                    silk_inner_prod16( x_ptr, x_ptr + n, subfr_length - n, arch ), rshifts );
            }
        }
    } else {
        for( s = 0; s < nb_subfr; s++ ) {
            int i;
            opus_int32 d;
            x_ptr = x + s * subfr_length;
            celt_pitch_xcorr(x_ptr, x_ptr + 1, xcorr, subfr_length - D, D, arch );
            for( n = 1; n < D + 1; n++ ) {
               for ( i = n + subfr_length - D, d = 0; i < subfr_length; i++ )
                  d = MAC16_16( d, x_ptr[ i ], x_ptr[ i - n ] );
               xcorr[ n - 1 ] += d;
            }
            for( n = 1; n < D + 1; n++ ) {
                C_first_row[ n - 1 ] += silk_LSHIFT32( xcorr[ n - 1 ], -rshifts );
            }
        }
    }
    silk_memcpy( C_last_row, C_first_row, SILK_MAX_ORDER_LPC * sizeof( opus_int32 ) );

    /* Initialize */
    CAb[ 0 ] = CAf[ 0 ] = C0 + silk_SMMUL( SILK_FIX_CONST( FIND_LPC_COND_FAC, 32 ), C0 ) + 1;                                /* Q(-rshifts) */

    invGain_Q30 = (opus_int32)1 << 30;
    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) */
        /* Update C * Af */
        /* Update C * flipud(Af) (stored in reversed order) */
        if( rshifts > -2 ) {
            for( s = 0; s < nb_subfr; s++ ) {
                x_ptr = x + s * subfr_length;
                x1  = -silk_LSHIFT32( (opus_int32)x_ptr[ n ],                    16 - rshifts );        /* Q(16-rshifts) */
                x2  = -silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n - 1 ], 16 - rshifts );        /* Q(16-rshifts) */
                tmp1 = silk_LSHIFT32( (opus_int32)x_ptr[ n ],                    QA - 16 );             /* Q(QA-16) */
                tmp2 = silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n - 1 ], QA - 16 );             /* Q(QA-16) */
                for( k = 0; k < n; k++ ) {
                    C_first_row[ k ] = silk_SMLAWB( C_first_row[ k ], x1, x_ptr[ n - k - 1 ]            ); /* Q( -rshifts ) */
                    C_last_row[ k ]  = silk_SMLAWB( C_last_row[ k ],  x2, x_ptr[ subfr_length - n + k ] ); /* Q( -rshifts ) */
                    Atmp_QA = Af_QA[ k ];
                    tmp1 = silk_SMLAWB( tmp1, Atmp_QA, x_ptr[ n - k - 1 ]            );                 /* Q(QA-16) */
                    tmp2 = silk_SMLAWB( tmp2, Atmp_QA, x_ptr[ subfr_length - n + k ] );                 /* Q(QA-16) */
                }
                tmp1 = silk_LSHIFT32( -tmp1, 32 - QA - rshifts );                                       /* Q(16-rshifts) */
                tmp2 = silk_LSHIFT32( -tmp2, 32 - QA - rshifts );                                       /* Q(16-rshifts) */
                for( k = 0; k <= n; k++ ) {
                    CAf[ k ] = silk_SMLAWB( CAf[ k ], tmp1, x_ptr[ n - k ]                    );        /* Q( -rshift ) */
                    CAb[ k ] = silk_SMLAWB( CAb[ k ], tmp2, x_ptr[ subfr_length - n + k - 1 ] );        /* Q( -rshift ) */
                }
            }
        } else {
            for( s = 0; s < nb_subfr; s++ ) {
                x_ptr = x + s * subfr_length;
                x1  = -silk_LSHIFT32( (opus_int32)x_ptr[ n ],                    -rshifts );            /* Q( -rshifts ) */
                x2  = -silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n - 1 ], -rshifts );            /* Q( -rshifts ) */
                tmp1 = silk_LSHIFT32( (opus_int32)x_ptr[ n ],                    17 );                  /* Q17 */
                tmp2 = silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n - 1 ], 17 );                  /* Q17 */
                for( k = 0; k < n; k++ ) {
                    C_first_row[ k ] = silk_MLA( C_first_row[ k ], x1, x_ptr[ n - k - 1 ]            ); /* Q( -rshifts ) */
                    C_last_row[ k ]  = silk_MLA( C_last_row[ k ],  x2, x_ptr[ subfr_length - n + k ] ); /* Q( -rshifts ) */
                    Atmp1 = silk_RSHIFT_ROUND( Af_QA[ k ], QA - 17 );                                   /* Q17 */
                    /* We sometimes get overflows in the multiplications (even beyond +/- 2^32),
                       but they cancel each other and the real result seems to always fit in a 32-bit
                       signed integer. This was determined experimentally, not theoretically (unfortunately). */
                    tmp1 = silk_MLA_ovflw( tmp1, x_ptr[ n - k - 1 ],            Atmp1 );                      /* Q17 */
                    tmp2 = silk_MLA_ovflw( tmp2, x_ptr[ subfr_length - n + k ], Atmp1 );                      /* Q17 */
                }
                tmp1 = -tmp1;                                                                           /* Q17 */
                tmp2 = -tmp2;                                                                           /* Q17 */
                for( k = 0; k <= n; k++ ) {
                    CAf[ k ] = silk_SMLAWW( CAf[ k ], tmp1,
                        silk_LSHIFT32( (opus_int32)x_ptr[ n - k ], -rshifts - 1 ) );                    /* Q( -rshift ) */
                    CAb[ k ] = silk_SMLAWW( CAb[ k ], tmp2,
                        silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n + k - 1 ], -rshifts - 1 ) ); /* Q( -rshift ) */
                }
            }
        }

        /* Calculate nominator and denominator for the next order reflection (parcor) coefficient */
        tmp1 = C_first_row[ n ];                                                                        /* Q( -rshifts ) */
        tmp2 = C_last_row[ n ];                                                                         /* Q( -rshifts ) */
        num  = 0;                                                                                       /* Q( -rshifts ) */
        nrg  = silk_ADD32( CAb[ 0 ], CAf[ 0 ] );                                                        /* Q( 1-rshifts ) */
        for( k = 0; k < n; k++ ) {
            Atmp_QA = Af_QA[ k ];
            lz = silk_CLZ32( silk_abs( Atmp_QA ) ) - 1;
            lz = silk_min( 32 - QA, lz );
            Atmp1 = silk_LSHIFT32( Atmp_QA, lz );                                                       /* Q( QA + lz ) */

            tmp1 = silk_ADD_LSHIFT32( tmp1, silk_SMMUL( C_last_row[  n - k - 1 ], Atmp1 ), 32 - QA - lz );  /* Q( -rshifts ) */
            tmp2 = silk_ADD_LSHIFT32( tmp2, silk_SMMUL( C_first_row[ n - k - 1 ], Atmp1 ), 32 - QA - lz );  /* Q( -rshifts ) */
            num  = silk_ADD_LSHIFT32( num,  silk_SMMUL( CAb[ n - k ],             Atmp1 ), 32 - QA - lz );  /* Q( -rshifts ) */
            nrg  = silk_ADD_LSHIFT32( nrg,  silk_SMMUL( silk_ADD32( CAb[ k + 1 ], CAf[ k + 1 ] ),
                                                                                Atmp1 ), 32 - QA - lz );    /* Q( 1-rshifts ) */
        }
        CAf[ n + 1 ] = tmp1;                                                                            /* Q( -rshifts ) */
        CAb[ n + 1 ] = tmp2;                                                                            /* Q( -rshifts ) */
        num = silk_ADD32( num, tmp2 );                                                                  /* Q( -rshifts ) */
        num = silk_LSHIFT32( -num, 1 );                                                                 /* Q( 1-rshifts ) */

        /* Calculate the next order reflection (parcor) coefficient */
        if( silk_abs( num ) < nrg ) {
            rc_Q31 = silk_DIV32_varQ( num, nrg, 31 );
        } else {
            rc_Q31 = ( num > 0 ) ? silk_int32_MAX : silk_int32_MIN;
        }

        /* Update inverse prediction gain */
        tmp1 = ( (opus_int32)1 << 30 ) - silk_SMMUL( rc_Q31, rc_Q31 );
        tmp1 = silk_LSHIFT( silk_SMMUL( invGain_Q30, tmp1 ), 2 );
        if( tmp1 <= minInvGain_Q30 ) {
            /* Max prediction gain exceeded; set reflection coefficient such that max prediction gain is exactly hit */
            tmp2 = ( (opus_int32)1 << 30 ) - silk_DIV32_varQ( minInvGain_Q30, invGain_Q30, 30 );            /* Q30 */
            rc_Q31 = silk_SQRT_APPROX( tmp2 );                                                  /* Q15 */
            if( rc_Q31 > 0 ) {
                /* Newton-Raphson iteration */
                rc_Q31 = silk_RSHIFT32( rc_Q31 + silk_DIV32( tmp2, rc_Q31 ), 1 );                       /* Q15 */
                rc_Q31 = silk_LSHIFT32( rc_Q31, 16 );                                                   /* Q31 */
                if( num < 0 ) {
                    /* Ensure adjusted reflection coefficients has the original sign */
                    rc_Q31 = -rc_Q31;
                }
            }
            invGain_Q30 = minInvGain_Q30;
            reached_max_gain = 1;
        } else {
            invGain_Q30 = tmp1;
        }

        /* Update the AR coefficients */
        for( k = 0; k < (n + 1) >> 1; k++ ) {
            tmp1 = Af_QA[ k ];                                                                  /* QA */
            tmp2 = Af_QA[ n - k - 1 ];                                                          /* QA */
            Af_QA[ k ]         = silk_ADD_LSHIFT32( tmp1, silk_SMMUL( tmp2, rc_Q31 ), 1 );      /* QA */
            Af_QA[ n - k - 1 ] = silk_ADD_LSHIFT32( tmp2, silk_SMMUL( tmp1, rc_Q31 ), 1 );      /* QA */
        }
        Af_QA[ n ] = silk_RSHIFT32( rc_Q31, 31 - QA );                                          /* QA */

        if( reached_max_gain ) {
            /* Reached max prediction gain; set remaining coefficients to zero and exit loop */
            for( k = n + 1; k < D; k++ ) {
                Af_QA[ k ] = 0;
            }
            break;
        }

        /* Update C * Af and C * Ab */
        for( k = 0; k <= n + 1; k++ ) {
            tmp1 = CAf[ k ];                                                                    /* Q( -rshifts ) */
            tmp2 = CAb[ n - k + 1 ];                                                            /* Q( -rshifts ) */
            CAf[ k ]         = silk_ADD_LSHIFT32( tmp1, silk_SMMUL( tmp2, rc_Q31 ), 1 );        /* Q( -rshifts ) */
            CAb[ n - k + 1 ] = silk_ADD_LSHIFT32( tmp2, silk_SMMUL( tmp1, rc_Q31 ), 1 );        /* Q( -rshifts ) */
        }
    }

    if( reached_max_gain ) {
        for( k = 0; k < D; k++ ) {
            /* Scale coefficients */
            A_Q16[ k ] = -silk_RSHIFT_ROUND( Af_QA[ k ], QA - 16 );
        }
        /* Subtract energy of preceding samples from C0 */
        if( rshifts > 0 ) {
            for( s = 0; s < nb_subfr; s++ ) {
                x_ptr = x + s * subfr_length;
                C0 -= (opus_int32)silk_RSHIFT64( silk_inner_prod16( x_ptr, x_ptr, D, arch ), rshifts );
            }
        } else {
            for( s = 0; s < nb_subfr; s++ ) {
                x_ptr = x + s * subfr_length;
                C0 -= silk_LSHIFT32( silk_inner_prod_aligned( x_ptr, x_ptr, D, arch), -rshifts);
            }
        }
        /* Approximate residual energy */
        *res_nrg = silk_LSHIFT( silk_SMMUL( invGain_Q30, C0 ), 2 );
        *res_nrg_Q = -rshifts;
    } else {
        /* Return residual energy */
        nrg  = CAf[ 0 ];                                                                            /* Q( -rshifts ) */
        tmp1 = (opus_int32)1 << 16;                                                                             /* Q16 */
        for( k = 0; k < D; k++ ) {
            Atmp1 = silk_RSHIFT_ROUND( Af_QA[ k ], QA - 16 );                                       /* Q16 */
            nrg  = silk_SMLAWW( nrg, CAf[ k + 1 ], Atmp1 );                                         /* Q( -rshifts ) */
            tmp1 = silk_SMLAWW( tmp1, Atmp1, Atmp1 );                                               /* Q16 */
            A_Q16[ k ] = -Atmp1;
        }
        *res_nrg = silk_SMLAWW( nrg, silk_SMMUL( SILK_FIX_CONST( FIND_LPC_COND_FAC, 32 ), C0 ), -tmp1 );/* Q( -rshifts ) */
        *res_nrg_Q = -rshifts;
    }
}

Coverage Report

Created: 2026-02-14 07:28

Line	Count	Source
1		/***********************************************************************
2		Copyright (c) 2006-2011, Skype Limited. All rights reserved.
3		Redistribution and use in source and binary forms, with or without
4		modification, are permitted provided that the following conditions
5		are met:
6		- Redistributions of source code must retain the above copyright notice,
7		this list of conditions and the following disclaimer.
8		- Redistributions in binary form must reproduce the above copyright
9		notice, this list of conditions and the following disclaimer in the
10		documentation and/or other materials provided with the distribution.
11		- Neither the name of Internet Society, IETF or IETF Trust, nor the
12		names of specific contributors, may be used to endorse or promote
13		products derived from this software without specific prior written
14		permission.
15		THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
16		AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
17		IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
18		ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
19		LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
20		CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
21		SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
22		INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
23		CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
24		ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
25		POSSIBILITY OF SUCH DAMAGE.
26		***********************************************************************/
27
28		#ifdef HAVE_CONFIG_H
29		#include "config.h"
30		#endif
31
32		#include "SigProc_FIX.h"
33		#include "define.h"
34		#include "tuning_parameters.h"
35		#include "pitch.h"
36
37		#define MAX_FRAME_SIZE 384 /* subfr_length * nb_subfr = ( 0.005 * 16000 + 16 ) * 4 = 384 */
38
39	27.4M	#define QA 25
40	27.4M	#define N_BITS_HEAD_ROOM 3
41	54.1M	#define MIN_RSHIFTS -16
42	27.4M	#define MAX_RSHIFTS (32 - QA)
43
44		/* Compute reflection coefficients from input signal */
45		void silk_burg_modified_c(
46		opus_int32 res_nrg, / O Residual energy */
47		opus_int res_nrg_Q, / O Residual energy Q value */
48		opus_int32 A_Q16[], /* O Prediction coefficients (length order) */
49		const opus_int16 x[], /* I Input signal, length: nb_subfr * ( D + subfr_length ) */
50		const opus_int32 minInvGain_Q30, /* I Inverse of max prediction gain */
51		const opus_int subfr_length, /* I Input signal subframe length (incl. D preceding samples) */
52		const opus_int nb_subfr, /* I Number of subframes stacked in x */
53		const opus_int D, /* I Order */
54		int arch /* I Run-time architecture */
55		)
56	27.4M	{
57	27.4M	opus_int k, n, s, lz, rshifts, reached_max_gain;
58	27.4M	opus_int32 C0, num, nrg, rc_Q31, invGain_Q30, Atmp_QA, Atmp1, tmp1, tmp2, x1, x2;
59	27.4M	const opus_int16 *x_ptr;
60	27.4M	opus_int32 C_first_row[ SILK_MAX_ORDER_LPC ];
61	27.4M	opus_int32 C_last_row[ SILK_MAX_ORDER_LPC ];
62	27.4M	opus_int32 Af_QA[ SILK_MAX_ORDER_LPC ];
63	27.4M	opus_int32 CAf[ SILK_MAX_ORDER_LPC + 1 ];
64	27.4M	opus_int32 CAb[ SILK_MAX_ORDER_LPC + 1 ];
65	27.4M	opus_int32 xcorr[ SILK_MAX_ORDER_LPC ];
66	27.4M	opus_int64 C0_64;
67
68	27.4M	celt_assert( subfr_length * nb_subfr <= MAX_FRAME_SIZE );
69
70		/* Compute autocorrelations, added over subframes */
71	27.4M	C0_64 = silk_inner_prod16( x, x, subfr_length*nb_subfr, arch );
72	27.4M	lz = silk_CLZ64(C0_64);
73	27.4M	rshifts = 32 + 1 + N_BITS_HEAD_ROOM - lz;
74	27.4M	if (rshifts > MAX_RSHIFTS) rshifts = MAX_RSHIFTS;
75	27.4M	if (rshifts < MIN_RSHIFTS) rshifts = MIN_RSHIFTS;
76
77	27.4M	if (rshifts > 0) {
78	36.8k	C0 = (opus_int32)silk_RSHIFT64(C0_64, rshifts );
79	27.4M	} else {
80	27.4M	C0 = silk_LSHIFT32((opus_int32)C0_64, -rshifts );
81	27.4M	}
82
83	27.4M	CAb[ 0 ] = CAf[ 0 ] = C0 + silk_SMMUL( SILK_FIX_CONST( FIND_LPC_COND_FAC, 32 ), C0 ) + 1; /* Q(-rshifts) */
84	27.4M	silk_memset( C_first_row, 0, SILK_MAX_ORDER_LPC * sizeof( opus_int32 ) );
85	27.4M	if( rshifts > 0 ) {
86	170k	for( s = 0; s < nb_subfr; s++ ) {
87	134k	x_ptr = x + s * subfr_length;
88	1.56M	for( n = 1; n < D + 1; n++ ) {
89	1.43M	C_first_row[ n - 1 ] += (opus_int32)silk_RSHIFT64(
90	1.43M	silk_inner_prod16( x_ptr, x_ptr + n, subfr_length - n, arch ), rshifts );
91	1.43M	}
92	134k	}
93	27.4M	} else {
94	113M	for( s = 0; s < nb_subfr; s++ ) {
95	85.8M	int i;
96	85.8M	opus_int32 d;
97	85.8M	x_ptr = x + s * subfr_length;
98	85.8M	celt_pitch_xcorr(x_ptr, x_ptr + 1, xcorr, subfr_length - D, D, arch );
99	969M	for( n = 1; n < D + 1; n++ ) {
100	5.07G	for ( i = n + subfr_length - D, d = 0; i < subfr_length; i++ )
101	4.18G	d = MAC16_16( d, x_ptr[ i ], x_ptr[ i - n ] );
102	884M	xcorr[ n - 1 ] += d;
103	884M	}
104	969M	for( n = 1; n < D + 1; n++ ) {
105	884M	C_first_row[ n - 1 ] += silk_LSHIFT32( xcorr[ n - 1 ], -rshifts );
106	884M	}
107	85.8M	}
108	27.4M	}
109	27.4M	silk_memcpy( C_last_row, C_first_row, SILK_MAX_ORDER_LPC * sizeof( opus_int32 ) );
110
111		/* Initialize */
112	27.4M	CAb[ 0 ] = CAf[ 0 ] = C0 + silk_SMMUL( SILK_FIX_CONST( FIND_LPC_COND_FAC, 32 ), C0 ) + 1; /* Q(-rshifts) */
113
114	27.4M	invGain_Q30 = (opus_int32)1 << 30;
115	27.4M	reached_max_gain = 0;
116	301M	for( n = 0; n < D; n++ ) {
117		/* Update first row of correlation matrix (without first element) */
118		/* Update last row of correlation matrix (without last element, stored in reversed order) */
119		/* Update C * Af */
120		/* Update C * flipud(Af) (stored in reversed order) */
121	274M	if( rshifts > -2 ) {
122	4.56M	for( s = 0; s < nb_subfr; s++ ) {
123	3.56M	x_ptr = x + s * subfr_length;
124	3.56M	x1 = -silk_LSHIFT32( (opus_int32)x_ptr[ n ], 16 - rshifts ); /* Q(16-rshifts) */
125	3.56M	x2 = -silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n - 1 ], 16 - rshifts ); /* Q(16-rshifts) */
126	3.56M	tmp1 = silk_LSHIFT32( (opus_int32)x_ptr[ n ], QA - 16 ); /* Q(QA-16) */
127	3.56M	tmp2 = silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n - 1 ], QA - 16 ); /* Q(QA-16) */
128	21.2M	for( k = 0; k < n; k++ ) {
129	17.7M	C_first_row[ k ] = silk_SMLAWB( C_first_row[ k ], x1, x_ptr[ n - k - 1 ] ); /* Q( -rshifts ) */
130	17.7M	C_last_row[ k ] = silk_SMLAWB( C_last_row[ k ], x2, x_ptr[ subfr_length - n + k ] ); /* Q( -rshifts ) */
131	17.7M	Atmp_QA = Af_QA[ k ];
132	17.7M	tmp1 = silk_SMLAWB( tmp1, Atmp_QA, x_ptr[ n - k - 1 ] ); /* Q(QA-16) */
133	17.7M	tmp2 = silk_SMLAWB( tmp2, Atmp_QA, x_ptr[ subfr_length - n + k ] ); /* Q(QA-16) */
134	17.7M	}
135	3.56M	tmp1 = silk_LSHIFT32( -tmp1, 32 - QA - rshifts ); /* Q(16-rshifts) */
136	3.56M	tmp2 = silk_LSHIFT32( -tmp2, 32 - QA - rshifts ); /* Q(16-rshifts) */
137	24.8M	for( k = 0; k <= n; k++ ) {
138	21.2M	CAf[ k ] = silk_SMLAWB( CAf[ k ], tmp1, x_ptr[ n - k ] ); /* Q( -rshift ) */
139	21.2M	CAb[ k ] = silk_SMLAWB( CAb[ k ], tmp2, x_ptr[ subfr_length - n + k - 1 ] ); /* Q( -rshift ) */
140	21.2M	}
141	3.56M	}
142	273M	} else {
143	1.12G	for( s = 0; s < nb_subfr; s++ ) {
144	856M	x_ptr = x + s * subfr_length;
145	856M	x1 = -silk_LSHIFT32( (opus_int32)x_ptr[ n ], -rshifts ); /* Q( -rshifts ) */
146	856M	x2 = -silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n - 1 ], -rshifts ); /* Q( -rshifts ) */
147	856M	tmp1 = silk_LSHIFT32( (opus_int32)x_ptr[ n ], 17 ); /* Q17 */
148	856M	tmp2 = silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n - 1 ], 17 ); /* Q17 */
149	4.88G	for( k = 0; k < n; k++ ) {
150	4.02G	C_first_row[ k ] = silk_MLA( C_first_row[ k ], x1, x_ptr[ n - k - 1 ] ); /* Q( -rshifts ) */
151	4.02G	C_last_row[ k ] = silk_MLA( C_last_row[ k ], x2, x_ptr[ subfr_length - n + k ] ); /* Q( -rshifts ) */
152	4.02G	Atmp1 = silk_RSHIFT_ROUND( Af_QA[ k ], QA - 17 ); /* Q17 */
153		/* We sometimes get overflows in the multiplications (even beyond +/- 2^32),
154		but they cancel each other and the real result seems to always fit in a 32-bit
155		signed integer. This was determined experimentally, not theoretically (unfortunately). */
156	4.02G	tmp1 = silk_MLA_ovflw( tmp1, x_ptr[ n - k - 1 ], Atmp1 ); /* Q17 */
157	4.02G	tmp2 = silk_MLA_ovflw( tmp2, x_ptr[ subfr_length - n + k ], Atmp1 ); /* Q17 */
158	4.02G	}
159	856M	tmp1 = -tmp1; /* Q17 */
160	856M	tmp2 = -tmp2; /* Q17 */
161	5.73G	for( k = 0; k <= n; k++ ) {
162	4.88G	CAf[ k ] = silk_SMLAWW( CAf[ k ], tmp1,
163	4.88G	silk_LSHIFT32( (opus_int32)x_ptr[ n - k ], -rshifts - 1 ) ); /* Q( -rshift ) */
164	4.88G	CAb[ k ] = silk_SMLAWW( CAb[ k ], tmp2,
165	4.88G	silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n + k - 1 ], -rshifts - 1 ) ); /* Q( -rshift ) */
166	4.88G	}
167	856M	}
168	273M	}
169
170		/* Calculate nominator and denominator for the next order reflection (parcor) coefficient */
171	274M	tmp1 = C_first_row[ n ]; /* Q( -rshifts ) */
172	274M	tmp2 = C_last_row[ n ]; /* Q( -rshifts ) */
173	274M	num = 0; /* Q( -rshifts ) */
174	274M	nrg = silk_ADD32( CAb[ 0 ], CAf[ 0 ] ); /* Q( 1-rshifts ) */
175	1.56G	for( k = 0; k < n; k++ ) {
176	1.28G	Atmp_QA = Af_QA[ k ];
177	1.28G	lz = silk_CLZ32( silk_abs( Atmp_QA ) ) - 1;
178	1.28G	lz = silk_min( 32 - QA, lz );
179	1.28G	Atmp1 = silk_LSHIFT32( Atmp_QA, lz ); /* Q( QA + lz ) */
180
181	1.28G	tmp1 = silk_ADD_LSHIFT32( tmp1, silk_SMMUL( C_last_row[ n - k - 1 ], Atmp1 ), 32 - QA - lz ); /* Q( -rshifts ) */
182	1.28G	tmp2 = silk_ADD_LSHIFT32( tmp2, silk_SMMUL( C_first_row[ n - k - 1 ], Atmp1 ), 32 - QA - lz ); /* Q( -rshifts ) */
183	1.28G	num = silk_ADD_LSHIFT32( num, silk_SMMUL( CAb[ n - k ], Atmp1 ), 32 - QA - lz ); /* Q( -rshifts ) */
184	1.28G	nrg = silk_ADD_LSHIFT32( nrg, silk_SMMUL( silk_ADD32( CAb[ k + 1 ], CAf[ k + 1 ] ),
185	1.28G	Atmp1 ), 32 - QA - lz ); /* Q( 1-rshifts ) */
186	1.28G	}
187	274M	CAf[ n + 1 ] = tmp1; /* Q( -rshifts ) */
188	274M	CAb[ n + 1 ] = tmp2; /* Q( -rshifts ) */
189	274M	num = silk_ADD32( num, tmp2 ); /* Q( -rshifts ) */
190	274M	num = silk_LSHIFT32( -num, 1 ); /* Q( 1-rshifts ) */
191
192		/* Calculate the next order reflection (parcor) coefficient */
193	274M	if( silk_abs( num ) < nrg ) {
194	274M	rc_Q31 = silk_DIV32_varQ( num, nrg, 31 );
195	274M	} else {
196	1	rc_Q31 = ( num > 0 ) ? silk_int32_MAX : silk_int32_MIN;
197	1	}
198
199		/* Update inverse prediction gain */
200	274M	tmp1 = ( (opus_int32)1 << 30 ) - silk_SMMUL( rc_Q31, rc_Q31 );
201	274M	tmp1 = silk_LSHIFT( silk_SMMUL( invGain_Q30, tmp1 ), 2 );
202	274M	if( tmp1 <= minInvGain_Q30 ) {
203		/* Max prediction gain exceeded; set reflection coefficient such that max prediction gain is exactly hit */
204	794k	tmp2 = ( (opus_int32)1 << 30 ) - silk_DIV32_varQ( minInvGain_Q30, invGain_Q30, 30 ); /* Q30 */
205	794k	rc_Q31 = silk_SQRT_APPROX( tmp2 ); /* Q15 */
206	794k	if( rc_Q31 > 0 ) {
207		/* Newton-Raphson iteration */
208	794k	rc_Q31 = silk_RSHIFT32( rc_Q31 + silk_DIV32( tmp2, rc_Q31 ), 1 ); /* Q15 */
209	794k	rc_Q31 = silk_LSHIFT32( rc_Q31, 16 ); /* Q31 */
210	794k	if( num < 0 ) {
211		/* Ensure adjusted reflection coefficients has the original sign */
212	776k	rc_Q31 = -rc_Q31;
213	776k	}
214	794k	}
215	794k	invGain_Q30 = minInvGain_Q30;
216	794k	reached_max_gain = 1;
217	273M	} else {
218	273M	invGain_Q30 = tmp1;
219	273M	}
220
221		/* Update the AR coefficients */
222	988M	for( k = 0; k < (n + 1) >> 1; k++ ) {
223	713M	tmp1 = Af_QA[ k ]; /* QA */
224	713M	tmp2 = Af_QA[ n - k - 1 ]; /* QA */
225	713M	Af_QA[ k ] = silk_ADD_LSHIFT32( tmp1, silk_SMMUL( tmp2, rc_Q31 ), 1 ); /* QA */
226	713M	Af_QA[ n - k - 1 ] = silk_ADD_LSHIFT32( tmp2, silk_SMMUL( tmp1, rc_Q31 ), 1 ); /* QA */
227	713M	}
228	274M	Af_QA[ n ] = silk_RSHIFT32( rc_Q31, 31 - QA ); /* QA */
229
230	274M	if( reached_max_gain ) {
231		/* Reached max prediction gain; set remaining coefficients to zero and exit loop */
232	8.65M	for( k = n + 1; k < D; k++ ) {
233	7.86M	Af_QA[ k ] = 0;
234	7.86M	}
235	794k	break;
236	794k	}
237
238		/* Update C * Af and C * Ab */
239	2.11G	for( k = 0; k <= n + 1; k++ ) {
240	1.83G	tmp1 = CAf[ k ]; /* Q( -rshifts ) */
241	1.83G	tmp2 = CAb[ n - k + 1 ]; /* Q( -rshifts ) */
242	1.83G	CAf[ k ] = silk_ADD_LSHIFT32( tmp1, silk_SMMUL( tmp2, rc_Q31 ), 1 ); /* Q( -rshifts ) */
243	1.83G	CAb[ n - k + 1 ] = silk_ADD_LSHIFT32( tmp2, silk_SMMUL( tmp1, rc_Q31 ), 1 ); /* Q( -rshifts ) */
244	1.83G	}
245	273M	}
246
247	27.4M	if( reached_max_gain ) {
248	9.61M	for( k = 0; k < D; k++ ) {
249		/* Scale coefficients */
250	8.82M	A_Q16[ k ] = -silk_RSHIFT_ROUND( Af_QA[ k ], QA - 16 );
251	8.82M	}
252		/* Subtract energy of preceding samples from C0 */
253	794k	if( rshifts > 0 ) {
254	67.8k	for( s = 0; s < nb_subfr; s++ ) {
255	53.5k	x_ptr = x + s * subfr_length;
256	53.5k	C0 -= (opus_int32)silk_RSHIFT64( silk_inner_prod16( x_ptr, x_ptr, D, arch ), rshifts );
257	53.5k	}
258	779k	} else {
259	3.32M	for( s = 0; s < nb_subfr; s++ ) {
260	2.54M	x_ptr = x + s * subfr_length;
261	2.54M	C0 -= silk_LSHIFT32( silk_inner_prod_aligned( x_ptr, x_ptr, D, arch), -rshifts);
262	2.54M	}
263	779k	}
264		/* Approximate residual energy */
265	794k	*res_nrg = silk_LSHIFT( silk_SMMUL( invGain_Q30, C0 ), 2 );
266	794k	*res_nrg_Q = -rshifts;
267	26.6M	} else {
268		/* Return residual energy */
269	26.6M	nrg = CAf[ 0 ]; /* Q( -rshifts ) */
270	26.6M	tmp1 = (opus_int32)1 << 16; /* Q16 */
271	300M	for( k = 0; k < D; k++ ) {
272	273M	Atmp1 = silk_RSHIFT_ROUND( Af_QA[ k ], QA - 16 ); /* Q16 */
273	273M	nrg = silk_SMLAWW( nrg, CAf[ k + 1 ], Atmp1 ); /* Q( -rshifts ) */
274	273M	tmp1 = silk_SMLAWW( tmp1, Atmp1, Atmp1 ); /* Q16 */
275	273M	A_Q16[ k ] = -Atmp1;
276	273M	}
277	26.6M	res_nrg = silk_SMLAWW( nrg, silk_SMMUL( SILK_FIX_CONST( FIND_LPC_COND_FAC, 32 ), C0 ), -tmp1 );/ Q( -rshifts ) */
278	26.6M	*res_nrg_Q = -rshifts;
279	26.6M	}
280	27.4M	}