FEC_HQ_phase_ecu.c

        if ( ( *num_plocs > 0 ) && sub( *num_plocs, 3 ) < 0 )
        {
            one_peak_flag_mask = noise_fac; /* all zeroes  mask -> zero  */
        }
        if ( *num_plocs == 0 )
        {
            X[0] = 0;               /* reset DC if there are no  peaks */
            X[shr( Lprot, 1 )] = 0; /* also reset fs/2 if there are no peaks */
        }
    }

    i = 1;
    k = 0;
    im_ind = Lprot - 1;
    for ( m = 0; m < *num_plocs; m++ )
    {
        delta_corr_dn = DELTA_CORR;
        delta_corr_up = DELTA_CORR;

        if ( m > 0 )
        {
            delta_tmp = ( plocs[m] - plocs[m - 1] - 1 ) / 2;
            if ( delta_tmp < DELTA_CORR )
            {
                delta_corr_dn = delta_tmp;
            }
        }

        if ( m < *num_plocs - 1 )
        {
            delta_tmp = ( plocs[m + 1] - plocs[m] - 1 ) / 2;
            if ( delta_tmp < DELTA_CORR )
            {
                delta_corr_up = delta_tmp;
            }
        }

        /* Input Xph */
        while ( i < plocs[m] - delta_corr_dn )
        {
            *seed = own_random( seed );

            if ( *seed & 0x40 )
            {
                sin_F = sincos[*seed >> 8];
            }
            else
            {
                sin_F = -sincos[*seed >> 8];
            }

            if ( *seed & 0x80 )
            {
                cos_F = sincos[-( *seed >> 8 )];
            }
            else
            {
                cos_F = -sincos[-( *seed >> 8 )];
            }

            if ( element_mode == EVS_MONO )
            {
                tmp = ( X[i] * cos_F - X[im_ind] * sin_F );
                X[im_ind] = ( X[i] * sin_F + X[im_ind] * cos_F );
            }
            else
            {
                tmp = one_peak_flag_mask * ( X[i] * cos_F - X[im_ind] * sin_F );
                X[im_ind] = one_peak_flag_mask * ( X[i] * sin_F + X[im_ind] * cos_F );
            }

            if ( alpha[k] < 1.0f )
            {
                *seed = rand_phase( *seed, &sin_F, &cos_F );
                X[i] = alpha[k] * tmp + beta[k] * Xavg[k] * cos_F;
                X[im_ind] = alpha[k] * X[im_ind] + beta[k] * Xavg[k] * sin_F;
            }
            else
            {
                X[i] = mag_chg[k] * tmp;
                X[im_ind] *= mag_chg[k];
            }
            i++;
            im_ind--;
            if ( i >= ivas_gwlpr[k + 1] )
            {
                k++;
            }
        }

        e = plocs[m] + delta_corr_up;
        if ( e > Lprot / 2 - 1 )
        {
            e = Lprot / 2 - 1;
        }

        Xph = corr_phase[m];
#ifdef FIX_1002_DEC_PHASE_ECU_USAN_OF_PHASE
        /* BE , but with maintained dual mult of a constant,  bnut mask the  10 bits in Word32, truncation  */
#define WMC_TOOL_SKIP
        Xph_short = ( int16_t )( 0x000003ff & ( int32_t )( ( Xph * 512 ) / EVS_PI ) );
        MULT( 2 );  /* mult with constant, twice due to legacy precedence use */
        LOGIC( 1 ); /* L_and        */
        MISC( 1 );  /* extract_l   */
#undef WMC_TOOL_SKIP
#if 1
        int16_t Xph_short_orig = ( int16_t )( ( ( int32_t )( Xph * 512 / EVS_PI ) ) % 32768 ) & 0x03ff; /* BE: 2 mults, 1 modulo */     /* bad line(for WMC):  WMC-costly modulo % in use,  no rounding   */
        int16_t Xph_short_1mult_nonbe = ( int16_t )( 0x000003ff & ( int32_t )( ( Xph * ( 512.0 / EVS_PI ) ) ) );                        /* nonBE:  single mult, mask the  10 bits in Word32 ,  still no rounding */
        int16_t Xph_short_1mult_rnd_nonbe = ( int16_t )( 0x000003ff & ( ( 1L + ( ( int32_t )( Xph * ( 1024.0 / EVS_PI ) ) ) ) >> 1 ) ); /* nonBE:  single mult, add 1(round) shift1,  mask the  10 bits in Word32 , i.e. with rounding */

        assert( Xph_short == Xph_short_orig && "Xph_short == Xph_short_orig" );
#endif

#ifdef DEBUGGING

        int32_t Xph_long_new = ( int32_t )( Xph * ( 512 / EVS_PI ) ); /* mult by a single constant and float cast to int32_t */
        int32_t Xph_long_old = ( int32_t )( ( Xph * 512 ) / EVS_PI ); /* mult by 2 constants ! and float cast to int32_t , (Xph<<9) / pi */
        int32_t Xph_long = Xph_long_old;

        //assert( Xph_long_new == Xph_long_old && __func__ && " Xph_long_new != Xph_long_old" ); // mat differ due to order of division

        int32_t L_tmp1 = Xph_long % 32768;      /* counts as division:  6-25 cycles */
        int32_t L_tmp2 = Xph_long & 0x00007fff; /* L_and mask : 1 cycle */
       // assert( L_tmp1 == L_tmp2 && __func__ && " L_tmp1 != L_tmp2 " );


        Xph_long = ( Xph_long & 0x00007fff ); /*  L_and()  equivalent to  % 32768 */
        int16_t Xph_short_new1 = ( (int16_t) Xph_long ) & 0x03ff;
        int16_t Xph_short_new2 = ( int16_t )( Xph_long_old & 0x000003ff );

        //assert( Xph_short == Xph_short_new1 && __func__ && " Xph_short != Xph_short_new1" );
        //assert( Xph_short == Xph_short_new2 && __func__ && " Xph_short != Xph_short_new2" );

#if 0
            /* guarantee no wrap around to negative side of 32 bit signed integer value  */
#define WMC_TOOL_SKIP
            tmpf = Xph * ( 1.0f / PI2 );                           /* normalize         :  mult() with a constant                           */
            tmpf = ( 1024.0f * ( tmpf - floorf( tmpf ) ) ) + 0.5f; /* obtain fraction  :  floor(), sub(), mult(), rnd(), out:  [0...1024.5[ */
            int16_t Xph_short_flt_correct = ( (int16_t) tmpf ) & 0x03ff;               /* mask to [0..1023]:  s_and()                      , out:  [0...1023]   */

            MULT( 2 );
            MISC( 1 );  /* floor    */
            ADD( 2 );   /* sub, rnd */
            LOGIC( 1 ); /* s_and    */
#undef WMC_TOOL_SKIP
#endif 
#endif

#else
        Xph_short = ( int16_t )( ( ( int32_t )( Xph * 512 / EVS_PI ) ) % 32768 ) & 0x03ff;
#endif


        if ( Xph_short >= 512 )
        {
            sin_F = -sincos_t_ext[Xph_short - 512];
            if ( Xph_short < 768 )
            {
                cos_F = -sincos_t_ext[Xph_short - 512 + 256];
            }
            else
            {
                cos_F = sincos_t_ext[-Xph_short + 1024 + 256];
            }
        }
        else
        {
            sin_F = sincos_t_ext[Xph_short];
            if ( Xph_short < 256 )
            {
                cos_F = sincos_t_ext[Xph_short + 256];
            }
            else
            {
                cos_F = -sincos_t_ext[-Xph_short + 256 + 512];
            }
        }

        while ( i <= e )
        {
            mag_chg_local = mag_chg[k];

            if ( ph_dith != 0.0f )
            {
                /* Call phase randomization only when needed */
                Xph = corr_phase[m];
                *seed = own_random( seed );
                Xph += *seed * ph_dith * PHASE_DITH_SCALE; /* where ph_dith is  0..2PI,  or -2PI (in transient), bin phase scaling factor from trans_ana */

                if ( ph_dith > 0.0f )
                {
                    /* up to 6 dB additional att of peaks in non_transient longer bursts, (when peak phase is randomized) */
                    /* 0.5~= sqrt((float)pow(10.0,-6/10.0));    ph_dith= 0..2pi,--> scale=1.0 ...5 */
                    mag_chg_local *= 0.5f + ( 1.0f - ( 1.0f / PHASE_DITH ) * ph_dith ) * 0.5f;
                }

#ifdef FIX_1002_DEC_PHASE_ECU_USAN_OF_PHASE
#define WMC_TOOL_SKIP
                Xph_short = ( int16_t )( ( ( int32_t )( ( Xph * 512 ) / EVS_PI ) ) & 0x000003ff );
                MULT( 2 );  /* mult with constant, twice due to legacy precedence  */
                LOGIC( 1 ); /* L_and        */
                MISC( 1 );  /* extract_l   */
#undef WMC_TOOL_SKIP
#if 1
                int16_t Xph_short_orig_ph_dith = ( int16_t )( Xph * 512 / EVS_PI ) & 0x03ff; /* USAN bad cast,  bad dual mults */
                assert( Xph_short == Xph_short_orig_ph_dith ); 
#endif

#ifdef DEBUGGING
                /* NB! USAN undefined behaviour warning for float-> short cast adressed, no rounding required here as phase is scrambled anyway */
                /* applied also  EVS_MONO */
#endif


#else
                Xph_short = ( int16_t )( Xph * 512 / EVS_PI ) & 0x03ff;
#endif

                    if ( Xph_short >= 512 )
                {
                    sin_F = -sincos_t_ext[Xph_short - 512];
                    if ( Xph_short < 768 )
                    {
                        cos_F = -sincos_t_ext[Xph_short - 512 + 256];
                    }
                    else
                    {
                        cos_F = sincos_t_ext[-Xph_short + 1024 + 256];
                    }
                }
                else
                {
                    sin_F = sincos_t_ext[Xph_short];
                    if ( Xph_short < 256 )
                    {
                        cos_F = sincos_t_ext[Xph_short + 256];
                    }
                    else
                    {
                        cos_F = -sincos_t_ext[-Xph_short + 256 + 512];
                    }
                }
            }

            tmp = ( X[i] * cos_F - X[im_ind] * sin_F );
            X[im_ind] = ( X[i] * sin_F + X[im_ind] * cos_F );
            if ( alpha[k] < 1.0f )
            {
                alpha_local = mag_chg_local;
                beta_local = (float) ( beta_mute * sqrt( 1.0f - SQR( alpha_local ) ) );
                if ( k >= LGW32k - 1 )
                {
                    beta_local *= 0.1f;
                }
                else if ( k >= LGW16k - 1 )
                {
                    beta_local *= 0.5f;
                }

                *seed = rand_phase( *seed, &sin_F, &cos_F );
                X[i] = alpha_local * tmp + beta_local * Xavg[k] * cos_F;
                X[im_ind] = alpha_local * X[im_ind] + beta_local * Xavg[k] * sin_F;
            }
            else
            {
                X[i] = mag_chg_local * tmp;
                X[im_ind] *= mag_chg_local;
            }

            i++;
            im_ind--;
            if ( i >= ivas_gwlpr[k + 1] )
            {
                k++;
            }
        }
    }

    while ( i < Lprot / 2 )
    {
        *seed = own_random( seed );

        if ( *seed & 0x40 )
        {
            sin_F = sincos[*seed >> 8];
        }
        else
        {
            sin_F = -sincos[*seed >> 8];
        }

        if ( *seed & 0x80 )
        {
            cos_F = sincos[-( *seed >> 8 )];
        }
        else
        {
            cos_F = -sincos[-( *seed >> 8 )];
        }


        if ( element_mode == EVS_MONO )
        {
            tmp = ( X[i] * cos_F - X[im_ind] * sin_F );
            X[im_ind] = ( X[i] * sin_F + X[im_ind] * cos_F );
        }
        else
        {
            tmp = one_peak_flag_mask * ( X[i] * cos_F - X[im_ind] * sin_F );
            X[im_ind] = one_peak_flag_mask * ( X[i] * sin_F + X[im_ind] * cos_F );
        }

        if ( alpha[k] < 1.0f )
        {
            *seed = rand_phase( *seed, &sin_F, &cos_F );
            X[i] = alpha[k] * tmp + beta[k] * Xavg[k] * cos_F;
            X[im_ind] = alpha[k] * X[im_ind] + beta[k] * Xavg[k] * sin_F;
            im_ind--;
        }
        else
        {
            X[i] = mag_chg[k] * tmp;
            X[im_ind--] *= mag_chg[k];
        }
        i++;

        if ( i >= ivas_gwlpr[k + 1] )
        {
            k++;
        }
    }

    return;
}

