Fully ported prepare_cldfb_params() - still in draft

void ivas_reverb_interpolate_energies_fx(

    const Word16 input_table_size,

    const Word32 *pInput_fc,    // input in Q16

    const Word32 *pInput_ene_l, // input in Q28

    const Word32 *pInput_ene_r, // input in Q28

    const Word32 *pInput_ene_l, // input in Q28 //Changed to Q31

    const Word32 *pInput_ene_r, // input in Q28 //Changes to q31

    const Word16 output_table_size,

    const Word32 *pOutput_fc, // Q16

    Word32 *pOutput_ene_l,      // output in Q28

    input_idx = 0;

    input_idx_next = 0;

    move16();

    move16();

    FOR( output_idx = 0; output_idx < output_table_size; output_idx++ )

    {

        {

            input_idx = 0;

            move16();

            input_idx_next = 1;

            input_idx_next = 0;

            move16();

            rel_offset = 0;

            move32();

            /* otherwise use linear interpolation between 2 consecutive points in the input table */

            ELSE

            {

                WHILE( GT_32( pOutput_fc[output_idx], pInput_fc[input_idx_next] ) )

                WHILE( GT_32( pOutput_fc[output_idx], pInput_fc[input_idx + 1] ) )

                {

                    input_idx = add( input_idx, 1 );

                }

                input_idx_next = add( input_idx, 1 );

                // Rel_offset

                rel_offset = BASOP_Util_Divide3232_Scale( L_sub( pOutput_fc[output_idx], pInput_fc[input_idx] ), L_sub( pInput_fc[input_idx_next], pInput_fc[input_idx] ), &rel_offset_e ); // Q15

                rel_offset = L_shl_sat( rel_offset, rel_offset_e );

                rel_offset = BASOP_Util_Divide3232_Scale( L_sub( pOutput_fc[output_idx], pInput_fc[input_idx] ), L_sub( pInput_fc[input_idx + 1], pInput_fc[input_idx] ), &rel_offset_e ); // Q15

                rel_offset = L_shl_sat( rel_offset, add(16, rel_offset_e ) );

                rel_offset_e = 0;

                move16();

            }

        }

    const float dist = DEFAULT_SRC_DIST;

    const float dmx_gain_2 = 4.0f * EVS_PI * dist * dist / 0.001f;*/

    #define CLDFB_BAND_WIDTH (400) // TODO: Fixme!

    Word32 output_fc_fx[CLDFB_NO_CHANNELS_MAX];

    Word32 output_t60_fx[CLDFB_NO_CHANNELS_MAX];

    Word32 output_ene_fx[CLDFB_NO_CHANNELS_MAX];

    Word32 delay_diff_fx, ln_1e6_inverted_fx, L_tmp;

    const Word32 dmx_gain_2_fx = 1852986624; // Q16

    Word32 exp_argument_fx, exp_argument_e, tmp, tmp_flag;

    Word16 pow_exp;

    Word16 pow_exp, tmp_exp;

    Word32 tmp, exp_argument_fx;

    Word32 avg_pwr_left_fx[CLDFB_NO_CHANNELS_MAX];

    Word32 avg_pwr_right_fx[CLDFB_NO_CHANNELS_MAX];

    for ( idx = 0; idx < CLDFB_NO_CHANNELS_MAX; idx++ )

    FOR ( idx = 0; idx < CLDFB_NO_CHANNELS_MAX; idx++ )

    {

        output_fc_fx[idx] = (Word32) ( ( idx + 0.5f ) * ( MAX_SAMPLING_RATE / ( 2 * CLDFB_NO_CHANNELS_MAX ) ) ) * ONE_IN_Q16;

        //output_fc_fx[idx] = add( shr(CLDFB_BAND_WIDTH, 1), mult(CLDFB_BAND_WIDTH, idx));

    }

    // TODO: Replace all these dedicated interpolate functions with just one taking a single input vector and map it to a single output vector. Period.

    for ( idx = 0; idx < CLDFB_NO_CHANNELS_MAX; idx++ )

    {

        L_tmp = Mpy_32_32( output_t60_fx[idx], ln_1e6_inverted_fx );// exp = pRt60_e[bin_idx] + 0

        exp_argument_fx = BASOP_Util_Divide3232_Scale( delay_diff_fx, L_tmp, &exp_argument_e ); //exp_argument_fx in Q14

        //exp_argument_fx = L_shl_sat( exp_argument_fx, exp_argument_e ); //exp_argument_fx converted to Q16

        /* Test Segment Begin*/

        //exp_argument_fx = shl_sat( exp_argument_fx, exp_argument_e ); // Q22

        /* Test Segment End*/

        L_tmp = Mpy_32_32( output_t60_fx[idx], ln_1e6_inverted_fx );// L_tmp in Q26 // exp = pRt60_e[bin_idx] + 0

        //exp_argument_e = add( exp_argument_e, sub( 4, pRt60_e[idx] ) ); // Q27 -> e4 -> This step must be adjusted accordingly as per the above two steps

        exp_argument_fx = BASOP_Util_Divide3232_Scale_newton( delay_diff_fx, L_tmp, &tmp_exp );//exp_argument_L_fx in Q30 (not alwys, because tmp_exp values are not constant) //Scale() returns values in the highest possible precision

        exp_argument_fx = L_shr_sat( exp_argument_fx, sub(6,tmp_exp));//exp_argument_L_fx in Q26

        /* Limit exponent to approx +/-100 dB in case of incoherent value of delay_diff, to prevent overflow */

        // 23 in Q26

        // Replace by 23 in Q14, only for _fx. 23 in Q14 = 376832

        // Pseudocode

        // IF (GT_32(exp_argument_fx, 376832 ))

        /*{

            assign Q14 equivalent of Left_Shift(23552,5) to exp_argument_fx 

          }

          IF (LT_32(exp_argument_fx, 0))

          {

        IF( GT_32( exp_argument_fx, 1543503872 ) ) //23 in Q26

        {

                If(LT_32(negate(exp_argument_fx), 376832))

                {

                    assign Q14 equivalent of Left_Shift(-23552,5) to exp_argument_fx

                }

            exp_argument_fx = 1543503872;

        }

        IF( LT_32( exp_argument_fx, -1543503872 ) )

        {

              exp_argument_fx = -1543503872;

        }

        */

        // TODO: Convert the capping to the limits

        //tmp_flag = BASOP_Util_Cmp_Mant32Exp( L_deposit_h( exp_argument_fx ), exp_argument_e, 1543503872, 5 );

        //IF( tmp_flag > 0 )

        //{

        //    exp_argument_fx = 23552;

        //    move16();

        //    exp_argument_e = 5; //exponent is 5 because 5 is the offset from Q31 to Q26

        //    move16(); // Same value in Q14 = 184??

        //}

        //tmp_flag = BASOP_Util_Cmp_Mant32Exp( L_deposit_h( exp_argument_fx ), exp_argument_e, 0, 31 );

        //IF( tmp_flag < 0 )

        //{

        //    tmp_flag = BASOP_Util_Cmp_Mant32Exp( L_deposit_h( negate( exp_argument_fx ) ), exp_argument_e, 1543503872, 5 );

        //    IF( tmp_flag < 0 )

        //    {

        //        exp_argument_fx = -23552;

        //        move16();

        //        exp_argument_e = 5;

        //        move16();

        //    }

        //}

        Word16 tmp_exp;

        /* expf(exp_argument) -> pow(2, log2(e) * exp_argument) */

        tmp = mult( 23637, exp_argument_fx ); // Q13 : mult(Q14,Q14) results in Q13 // this step denotes x * 1/ln(2)

        //tmp_exp = add( exp_argument_e, 1 ); // exp_argument_e + 1

        //L_tmp = BASOP_util_Pow2( L_deposit_h( tmp ), tmp_exp, &pow_exp );

        L_tmp = BASOP_util_Pow2( tmp, 0, &pow_exp );

        /*  L_tmp = mantissa of the operation

            pow_exp = exponent of the operation

            operation returns 2 ^ ( tmp * 2^(tmp_exp))

        */

        L_tmp = Mpy_32_32( L_tmp, output_ene_fx[idx] );

//        tmp_exp = add( pow_exp, output_ene_e[idx] );

        output_ene_fx[idx] = L_tmp;

        move32();

//        output_ene_e[idx] = tmp_exp;

        move16();

    }

        tmp = Mpy_32_32( 96817114, exp_argument_fx );//Q21 -> multiplication of Q26 by Q26 // x * 1/ln(2) it is ensured that exp_argument_fx is always in Q26

    // interpolate_acoustic_data() is returning t60 values in Q31.

    // Return  T60 in fixed-point

    //The below loop is not needed once interpolate_data() returns only fx values. 

    for ( int i = 0; i < CLDFB_NO_CHANNELS_MAX; i++ )

    {

        // Assign t60 to output variable

        //pOutput_t60[i] = L_shl( output_t60_fx[i], sub( output_t60_e[i], 5 ) );

        // Assign energy to output variable

        //pOutput_ene[i] = L_shl( output_ene_fx[i], sub( output_ene_e[i], 0 ) );

        // pOutput_t60[i] = L_shl( output_t60_fx[i], output_t60_e[i] );

        // pOutput_t60[i] = output_t60_fx[i] * ( 1 << output_t60_e[i] ); // Will not work - mantissa is 16 bit

     }

        L_tmp = BASOP_util_Pow2( tmp, 10, &pow_exp );//Q28 -> not always, because pow_exp is not fixed

        // tmp is a 32-bit signed integer in Q21, therefore expoenent to the power function is 10

        //L_tmp = L_shl_sat( L_tmp, pow_exp );// Q28??

        //L_tmp = L_shl( L_tmp, add( 1, pow_exp ) );

        //L_tmp = L_shl( L_tmp, pow_exp );

        L_tmp = Mpy_32_32( L_tmp, output_ene_fx[idx] ); // Q27. Q28 multiplied by Q30 gives Q27 -> Finally, energy should be in Q31 ( this is mandatory for further computations)

        L_tmp = L_shl_sat( L_tmp, add( 1, pow_exp ) );//

    /*for ( int i = 0; i < CLDFB_NO_CHANNELS_MAX; i++ ) //Should not be required in fixed-point

    {

        pOutput_t60[i] = (float) fabs( me2f( output_t60_fx[i], output_t60_e[i] ) ); -->This operation must be recreated in fixed-point to return T60

        pOutput_ene[i] = (float) fabs( me2f( output_ene_fx[i], output_ene_e[i] ) );

    }*/

        pOutput_ene[idx] = L_tmp;

        move32();

        //Assign t60 to output pointer

        pOutput_t60[idx] = output_t60_fx[idx]; //Q26

    }

    ivas_reverb_set_energies( hHrtfStatistics->average_energy_l, hHrtfStatistics->average_energy_r, output_Fs, avg_pwr_left_fx, avg_pwr_right_fx );

    ivas_reverb_set_energies( hHrtfStatistics->average_energy_l, hHrtfStatistics->average_energy_r, output_Fs, avg_pwr_left_fx, avg_pwr_right_fx ); //This function returns avg powers for left and right in Q28 -> why?? how??

    for ( idx = 0; idx < CLDFB_NO_CHANNELS_MAX; idx++ )

    {

        // Step 1 : add avg_pwr_left and avg_pwr_right (assuming same Q factor) and store in a temp variable -> use L_add : 32bit addition with saturation

        // Step 2 : scale the tmep variable by 0.5 -> can this be done by right shifting by 1?

        tmp_ene = ( L_add( avg_pwr_left_fx[idx], avg_pwr_right_fx[idx] ) ) >> 1; // Or Madd32_32()??

        tmp_ene = ( L_add( avg_pwr_left_fx[idx], avg_pwr_right_fx[idx] ) ) >> 1; //Q28

        // Step 3 : multiply temp variable by gain, with appropriate scaling for Q31 and Q16 multiplication

        tmp_ene = Mpy_32_32( tmp_ene, dmx_gain_2_fx ); // Mpy_32_32_ss multiplies two 32 bit numberrs in 64-bit with saturation or L_mult -> outputs word16 // Qfactor of tmp_ene?? dmx_gain_2_fx is in Q16

        // incomplete implentation

        // Step 4 : multiply step 3 result with output energy

        pOutput_ene[idx] =Mpy_32_32(pOutput_ene[idx] ,tmp_ene);

        pOutput_ene[idx] =Mpy_32_32(pOutput_ene[idx] ,tmp_ene);//result seems to be correct, but in Q13. This should be scaled to Q31 for correct output.

        /*floating-point