/*
**
** PsychoC.cpp
**
** $VER: PsychoC.cpp 1.0 (01.06.98)
**
** Analyseur psycho-acoustique
**
** $Revision: 1.4 $
** $State: Exp $
** $Date: 1998/08/16 19:03:41 $
**
** $Log: PsychoC.cpp $
** Revision 1.4  1998/08/16 19:03:41  kakace
** Version Beta3+
**
** Revision 1.3  1998/07/27 14:13:05  kakace
** Validation de toutes les fonctions
**
** Revision 1.2  1998/07/27 02:15:41  kakace
** Validation des fonctions (Layer II)
**
*/


//----------------------------------------------------------------------------------------------------

/// Includes
#ifndef  _INCLUDE_ASSERT_H
#include <assert.h>
#endif

#ifndef  _INCLUDE_IOSTREAM_H
#include <iostream.h>
#endif

#ifndef  _INCLUDE_IOMANIP_H
#include <iomanip.h>
#endif

#ifndef  _INCLUDE_STDLIB_H
#include <stdlib.h>
#endif

#ifndef  _INCLUDE_MATH_H
#include <math.h>
#endif

#ifndef  _MEMORY_HPP
#include "Memory.hpp"
#endif

#ifndef  _CODER_CLASS_HPP
#include "CoderC.hpp"
#endif

#ifndef  _PSYCHO_CLASS_HPP
#include "PsychoC.hpp"
#endif

#ifndef  _PEGASECOND_HPP
#include "PegaseCond.hpp"
#endif
///

#define MINTHRESHOLD  60802371420160.0;


//----------------------------------------------------------------------------------------------------
//========================================== Classe PsychoC ==========================================
//----------------------------------------------------------------------------------------------------

/// PsychoC::PsychoC()
PsychoC::PsychoC(CoderC &roCoderC) : psy_roCoderC(roCoderC)
{
    psy_MemoryError = FFT_MemoryError;

    // Allouer tous les tampons nécessaires.

    psy_MemoryError |= !(psy_pdSpreading = (T_FCBCB *) AllocAlignedBuffer(sizeof(T_FCBCB)));
    psy_MemoryError |= !(psy_pdTmpSave = (T_F2HBLK32 *) AllocAlignedBuffer(sizeof(T_F2HBLK32)));
}
///
/// PsychoC::~PsychoC() (virtual)
PsychoC::~PsychoC()
{
    if (psy_pdSpreading)        FreeAlignedBuffer(psy_pdSpreading);
    if (psy_pdTmpSave)          FreeAlignedBuffer(psy_pdTmpSave);
}
///

/// PsychoC::CalcPermissibleNoise()
/****** Class PsychoC/CalcPermissibleNoise **********************************
*
*   NAME
*   PsychoC::CalcPermissibleNoise --
*
*   SYNOPSIS
*   PsychoC::CalcPermissibleNoise()
*
*   void PsychoC::CalcPermissibleNoise()
*
*   FUNCTION
*
*   EXAMPLE
*
*   NOTES
*   Fonction fortement modifiée pour optimiser la vitesse de traitement.
*   A contrôler.
*
*   BUGS
*   Fonction contrôlée conforme le 26.07.98
*
*****************************************************************************
*
*   Benchmarks (25 cycles par seconde) :
*       14.06.98 : 1'05.63
*
*   Estimation de temps de calcul (pour 31768 frames) :
*   · Log + Exp : 22 secondes
*   · Calcul initial "cb" et "ecb" : 30 secondes.
*/