/*--------------------------------------------------------------------------
 *  rec_wtda()
 *
 *  Windowing and TDA of reconstructed frame
 *--------------------------------------------------------------------------*/

static void rec_wtda(
    float *X,                   /* i/o: ECU frame / unwindowed ECU frame    */
    float *ecu_rec,             /* o  : Reconstructed frame in tda domain   */
    const int16_t output_frame, /* i  : Frame length                        */
    const int16_t Lprot,        /* i  : Prototype frame length              */
    const float old_dec[270],   /* i  : end of last decoded for OLA before tda and itda */
    const int16_t element_mode, /* i  : IVAS element mode                   */
    const int16_t *num_p,       /* i  : Number of peaks                     */
    const int16_t *plocs        /* i  : Peak locations                      */
)
{
    int16_t timesh;
    int16_t xf_len;
    int16_t i;
    float *p_ecu;
    float g;
    float tbl_delta;
    float xsubst_[2 * L_FRAME48k];
    const float *w_hamm;
    float *pX_start, *pX_end;
    float tmp;
    int16_t hamm_len2;
    float *pNew;
    const float *pOldW, *pNewW;
    float xfwin[NS2SA( L_FRAME48k * FRAMES_PER_SEC, N_ZERO_MDCT_NS - ( 2 * FRAME_SIZE_NS - L_PROT_NS ) / 2 )];
    const float *pOld;
    int16_t copy_len;
    int16_t ola_len;

    copy_len = NS2SA( output_frame * FRAMES_PER_SEC, ( 2 * FRAME_SIZE_NS - L_PROT_NS ) / 2 );                 /* prototype fill on each side of xsubst to fill MDCT Frame */
    ola_len = NS2SA( output_frame * FRAMES_PER_SEC, N_ZERO_MDCT_NS - ( 2 * FRAME_SIZE_NS - L_PROT_NS ) / 2 ); /* remaining lengt of LA_ZEROS to overlap add decoded with xsubst */

    if ( output_frame == L_FRAME48k )
    {
        w_hamm = w_hamm_sana48k_2;
        hamm_len2 = L_PROT_HAMM_LEN2_48k;
    }
    else if ( output_frame == L_FRAME32k )
    {
        w_hamm = w_hamm_sana32k_2;
        hamm_len2 = L_PROT_HAMM_LEN2_32k;
    }
    else
    {
        w_hamm = w_hamm_sana16k_2;
        hamm_len2 = L_PROT_HAMM_LEN2_16k;
    }

    if ( element_mode != EVS_MONO && *num_p > 0 && plocs[0] > 3 )
    {
        /* Perform inverse windowing of hammrect */
        pX_start = X;
        pX_end = X + Lprot - 1;
        for ( i = 0; i < hamm_len2; i++ )
        {
            tmp = 1.0f / *w_hamm;
            *pX_start *= tmp;
            *pX_end *= tmp;
            pX_start++;
            pX_end--;
            w_hamm++;
        }
    }

    /* extract reconstructed frame with aldo window */
    timesh = NS2SA( output_frame * FRAMES_PER_SEC, N_ZERO_MDCT_NS ) - ( 2 * output_frame - Lprot ) / 2;

    set_f( xsubst_, 0.0f, 2 * output_frame - Lprot + timesh );
    mvr2r( X, xsubst_ + 2 * output_frame - Lprot + timesh, Lprot - timesh );

    /* Copy and OLA look ahead zero part of MDCT window from decoded signal  */
    if ( element_mode != EVS_MONO )
    {
        mvr2r( old_dec, xsubst_ + NS2SA( output_frame * FRAMES_PER_SEC, N_ZERO_MDCT_NS ), copy_len ); /* also need to scale to Q0 ?? */
        pOld = old_dec + copy_len;
        pNew = xsubst_ + copy_len + NS2SA( output_frame * FRAMES_PER_SEC, N_ZERO_MDCT_NS );
        sinq( EVS_PI / ( ola_len * 2 ), 0.0f, ola_len, xfwin );
        v_mult( xfwin, xfwin, xfwin, ola_len ); /* xfwin = sin^2 of 0..pi/4 */
        pOldW = xfwin + ola_len - 1;
        pNewW = xfwin;
        for ( i = 0; i < ola_len; i++ )
        {
            *pNew = *pOld * *pOldW + *pNew * *pNewW;
            pOld += 1;
            pNew += 1;
            pOldW -= 1;
            pNewW += 1;
        }
    }
    else
    {
        /* Smoothen onset of ECU frame */
        xf_len = ( int16_t )( (float) output_frame * N_ZERO_MDCT_NS / FRAME_SIZE_NS ) - ( output_frame - Lprot / 2 );
        p_ecu = xsubst_ + 2 * output_frame - Lprot + timesh;
        tbl_delta = 64.f / xf_len; /* 64 samples = 1/4 cycle in sincos_t */
        for ( i = 0; i < xf_len; i++, p_ecu++ )
        {
            g = sincos_t[( ( int16_t )( i * tbl_delta ) )];
            g *= g;
            *p_ecu = g * ( *p_ecu );
        }
    }

    /* Apply TDA and windowing to ECU frame */
    wtda( xsubst_ + output_frame, ecu_rec, NULL, ALDO_WINDOW, ALDO_WINDOW, output_frame );

    return;
}


