/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-05-16 07:41

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	50.3M	#define QA 25
40	50.3M	#define N_BITS_HEAD_ROOM 3
41	99.3M	#define MIN_RSHIFTS -16
42	50.3M	#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	50.3M	{
57	50.3M	opus_int k, n, s, lz, rshifts, reached_max_gain;
58	50.3M	opus_int32 C0, num, nrg, rc_Q31, invGain_Q30, Atmp_QA, Atmp1, tmp1, tmp2, x1, x2;
59	50.3M	const opus_int16 *x_ptr;
60	50.3M	opus_int32 C_first_row[ SILK_MAX_ORDER_LPC ];
61	50.3M	opus_int32 C_last_row[ SILK_MAX_ORDER_LPC ];
62	50.3M	opus_int32 Af_QA[ SILK_MAX_ORDER_LPC ];
63	50.3M	opus_int32 CAf[ SILK_MAX_ORDER_LPC + 1 ];
64	50.3M	opus_int32 CAb[ SILK_MAX_ORDER_LPC + 1 ];
65	50.3M	opus_int32 xcorr[ SILK_MAX_ORDER_LPC ];
66	50.3M	opus_int64 C0_64;
67
68	50.3M	celt_assert( subfr_length * nb_subfr <= MAX_FRAME_SIZE );
69
70		/* Compute autocorrelations, added over subframes */
71	50.3M	C0_64 = silk_inner_prod16( x, x, subfr_length*nb_subfr, arch );
72	50.3M	lz = silk_CLZ64(C0_64);
73	50.3M	rshifts = 32 + 1 + N_BITS_HEAD_ROOM - lz;
74	50.3M	if (rshifts > MAX_RSHIFTS) rshifts = MAX_RSHIFTS;
75	50.3M	if (rshifts < MIN_RSHIFTS) rshifts = MIN_RSHIFTS;
76
77	50.3M	if (rshifts > 0) {
78	71.2k	C0 = (opus_int32)silk_RSHIFT64(C0_64, rshifts );
79	50.2M	} else {
80	50.2M	C0 = silk_LSHIFT32((opus_int32)C0_64, -rshifts );
81	50.2M	}
82
83	50.3M	CAb[ 0 ] = CAf[ 0 ] = C0 + silk_SMMUL( SILK_FIX_CONST( FIND_LPC_COND_FAC, 32 ), C0 ) + 1; /* Q(-rshifts) */
84	50.3M	silk_memset( C_first_row, 0, SILK_MAX_ORDER_LPC * sizeof( opus_int32 ) );
85	50.3M	if( rshifts > 0 ) {
86	335k	for( s = 0; s < nb_subfr; s++ ) {
87	264k	x_ptr = x + s * subfr_length;
88	3.16M	for( n = 1; n < D + 1; n++ ) {
89	2.90M	C_first_row[ n - 1 ] += (opus_int32)silk_RSHIFT64(
90	2.90M	silk_inner_prod16( x_ptr, x_ptr + n, subfr_length - n, arch ), rshifts );
91	2.90M	}
92	264k	}
93	50.2M	} else {
94	210M	for( s = 0; s < nb_subfr; s++ ) {
95	160M	int i;
96	160M	opus_int32 d;
97	160M	x_ptr = x + s * subfr_length;
98	160M	celt_pitch_xcorr(x_ptr, x_ptr + 1, xcorr, subfr_length - D, D, arch );
99	1.85G	for( n = 1; n < D + 1; n++ ) {
100	10.0G	for ( i = n + subfr_length - D, d = 0; i < subfr_length; i++ )
101	8.33G	d = MAC16_16( d, x_ptr[ i ], x_ptr[ i - n ] );
102	1.69G	xcorr[ n - 1 ] += d;
103	1.69G	}
104	1.85G	for( n = 1; n < D + 1; n++ ) {
105	1.69G	C_first_row[ n - 1 ] += silk_LSHIFT32( xcorr[ n - 1 ], -rshifts );
106	1.69G	}
107	160M	}
108	50.2M	}
109	50.3M	silk_memcpy( C_last_row, C_first_row, SILK_MAX_ORDER_LPC * sizeof( opus_int32 ) );
110
111		/* Initialize */
112	50.3M	CAb[ 0 ] = CAf[ 0 ] = C0 + silk_SMMUL( SILK_FIX_CONST( FIND_LPC_COND_FAC, 32 ), C0 ) + 1; /* Q(-rshifts) */
113
114	50.3M	invGain_Q30 = (opus_int32)1 << 30;
115	50.3M	reached_max_gain = 0;
116	561M	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	513M	if( rshifts > -2 ) {
122	8.64M	for( s = 0; s < nb_subfr; s++ ) {
123	6.75M	x_ptr = x + s * subfr_length;
124	6.75M	x1 = -silk_LSHIFT32( (opus_int32)x_ptr[ n ], 16 - rshifts ); /* Q(16-rshifts) */
125	6.75M	x2 = -silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n - 1 ], 16 - rshifts ); /* Q(16-rshifts) */
126	6.75M	tmp1 = silk_LSHIFT32( (opus_int32)x_ptr[ n ], QA - 16 ); /* Q(QA-16) */
127	6.75M	tmp2 = silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n - 1 ], QA - 16 ); /* Q(QA-16) */
128	40.1M	for( k = 0; k < n; k++ ) {
129	33.4M	C_first_row[ k ] = silk_SMLAWB( C_first_row[ k ], x1, x_ptr[ n - k - 1 ] ); /* Q( -rshifts ) */
130	33.4M	C_last_row[ k ] = silk_SMLAWB( C_last_row[ k ], x2, x_ptr[ subfr_length - n + k ] ); /* Q( -rshifts ) */
131	33.4M	Atmp_QA = Af_QA[ k ];
132	33.4M	tmp1 = silk_SMLAWB( tmp1, Atmp_QA, x_ptr[ n - k - 1 ] ); /* Q(QA-16) */
133	33.4M	tmp2 = silk_SMLAWB( tmp2, Atmp_QA, x_ptr[ subfr_length - n + k ] ); /* Q(QA-16) */
134	33.4M	}
135	6.75M	tmp1 = silk_LSHIFT32( -tmp1, 32 - QA - rshifts ); /* Q(16-rshifts) */
136	6.75M	tmp2 = silk_LSHIFT32( -tmp2, 32 - QA - rshifts ); /* Q(16-rshifts) */
137	46.9M	for( k = 0; k <= n; k++ ) {
138	40.1M	CAf[ k ] = silk_SMLAWB( CAf[ k ], tmp1, x_ptr[ n - k ] ); /* Q( -rshift ) */
139	40.1M	CAb[ k ] = silk_SMLAWB( CAb[ k ], tmp2, x_ptr[ subfr_length - n + k - 1 ] ); /* Q( -rshift ) */
140	40.1M	}
141	6.75M	}
142	511M	} else {
143	2.14G	for( s = 0; s < nb_subfr; s++ ) {
144	1.63G	x_ptr = x + s * subfr_length;
145	1.63G	x1 = -silk_LSHIFT32( (opus_int32)x_ptr[ n ], -rshifts ); /* Q( -rshifts ) */
146	1.63G	x2 = -silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n - 1 ], -rshifts ); /* Q( -rshifts ) */
147	1.63G	tmp1 = silk_LSHIFT32( (opus_int32)x_ptr[ n ], 17 ); /* Q17 */
148	1.63G	tmp2 = silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n - 1 ], 17 ); /* Q17 */
149	9.59G	for( k = 0; k < n; k++ ) {
150	7.96G	C_first_row[ k ] = silk_MLA( C_first_row[ k ], x1, x_ptr[ n - k - 1 ] ); /* Q( -rshifts ) */
151	7.96G	C_last_row[ k ] = silk_MLA( C_last_row[ k ], x2, x_ptr[ subfr_length - n + k ] ); /* Q( -rshifts ) */
152	7.96G	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	7.96G	tmp1 = silk_MLA_ovflw( tmp1, x_ptr[ n - k - 1 ], Atmp1 ); /* Q17 */
157	7.96G	tmp2 = silk_MLA_ovflw( tmp2, x_ptr[ subfr_length - n + k ], Atmp1 ); /* Q17 */
158	7.96G	}
159	1.63G	tmp1 = -tmp1; /* Q17 */
160	1.63G	tmp2 = -tmp2; /* Q17 */
161	11.2G	for( k = 0; k <= n; k++ ) {
162	9.59G	CAf[ k ] = silk_SMLAWW( CAf[ k ], tmp1,
163	9.59G	silk_LSHIFT32( (opus_int32)x_ptr[ n - k ], -rshifts - 1 ) ); /* Q( -rshift ) */
164	9.59G	CAb[ k ] = silk_SMLAWW( CAb[ k ], tmp2,
165	9.59G	silk_LSHIFT32( (opus_int32)x_ptr[ subfr_length - n + k - 1 ], -rshifts - 1 ) ); /* Q( -rshift ) */
166	9.59G	}
167	1.63G	}
168	511M	}
169
170		/* Calculate nominator and denominator for the next order reflection (parcor) coefficient */
171	513M	tmp1 = C_first_row[ n ]; /* Q( -rshifts ) */
172	513M	tmp2 = C_last_row[ n ]; /* Q( -rshifts ) */
173	513M	num = 0; /* Q( -rshifts ) */
174	513M	nrg = silk_ADD32( CAb[ 0 ], CAf[ 0 ] ); /* Q( 1-rshifts ) */
175	3.01G	for( k = 0; k < n; k++ ) {
176	2.50G	Atmp_QA = Af_QA[ k ];
177	2.50G	lz = silk_CLZ32( silk_abs( Atmp_QA ) ) - 1;
178	2.50G	lz = silk_min( 32 - QA, lz );
179	2.50G	Atmp1 = silk_LSHIFT32( Atmp_QA, lz ); /* Q( QA + lz ) */
180
181	2.50G	tmp1 = silk_ADD_LSHIFT32( tmp1, silk_SMMUL( C_last_row[ n - k - 1 ], Atmp1 ), 32 - QA - lz ); /* Q( -rshifts ) */
182	2.50G	tmp2 = silk_ADD_LSHIFT32( tmp2, silk_SMMUL( C_first_row[ n - k - 1 ], Atmp1 ), 32 - QA - lz ); /* Q( -rshifts ) */
183	2.50G	num = silk_ADD_LSHIFT32( num, silk_SMMUL( CAb[ n - k ], Atmp1 ), 32 - QA - lz ); /* Q( -rshifts ) */
184	2.50G	nrg = silk_ADD_LSHIFT32( nrg, silk_SMMUL( silk_ADD32( CAb[ k + 1 ], CAf[ k + 1 ] ),
185	2.50G	Atmp1 ), 32 - QA - lz ); /* Q( 1-rshifts ) */
186	2.50G	}
187	513M	CAf[ n + 1 ] = tmp1; /* Q( -rshifts ) */
188	513M	CAb[ n + 1 ] = tmp2; /* Q( -rshifts ) */
189	513M	num = silk_ADD32( num, tmp2 ); /* Q( -rshifts ) */
190	513M	num = silk_LSHIFT32( -num, 1 ); /* Q( 1-rshifts ) */
191
192		/* Calculate the next order reflection (parcor) coefficient */
193	513M	if( silk_abs( num ) < nrg ) {
194	513M	rc_Q31 = silk_DIV32_varQ( num, nrg, 31 );
195	513M	} else {
196	1	rc_Q31 = ( num > 0 ) ? silk_int32_MAX : silk_int32_MIN;
197	1	}
198
199		/* Update inverse prediction gain */
200	513M	tmp1 = ( (opus_int32)1 << 30 ) - silk_SMMUL( rc_Q31, rc_Q31 );
201	513M	tmp1 = silk_LSHIFT( silk_SMMUL( invGain_Q30, tmp1 ), 2 );
202	513M	if( tmp1 <= minInvGain_Q30 ) {
203		/* Max prediction gain exceeded; set reflection coefficient such that max prediction gain is exactly hit */
204	1.71M	tmp2 = ( (opus_int32)1 << 30 ) - silk_DIV32_varQ( minInvGain_Q30, invGain_Q30, 30 ); /* Q30 */
205	1.71M	rc_Q31 = silk_SQRT_APPROX( tmp2 ); /* Q15 */
206	1.71M	if( rc_Q31 > 0 ) {
207		/* Newton-Raphson iteration */
208	1.71M	rc_Q31 = silk_RSHIFT32( rc_Q31 + silk_DIV32( tmp2, rc_Q31 ), 1 ); /* Q15 */
209	1.71M	rc_Q31 = silk_LSHIFT32( rc_Q31, 16 ); /* Q31 */
210	1.71M	if( num < 0 ) {
211		/* Ensure adjusted reflection coefficients has the original sign */
212	1.67M	rc_Q31 = -rc_Q31;
213	1.67M	}
214	1.71M	}
215	1.71M	invGain_Q30 = minInvGain_Q30;
216	1.71M	reached_max_gain = 1;
217	511M	} else {
218	511M	invGain_Q30 = tmp1;
219	511M	}
220
221		/* Update the AR coefficients */
222	1.89G	for( k = 0; k < (n + 1) >> 1; k++ ) {
223	1.37G	tmp1 = Af_QA[ k ]; /* QA */
224	1.37G	tmp2 = Af_QA[ n - k - 1 ]; /* QA */
225	1.37G	Af_QA[ k ] = silk_ADD_LSHIFT32( tmp1, silk_SMMUL( tmp2, rc_Q31 ), 1 ); /* QA */
226	1.37G	Af_QA[ n - k - 1 ] = silk_ADD_LSHIFT32( tmp2, silk_SMMUL( tmp1, rc_Q31 ), 1 ); /* QA */
227	1.37G	}
228	513M	Af_QA[ n ] = silk_RSHIFT32( rc_Q31, 31 - QA ); /* QA */
229
230	513M	if( reached_max_gain ) {
231		/* Reached max prediction gain; set remaining coefficients to zero and exit loop */
232	19.2M	for( k = n + 1; k < D; k++ ) {
233	17.5M	Af_QA[ k ] = 0;
234	17.5M	}
235	1.71M	break;
236	1.71M	}
237
238		/* Update C * Af and C * Ab */
239	4.03G	for( k = 0; k <= n + 1; k++ ) {
240	3.52G	tmp1 = CAf[ k ]; /* Q( -rshifts ) */
241	3.52G	tmp2 = CAb[ n - k + 1 ]; /* Q( -rshifts ) */
242	3.52G	CAf[ k ] = silk_ADD_LSHIFT32( tmp1, silk_SMMUL( tmp2, rc_Q31 ), 1 ); /* Q( -rshifts ) */
243	3.52G	CAb[ n - k + 1 ] = silk_ADD_LSHIFT32( tmp2, silk_SMMUL( tmp1, rc_Q31 ), 1 ); /* Q( -rshifts ) */
244	3.52G	}
245	511M	}
246
247	50.3M	if( reached_max_gain ) {
248	21.3M	for( k = 0; k < D; k++ ) {
249		/* Scale coefficients */
250	19.6M	A_Q16[ k ] = -silk_RSHIFT_ROUND( Af_QA[ k ], QA - 16 );
251	19.6M	}
252		/* Subtract energy of preceding samples from C0 */
253	1.71M	if( rshifts > 0 ) {
254	150k	for( s = 0; s < nb_subfr; s++ ) {
255	119k	x_ptr = x + s * subfr_length;
256	119k	C0 -= (opus_int32)silk_RSHIFT64( silk_inner_prod16( x_ptr, x_ptr, D, arch ), rshifts );
257	119k	}
258	1.68M	} else {
259	7.12M	for( s = 0; s < nb_subfr; s++ ) {
260	5.44M	x_ptr = x + s * subfr_length;
261	5.44M	C0 -= silk_LSHIFT32( silk_inner_prod_aligned( x_ptr, x_ptr, D, arch), -rshifts);
262	5.44M	}
263	1.68M	}
264		/* Approximate residual energy */
265	1.71M	*res_nrg = silk_LSHIFT( silk_SMMUL( invGain_Q30, C0 ), 2 );
266	1.71M	*res_nrg_Q = -rshifts;
267	48.6M	} else {
268		/* Return residual energy */
269	48.6M	nrg = CAf[ 0 ]; /* Q( -rshifts ) */
270	48.6M	tmp1 = (opus_int32)1 << 16; /* Q16 */
271	559M	for( k = 0; k < D; k++ ) {
272	511M	Atmp1 = silk_RSHIFT_ROUND( Af_QA[ k ], QA - 16 ); /* Q16 */
273	511M	nrg = silk_SMLAWW( nrg, CAf[ k + 1 ], Atmp1 ); /* Q( -rshifts ) */
274	511M	tmp1 = silk_SMLAWW( tmp1, Atmp1, Atmp1 ); /* Q16 */
275	511M	A_Q16[ k ] = -Atmp1;
276	511M	}
277	48.6M	res_nrg = silk_SMLAWW( nrg, silk_SMMUL( SILK_FIX_CONST( FIND_LPC_COND_FAC, 32 ), C0 ), -tmp1 );/ Q( -rshifts ) */
278	48.6M	*res_nrg_Q = -rshifts;
279	48.6M	}
280	50.3M	}