void PsychoC::CalcPermissibleNoise()
{
    static float bmax[27] = {20.0, 20.0, 20.0, 20.0, 20.0, 17.0, 15.0,
                             10.0,  7.0,  4.4,  4.5,  4.5,  4.5,  4.5,
                              4.5,  4.5,  4.5,  4.5,  4.5,  4.5,  4.5,
                              4.5,  4.5,  4.5,  3.5,  3.5,  3.5
                            };

    // bmax_exp[n] = exp(-bmax[n] * LN_TO_LOG10)

    static REAL bmax_exp[27] = {0.01,             0.01,             0.01,
                                0.01,             0.01,             0.01995262314969,
                                0.03162277660168, 0.1,              0.19952623149689,
                                0.36307805477010, 0.35481338923358, 0.35481338923358,
                                0.35481338923358, 0.35481338923358, 0.35481338923358,
                                0.35481338923358, 0.35481338923358, 0.35481338923358,
                                0.35481338923358, 0.35481338923358, 0.35481338923358,
                                0.35481338923358, 0.35481338923358, 0.35481338923358,
                                0.44668359215096, 0.44668359215096, 0.44668359215096
                               };

//----------------------------------------------------------------------------------------------------

#if 1

    register int j, k, l;
    register REAL cb, ecb, temp;
    register REAL *pSpreading;

    register REAL  *pdTMN = psy_adToneMaskingNoise;
    register REAL  *pdNB  = psy_adNB;
    register BYTE  *pbCB  = &psy_abCBands[0][0];

    // NB : CBANDS a été remplacé par "psy_wCBands" pour tenir compte du fait qu'il n'y a pas
    //      63 bandes lorsqu'on encode en 32 KHz.

    for (j = 0; j <= psy_wCBands; j++)
    {
        REAL *pGrp;
        cb = ecb = 0.0;

        // La valeur "Spreading" est non nulle pour des indices bien particuliers. C'est pourquoi
        // le tableau "psy_abCBands" fournit l'indice de début et l'indice de fin des valeurs non
        // nulles, ce qui élimine un test inutile.

        l = pbCB[64];
        k = *pbCB++;
        pGrp       = &psy_adGroupedValues[k][0];
        pSpreading = &(*psy_pdSpreading)[j][k];

        for (; k <= l; k++)
        {
            temp = *pSpreading++;
            ecb += temp * *pGrp++;      // = grouped_e[]
            cb  += temp * *pGrp++;      // = grouped_c[]
        }

        if (ecb)
        {
            cb = cb / ecb;
            if      (cb < 0.05) cb = LOG10_005 + LOG10_2;   // = log10(0.05)
            else if (cb > 0.5)  cb = LOG10_05  + LOG10_2;   // = log10(0.5)
            else                cb = log10(cb) + LOG10_2;
        }
        else                    cb = LOG10_005 + LOG10_2;   // = log10(0.05)

        // Partie fortement modifiée par rapport à la version originale
        // (Les tables ont été modifiées pour prendre en compte encore plus
        //  de données précalculées).

        cb = 5.5 - cb * *pdTMN++;       // psy_adToneMaskingNoise[j];
        k  = psy_abCenterBarkVal[j];

        if (cb < bmax[k]) cb = bmax_exp[k];
        else              cb = exp(cb * (-LN_TO_LOG10));

        temp = psy_adRNorm[j];

        if (temp)   *pdNB++ = ecb * cb * temp;
        else        *pdNB++ = temp;                     // = 0.0
    }

#else   // Original code

    int j, k, l;
    REAL cb, ecb, temp;
    REAL *pSpreading;

    // NB : CBANDS a été remplacé par "psy_wCBands" pour tenir compte du fait qu'il n'y a pas
    //      63 bandes lorsqu'on encode en 32 KHz.

    for (j = 0; j <= psy_wCBands; j++)
    {
        REAL *pGrp;
        cb = ecb = 0.0;

        // La valeur "Spreading" est non nulle pour des indices bien particuliers. C'est pourquoi
        // le tableau "psy_abCBands" fournit l'indice de début et l'indice de fin des valeurs non
        // nulles, ce qui élimine un test inutile.

        k = psy_abCBands[0][j];
        l = psy_abCBands[1][j];
        pGrp       = &psy_adGroupedValues[k][0];
        pSpreading = &(*psy_pdSpreading)[j][k];

        for (; k <= l; k++)
        {
            temp = *pSpreading++;
            ecb += temp * *pGrp++;      // = grouped_e[]
            cb  += temp * *pGrp++;      // = grouped_c[]
        }

        if (ecb != 0.0)
        {
            cb = cb / ecb;
            if      (cb < 0.05) cb = LOG10_005 + LOG10_2;   // = log10(0.05)
            else if (cb > 0.5)  cb = LOG10_05 + LOG10_2;    // = log10(0.5)
            else                cb = log10(cb) + LOG10_2;
        }
        else                    cb = LOG10_005 + LOG10_2;   // = log10(0.05)

        // Partie fortement modifiée par rapport à la version originale
        // (Les tables ont été modifiées pour prendre en compte encore plus
        //  de données précalculées).

        cb = 5.5 - cb * psy_adToneMaskingNoise[j];
        k  = psy_abCenterBarkVal[j];

        if (cb < bmax[k]) cb = bmax_exp[k];
        else              cb = exp(cb * (-LN_TO_LOG10));

        temp = psy_adRNorm[j];

        if (temp)   psy_adNB[j] = ecb * cb * temp;
        else        psy_adNB[j] = temp;                     // = 0.0
    }
#endif
}
///
/// PsychoC::CalcGroupedValues()
/****** Class PsychoC/CalcGroupedValues *************************************
*
*   NAME
*   PsychoC::CalcGroupedValues --
*
*   SYNOPSIS
*   PsychoC::CalcGroupedValues(wChannel, wNew)
*
*   void PsychoC::CalcGroupedValues(UWORD, UWORD);
*
*   FUNCTION
*   Calcule des vecteurs complexes à partir des données fournies par la FFT,
*   et les cumule par bandes de fréquences.
*
*   INPUTS
*   wChannel - Numéro de canal.
*   wNew     - Indice courant
*
*   NOTES
*
*   BUGS
*   Cette fonction n'est pas strictement conforme au source originale : L'
*   énergie du complexe C[512] est ici bornée (borne inférieure = 0.0005).
*
*   Fonction validée le 26/07/98
*
*   SEE ALSO
*
*****************************************************************************
*
*   Benchmarks (3 cycles par seconde) :
*       21.06.98 : 2'47.80
*
*   NB : La clarté du code souffre un peu de l'entrelacement des calculs
*        destinés à occuper le FPU pendant que le cache termine (en principe)
*        la lecture d'une ligne.
*/