/*--------------------------------------------------------------------------
 *  rec_frame()
 *
 *  Frame reconstruction
 *--------------------------------------------------------------------------*/

static void rec_frame(
    float *X,                   /* i/o: FFT spectrum / IFFT of spectrum       */
    float *ecu_rec,             /* o  : Reconstructed frame in tda domain     */
    const int16_t output_frame, /* i  : Frame length                          */
    const float *old_dec,       /* i  : end of last decoded for OLA before tda and itda */
    const int16_t element_mode, /* i  : IVAS element mode                     */
    const int16_t *num_p,       /* i  : Number of peaks                       */
    const int16_t *plocs        /* i  : Peak locations                        */
)
{
    int16_t Lprot, LprotLog2 = 0;

    Lprot = 2 * output_frame * L_PROT32k / 1280;

    if ( output_frame == L_FRAME48k )
    {
        LprotLog2 = 9;
    }
    else if ( output_frame == L_FRAME32k )
    {
        LprotLog2 = 10;
    }
    else
    {
        LprotLog2 = 9;
    }

    /* extend spectrum and IDFT */
    if ( output_frame == L_FRAME48k )
    {
        ifft3( X, X, Lprot );
    }
    else
    {
        ifft_rel( X, Lprot, LprotLog2 );
    }

    rec_wtda( X, ecu_rec, output_frame, Lprot, old_dec, element_mode, num_p, plocs );

    return;
}


