home *** CD-ROM | disk | FTP | other *** search
/ Magazyn Amiga 14 / MA_Cover_14.iso / source / c / pegase_src / psychoc.cpp < prev    next >
Encoding:
C/C++ Source or Header  |  1999-06-17  |  32.3 KB  |  1,133 lines
  1. /*
  2. **
  3. ** PsychoC.cpp
  4. **
  5. ** $VER: PsychoC.cpp 1.0 (01.06.98)
  6. **
  7. ** Analyseur psycho-acoustique
  8. **
  9. ** $Revision: 1.4 $
  10. ** $State: Exp $
  11. ** $Date: 1998/08/16 19:03:41 $
  12. **
  13. ** $Log: PsychoC.cpp $
  14. ** Revision 1.4  1998/08/16 19:03:41  kakace
  15. ** Version Beta3+
  16. **
  17. ** Revision 1.3  1998/07/27 14:13:05  kakace
  18. ** Validation de toutes les fonctions
  19. **
  20. ** Revision 1.2  1998/07/27 02:15:41  kakace
  21. ** Validation des fonctions (Layer II)
  22. **
  23. */
  24.  
  25.  
  26. //----------------------------------------------------------------------------------------------------
  27.  
  28. /// Includes
  29. #ifndef  _INCLUDE_ASSERT_H
  30. #include <assert.h>
  31. #endif
  32.  
  33. #ifndef  _INCLUDE_IOSTREAM_H
  34. #include <iostream.h>
  35. #endif
  36.  
  37. #ifndef  _INCLUDE_IOMANIP_H
  38. #include <iomanip.h>
  39. #endif
  40.  
  41. #ifndef  _INCLUDE_STDLIB_H
  42. #include <stdlib.h>
  43. #endif
  44.  
  45. #ifndef  _INCLUDE_MATH_H
  46. #include <math.h>
  47. #endif
  48.  
  49. #ifndef  _MEMORY_HPP
  50. #include "Memory.hpp"
  51. #endif
  52.  
  53. #ifndef  _CODER_CLASS_HPP
  54. #include "CoderC.hpp"
  55. #endif
  56.  
  57. #ifndef  _PSYCHO_CLASS_HPP
  58. #include "PsychoC.hpp"
  59. #endif
  60.  
  61. #ifndef  _PEGASECOND_HPP
  62. #include "PegaseCond.hpp"
  63. #endif
  64. ///
  65.  
  66. #define MINTHRESHOLD  60802371420160.0;
  67.  
  68.  
  69. //----------------------------------------------------------------------------------------------------
  70. //========================================== Classe PsychoC ==========================================
  71. //----------------------------------------------------------------------------------------------------
  72.  
  73. /// PsychoC::PsychoC()
  74. PsychoC::PsychoC(CoderC &roCoderC) : psy_roCoderC(roCoderC)
  75. {
  76.     psy_MemoryError = FFT_MemoryError;
  77.  
  78.     // Allouer tous les tampons nécessaires.
  79.  
  80.     psy_MemoryError |= !(psy_pdSpreading = (T_FCBCB *) AllocAlignedBuffer(sizeof(T_FCBCB)));
  81.     psy_MemoryError |= !(psy_pdTmpSave = (T_F2HBLK32 *) AllocAlignedBuffer(sizeof(T_F2HBLK32)));
  82. }
  83. ///
  84. /// PsychoC::~PsychoC() (virtual)
  85. PsychoC::~PsychoC()
  86. {
  87.     if (psy_pdSpreading)        FreeAlignedBuffer(psy_pdSpreading);
  88.     if (psy_pdTmpSave)          FreeAlignedBuffer(psy_pdTmpSave);
  89. }
  90. ///
  91.  
  92. /// PsychoC::CalcPermissibleNoise()
  93. /****** Class PsychoC/CalcPermissibleNoise **********************************
  94. *
  95. *   NAME
  96. *   PsychoC::CalcPermissibleNoise --
  97. *
  98. *   SYNOPSIS
  99. *   PsychoC::CalcPermissibleNoise()
  100. *
  101. *   void PsychoC::CalcPermissibleNoise()
  102. *
  103. *   FUNCTION
  104. *
  105. *   EXAMPLE
  106. *
  107. *   NOTES
  108. *   Fonction fortement modifiée pour optimiser la vitesse de traitement.
  109. *   A contrôler.
  110. *
  111. *   BUGS
  112. *   Fonction contrôlée conforme le 26.07.98
  113. *
  114. *****************************************************************************
  115. *
  116. *   Benchmarks (25 cycles par seconde) :
  117. *       14.06.98 : 1'05.63
  118. *
  119. *   Estimation de temps de calcul (pour 31768 frames) :
  120. *   · Log + Exp : 22 secondes
  121. *   · Calcul initial "cb" et "ecb" : 30 secondes.
  122. */
  123.  
  124.  
  125. void PsychoC::CalcPermissibleNoise()
  126. {
  127.     static float bmax[27] = {20.0, 20.0, 20.0, 20.0, 20.0, 17.0, 15.0,
  128.                              10.0,  7.0,  4.4,  4.5,  4.5,  4.5,  4.5,
  129.                               4.5,  4.5,  4.5,  4.5,  4.5,  4.5,  4.5,
  130.                               4.5,  4.5,  4.5,  3.5,  3.5,  3.5
  131.                             };
  132.  
  133.     // bmax_exp[n] = exp(-bmax[n] * LN_TO_LOG10)
  134.  
  135.     static REAL bmax_exp[27] = {0.01,             0.01,             0.01,
  136.                                 0.01,             0.01,             0.01995262314969,
  137.                                 0.03162277660168, 0.1,              0.19952623149689,
  138.                                 0.36307805477010, 0.35481338923358, 0.35481338923358,
  139.                                 0.35481338923358, 0.35481338923358, 0.35481338923358,
  140.                                 0.35481338923358, 0.35481338923358, 0.35481338923358,
  141.                                 0.35481338923358, 0.35481338923358, 0.35481338923358,
  142.                                 0.35481338923358, 0.35481338923358, 0.35481338923358,
  143.                                 0.44668359215096, 0.44668359215096, 0.44668359215096
  144.                                };
  145.  
  146. //----------------------------------------------------------------------------------------------------
  147.  
  148. #if 1
  149.  
  150.     register int j, k, l;
  151.     register REAL cb, ecb, temp;
  152.     register REAL *pSpreading;
  153.  
  154.     register REAL  *pdTMN = psy_adToneMaskingNoise;
  155.     register REAL  *pdNB  = psy_adNB;
  156.     register BYTE  *pbCB  = &psy_abCBands[0][0];
  157.  
  158.     // NB : CBANDS a été remplacé par "psy_wCBands" pour tenir compte du fait qu'il n'y a pas
  159.     //      63 bandes lorsqu'on encode en 32 KHz.
  160.  
  161.     for (j = 0; j <= psy_wCBands; j++)
  162.     {
  163.         REAL *pGrp;
  164.         cb = ecb = 0.0;
  165.  
  166.         // La valeur "Spreading" est non nulle pour des indices bien particuliers. C'est pourquoi
  167.         // le tableau "psy_abCBands" fournit l'indice de début et l'indice de fin des valeurs non
  168.         // nulles, ce qui élimine un test inutile.
  169.  
  170.         l = pbCB[64];
  171.         k = *pbCB++;
  172.         pGrp       = &psy_adGroupedValues[k][0];
  173.         pSpreading = &(*psy_pdSpreading)[j][k];
  174.  
  175.         for (; k <= l; k++)
  176.         {
  177.             temp = *pSpreading++;
  178.             ecb += temp * *pGrp++;      // = grouped_e[]
  179.             cb  += temp * *pGrp++;      // = grouped_c[]
  180.         }
  181.  
  182.         if (ecb)
  183.         {
  184.             cb = cb / ecb;
  185.             if      (cb < 0.05) cb = LOG10_005 + LOG10_2;   // = log10(0.05)
  186.             else if (cb > 0.5)  cb = LOG10_05  + LOG10_2;   // = log10(0.5)
  187.             else                cb = log10(cb) + LOG10_2;
  188.         }
  189.         else                    cb = LOG10_005 + LOG10_2;   // = log10(0.05)
  190.  
  191.         // Partie fortement modifiée par rapport à la version originale
  192.         // (Les tables ont été modifiées pour prendre en compte encore plus
  193.         //  de données précalculées).
  194.  
  195.         cb = 5.5 - cb * *pdTMN++;       // psy_adToneMaskingNoise[j];
  196.         k  = psy_abCenterBarkVal[j];
  197.  
  198.         if (cb < bmax[k]) cb = bmax_exp[k];
  199.         else              cb = exp(cb * (-LN_TO_LOG10));
  200.  
  201.         temp = psy_adRNorm[j];
  202.  
  203.         if (temp)   *pdNB++ = ecb * cb * temp;
  204.         else        *pdNB++ = temp;                     // = 0.0
  205.     }
  206.  
  207. #else   // Original code
  208.  
  209.     int j, k, l;
  210.     REAL cb, ecb, temp;
  211.     REAL *pSpreading;
  212.  
  213.     // NB : CBANDS a été remplacé par "psy_wCBands" pour tenir compte du fait qu'il n'y a pas
  214.     //      63 bandes lorsqu'on encode en 32 KHz.
  215.  
  216.     for (j = 0; j <= psy_wCBands; j++)
  217.     {
  218.         REAL *pGrp;
  219.         cb = ecb = 0.0;
  220.  
  221.         // La valeur "Spreading" est non nulle pour des indices bien particuliers. C'est pourquoi
  222.         // le tableau "psy_abCBands" fournit l'indice de début et l'indice de fin des valeurs non
  223.         // nulles, ce qui élimine un test inutile.
  224.  
  225.         k = psy_abCBands[0][j];
  226.         l = psy_abCBands[1][j];
  227.         pGrp       = &psy_adGroupedValues[k][0];
  228.         pSpreading = &(*psy_pdSpreading)[j][k];
  229.  
  230.         for (; k <= l; k++)
  231.         {
  232.             temp = *pSpreading++;
  233.             ecb += temp * *pGrp++;      // = grouped_e[]
  234.             cb  += temp * *pGrp++;      // = grouped_c[]
  235.         }
  236.  
  237.         if (ecb != 0.0)
  238.         {
  239.             cb = cb / ecb;
  240.             if      (cb < 0.05) cb = LOG10_005 + LOG10_2;   // = log10(0.05)
  241.             else if (cb > 0.5)  cb = LOG10_05 + LOG10_2;    // = log10(0.5)
  242.             else                cb = log10(cb) + LOG10_2;
  243.         }
  244.         else                    cb = LOG10_005 + LOG10_2;   // = log10(0.05)
  245.  
  246.         // Partie fortement modifiée par rapport à la version originale
  247.         // (Les tables ont été modifiées pour prendre en compte encore plus
  248.         //  de données précalculées).
  249.  
  250.         cb = 5.5 - cb * psy_adToneMaskingNoise[j];
  251.         k  = psy_abCenterBarkVal[j];
  252.  
  253.         if (cb < bmax[k]) cb = bmax_exp[k];
  254.         else              cb = exp(cb * (-LN_TO_LOG10));
  255.  
  256.         temp = psy_adRNorm[j];
  257.  
  258.         if (temp)   psy_adNB[j] = ecb * cb * temp;
  259.         else        psy_adNB[j] = temp;                     // = 0.0
  260.     }
  261. #endif
  262. }
  263. ///
  264. /// PsychoC::CalcGroupedValues()
  265. /****** Class PsychoC/CalcGroupedValues *************************************
  266. *
  267. *   NAME
  268. *   PsychoC::CalcGroupedValues --
  269. *
  270. *   SYNOPSIS
  271. *   PsychoC::CalcGroupedValues(wChannel, wNew)
  272. *
  273. *   void PsychoC::CalcGroupedValues(UWORD, UWORD);
  274. *
  275. *   FUNCTION
  276. *   Calcule des vecteurs complexes à partir des données fournies par la FFT,
  277. *   et les cumule par bandes de fréquences.
  278. *
  279. *   INPUTS
  280. *   wChannel - Numéro de canal.
  281. *   wNew     - Indice courant
  282. *
  283. *   NOTES
  284. *
  285. *   BUGS
  286. *   Cette fonction n'est pas strictement conforme au source originale : L'
  287. *   énergie du complexe C[512] est ici bornée (borne inférieure = 0.0005).
  288. *
  289. *   Fonction validée le 26/07/98
  290. *
  291. *   SEE ALSO
  292. *
  293. *****************************************************************************
  294. *
  295. *   Benchmarks (3 cycles par seconde) :
  296. *       21.06.98 : 2'47.80
  297. *
  298. *   NB : La clarté du code souffre un peu de l'entrelacement des calculs
  299. *        destinés à occuper le FPU pendant que le cache termine (en principe)
  300. *        la lecture d'une ligne.
  301. */
  302.  
  303. #if _USE_ASM_FCT_ == 1
  304.  
  305. extern "C"{
  306. void CalcGroupedValues(UWORD wNew, REAL *GroupedValues, T_FFTDATAS *FFTOut, REAL *Energy2, T_FFTDATAS *tmp, UBYTE *Part);
  307. }
  308.  
  309. inline void PsychoC::CalcGroupedValues(UWORD wChannel, UWORD wNew)
  310. {
  311.     ::CalcGroupedValues(wNew, (REAL *) psy_adGroupedValues, &(*FFT_pdWinSmplShuffled)[wChannel][0][0],
  312.                         psy_adEnergy2, &(*psy_pdTmpSave)[wChannel][0][0][0], psy_abPartition);
  313. }
  314.  
  315. #else
  316.  
  317. void PsychoC::CalcGroupedValues(UWORD wChannel, UWORD wNew)
  318. {
  319. #define pOldest  pNew                       // Synonyme.
  320.  
  321.     const ULONG COSINE = 0;                 // Index d'accès aux données.
  322.     const ULONG SINE   = 2;
  323.     const ULONG ENERGY = 4;
  324.     const ULONG NEXT_OFFSET = 6;            // {COSINE, SINE, ENERGY} × {NEW, OLD}
  325.  
  326.     UWORD wOld = 1 - wNew;
  327.  
  328.     int j, part, idx;
  329.  
  330.     double cos_phi, cos_2phi;
  331.     double sin_phi, sin_2phi;
  332.     double energy, r, r_prime;
  333.  
  334.     T_FFTDATAS *p;
  335.     T_FFTDATAS *pNew, *pOld;
  336.     REAL *Grp;
  337.  
  338.     // Initialiser les tableaux destinés à accueillir les sommes intermédiaires.
  339.  
  340.     {
  341.         REAL *p = &psy_adGroupedValues[0][0];
  342.  
  343.         for (j = 0; j <= psy_wCBands; j++)
  344.         {
  345.             *p++ = 0.0;                 // Grouped_e[j]
  346.             *p++ = 0.0;                 // Grouped_c[j]
  347.         }
  348.     }
  349.  
  350.     p       = &(*FFT_pdWinSmplShuffled)[wChannel][0][0];
  351.     pOld    = &(*psy_pdTmpSave)[wChannel][0][0][wOld];
  352.     pNew    = &(*psy_pdTmpSave)[wChannel][0][0][wNew];
  353.  
  354.     for (j = 0; j < HBLKSIZE; j++)
  355.     {
  356.         // Initialiser les variables intermédiaires.
  357.         // Les tables pOld et pNew se présentent sous la forme :
  358.         //      Cos_new, Cos_old, Sin_new, Sin_old, Nrj_new, Nrj_old, Cos_new...
  359.         // NB : Cette partie s'efforce d'effectuer un maximum de calculs entre chaque
  360.         //      lecture dans une table de façon à ne pas être bloqué par un cache en
  361.         //      train de lire une ligne.
  362.  
  363.         cos_phi  = pOld[COSINE];                        // = cos(a)
  364.         r        = cos_phi + cos_phi;                   // = 2.cos(a)
  365.         r_prime  = r * cos_phi;                         // = 2.cos²(a)
  366.         cos_phi  = pOldest[COSINE];                     // = cos(b)
  367.  
  368.         sin_2phi = pOld[SINE] * r;                      // = sin(2a) = 2.cos(a).sin(a)
  369.         r        = r_prime - 1.0;                       // = cos(2a) = 2.cos²(a) - 1
  370.         energy   = cos_phi * r;                         // Commencer le calcul de cos(2a-b) ( = cos(2a).cos(b))
  371.         sin_phi  = pOldest[SINE];                       // = sin(b)
  372.  
  373.         r_prime   = pOld[ENERGY];
  374.         energy  += sin_phi * sin_2phi;                  // cos(2a-b) = cos(2a).cos(b) + sin(2a).sin(b)
  375.         sin_2phi = sin_2phi * cos_phi - r * sin_phi;    // sin(2a-b) = sin(2a).cos(b) - cos(2a).sin(b)
  376.         r_prime += r_prime - pOldest[ENERGY];
  377.  
  378.         cos_2phi = energy;                              // HINT : "energy" est dans un registre, ce qui
  379.                                                         // n'est pas le cas de "cos_2phi"
  380.  
  381.         // Récupérer les données en sortie de FFT (elle fournie gracieusement les
  382.         // sin() et cos() voulus : inutile de passer par un angle intermédiaire).
  383.  
  384.         cos_phi = *p++;                                 // (*FFT_pdWinSmplShuffled)[wChannel][j][0];
  385.         energy  = cos_phi * cos_phi;
  386.         r       = 0.0005;
  387.         sin_phi = *p++;                                 // (*FFT_pdWinSmplShuffled)[wChannel][j][1];
  388.         energy += sin_phi * sin_phi;
  389.  
  390.         // Donner une valeur minimum à l'énergie.
  391.         // NOTE : Dans le source d'origine, l'énergie n'est pas bornée (borne inférieure)
  392.         //        pour le complexe C[512].
  393.  
  394.         if (energy <= r)
  395.         {
  396.             energy  = r;                                // energy = 0.0005
  397.             r       = SQRT_00005;                       // = sqrt(0.0005) = sqrt(5) / 100
  398.             cos_phi = r;                                // => phi = 0
  399.             sin_phi = 0.0;
  400.         }
  401.         else
  402.         {
  403.             r = sqrt(energy);
  404.         }
  405.  
  406.         // Préparer la somme des énergies pour TranslateThreshold().
  407.         // Par la même occasion, on occupe l'IU pendant que le FPU calcule SQRT().
  408.  
  409.         part = psy_abPartition[j];
  410.         Grp  = &psy_adGroupedValues[part][0];
  411.  
  412.         idx  = j >> 4;
  413.         part = j & 15;
  414.  
  415.         if (part)           psy_adEnergy2[idx] += energy;   // part != 0 : Ajouter l'énergie à la bande courante.
  416.         else
  417.         {
  418.             if (idx != 32)  psy_adEnergy2[idx]  = energy;   // part == 0 && idx != 32 : Initialiser la somme.
  419.             if (idx--)      psy_adEnergy2[idx] += energy;   // part == 0 && idx != 0  : Ajouter à la bande précédente.
  420.         }
  421.  
  422.         // Mémoriser les valeurs courantes pour la prochaine passe.
  423.  
  424.         pNew[ENERGY] = r;                               // r != 0 puisque l'énergie est bornée.
  425.  
  426.         r = 1.0 / r;
  427.         pNew[COSINE] = cos_phi * r;
  428.         pNew[SINE]   = sin_phi * r;
  429.  
  430.         r = pNew[ENERGY] + fabs(r_prime);               // r > 0 puisque l'energie est bornée.
  431.  
  432.         cos_phi -= r_prime * cos_2phi;
  433.         sin_phi -= r_prime * sin_2phi;
  434.         cos_phi  = cos_phi * cos_phi + sin_phi * sin_phi;
  435.               // = energy + r'² - 2.r.r'.cos(phi + phi')
  436.  
  437.         Grp[0] += energy;
  438.         Grp[1] += energy * sqrt(cos_phi) / r;
  439.  
  440.         pNew += NEXT_OFFSET;                            // Ajuster les pointeurs.
  441.         pOld += NEXT_OFFSET;                            // (pOldest = pNew)
  442.     }
  443. }
  444.  
  445. #endif  // _USE_ASM_FCT_
  446. ///
  447. /// PsychoC::InitPsychoAnalyzer()
  448. /****** Class PsychoC/InitPsychoAnalyzer ************************************
  449. *
  450. *   NAME
  451. *   PsychoC::InitPsychoAnalyzer -- Initialiser l'analyseur psycho-acoustique.
  452. *
  453. *   SYNOPSIS
  454. *   result = PsychoC::InitPsychoAnalyzer(iSBLimit)
  455. *
  456. *   bool PsychoC::InitPsychoAnalyzer(int);
  457. *
  458. *   FUNCTION
  459. *   Initialise l'analyseur psycho-acoustique en préparant les différentes
  460. *   tables nécessaires à son fonctionnement.
  461. *
  462. *   INPUTS
  463. *   iSBLimit - Nombre de sous-bandes.
  464. *
  465. *   RESULT
  466. *   result = FALSE en cas d'erreur (tampons non alloués par exemple).
  467. *
  468. *   NOTES
  469. *   Une table SinCos() devrait aussi être initialisée afin d'améliorer la
  470. *   vitesse des FFT lorsque le CPU utilisé est plus handicapé par une capa-
  471. *   cité de traitement limitée que par ses accès mémoire (68020/68030).
  472. *
  473. *   BUGS
  474. *   Il n'existe pas de table AbsoluteThreshold pour le MPEG-2 !!!
  475. *   Fonction contrôlée conforme le 27/07/98
  476. *
  477. *   SEE ALSO
  478. *
  479. *****************************************************************************
  480. *
  481. */
  482.  
  483.  
  484. int PsychoC::InitPsychoAnalyzer(int iSBLimit)
  485. {
  486.     extern const REAL a_dAbsoluteThresholds[3][HBLKSIZE];
  487.  
  488.     REAL adCenterBarkValue[CBANDS];         // (cbval[])
  489.  
  490.     int Result = TRUE;
  491.     int i, j;
  492.  
  493.     // Initialiser les paramètres de l'analyseur psycho-acoustique selon le type de Layer.
  494.  
  495.     psy_wNew     = 0;
  496.     psy_wSBLimit = iSBLimit;
  497.  
  498.     // Adresser les tables des Thresholds en fonction de la fréquence d'échantillonage
  499.  
  500.     switch (psy_roCoderC.GetSampleFreq())
  501.     {
  502.         case 32000: psy_pdAbsThreshold = &a_dAbsoluteThresholds[0][0]; break;
  503.         case 44100: psy_pdAbsThreshold = &a_dAbsoluteThresholds[1][0]; break;
  504.         case 48000: psy_pdAbsThreshold = &a_dAbsoluteThresholds[2][0]; break;
  505.         default:
  506.             Result = FALSE;
  507.             break;
  508.     }
  509.  
  510.     // Initialiser les données pour le FFT.
  511.  
  512.     {
  513.         static int crit_band[27] = {    0,   100,   200,   300,   400,   510,   630,   770,
  514.                                       920,  1080,  1270,  1480,  1720,  2000,  2320,  2700,
  515.                                      3150,  3700,  4400,  5300,  6400,  7700,  9500, 12000,
  516.                                     15500, 25000, 30000
  517.                                    };
  518.  
  519.         const REAL INV_BLKSIZE = 1.0 / BLKSIZE;
  520.         REAL  t, cb, cb_pred;
  521.         REAL  bval_lo;
  522.         int   part, idx;
  523.  
  524.         REAL  freq_mult = psy_roCoderC.GetSampleFreq() * INV_BLKSIZE;
  525.  
  526.         idx = 1;
  527.         t = freq_mult;
  528.         cb = crit_band[1];
  529.         cb_pred = 0.0;
  530.  
  531.         bval_lo = 0.0;                              // Fréquence de départ de la bande courante.
  532.         adCenterBarkValue[0] = 0.0;
  533.         part = 0;                                   // Numéro d'index de la bande courante.
  534.         j = 1;                                      // Nombre de fréquences dans la bande courante.
  535.  
  536.         for (i = 1; i < HBLKSIZE; i++, t += freq_mult)
  537.         {
  538.             REAL bval;
  539.  
  540.             if (t > cb)
  541.             {
  542.                 idx++;
  543.                 cb_pred = cb;
  544.                 cb = crit_band[idx];
  545.             }
  546.  
  547.             bval = idx - 1 + (t - cb_pred) / (cb - cb_pred);
  548.  
  549.             if ( (bval - bval_lo) > CB_FRACTION)
  550.             {
  551.                 adCenterBarkValue[part] /= j;       // Ajuster les infos de la bande précédente.
  552.                 psy_abNumLines[part] = j;
  553.  
  554.                 psy_abPartition[i] = ++part;        // Initialiser la bande suivante.
  555.                 adCenterBarkValue[part] = bval_lo = bval;
  556.                 j = 1;
  557.             }
  558.             else
  559.             {
  560.                 psy_abPartition[i] = part;          // Ajouter la fréquence dans la bande courante.
  561.                 adCenterBarkValue[part] += bval;
  562.                 j++;
  563.             }
  564.         }
  565.  
  566.         // Ajuster les informations de la dernière bande.
  567.  
  568.         psy_wCBands = part = psy_abPartition[HBLKSIZE-1];
  569.         psy_abNumLines[part] = j;
  570.         adCenterBarkValue[part] /= j;
  571.     }
  572.  
  573.     // Calculer la fonction de "spreading"
  574.  
  575.     {
  576.         const REAL cst = 17.5 * 1.106650803 - 3.555;      // = 17.5 * sqrt(1.224676) - 3.555 = 15.811389
  577.         REAL t1, t2, t3;
  578.  
  579.         // NB : CBANDS a été remplacé par "psy_wCBands" pour tenir compte du fait qu'il n'y a pas
  580.         //      63 bandes lorsqu'on encode en 32 KHz.
  581.  
  582.         for (j = 0; j <= psy_wCBands; j++)
  583.         {
  584.             psy_abCBands[0][j] = psy_abCBands[1][j] = -1;
  585.  
  586.             for (i = 0; i <= psy_wCBands; i++)
  587.             {
  588.                 t1 = (adCenterBarkValue[j] - adCenterBarkValue[i]) * 1.05;
  589.  
  590.                 if (t1 >= 0.5 && t1 <= 2.5)
  591.                 {
  592.                     t2 = t1 - 0.5;
  593.                     t2 = 8.0 * t2 * (t2 - 2.0);     // = 8.(t2.t2 - 2.t2)
  594.                 }
  595.                 else t2 = 0.0;
  596.  
  597.                 t1 += 0.474;
  598.                 t3 = cst + 7.5 * t1 - 17.5 * sqrt(1.0 + t1 * t1);
  599.  
  600.                 if (t3 <= -100.0) (*psy_pdSpreading)[j][i] = 0.0;
  601.                 else
  602.                 {
  603.                     if (psy_abCBands[0][j] < 0) psy_abCBands[0][j] = i;
  604.                     psy_abCBands[1][j] = i;
  605.  
  606.                     t3 = (t2 + t3) * LN_TO_LOG10;
  607.                     (*psy_pdSpreading)[j][i] = exp(t3);
  608.                 }
  609.             }
  610.         }
  611.     }
  612.  
  613.     // Calculer les valeurs "Tone Masking Noise"
  614.  
  615.     {
  616.         REAL t1;
  617.  
  618.         // NB : CBANDS a été remplacé par "psy_wCBands" pour tenir compte du fait qu'il n'y a pas
  619.         //      63 bandes lorsqu'on encode en 32 KHz.
  620.  
  621.         for (j = 0; j <= psy_wCBands; j++)
  622.         {
  623.             t1 = adCenterBarkValue[j];
  624.             psy_abCenterBarkVal[j] = (UBYTE) ceil(t1);
  625.  
  626.             // Valeur précalculée pour CalcPersmissibleNoise()
  627.             psy_adToneMaskingNoise[j] = (t1 > 9.0) ? (t1 + 10.0) : 19.0;
  628.  
  629.             t1 = (*psy_pdSpreading)[j][0];
  630.  
  631.             for (i = 1; i <= psy_wCBands; i++)
  632.             {
  633.                 t1 += (*psy_pdSpreading)[j][i];
  634.             }
  635.  
  636.             // Précalculer la valeur utilisée par CalcPermissibleNoise().
  637.  
  638.             t1 *= psy_abNumLines[j];
  639.             if (t1 != 0.0) t1 = 1.0 / t1;
  640.  
  641.             psy_adRNorm[j] = t1;
  642.         }
  643.     }
  644.  
  645.     // Initialiser la table TmpSave[]
  646.  
  647.     for (j = 0; j < HBLKSIZE; j++)
  648.     {
  649.         (*psy_pdTmpSave)[0][j][0][0] = (*psy_pdTmpSave)[0][j][0][1] = 1;
  650.         (*psy_pdTmpSave)[1][j][0][0] = (*psy_pdTmpSave)[1][j][0][1] = 1;
  651.     }
  652.  
  653.     return Result;
  654. }
  655. ///
  656.  
  657. //----------------------------------------------------------------------------------------------------
  658. //======================================== Classe Psycho_L1C =========================================
  659. //----------------------------------------------------------------------------------------------------
  660.  
  661. /// Psycho_L1C::TranslateThreshold()  (virtual)
  662. /****** Class Psycho_L1C/TranslateThreshold *********************************
  663. *
  664. *   NAME
  665. *   Psycho_L1C::TranslateThreshold --
  666. *
  667. *   SYNOPSIS
  668. *   Psycho_L1C::TranslateThreshold(wChannel)
  669. *
  670. *   void Psycho_L1C::TranslateThreshold(UWORD);
  671. *
  672. *   FUNCTION
  673. *
  674. *   INPUTS
  675. *   wChannel - Numéro de canal.
  676. *
  677. *   NOTES
  678. *
  679. *   BUGS
  680. *   Fonction contrôlée conforme le 27/07/98
  681. *
  682. *****************************************************************************
  683. *
  684. */
  685.  
  686. /*
  687. inline REAL Psycho_L1C::GetNoiseLevel(REAL dNoiseBand, REAL dThreshold, REAL *p_LThr)
  688. {
  689.     if (dNoiseBand < dThreshold) dNoiseBand = dThreshold;
  690.  
  691.     dThreshold = *p_LThr;
  692.     *p_LThr = LXMIN * dNoiseBand;       // = 32 * dNoiseBand
  693.  
  694.     if (dNoiseBand < dThreshold) dThreshold = dNoiseBand;
  695.     dNoiseBand *= 0.00316;
  696.  
  697.     if (dNoiseBand < dThreshold) dNoiseBand = dThreshold;
  698.  
  699.     return dNoiseBand;
  700. }
  701. */
  702.  
  703. void Psycho_L1C::TranslateThreshold(UWORD wChannel)
  704. {
  705. #if 1
  706.  
  707.     register int j, k;
  708.     register REAL dNoiseBand, dThreshold, minthres;
  709.     register ULONG   SBLimit = psy_wSBLimit;
  710.  
  711.     register const REAL *pThreshold  = psy_pdAbsThreshold;
  712.     register REAL       *pLThr       = &(*psy_pdLThr)[wChannel][0];
  713.     register REAL       *pNB         = psy_adNB;
  714.     register UBYTE      *pbPartition = &psy_abPartition[1];
  715.     register REAL       *pdEnergy2   = psy_adEnergy2;
  716.     register REAL       *pdSNR       = psy_adSNRTemp;
  717.  
  718.     dNoiseBand = *pNB;
  719.     dThreshold = *pThreshold++;
  720.  
  721.     if (dNoiseBand < dThreshold) dNoiseBand = dThreshold;
  722.  
  723.     dThreshold = *pLThr;
  724.     *pLThr++ = LXMIN * dNoiseBand;          // = 32 * dNoiseBand
  725.  
  726.     if (dNoiseBand < dThreshold) dThreshold = dNoiseBand;
  727.     dNoiseBand *= 0.00316;
  728.  
  729.     if (dNoiseBand < dThreshold) dNoiseBand = dThreshold;
  730.  
  731.     for (j = 0; j < 13; j++)
  732.     {
  733.         minthres = MINTHRESHOLD;
  734.         if (minthres > dNoiseBand) minthres = dNoiseBand;
  735.  
  736.         for (k = 16; k > 0; k--)
  737.         {
  738.             dNoiseBand = pNB[ *pbPartition++ ];
  739.             dThreshold = *pThreshold++;
  740.  
  741.             if (dNoiseBand < dThreshold) dNoiseBand = dThreshold;
  742.  
  743.             dThreshold = *pLThr;
  744.             *pLThr++ = LXMIN * dNoiseBand;          // = 32 * dNoiseBand
  745.  
  746.             if (dNoiseBand < dThreshold) dThreshold = dNoiseBand;
  747.             dNoiseBand *= 0.00316;
  748.  
  749.             if (dNoiseBand < dThreshold) dNoiseBand = dThreshold;
  750.  
  751.             if (minthres > dNoiseBand) minthres = dNoiseBand;
  752.         }
  753.         *pdSNR++ = *pdEnergy2++ / (minthres * 17.0);
  754.     }
  755.  
  756.     for (j = 13; j < SBLimit; j++)
  757.     {
  758.         minthres = dNoiseBand;
  759.  
  760.         for (k = 16; k > 0; k--)
  761.         {
  762.             dNoiseBand = pNB[ *pbPartition++ ];
  763.             dThreshold = *pThreshold++;
  764.  
  765.             if (dNoiseBand < dThreshold) dNoiseBand = dThreshold;
  766.  
  767.             dThreshold = *pLThr;
  768.             *pLThr++ = LXMIN * dNoiseBand;          // = 32 * dNoiseBand
  769.  
  770.             if (dNoiseBand < dThreshold) dThreshold = dNoiseBand;
  771.             dNoiseBand *= 0.00316;
  772.  
  773.             if (dNoiseBand < dThreshold) dNoiseBand = dThreshold;
  774.  
  775.             minthres += dNoiseBand;
  776.         }
  777.         *pdSNR++ = *pdEnergy2++ / minthres;
  778.     }
  779.  
  780. #else   // Original code.
  781.  
  782.     int j, k, idx;
  783.     REAL dNoiseBand, dThreshold, minthres;
  784.  
  785.     const REAL *pThreshold = psy_pdAbsThreshold;
  786.     REAL       *pLThr      = &(*psy_pdLThr)[wChannel][0];
  787.     REAL       *pNB        = psy_adNB;
  788.  
  789.     dNoiseBand = GetNoiseLevel( pNB[0], pThreshold[0], pLThr);
  790.  
  791.     for (j = 0; j < 13; j++)
  792.     {
  793.         minthres = MINTHRESHOLD;
  794.         if (minthres > dNoiseBand) minthres = dNoiseBand;
  795.  
  796.         for (k = 1, idx = (j * 16) + 1; k < 17; k++, idx++)
  797.         {
  798.             dNoiseBand = pNB[ psy_abPartition[idx] ];
  799.             dThreshold = pThreshold[idx];
  800.  
  801.             dNoiseBand = GetNoiseLevel(dNoiseBand, dThreshold, &pLThr[idx]);
  802.  
  803.             if (minthres > dNoiseBand) minthres = dNoiseBand;
  804.         }
  805.         psy_adSNRTemp[j] = psy_adEnergy2[j] / (minthres * 17.0);
  806.     }
  807.  
  808.     for (j = 13; j < psy_wSBLimit; j++)
  809.     {
  810.         minthres = dNoiseBand;
  811.  
  812.         for (k = 1, idx = (j * 16) + 1; k < 17; k++, idx++)
  813.         {
  814.             dNoiseBand = pNB[ psy_abPartition[idx] ];
  815.             dThreshold = pThreshold[idx];
  816.  
  817.             dNoiseBand = GetNoiseLevel(dNoiseBand, dThreshold, &pLThr[idx]);
  818.             minthres += dNoiseBand;
  819.         }
  820.         psy_adSNRTemp[j] = psy_adEnergy2[j] / minthres;
  821.     }
  822. #endif
  823. }
  824. ///
  825. /// Psycho_L1C::PsychoAnal()
  826. /****** Class Psycho_L1C/PsychoAnal *****************************************
  827. *
  828. *   NAME
  829. *   Psycho_L1C::PsychoAnal -- Fonction principale.
  830. *
  831. *   SYNOPSIS
  832. *   Psycho_L1C::PsychoAnal(LTMin[2][32])
  833. *
  834. *   void Psycho_L1C::PsychoAnal(LONG);
  835. *
  836. *   FUNCTION
  837. *   Effectue l'analyse psycho-acoustique.
  838. *
  839. *   NOTES
  840. *
  841. *   BUGS
  842. *
  843. *****************************************************************************
  844. *
  845. */
  846.  
  847.  
  848. void Psycho_L1C::PsychoAnal(LONG LTMin[MAX_CH][SBLIMIT])
  849. {
  850.     const ULONG FLUSH = 384;
  851.  
  852.     ULONG j, k;
  853.     UWORD wNumChannels = psy_roCoderC.GetNumChannels();
  854.     WORD  *pSamples    = psy_roCoderC.GetSamplePtr();
  855.  
  856.     // Appliquer le filtre d'analyse (HANN window).
  857.  
  858.     SetupWinSamples_1024(pSamples, wNumChannels);
  859.  
  860.     // Calculer la transformée de Fourrier.
  861.  
  862.     if (wNumChannels == 2)
  863.     {
  864.         pSamples += FLUSH * 2;
  865.         ComputeFastFFT_Stereo(BLKSIZE);
  866.     }
  867.     else
  868.     {
  869.         pSamples += FLUSH;
  870.         ComputeFastFFT_Mono(BLKSIZE);
  871.     }
  872.  
  873.     // Exploiter les résultats de la transformée de Fourrier pour chaque canal.
  874.  
  875.     for (k = 0, psy_wNew = 1 - psy_wNew; k < wNumChannels; k++)
  876.     {
  877.         CalcGroupedValues(k, psy_wNew);
  878.         CalcPermissibleNoise();
  879.         TranslateThreshold(k);
  880.  
  881.         for (j = 0; j < psy_wSBLimit; j++)
  882.         {
  883.         #ifndef USE_OLD_MATHLIB
  884.             LTMin[k][j] = (LONG) ( (LOG10_TO_LN * 1000) * log(psy_adSNRTemp[j]) );
  885.         #else
  886.             LTMin[k][j] = (LONG) floor( (LOG10_TO_LN * 1000) * log(psy_adSNRTemp[j]) );
  887.         #endif
  888.         }
  889.     }
  890. }
  891. ///
  892. /// Psycho_L1C::InitPsychoAnalyzer()
  893. //
  894. //  Initialiser la table "psy_pdLThr", uniquement utilisée avec le Layer I
  895.  
  896. int Psycho_L1C::InitPsychoAnalyzer(int iSBLimit)
  897. {
  898.     // Initialiser la table des Threshold mini
  899.  
  900.     REAL *p1 = &(*psy_pdLThr)[0][0];
  901.     REAL *p2 = &(*psy_pdLThr)[1][0];
  902.     REAL d   = MINTHRESHOLD;
  903.  
  904.     for (int i = 0; i < HBLKSIZE; i++)
  905.     {
  906.         *p1++ = *p2++ = d;
  907.     }
  908.  
  909.     return PsychoC::InitPsychoAnalyzer(iSBLimit);          // Initialisation générale.
  910. }
  911. ///
  912.  
  913. //----------------------------------------------------------------------------------------------------
  914. //======================================== Classe Psycho_L2C =========================================
  915. //----------------------------------------------------------------------------------------------------
  916.  
  917. /// Psycho_L2C::TranslateThreshold()  (virtual)
  918. /****** Class Psycho_L2C/TranslateThreshold *********************************
  919. *
  920. *   NAME
  921. *   Psycho_L2C::TranslateThreshold --
  922. *
  923. *   SYNOPSIS
  924. *   Psycho_L2C::TranslateThreshold(pSNR)
  925. *
  926. *   void Psycho_L2C::TranslateThreshold(REAL *);
  927. *
  928. *   FUNCTION
  929. *
  930. *   INPUTS
  931. *
  932. *   NOTES
  933. *
  934. *   BUGS
  935. *
  936. *****************************************************************************
  937. *
  938. */
  939.  
  940.  
  941. void Psycho_L2C::TranslateThreshold(REAL *pSNR)
  942. {
  943. #if 1
  944.  
  945.     register int j, k;
  946.     register REAL dNoiseBand, dThreshold, minthres;
  947.     register ULONG      SBLimit      = psy_wSBLimit;
  948.  
  949.     register const REAL *pThreshold  = psy_pdAbsThreshold;
  950.     register REAL       *pNB         = psy_adNB;
  951.     register UBYTE      *pbPartition = &psy_abPartition[1];
  952.     register REAL       *pdEnergy2   = psy_adEnergy2;
  953.  
  954.     dNoiseBand = *pNB;
  955.     dThreshold = *pThreshold++;
  956.  
  957.     if (dNoiseBand < dThreshold) dNoiseBand = dThreshold;
  958.  
  959.     for (j = 0; j < 13; j++)
  960.     {
  961.         minthres = MINTHRESHOLD;
  962.         if (minthres > dNoiseBand) minthres = dNoiseBand;
  963.  
  964.         for (k = 16; k > 0; k--)
  965.         {
  966.             dNoiseBand = pNB[ *pbPartition++ ];
  967.             dThreshold = *pThreshold++;
  968.  
  969.             if (dNoiseBand < dThreshold) dNoiseBand = dThreshold;
  970.  
  971.             if (minthres > dNoiseBand) minthres = dNoiseBand;
  972.         }
  973.         *pSNR++ = *pdEnergy2++ / (minthres * 17.0);
  974.     }
  975.  
  976.     for ( ; j < SBLimit; j++)
  977.     {
  978.         minthres = dNoiseBand;
  979.  
  980.         for (k = 16; k > 0; k--)
  981.         {
  982.             dNoiseBand = pNB[ *pbPartition++ ];
  983.             dThreshold = *pThreshold++;
  984.  
  985.             if (dNoiseBand < dThreshold) dNoiseBand = dThreshold;
  986.  
  987.             minthres += dNoiseBand;
  988.         }
  989.         *pSNR++ = *pdEnergy2++ / minthres;
  990.     }
  991.  
  992. #else   // Original code.
  993.  
  994.     int j, k, idx;
  995.     REAL dNoiseBand, dThreshold, minthres;
  996.  
  997.     const REAL *pThreshold = psy_pdAbsThreshold;
  998.     REAL       *pNB        = psy_adNB;
  999.  
  1000.     dNoiseBand = pNB[0];
  1001.     dThreshold = pThreshold[0];
  1002.  
  1003.     if (dNoiseBand < dThreshold) dNoiseBand = dThreshold;
  1004.  
  1005.     for (j = 0; j < 13; j++)
  1006.     {
  1007.         minthres = MINTHRESHOLD;
  1008.         if (minthres > dNoiseBand) minthres = dNoiseBand;
  1009.  
  1010.         for (k = 1, idx = (j * 16) + 1; k < 17; k++, idx++)
  1011.         {
  1012.             dNoiseBand = pNB[ psy_abPartition[idx] ];
  1013.             dThreshold = pThreshold[idx];
  1014.  
  1015.             if (dNoiseBand < dThreshold) dNoiseBand = dThreshold;
  1016.  
  1017.             if (minthres > dNoiseBand) minthres = dNoiseBand;
  1018.         }
  1019.         *pSNR++ = psy_adEnergy2[j] / (minthres * 17.0);
  1020.     }
  1021.  
  1022.     for (j = 13; j < psy_wSBLimit; j++)
  1023.     {
  1024.         minthres = dNoiseBand;
  1025.  
  1026.         for (k = 1, idx = (j * 16) + 1; k < 17; k++, idx++)
  1027.         {
  1028.             dNoiseBand = pNB[ psy_abPartition[idx] ];
  1029.             dThreshold = pThreshold[idx];
  1030.  
  1031.             if (dNoiseBand < dThreshold) dNoiseBand = dThreshold;
  1032.  
  1033.             minthres += dNoiseBand;
  1034.         }
  1035.         *pSNR++ = psy_adEnergy2[j] / minthres;
  1036.     }
  1037. #endif
  1038. }
  1039. ///
  1040. /// Psycho_L2C::PsychoAnal()
  1041. /****** Class Psycho_L2C/PsychoAnal *****************************************
  1042. *
  1043. *   NAME
  1044. *   Psycho_L2C::PsychoAnal -- Fonction principale.
  1045. *
  1046. *   SYNOPSIS
  1047. *   Psycho_L2C::PsychoAnal(LTMin[2][32])
  1048. *
  1049. *   void Psycho_L2C::PsychoAnal(LONG);
  1050. *
  1051. *   FUNCTION
  1052. *   Effectue l'analyse psycho-acoustique.
  1053. *
  1054. *   NOTES
  1055. *
  1056. *   BUGS
  1057. *
  1058. *****************************************************************************
  1059. *
  1060. */
  1061.  
  1062.  
  1063. void Psycho_L2C::PsychoAnal(LONG LTMin[MAX_CH][SBLIMIT])
  1064. {
  1065.     const ULONG FLUSH = 576;
  1066.  
  1067.     ULONG i, j, k;
  1068.     UWORD wNumChannels = psy_roCoderC.GetNumChannels();
  1069.     WORD  *pSamples    = psy_roCoderC.GetSamplePtr();
  1070.  
  1071.     for (i = 0; i < 2; i++)
  1072.     {
  1073.         // Appliquer le filtre d'analyse (HANN window).
  1074.  
  1075.         SetupWinSamples_1024(pSamples, wNumChannels);
  1076.  
  1077.         // Calculer la transformée de Fourrier.
  1078.  
  1079.         if (wNumChannels == 2)
  1080.         {
  1081.             pSamples += FLUSH * 2;
  1082.             ComputeFastFFT_Stereo(BLKSIZE);
  1083.         }
  1084.         else
  1085.         {
  1086.             pSamples += FLUSH;
  1087.             ComputeFastFFT_Mono(BLKSIZE);
  1088.         }
  1089.  
  1090.         // Traiter le résultat de la transformée de Fourrier pour chaque canal.
  1091.  
  1092.         for (k = 0; k < wNumChannels; k++)
  1093.         {
  1094.             REAL *pSNR = &psy_adSNRTemp[k][0][0];
  1095.  
  1096.             psy_wNew = 1 - psy_wNew;
  1097.  
  1098.             CalcGroupedValues(k, psy_wNew);
  1099.             CalcPermissibleNoise();
  1100.  
  1101.             if (i)
  1102.             {
  1103.                 LONG *pLTMin = <Min[k][0];
  1104.  
  1105.                 TranslateThreshold(&pSNR[SBLIMIT]);
  1106.  
  1107.                 for (j = 0; j < psy_wSBLimit; j++)
  1108.                 {
  1109.                     REAL snr_1 = pSNR[SBLIMIT];
  1110.                     REAL snr_0 = *pSNR++;
  1111.  
  1112.                     if (snr_1 > snr_0) snr_0 = snr_1;
  1113.  
  1114.                 #ifndef USE_OLD_MATHLIB
  1115.                     *pLTMin++ = (LONG) ( (LOG10_TO_LN * 1000) * log(snr_0) );
  1116.                 #else
  1117.                     *pLTMin++ = (LONG) floor( (LOG10_TO_LN * 1000) * log(snr_0) );
  1118.                 #endif
  1119.                 }
  1120.             }
  1121.             else TranslateThreshold(pSNR);
  1122.         }
  1123.     }
  1124. }
  1125. ///
  1126.  
  1127. #if     _USE_ASM_FCT_ == 1
  1128. #if     _CPUMODEL_ == 68020
  1129. #error  "Can't use ASM optimized version (No FPU available)!"
  1130. #endif
  1131. #endif
  1132.  
  1133.