#if _USE_ASM_FCT_ == 1

extern "C"{
void CalcGroupedValues(UWORD wNew, REAL *GroupedValues, T_FFTDATAS *FFTOut, REAL *Energy2, T_FFTDATAS *tmp, UBYTE *Part);
}

inline void PsychoC::CalcGroupedValues(UWORD wChannel, UWORD wNew)
{
    ::CalcGroupedValues(wNew, (REAL *) psy_adGroupedValues, &(*FFT_pdWinSmplShuffled)[wChannel][0][0],
                        psy_adEnergy2, &(*psy_pdTmpSave)[wChannel][0][0][0], psy_abPartition);
}

#else

void PsychoC::CalcGroupedValues(UWORD wChannel, UWORD wNew)
{
#define pOldest  pNew                       // Synonyme.

    const ULONG COSINE = 0;                 // Index d'accès aux données.
    const ULONG SINE   = 2;
    const ULONG ENERGY = 4;
    const ULONG NEXT_OFFSET = 6;            // {COSINE, SINE, ENERGY} × {NEW, OLD}

    UWORD wOld = 1 - wNew;

    int j, part, idx;

    double cos_phi, cos_2phi;
    double sin_phi, sin_2phi;
    double energy, r, r_prime;

    T_FFTDATAS *p;
    T_FFTDATAS *pNew, *pOld;
    REAL *Grp;

    // Initialiser les tableaux destinés à accueillir les sommes intermédiaires.

    {
        REAL *p = &psy_adGroupedValues[0][0];

        for (j = 0; j <= psy_wCBands; j++)
        {
            *p++ = 0.0;                 // Grouped_e[j]
            *p++ = 0.0;                 // Grouped_c[j]
        }
    }

    p       = &(*FFT_pdWinSmplShuffled)[wChannel][0][0];
    pOld    = &(*psy_pdTmpSave)[wChannel][0][0][wOld];
    pNew    = &(*psy_pdTmpSave)[wChannel][0][0][wNew];

    for (j = 0; j < HBLKSIZE; j++)
    {
        // Initialiser les variables intermédiaires.
        // Les tables pOld et pNew se présentent sous la forme :
        //      Cos_new, Cos_old, Sin_new, Sin_old, Nrj_new, Nrj_old, Cos_new...
        // NB : Cette partie s'efforce d'effectuer un maximum de calculs entre chaque
        //      lecture dans une table de façon à ne pas être bloqué par un cache en
        //      train de lire une ligne.

        cos_phi  = pOld[COSINE];                        // = cos(a)
        r        = cos_phi + cos_phi;                   // = 2.cos(a)
        r_prime  = r * cos_phi;                         // = 2.cos²(a)
        cos_phi  = pOldest[COSINE];                     // = cos(b)

        sin_2phi = pOld[SINE] * r;                      // = sin(2a) = 2.cos(a).sin(a)
        r        = r_prime - 1.0;                       // = cos(2a) = 2.cos²(a) - 1
        energy   = cos_phi * r;                         // Commencer le calcul de cos(2a-b) ( = cos(2a).cos(b))
        sin_phi  = pOldest[SINE];                       // = sin(b)

        r_prime   = pOld[ENERGY];
        energy  += sin_phi * sin_2phi;                  // cos(2a-b) = cos(2a).cos(b) + sin(2a).sin(b)
        sin_2phi = sin_2phi * cos_phi - r * sin_phi;    // sin(2a-b) = sin(2a).cos(b) - cos(2a).sin(b)
        r_prime += r_prime - pOldest[ENERGY];

        cos_2phi = energy;                              // HINT : "energy" est dans un registre, ce qui
                                                        // n'est pas le cas de "cos_2phi"

        // Récupérer les données en sortie de FFT (elle fournie gracieusement les
        // sin() et cos() voulus : inutile de passer par un angle intermédiaire).

        cos_phi = *p++;                                 // (*FFT_pdWinSmplShuffled)[wChannel][j][0];
        energy  = cos_phi * cos_phi;
        r       = 0.0005;
        sin_phi = *p++;                                 // (*FFT_pdWinSmplShuffled)[wChannel][j][1];
        energy += sin_phi * sin_phi;

        // Donner une valeur minimum à l'énergie.
        // NOTE : Dans le source d'origine, l'énergie n'est pas bornée (borne inférieure)
        //        pour le complexe C[512].

        if (energy <= r)
        {
            energy  = r;                                // energy = 0.0005
            r       = SQRT_00005;                       // = sqrt(0.0005) = sqrt(5) / 100
            cos_phi = r;                                // => phi = 0
            sin_phi = 0.0;
        }
        else
        {
            r = sqrt(energy);
        }

        // Préparer la somme des énergies pour TranslateThreshold().
        // Par la même occasion, on occupe l'IU pendant que le FPU calcule SQRT().

        part = psy_abPartition[j];
        Grp  = &psy_adGroupedValues[part][0];

        idx  = j >> 4;
        part = j & 15;

        if (part)           psy_adEnergy2[idx] += energy;   // part != 0 : Ajouter l'énergie à la bande courante.
        else
        {
            if (idx != 32)  psy_adEnergy2[idx]  = energy;   // part == 0 && idx != 32 : Initialiser la somme.
            if (idx--)      psy_adEnergy2[idx] += energy;   // part == 0 && idx != 0  : Ajouter à la bande précédente.
        }

        // Mémoriser les valeurs courantes pour la prochaine passe.

        pNew[ENERGY] = r;                               // r != 0 puisque l'énergie est bornée.

        r = 1.0 / r;
        pNew[COSINE] = cos_phi * r;
        pNew[SINE]   = sin_phi * r;

        r = pNew[ENERGY] + fabs(r_prime);               // r > 0 puisque l'energie est bornée.

        cos_phi -= r_prime * cos_2phi;
        sin_phi -= r_prime * sin_2phi;
        cos_phi  = cos_phi * cos_phi + sin_phi * sin_phi;
              // = energy + r'² - 2.r.r'.cos(phi + phi')

        Grp[0] += energy;
        Grp[1] += energy * sqrt(cos_phi) / r;

        pNew += NEXT_OFFSET;                            // Ajuster les pointeurs.
        pOld += NEXT_OFFSET;                            // (pOldest = pNew)
    }
}