/*--------------------------------------------------------------------------
 *  fir_dwn()
 *
 *  FIR downsampling filter
 *--------------------------------------------------------------------------*/

static void fir_dwn(
    const float x[],         /* i  : input vector                              */
    const float h[],         /* i  : impulse response of the FIR filter        */
    float y[],               /* o  : output vector (result of filtering)       */
    const int16_t L,         /* i  : input vector size                         */
    const int16_t K,         /* i  : order of the FIR filter (K+1 coefs.)      */
    const int16_t decimation /* i  : decimation                                */
)
{
    float s;
    int16_t i, j, k, mmax;

    k = 0;

    /* do the filtering */
    for ( i = K / 2; i < L; i += decimation )
    {
        s = x[i] * h[0];

        if ( i < K )
        {
            mmax = i;
        }
        else
        {
            mmax = K;
        }

        for ( j = 1; j <= mmax; j++ )
        {
            s += h[j] * x[i - j];
        }

        y[k] = s;
        k++;
    }

    for ( ; i < L + K / 2; i += decimation )
    {
        s = 0;

        for ( j = i - L + 1; j <= K; j++ )
        {
            s += h[j] * x[i - j];
        }

        y[k] = s;
        k++;
    }

    return;
}


/*--------------------------------------------------------------------------
 *  fec_ecu_pitch()
 *
 *  Pitch/correlation analysis and adaptive analysis frame length calculation
 *--------------------------------------------------------------------------*/