#endif  // _USE_ASM_FCT_
///
/// PsychoC::InitPsychoAnalyzer()
/****** Class PsychoC/InitPsychoAnalyzer ************************************
*
*   NAME
*   PsychoC::InitPsychoAnalyzer -- Initialiser l'analyseur psycho-acoustique.
*
*   SYNOPSIS
*   result = PsychoC::InitPsychoAnalyzer(iSBLimit)
*
*   bool PsychoC::InitPsychoAnalyzer(int);
*
*   FUNCTION
*   Initialise l'analyseur psycho-acoustique en préparant les différentes
*   tables nécessaires à son fonctionnement.
*
*   INPUTS
*   iSBLimit - Nombre de sous-bandes.
*
*   RESULT
*   result = FALSE en cas d'erreur (tampons non alloués par exemple).
*
*   NOTES
*   Une table SinCos() devrait aussi être initialisée afin d'améliorer la
*   vitesse des FFT lorsque le CPU utilisé est plus handicapé par une capa-
*   cité de traitement limitée que par ses accès mémoire (68020/68030).
*
*   BUGS
*   Il n'existe pas de table AbsoluteThreshold pour le MPEG-2 !!!
*   Fonction contrôlée conforme le 27/07/98
*
*   SEE ALSO
*
*****************************************************************************
*
*/


int PsychoC::InitPsychoAnalyzer(int iSBLimit)
{
    extern const REAL a_dAbsoluteThresholds[3][HBLKSIZE];

    REAL adCenterBarkValue[CBANDS];         // (cbval[])

    int Result = TRUE;
    int i, j;

    // Initialiser les paramètres de l'analyseur psycho-acoustique selon le type de Layer.

    psy_wNew     = 0;
    psy_wSBLimit = iSBLimit;

    // Adresser les tables des Thresholds en fonction de la fréquence d'échantillonage

    switch (psy_roCoderC.GetSampleFreq())
    {
        case 32000: psy_pdAbsThreshold = &a_dAbsoluteThresholds[0][0]; break;
        case 44100: psy_pdAbsThreshold = &a_dAbsoluteThresholds[1][0]; break;
        case 48000: psy_pdAbsThreshold = &a_dAbsoluteThresholds[2][0]; break;
        default:
            Result = FALSE;
            break;
    }

    // Initialiser les données pour le FFT.

    {
        static int crit_band[27] = {    0,   100,   200,   300,   400,   510,   630,   770,
                                      920,  1080,  1270,  1480,  1720,  2000,  2320,  2700,
                                     3150,  3700,  4400,  5300,  6400,  7700,  9500, 12000,
                                    15500, 25000, 30000
                                   };

        const REAL INV_BLKSIZE = 1.0 / BLKSIZE;
        REAL  t, cb, cb_pred;
        REAL  bval_lo;
        int   part, idx;

        REAL  freq_mult = psy_roCoderC.GetSampleFreq() * INV_BLKSIZE;

        idx = 1;
        t = freq_mult;
        cb = crit_band[1];
        cb_pred = 0.0;

        bval_lo = 0.0;                              // Fréquence de départ de la bande courante.
        adCenterBarkValue[0] = 0.0;
        part = 0;                                   // Numéro d'index de la bande courante.
        j = 1;                                      // Nombre de fréquences dans la bande courante.

        for (i = 1; i < HBLKSIZE; i++, t += freq_mult)
        {
            REAL bval;

            if (t > cb)
            {
                idx++;
                cb_pred = cb;
                cb = crit_band[idx];
            }

            bval = idx - 1 + (t - cb_pred) / (cb - cb_pred);

            if ( (bval - bval_lo) > CB_FRACTION)
            {
                adCenterBarkValue[part] /= j;       // Ajuster les infos de la bande précédente.
                psy_abNumLines[part] = j;

                psy_abPartition[i] = ++part;        // Initialiser la bande suivante.
                adCenterBarkValue[part] = bval_lo = bval;
                j = 1;
            }
            else
            {
                psy_abPartition[i] = part;          // Ajouter la fréquence dans la bande courante.
                adCenterBarkValue[part] += bval;
                j++;
            }
        }

        // Ajuster les informations de la dernière bande.

        psy_wCBands = part = psy_abPartition[HBLKSIZE-1];
        psy_abNumLines[part] = j;
        adCenterBarkValue[part] /= j;
    }

    // Calculer la fonction de "spreading"

    {
        const REAL cst = 17.5 * 1.106650803 - 3.555;      // = 17.5 * sqrt(1.224676) - 3.555 = 15.811389
        REAL t1, t2, t3;

        // NB : CBANDS a été remplacé par "psy_wCBands" pour tenir compte du fait qu'il n'y a pas
        //      63 bandes lorsqu'on encode en 32 KHz.

        for (j = 0; j <= psy_wCBands; j++)
        {
            psy_abCBands[0][j] = psy_abCBands[1][j] = -1;

            for (i = 0; i <= psy_wCBands; i++)
            {
                t1 = (adCenterBarkValue[j] - adCenterBarkValue[i]) * 1.05;

                if (t1 >= 0.5 && t1 <= 2.5)
                {
                    t2 = t1 - 0.5;
                    t2 = 8.0 * t2 * (t2 - 2.0);     // = 8.(t2.t2 - 2.t2)
                }
                else t2 = 0.0;

                t1 += 0.474;
                t3 = cst + 7.5 * t1 - 17.5 * sqrt(1.0 + t1 * t1);

                if (t3 <= -100.0) (*psy_pdSpreading)[j][i] = 0.0;
                else
                {
                    if (psy_abCBands[0][j] < 0) psy_abCBands[0][j] = i;
                    psy_abCBands[1][j] = i;

                    t3 = (t2 + t3) * LN_TO_LOG10;
                    (*psy_pdSpreading)[j][i] = exp(t3);
                }
            }
        }
    }

    // Calculer les valeurs "Tone Masking Noise"

    {
        REAL t1;

        // NB : CBANDS a été remplacé par "psy_wCBands" pour tenir compte du fait qu'il n'y a pas
        //      63 bandes lorsqu'on encode en 32 KHz.

        for (j = 0; j <= psy_wCBands; j++)
        {
            t1 = adCenterBarkValue[j];
            psy_abCenterBarkVal[j] = (UBYTE) ceil(t1);

            // Valeur précalculée pour CalcPersmissibleNoise()
            psy_adToneMaskingNoise[j] = (t1 > 9.0) ? (t1 + 10.0) : 19.0;

            t1 = (*psy_pdSpreading)[j][0];

            for (i = 1; i <= psy_wCBands; i++)
            {
                t1 += (*psy_pdSpreading)[j][i];
            }

            // Précalculer la valeur utilisée par CalcPermissibleNoise().

            t1 *= psy_abNumLines[j];
            if (t1 != 0.0) t1 = 1.0 / t1;

            psy_adRNorm[j] = t1;
        }
    }

    // Initialiser la table TmpSave[]

    for (j = 0; j < HBLKSIZE; j++)
    {
        (*psy_pdTmpSave)[0][j][0][0] = (*psy_pdTmpSave)[0][j][0][1] = 1;
        (*psy_pdTmpSave)[1][j][0][0] = (*psy_pdTmpSave)[1][j][0][1] = 1;
    }

    return Result;
}
///

//----------------------------------------------------------------------------------------------------
//======================================== Classe Psycho_L1C =========================================
//----------------------------------------------------------------------------------------------------