static void fec_ecu_pitch(
    const float *prevsynth, /* i  : previous synthesis            */
    float *prevsynth_LP,    /* o  : down-sampled synthesis        */
    const int16_t L,        /* i  : Output frame length           */
    int16_t *N,             /* o  : Analysis frame length in 8kHz */
    float *min_corr,        /* o  : correlation                   */
    int16_t *decimatefator, /* o  : Decimation factor             */
    const int16_t HqVoicing /* i  : Hq voicing flag               */
)
{
    int16_t i, filt_size;
    float accA, accB, accC, Ryy;
    int16_t delay_ind, k, cb_start, cb_end, tmp_short, Lon20;
    float Asr_LP[FEC_DCIM_FILT_SIZE_MAX + 1];

    switch ( L )
    {
        case L_FRAME48k:
            *decimatefator = 6;
            filt_size = 60;
            mvr2r( Asr_LP48, Asr_LP, filt_size + 1 );
            break;
        case L_FRAME32k:
            *decimatefator = 4;
            filt_size = 40;
            mvr2r( Asr_LP32, Asr_LP, filt_size + 1 );
            break;
        case L_FRAME16k:
            *decimatefator = 2;
            filt_size = 20;
            mvr2r( Asr_LP16, Asr_LP, filt_size + 1 );
            break;
        default:
            *decimatefator = 2;
            filt_size = 40;
            mvr2r( Asr_LP16, Asr_LP, filt_size + 1 );
            break;
    }

    /* We need to inverse the ALDO window */

    /* Resampling to work at 8Khz */
    fir_dwn( prevsynth, Asr_LP, prevsynth_LP, 2 * L, filt_size, *decimatefator ); /* resampling without delay */

    Lon20 = ( int16_t )( ( L / 20 ) / *decimatefator );

    /* Correlation analysis */
    *min_corr = 0;
    accC = 0;
    for ( k = 0; k < 6 * Lon20; k++ )
    {
        accC += prevsynth_LP[34 * Lon20 + k] * prevsynth_LP[34 * Lon20 + k];
    }

    if ( HqVoicing == 1 )
    {
        cb_start = 0;
        cb_end = 33 * Lon20;
    }
    else
    {
        cb_start = 0;
        cb_end = 28 * Lon20;
    }

    tmp_short = 34 * Lon20;
    accB = 0;
    delay_ind = cb_start;
    for ( i = cb_start; i < cb_end; i++ ) /* cb_end = 35 let 6 ms min of loop size */
    {
        accA = 0;
        if ( i == cb_start )
        {
            accB = 0;
            for ( k = 0; k < 6 * Lon20; k++ )
            {
                accA += prevsynth_LP[i + k] * prevsynth_LP[tmp_short + k];
                accB += prevsynth_LP[i + k] * prevsynth_LP[i + k];
            }
        }
        else
        {
            accB = accB - prevsynth_LP[i - 1] * prevsynth_LP[i - 1] + prevsynth_LP[i + 6 * Lon20 - 1] * prevsynth_LP[i + 6 * Lon20 - 1];
            for ( k = 0; k < 6 * Lon20; k++ )
            {
                accA += prevsynth_LP[i + k] * prevsynth_LP[tmp_short + k];
            }
        }

        /* tests to avoid division by zero */
        if ( accB <= 0 )
        {
            accB = 1.0f;
        }

        if ( accC <= 0 )
        {
            accC = 1.0f;
        }

        if ( accB * accC <= 0 )
        {
            Ryy = accA;
        }
        else
        {
            Ryy = accA / (float) sqrt( ( accB * accC ) );
        }

        if ( Ryy > *min_corr )
        {
            *min_corr = Ryy;
            delay_ind = i;
        }

        if ( HqVoicing == 0 && *min_corr > 0.95f )
        {
            break;
        }
    }

    *N = 40 * Lon20 - delay_ind - 6 * Lon20;

    return;
}