/// Psycho_L1C::TranslateThreshold()  (virtual)
/****** Class Psycho_L1C/TranslateThreshold *********************************
*
*   NAME
*   Psycho_L1C::TranslateThreshold --
*
*   SYNOPSIS
*   Psycho_L1C::TranslateThreshold(wChannel)
*
*   void Psycho_L1C::TranslateThreshold(UWORD);
*
*   FUNCTION
*
*   INPUTS
*   wChannel - Numéro de canal.
*
*   NOTES
*
*   BUGS
*   Fonction contrôlée conforme le 27/07/98
*
*****************************************************************************
*
*/

/*
inline REAL Psycho_L1C::GetNoiseLevel(REAL dNoiseBand, REAL dThreshold, REAL *p_LThr)
{
    if (dNoiseBand < dThreshold) dNoiseBand = dThreshold;

    dThreshold = *p_LThr;
    *p_LThr = LXMIN * dNoiseBand;       // = 32 * dNoiseBand

    if (dNoiseBand < dThreshold) dThreshold = dNoiseBand;
    dNoiseBand *= 0.00316;

    if (dNoiseBand < dThreshold) dNoiseBand = dThreshold;

    return dNoiseBand;
}
*/

void Psycho_L1C::TranslateThreshold(UWORD wChannel)
{
#if 1

    register int j, k;
    register REAL dNoiseBand, dThreshold, minthres;
    register ULONG   SBLimit = psy_wSBLimit;

    register const REAL *pThreshold  = psy_pdAbsThreshold;
    register REAL       *pLThr       = &(*psy_pdLThr)[wChannel][0];
    register REAL       *pNB         = psy_adNB;
    register UBYTE      *pbPartition = &psy_abPartition[1];
    register REAL       *pdEnergy2   = psy_adEnergy2;
    register REAL       *pdSNR       = psy_adSNRTemp;

    dNoiseBand = *pNB;
    dThreshold = *pThreshold++;

    if (dNoiseBand < dThreshold) dNoiseBand = dThreshold;

    dThreshold = *pLThr;
    *pLThr++ = LXMIN * dNoiseBand;          // = 32 * dNoiseBand

    if (dNoiseBand < dThreshold) dThreshold = dNoiseBand;
    dNoiseBand *= 0.00316;

    if (dNoiseBand < dThreshold) dNoiseBand = dThreshold;

    for (j = 0; j < 13; j++)
    {
        minthres = MINTHRESHOLD;
        if (minthres > dNoiseBand) minthres = dNoiseBand;

        for (k = 16; k > 0; k--)
        {
            dNoiseBand = pNB[ *pbPartition++ ];
            dThreshold = *pThreshold++;

            if (dNoiseBand < dThreshold) dNoiseBand = dThreshold;

            dThreshold = *pLThr;
            *pLThr++ = LXMIN * dNoiseBand;          // = 32 * dNoiseBand

            if (dNoiseBand < dThreshold) dThreshold = dNoiseBand;
            dNoiseBand *= 0.00316;

            if (dNoiseBand < dThreshold) dNoiseBand = dThreshold;

            if (minthres > dNoiseBand) minthres = dNoiseBand;
        }
        *pdSNR++ = *pdEnergy2++ / (minthres * 17.0);
    }

    for (j = 13; j < SBLimit; j++)
    {
        minthres = dNoiseBand;

        for (k = 16; k > 0; k--)
        {
            dNoiseBand = pNB[ *pbPartition++ ];
            dThreshold = *pThreshold++;

            if (dNoiseBand < dThreshold) dNoiseBand = dThreshold;

            dThreshold = *pLThr;
            *pLThr++ = LXMIN * dNoiseBand;          // = 32 * dNoiseBand

            if (dNoiseBand < dThreshold) dThreshold = dNoiseBand;
            dNoiseBand *= 0.00316;

            if (dNoiseBand < dThreshold) dNoiseBand = dThreshold;

            minthres += dNoiseBand;
        }
        *pdSNR++ = *pdEnergy2++ / minthres;
    }

#else   // Original code.

    int j, k, idx;
    REAL dNoiseBand, dThreshold, minthres;

    const REAL *pThreshold = psy_pdAbsThreshold;
    REAL       *pLThr      = &(*psy_pdLThr)[wChannel][0];
    REAL       *pNB        = psy_adNB;

    dNoiseBand = GetNoiseLevel( pNB[0], pThreshold[0], pLThr);

    for (j = 0; j < 13; j++)
    {
        minthres = MINTHRESHOLD;
        if (minthres > dNoiseBand) minthres = dNoiseBand;

        for (k = 1, idx = (j * 16) + 1; k < 17; k++, idx++)
        {
            dNoiseBand = pNB[ psy_abPartition[idx] ];
            dThreshold = pThreshold[idx];

            dNoiseBand = GetNoiseLevel(dNoiseBand, dThreshold, &pLThr[idx]);

            if (minthres > dNoiseBand) minthres = dNoiseBand;
        }
        psy_adSNRTemp[j] = psy_adEnergy2[j] / (minthres * 17.0);
    }

    for (j = 13; j < psy_wSBLimit; j++)
    {
        minthres = dNoiseBand;

        for (k = 1, idx = (j * 16) + 1; k < 17; k++, idx++)
        {
            dNoiseBand = pNB[ psy_abPartition[idx] ];
            dThreshold = pThreshold[idx];

            dNoiseBand = GetNoiseLevel(dNoiseBand, dThreshold, &pLThr[idx]);
            minthres += dNoiseBand;
        }
        psy_adSNRTemp[j] = psy_adEnergy2[j] / minthres;
    }
#endif
}
///
/// Psycho_L1C::PsychoAnal()
/****** Class Psycho_L1C/PsychoAnal *****************************************
*
*   NAME
*   Psycho_L1C::PsychoAnal -- Fonction principale.
*
*   SYNOPSIS
*   Psycho_L1C::PsychoAnal(LTMin[2][32])
*
*   void Psycho_L1C::PsychoAnal(LONG);
*
*   FUNCTION
*   Effectue l'analyse psycho-acoustique.
*
*   NOTES
*
*   BUGS
*
*****************************************************************************
*
*/