/*--------------------------------------------------------------------------
 *  fec_ecu_dft()
 *
 *  DFT analysis on adaptive frame length. Analysis frame stretched to
 *  next power of 2 using linear interpolation.
 *--------------------------------------------------------------------------*/

static void fec_ecu_dft(
    const float *prevsynth_LP, /* i  : Downsampled past synthesis (2*160 samples) */
    const int16_t N,           /* i  : Analysis frame length in 8 kHz (corr. max) */
    float *Tfr,                /* o  : DFT coefficients, real part                */
    float *Tfi,                /* o  : DFT coefficients, imag part                */
    float *sum_Tf_abs,         /* o  : Sum of magnitude spectrum                  */
    float *Tf_abs,             /* o  : Magnitude spectrum                         */
    int16_t *Nfft,             /* o  : DFT analysis length 2^(nextpow2(N))        */
    const int16_t element_mode /* i  : IVAS element mode */
)
{
    float target[2 * L_FRAME8k], tmp;
    int16_t i, Lon20, tmp_short, N_LP, k;
    int16_t alignment_point;

    Lon20 = (int16_t) 160 / 20;
    if ( element_mode == EVS_MONO )
    {
        alignment_point = 2 * 160 - 3 * Lon20;
    }
    else
    {
        alignment_point = 2 * 160;
    }


    for ( i = 0; i < N; i++ )
    {
        target[i] = prevsynth_LP[alignment_point - N + i];
    }

    /* DFT  */
    *sum_Tf_abs = 0;
    *Nfft = (int16_t) pow( 2, (int16_t) ceil( log( N ) / log( 2 ) ) );
    tmp = ( (float) N - 1.0f ) / ( (float) *Nfft - 1.0f );

    set_f( Tfr, 0.0f, *Nfft );
    set_f( Tfi, 0.0f, *Nfft );
    Tfr[0] = target[0];
    Tfr[*Nfft - 1] = target[N - 1];
    for ( i = 1; i < *Nfft - 1; i++ ) /* interpolation for FFT */
    {
        tmp_short = (int16_t) floorf( i * tmp );
        Tfr[i] = target[tmp_short] + ( (float) i * tmp - ( (float) tmp_short ) ) * ( target[tmp_short + 1] - target[tmp_short] );
    }

    DoRTFTn( Tfr, Tfi, *Nfft );
    N_LP = (int16_t) ceil( *Nfft / 2 );

    for ( k = 0; k < N_LP; k++ )
    {
        Tf_abs[k] = (float) sqrt( Tfr[k] * Tfr[k] + Tfi[k] * Tfi[k] );
        *sum_Tf_abs += Tf_abs[k];
    }

    return;
}

/*--------------------------------------------------------------------------*
 * singenerator()
 *
 * fast cosinus generator Amp*cos(2*pi*freq+phi)
 *--------------------------------------------------------------------------*/


static void singenerator(
    const int16_t L,     /* i  : size of output */
    const float cosfreq, /* i  : cosine of 1-sample dephasing at the given frequency */
    const float sinfreq, /* i  : sine   of 1-sample dephasing at the given frequency */
    const float a_re,    /* i  : real part of complex spectral coefficient at the given frequency */
    const float a_im,    /* i  : imag part of complex spectral coefficient at the given frequency */
    float xx[]           /* o  : output vector */
)
{

    float *ptr, L_C0, L_S0, L_C1, L_S1;
    float C0, S0, C1, S1;
    int16_t i;

    L_C0 = a_re;
    S0 = a_im;

    ptr = xx;

    *ptr = *ptr + L_C0;
    ptr++;

    for ( i = 0; i < L / 2 - 1; i++ )
    {
        C0 = L_C0;
        L_C1 = C0 * cosfreq;
        L_C1 = L_C1 - S0 * sinfreq;
        L_S1 = C0 * sinfreq;
        S1 = L_S1 + S0 * cosfreq;
        *ptr = *ptr + L_C1;
        ptr++;

        C1 = L_C1;
        L_C0 = C1 * cosfreq;
        L_C0 = L_C0 - S1 * sinfreq;
        L_S0 = C1 * sinfreq;
        S0 = L_S0 + S1 * cosfreq;
        *ptr = *ptr + L_C0;
        ptr++;
    }

    C0 = L_C0;
    L_C1 = C0 * cosfreq;
    L_C1 = L_C1 - S0 * sinfreq;
    *ptr = *ptr + L_C1;
    ptr++;

    return;
}


/*--------------------------------------------------------------------------
 *  sinusoidal_synthesis()
 *
 *  ECU frame sinusoid generation
 *--------------------------------------------------------------------------*/