void Psycho_L1C::PsychoAnal(LONG LTMin[MAX_CH][SBLIMIT])
{
    const ULONG FLUSH = 384;

    ULONG j, k;
    UWORD wNumChannels = psy_roCoderC.GetNumChannels();
    WORD  *pSamples    = psy_roCoderC.GetSamplePtr();

    // Appliquer le filtre d'analyse (HANN window).

    SetupWinSamples_1024(pSamples, wNumChannels);

    // Calculer la transformée de Fourrier.

    if (wNumChannels == 2)
    {
        pSamples += FLUSH * 2;
        ComputeFastFFT_Stereo(BLKSIZE);
    }
    else
    {
        pSamples += FLUSH;
        ComputeFastFFT_Mono(BLKSIZE);
    }

    // Exploiter les résultats de la transformée de Fourrier pour chaque canal.

    for (k = 0, psy_wNew = 1 - psy_wNew; k < wNumChannels; k++)
    {
        CalcGroupedValues(k, psy_wNew);
        CalcPermissibleNoise();
        TranslateThreshold(k);

        for (j = 0; j < psy_wSBLimit; j++)
        {
        #ifndef USE_OLD_MATHLIB
            LTMin[k][j] = (LONG) ( (LOG10_TO_LN * 1000) * log(psy_adSNRTemp[j]) );
        #else
            LTMin[k][j] = (LONG) floor( (LOG10_TO_LN * 1000) * log(psy_adSNRTemp[j]) );
        #endif
        }
    }
}
///
/// Psycho_L1C::InitPsychoAnalyzer()
//
//  Initialiser la table "psy_pdLThr", uniquement utilisée avec le Layer I

int Psycho_L1C::InitPsychoAnalyzer(int iSBLimit)
{
    // Initialiser la table des Threshold mini

    REAL *p1 = &(*psy_pdLThr)[0][0];
    REAL *p2 = &(*psy_pdLThr)[1][0];
    REAL d   = MINTHRESHOLD;

    for (int i = 0; i < HBLKSIZE; i++)
    {
        *p1++ = *p2++ = d;
    }

    return PsychoC::InitPsychoAnalyzer(iSBLimit);          // Initialisation générale.
}
///

//----------------------------------------------------------------------------------------------------
//======================================== Classe Psycho_L2C =========================================
//----------------------------------------------------------------------------------------------------

/// Psycho_L2C::TranslateThreshold()  (virtual)
/****** Class Psycho_L2C/TranslateThreshold *********************************
*
*   NAME
*   Psycho_L2C::TranslateThreshold --
*
*   SYNOPSIS
*   Psycho_L2C::TranslateThreshold(pSNR)
*
*   void Psycho_L2C::TranslateThreshold(REAL *);
*
*   FUNCTION
*
*   INPUTS
*
*   NOTES
*
*   BUGS
*
*****************************************************************************
*
*/


void Psycho_L2C::TranslateThreshold(REAL *pSNR)
{
#if 1

    register int j, k;
    register REAL dNoiseBand, dThreshold, minthres;
    register ULONG      SBLimit      = psy_wSBLimit;

    register const REAL *pThreshold  = psy_pdAbsThreshold;
    register REAL       *pNB         = psy_adNB;
    register UBYTE      *pbPartition = &psy_abPartition[1];
    register REAL       *pdEnergy2   = psy_adEnergy2;

    dNoiseBand = *pNB;
    dThreshold = *pThreshold++;

    if (dNoiseBand < dThreshold) dNoiseBand = dThreshold;

    for (j = 0; j < 13; j++)
    {
        minthres = MINTHRESHOLD;
        if (minthres > dNoiseBand) minthres = dNoiseBand;

        for (k = 16; k > 0; k--)
        {
            dNoiseBand = pNB[ *pbPartition++ ];
            dThreshold = *pThreshold++;

            if (dNoiseBand < dThreshold) dNoiseBand = dThreshold;

            if (minthres > dNoiseBand) minthres = dNoiseBand;
        }
        *pSNR++ = *pdEnergy2++ / (minthres * 17.0);
    }

    for ( ; j < SBLimit; j++)
    {
        minthres = dNoiseBand;

        for (k = 16; k > 0; k--)
        {
            dNoiseBand = pNB[ *pbPartition++ ];
            dThreshold = *pThreshold++;

            if (dNoiseBand < dThreshold) dNoiseBand = dThreshold;

            minthres += dNoiseBand;
        }
        *pSNR++ = *pdEnergy2++ / minthres;
    }

#else   // Original code.

    int j, k, idx;
    REAL dNoiseBand, dThreshold, minthres;

    const REAL *pThreshold = psy_pdAbsThreshold;
    REAL       *pNB        = psy_adNB;

    dNoiseBand = pNB[0];
    dThreshold = pThreshold[0];

    if (dNoiseBand < dThreshold) dNoiseBand = dThreshold;

    for (j = 0; j < 13; j++)
    {
        minthres = MINTHRESHOLD;
        if (minthres > dNoiseBand) minthres = dNoiseBand;

        for (k = 1, idx = (j * 16) + 1; k < 17; k++, idx++)
        {
            dNoiseBand = pNB[ psy_abPartition[idx] ];
            dThreshold = pThreshold[idx];

            if (dNoiseBand < dThreshold) dNoiseBand = dThreshold;

            if (minthres > dNoiseBand) minthres = dNoiseBand;
        }
        *pSNR++ = psy_adEnergy2[j] / (minthres * 17.0);
    }

    for (j = 13; j < psy_wSBLimit; j++)
    {
        minthres = dNoiseBand;

        for (k = 1, idx = (j * 16) + 1; k < 17; k++, idx++)
        {
            dNoiseBand = pNB[ psy_abPartition[idx] ];
            dThreshold = pThreshold[idx];

            if (dNoiseBand < dThreshold) dNoiseBand = dThreshold;

            minthres += dNoiseBand;
        }
        *pSNR++ = psy_adEnergy2[j] / minthres;
    }
#endif
}
///
/// Psycho_L2C::PsychoAnal()
/****** Class Psycho_L2C/PsychoAnal *****************************************
*
*   NAME
*   Psycho_L2C::PsychoAnal -- Fonction principale.
*
*   SYNOPSIS
*   Psycho_L2C::PsychoAnal(LTMin[2][32])
*
*   void Psycho_L2C::PsychoAnal(LONG);
*
*   FUNCTION
*   Effectue l'analyse psycho-acoustique.
*
*   NOTES
*
*   BUGS
*
*****************************************************************************
*
*/


void Psycho_L2C::PsychoAnal(LONG LTMin[MAX_CH][SBLIMIT])
{
    const ULONG FLUSH = 576;

    ULONG i, j, k;
    UWORD wNumChannels = psy_roCoderC.GetNumChannels();
    WORD  *pSamples    = psy_roCoderC.GetSamplePtr();

    for (i = 0; i < 2; i++)
    {
        // Appliquer le filtre d'analyse (HANN window).

        SetupWinSamples_1024(pSamples, wNumChannels);

        // Calculer la transformée de Fourrier.

        if (wNumChannels == 2)
        {
            pSamples += FLUSH * 2;
            ComputeFastFFT_Stereo(BLKSIZE);
        }
        else
        {
            pSamples += FLUSH;
            ComputeFastFFT_Mono(BLKSIZE);
        }

        // Traiter le résultat de la transformée de Fourrier pour chaque canal.

        for (k = 0; k < wNumChannels; k++)
        {
            REAL *pSNR = &psy_adSNRTemp[k][0][0];

            psy_wNew = 1 - psy_wNew;

            CalcGroupedValues(k, psy_wNew);
            CalcPermissibleNoise();

            if (i)
            {
                LONG *pLTMin = &LTMin[k][0];

                TranslateThreshold(&pSNR[SBLIMIT]);

                for (j = 0; j < psy_wSBLimit; j++)
                {
                    REAL snr_1 = pSNR[SBLIMIT];
                    REAL snr_0 = *pSNR++;

                    if (snr_1 > snr_0) snr_0 = snr_1;

                #ifndef USE_OLD_MATHLIB
                    *pLTMin++ = (LONG) ( (LOG10_TO_LN * 1000) * log(snr_0) );
                #else
                    *pLTMin++ = (LONG) floor( (LOG10_TO_LN * 1000) * log(snr_0) );
                #endif
                }
            }
            else TranslateThreshold(pSNR);
        }
    }
}
///

#if     _USE_ASM_FCT_ == 1
#if     _CPUMODEL_ == 68020
#error  "Can't use ASM optimized version (No FPU available)!"
#endif
#endif