static void sinusoidal_synthesis(
    const float *Tfr,              /* i  : DFT coefficients, real part                */
    const float *Tfi,              /* i  : DFT coefficients, imag part                */
    float *Tf_abs,                 /* i  : Magnitude spectrum                         */
    const int16_t N,               /* i  : Analysis frame length in 8 kHz (corr. max) */
    const int16_t L,               /* i  : Output frame length                        */
    const int16_t decimate_factor, /* i  : Downsampling factor                        */
    const int16_t Nfft,            /* i  : DFT analysis length 2^(nextpow2(N))        */
    const float sum_Tf_abs,        /* i  : Sum of magnitude spectrum                  */
    float *synthesis,              /* o  : ECU sinusoidal synthesis                   */
    const int16_t HqVoicing        /* i  : HQ voicing flag                            */
)
{
    int16_t i, k, nb_pulses, indmax = 0, nb_pulses_final;
    int16_t pulses[FEC_MAX / 2];
    float a_re[FEC_NB_PULSE_MAX], a_im[FEC_NB_PULSE_MAX], cosfreq, sinfreq;
    float freq[FEC_NB_PULSE_MAX], tmp;
    float mmax, cumsum;
    int16_t Lon20 = 8;

    /* peak selection  */
    int16_t PL, cpt;
    float old, new_s;
    int16_t *p_pulses;
    int16_t glued;

    p_pulses = pulses;
    nb_pulses = 0;
    new_s = Tf_abs[1];
    glued = 1;
    cpt = 0;
    old = 0;

    PL = 0;
    if ( N > Lon20 * 10 || HqVoicing )
    {
        PL = 1;
    }
    while ( cpt <= N / 2 - 1 - 2 )
    {
        if ( Tf_abs[cpt] > old && Tf_abs[cpt] > new_s )
        {
            glued = cpt;

            for ( i = glued; i < cpt + PL + 1; i++ )
            {
                *p_pulses++ = i;
                nb_pulses++;
            }
            old = Tf_abs[cpt + PL];
            new_s = Tf_abs[cpt + 2 + PL];
            cpt = cpt + PL + 1;
            glued = 1;
        }
        else
        {
            old = Tf_abs[cpt];
            new_s = Tf_abs[cpt + 2];
            cpt++;
            glued = 0;
        }
    }


    nb_pulses_final = 0;

    /* peak selection : keep the more energetics (max 20) */
    tmp = 1.0f / (float) ( Nfft / 2 );
    cumsum = 0;
    for ( i = 0; i < min( FEC_NB_PULSE_MAX, nb_pulses ); i++ )
    {
        mmax = 0;
        for ( k = 0; k < nb_pulses; k++ )
        {
            if ( Tf_abs[pulses[k]] > mmax )
            {
                mmax = Tf_abs[pulses[k]];
                indmax = pulses[k];
            }
        }
        cumsum += Tf_abs[indmax];

        if ( HqVoicing || cumsum < sum_Tf_abs * 0.7f )
        {
            a_re[i] = Tfr[indmax] * tmp;
            a_im[i] = Tfi[indmax] * tmp;
            freq[i] = (float) indmax * 2 / ( (float) N );
            nb_pulses_final++;
            Tf_abs[indmax] = -1;
        }
        else
        {
            a_re[i] = Tfr[indmax] * tmp;
            a_im[i] = Tfi[indmax] * tmp;
            freq[i] = (float) indmax * 2 / ( (float) N );
            nb_pulses_final++;
            break;
        }
    }

    nb_pulses = nb_pulses_final;

    /* sinusoidal synthesis */
    set_f( synthesis, 0.0f, 40 * L / 20 );

    if ( HqVoicing )
    {
        k = 40 * L / 20;
    }
    else
    {
        k = 40 * L / 20;
    }

    for ( i = 0; i < nb_pulses; i++ )
    {
        cosfreq = (float) cos( EVS_PI * freq[i] / (float) decimate_factor );
        sinfreq = (float) sin( EVS_PI * freq[i] / (float) decimate_factor );
        singenerator( k, cosfreq, sinfreq, a_re[i], a_im[i], synthesis );
    }

    return;
}

/*--------------------------------------------------------------------------
 *  fec_noise_filling()
 *
 *  Adds noise component to ECU frame by
 *  - subtracting the sinusoidal synthesis from the analysis frame
 *  - copying noise segments at random indices and adding them together
 *    with an overlap-add operation
 *  Copying the beginning of the frame from the past synthesis and aligning
 *  it to be inserted into wtda
 *--------------------------------------------------------------------------*/

static void fec_noise_filling(
    const float *prevsynth,     /* i  : Past synthesis buffer (length 2*L)         */
    float *synthesis,           /* i/o: Sinusoidal ECU / Sinusoidal ECU + noise    */
    const int16_t L,            /* i  : Output frame length                        */
    const int16_t N,            /* i  : Analysis frame length in 8 kHz (corr. max) */
    const int16_t HqVoicing,    /* i  : HQ voicing flag                            */
    float *gapsynth,            /* o  : Buffer for crossfade to good frame         */
    int16_t *ni_seed_forfec,    /* i/o: Random seed for picking noise segments     */
    const int16_t element_mode, /* i  : IVAS element mode                          */
    const float *old_out )
{
    float SS[L_FRAME48k / 2], tmp;
    int16_t Rnd_N_noise;
    int16_t k, kk, i;
    int16_t N_noise;
    float noisevect[34 * L_FRAME48k / 20];
    const float *p_mdct_ola;
    int16_t alignment_point;

    if ( element_mode == EVS_MONO )
    {
        alignment_point = 2 * L - 3 * L / 20;
    }
    else
    {
        alignment_point = 2 * L;
    }
    mvr2r( prevsynth + alignment_point - N, noisevect, N );


    /* Noise addition on full band  */
    /* residual  */
    if ( N < L )
    {
        N_noise = ( N / 2 );
        for ( k = 0; k < N; k++ )
        {
            noisevect[k] = ( noisevect[k] - synthesis[k] );
        }
    }
    else
    {
        N_noise = L / 2;