/* ** ** FastFFT.cpp ** ** $VER: FastFFT.cpp 1.0 (11.07.98) ** ** Calcul de la transformée de Fourier rapide ** ** $Revision: 1.4 $ ** $State: Exp $ ** $Date: 1998/08/16 19:03:44 $ ** ** $Log: FastFFTC.cpp $ ** Revision 1.4 1998/08/16 19:03:44 kakace ** Version Beta3+ ** ** Revision 1.3 1998/08/09 19:21:36 kakace ** Correction d'un bug dans la routine FFT (avec tableà ** ** Revision 1.2 1998/08/02 15:17:10 kakace ** Fonctionnement OK (Layer 1/2 + Stereo/JStereo) ** */ #ifndef _INCLUDE_MATH_H #include #endif #ifndef _MEMORY_HPP #include "Memory.hpp" #endif #ifndef _FASTFFT_CLASS_HPP #include "FastFFTC.hpp" #endif #ifndef _PEGASECOND_HPP #include "PegaseCond.hpp" #endif //---------------------------------------------------------------------------------------------------- // Fonctions "inline" /// FastFFTC::GetShuffleIndex_?() /****** Class FastFFTC/GetShuffleIndex_10 *********************************** * * NAME * FastFFTC::GetShuffleIndex_10 -- Obtenir un index sur 10 bits. * FastFFTC::GetShuffleIndex_8 -- Obtenir un index sur 8 bits. * * SYNOPSIS * Index = FastFFTC::GetShuffleIndex_10(input) * Index = FastFFTC::GetShuffleIndex_8(input) * * inline ULONG FastFFTC::GetShuffleIndex_10(ULONG); * inline ULONG FastFFTC::GetShuffleIndex_8(ULONG); * * FUNCTION * Retourne la valeur transmise en argument après avoir inversé le champs * binaire (bitreverse). La longueur du champs de bits doit correspondre à * longueur indiquée dans le nom de la fonction. * * INPUTS * input - Valeur à traiter. * * RESULT * Index - Index obtenu après inversion du masque binaire. * ***************************************************************************** * */ static UBYTE ShuffleArray[] = {0, 8, 4, 12, 2, 10, 6, 14, 1, 9, 5, 13, 3, 11, 7, 15}; inline ULONG FastFFTC::GetShuffleIndex_10(ULONG input) { ULONG k; ULONG bitrev = 0; for (k = 0; k < 3; k++) { bitrev = (bitrev * 16) | (::ShuffleArray[input & 15]); input >>= 4; } return (bitrev >> 2); } inline ULONG FastFFTC::GetShuffleIndex_8(ULONG input) { return ( (::ShuffleArray[input & 15] * 16) | (::ShuffleArray[input >> 4]) ); } /// // Constructeur et destructeur. /// FastFFTC::FastFFTC FastFFTC::FastFFTC() { int i; #if _USE_SINCOS_TABLE_ == 1 // Initialiser la table SinCos (68020/68030) REAL cos_alpha; REAL f = - (PI / 512.0); for (i = 0; i < SINCOSSIZE / 3; i++) { cos_alpha = cos(f * i); FFT_dSinCos[i] = cos_alpha; FFT_dSinCos[SINCOSSIZE * 2 / 3 - i] = - cos_alpha; FFT_dSinCos[SINCOSSIZE * 2 / 3 + i] = - cos_alpha; } FFT_dSinCos[SINCOSSIZE / 3] = 0; #else // Initialiser la table SinCos réduite (68040/68060) REAL f = - (PI / 512.0); for (i = 9; i >= 0; i--) { FFT_dSinCos[0][i] = cos(f); FFT_dSinCos[1][i] = sin(f); f += f; } #endif FFT_MemoryError = FALSE; FFT_MemoryError |= !(FFT_pdWinSmplShuffled = (T_FFTBUFF *) AllocAlignedBuffer(sizeof(T_FFTBUFF))); FFT_MemoryError |= !(FFT_pdHannWindowShuffled = (T_HANN *) AllocAlignedBuffer(sizeof(T_HANN))); // Initialiser la table HannWindow. // Cette table est réordonnée pour faciliter l'initialisation des tampons d'entrée // des FFTs. { REAL cos_alpha, sin_alpha; REAL cos_delta, sin_delta; REAL f; f = - (PI / 1024.0); cos_alpha = cos(f); sin_alpha = sin(f); f = cos_alpha + cos_alpha; // = 2.cos(a) cos_delta = f * cos_alpha - 1.0; // cos(2a) = 2.cos²(a) - 1 sin_delta = - f * sin_alpha; // sin(2a) = 2.sin(a).cos(a) for (i = 0; i < BLKSIZE; i++) { f = 0.5 * (1.0 - cos_alpha); (*FFT_pdHannWindowShuffled)[GetShuffleIndex_10(i)] = f; f = cos_alpha; cos_alpha = f * cos_delta - sin_alpha * sin_delta; sin_alpha = f * sin_delta + sin_alpha * cos_delta; } } } /// /// FastFFTC::~FastFFTC() FastFFTC::~FastFFTC() { if (FFT_pdWinSmplShuffled) FreeAlignedBuffer(FFT_pdWinSmplShuffled); if (FFT_pdHannWindowShuffled) FreeAlignedBuffer(FFT_pdHannWindowShuffled); } /// // Fonctions privées. /// FastFFTC::FastFFT() /****** Class FastFFTC/FastFFT ********************************************** * * NAME * FastFFTC::FastFFT -- Calculer la transformée de Fourier. * * SYNOPSIS * FastFFTC::FastFFT(pdFFTIn, iFFTSize) * * void FastFFTC::FastFFT(REAL *, UWORD); * * FUNCTION * Calcule la transformée de Fourier sur un ensemble d'échantillons. * * INPUTS * pdFFTIn - Tableau contenant les échantillons (série de nombres complexes * dans l'ordre : [partie réelle][partie imaginaire]. * iFFTSize - Nombre d'échantillons à traiter. * * NOTES * Cette fonction a été fortement optimisée en tirant partie de quelques * angles particuliers. Elle est capable de calculer une FFT sur un ensemble * de 1024 échantillons au maximum. * * BUGS * Fonction contrôlée conforme le 20.06.98 * * SEE ALSO * FastFFTC::ComputeFastFFT_Stereo() * ***************************************************************************** * * Benchmark (3148 cycles par seconde (global)) : * 21.06.98 : 1'45.53 * * 06.07.98 : * - Utilisation d'une table SinCos[] en lieu et place des calculs * trigonométriques (si _USE_SINCOS_TABLE_ == 1). */ #if _USE_ASM_FCT_ == 0 void FastFFTC::FastFFT(T_FFTDATAS *pdFFTIn, UWORD iFFTSize) { #define MAX_FFT_BLOCK 256 ULONG i, j, delta; LONG k; ULONG BlockSize, BlockEnd, n; ULONG iBandSize = (iFFTSize < MAX_FFT_BLOCK) ? iFFTSize : MAX_FFT_BLOCK; double tr, ti; T_FFTDATAS *pd_J, *pd_K, *pd_J2, *pd_K2; #if _USE_SINCOS_TABLE_ == 1 ULONG step; float *Cosin = &FFT_dSinCos[0]; float *Sinus = &FFT_dSinCos[SINCOSSIZE / 3]; #endif // Appliquer la FFT sur le tampon d'entrée. // Les passes 0 à 7 sont calculées ici, par colonnes de "iBandSize" échantillons // (pour tirer un meilleur parti des caches). // Les passes 8 et 9 sont quant à elles calculées à part. // On commence par la fin pour profiter du fait que le tampon d'entrée vient d'être // initialisé, et donc que les derniers échantillons sont toujours dans le cache. for (k = iFFTSize - iBandSize; k >= 0; k -= iBandSize) { BlockEnd = 1; #if _USE_SINCOS_TABLE_ == 0 delta = 0; #else step = 512; #endif for (BlockSize = 2; BlockSize <= iBandSize; BlockSize += BlockSize) { #if _USE_SINCOS_TABLE_ == 1 delta = 0; #endif n = BlockSize + BlockSize; pd_J = pd_J2 = &pdFFTIn[k + k]; pd_K = pd_K2 = &pd_J2[BlockSize]; // Première partie : ar = 1, ai = 0 (commune à toutes les passes) for (i = 0; i < iBandSize; i += BlockSize) // Boucle exécutée 1020 fois { tr = pd_K[0]; ti = pd_J[0]; pd_K[0] = ti - tr; pd_J[0] = ti + tr; tr = pd_J[1]; ti = pd_K[1]; pd_K[1] = tr - ti; pd_J[1] = tr + ti; pd_K = &pd_K[n]; // <==> pd_J += n pd_J = &pd_J[n]; // <==> pd_J += n } pd_J2 += 2; pd_K2 += 2; // Deuxième partie : ar et ai quelconques. if (BlockEnd > 4) { double ar, ai; #if _USE_SINCOS_TABLE_ == 0 double cos_alpha, sin_alpha; ar = cos_alpha = FFT_dSinCos[0][delta]; // = cos(2.PI/BlockSize) ai = sin_alpha = FFT_dSinCos[1][delta]; // = sin(2.PI/BlockSize) #endif for (j = 1; j < BlockEnd; j++) { pd_J = pd_J2; pd_K = pd_K2; #if _USE_SINCOS_TABLE_ == 1 delta += step; ar = Cosin[delta]; ai = Sinus[delta]; #endif for (i = 0; i < iBandSize; i += BlockSize) { ti = pd_K[0]; tr = ar * ti - ai * pd_K[1]; ti = ar * pd_K[1] + ai * ti; pd_K[0] = pd_J[0] - tr; pd_J[0] = pd_J[0] + tr; tr = pd_J[1]; pd_K[1] = tr - ti; pd_J[1] = tr + ti; pd_K = &pd_K[n]; // <==> pd_J += n pd_J = &pd_J[n]; // <==> pd_J += n } pd_J2 += 2; pd_K2 += 2; #if _USE_SINCOS_TABLE_ == 0 tr = ar; ar = ar * cos_alpha - ai * sin_alpha; ai = tr * sin_alpha + ai * cos_alpha; #endif } } else { // Troisième partie : ar = 0, ai = -1. pd_J = pd_J2; pd_K = pd_K2; if (BlockEnd == 2) { for (i = 0; i < iBandSize; i += 4) { double tz; ti = pd_K[1]; // ai = -1 tr = pd_K[0]; // ar = 0 tz = pd_J[0]; pd_K[0] = tz - ti; pd_J[0] = tz + ti; tz = pd_J[1]; pd_K[1] = tz + tr; pd_J[1] = tz - tr; pd_K += 8; pd_J += 8; } } // Quatrième partie : ar = {1; 0.707; 0; -0.707}, ai = {0; -0.707; -1; -0.707} if (BlockEnd == 4) { double ar, tz; #if _USE_SINCOS_TABLE_ == 0 ar = FFT_dSinCos[0][delta]; // = cos(2.PI/8) = sqrt(2) / 2 #else ar = Cosin[SINCOSSIZE / 6]; #endif for (i = 0; i < iBandSize; i += 8) { ti = pd_K[0]; tz = pd_K[1]; tr = ar * (ti + tz); // ar = -ai = sqrt(2) / 2 ti = ar * (tz - ti); tz = pd_J[0]; pd_K[0] = tz - tr; pd_J[0] = tz + tr; tz = pd_J[1]; pd_K[1] = tz - ti; pd_J[1] = tz + ti; tr = pd_K[2]; // ar = 0 ti = pd_K[3]; // ai = -1 tz = pd_J[2]; pd_K[2] = tz - ti; pd_J[2] = tz + ti; tz = pd_J[3]; pd_K[3] = tz + tr; pd_J[3] = tz - tr; ti = pd_K[4]; tz = pd_K[5]; tr = ar * (ti - tz); // ar = ai = - sqrt(2) / 2 ti = ar * (tz + ti); tz = pd_J[4]; pd_K[4] = tz + tr; pd_J[4] = tz - tr; tz = pd_J[5]; pd_K[5] = tz + ti; pd_J[5] = tz - ti; pd_K += 16; pd_J += 16; } } } BlockEnd = BlockSize; #if _USE_SINCOS_TABLE_ == 0 ++delta; #else step >>= 1; #endif } } // ================================================================ // Dernière passe : BlockSize = 512, BlockEnd = 256 (passes 8 et 9) // ================================================================ for ( ; BlockSize <= iFFTSize; BlockSize += BlockSize) { #if _USE_SINCOS_TABLE_ == 0 double cos_alpha, sin_alpha, ar, ai; cos_alpha = FFT_dSinCos[0][delta]; sin_alpha = FFT_dSinCos[1][delta]; #else double ar, ai; delta = 0; #endif pd_J = &pdFFTIn[0]; pd_K = &pdFFTIn[BlockSize]; ar = 1.0; ai = 0.0; for (j = 0; j < BlockEnd; j++) { if (BlockEnd == 256 && iFFTSize == 1024) { ti = pd_K[1024]; tr = ar * ti - ai * pd_K[1025]; ti = ar * pd_K[1025] + ai * ti; pd_K[1025] = pd_J[1025] - ti; pd_J[1025] = pd_J[1025] + ti; ti = pd_J[1024]; pd_K[1024] = ti - tr; pd_J[1024] = ti + tr; } ti = pd_K[0]; tr = ar * ti - ai * pd_K[1]; ti = ar * pd_K[1] + ai * ti; pd_K[0] = pd_J[0] - tr; pd_J[0] = pd_J[0] + tr; tr = pd_J[1]; pd_K[1] = tr - ti; pd_J[1] = tr + ti; pd_J += 2; pd_K += 2; #if _USE_SINCOS_TABLE_ == 0 tr = ar; ar = ar * cos_alpha - ai * sin_alpha; ai = tr * sin_alpha + ai * cos_alpha; #else delta += step; ar = Cosin[delta]; ai = Sinus[delta]; #endif } BlockEnd = BlockSize; #if _USE_SINCOS_TABLE_ == 0 ++delta; #else step >>= 1; #endif } } #endif // _USE_ASM_FCT_ /// /// FastFFTC::FixFFT_Stereo() /****** Class FastFFTC/FixFFT_Stereo **************************************** * * NAME * FastFFTC::FixFFT_Stereo -- Calcule une FFT. * * SYNOPSIS * FastFFTC::FixFFT_Stereo(pd_FFTIn, iFFTSize) * * void FastFFTC::FixFFT_Stereo(REAL *, UWORD); * * FUNCTION * Calcule la FFT d'un ensemble d'échantillons stéréo. * * INPUTS * pd_FFTIn - Tableau de sortie de la FFT. * iFFTSize - Nombre d'échantillons à traiter. * * NOTES * La routine FFT effectue son calcul sur les deux canaux simultanément (les * échantillons ayant été placés dans l'ordre correct auparavant). * Dès lors, il suffit de réordonner les valeurs fournies par la FFT en * profitant des symétries pour obtenir la véritable FFT de chacun des deux * canaux. * * BUGS * Fonction contrôlée conforme le 20.06.98 * * SEE ALSO * FastFFTC::FastFFT() * ***************************************************************************** * * Benchmark (3 cycles par seconde) : * 21.06.98 : 0'19.61 */ #if _USE_ASM_FCT_ == 0 void FastFFTC::FixFFT_Stereo(T_FFTDATAS *pd_FFTIn, UWORD iFFTSize) { T_FFTDATAS *pd_FFTOut = &pd_FFTIn[BLKSIZE * 2]; // &(*FFT_pdWinSmplShuffled)[1][0][0]; T_FFTDATAS *pd_FFTIn2 = &pd_FFTIn[iFFTSize * 2]; // &(*FFT_pdWinSmplShuffled)[0][iFFTSize][0]; double dOneHalf = 0.5; ULONG j; ULONG iFFTLen = iFFTSize >> 1; // Décoder les données en sortie de FFT pour obtenir la table finale. pd_FFTOut[0] = pd_FFTIn[1]; pd_FFTIn[1] = pd_FFTOut[1] = 0.0; pd_FFTIn += 2; pd_FFTOut += 2; for (j = iFFTLen; j > 0; j--) { pd_FFTIn2 += (-2); *pd_FFTOut++ = dOneHalf * (pd_FFTIn[1] + pd_FFTIn2[1]); *pd_FFTOut++ = dOneHalf * (pd_FFTIn2[0] - pd_FFTIn[0]); *pd_FFTIn++ = dOneHalf * (*pd_FFTIn + pd_FFTIn2[0]); *pd_FFTIn++ = dOneHalf * (*pd_FFTIn - pd_FFTIn2[1]); } } #endif // _USE_ASM_FCT_ /// /// FastFFTC::FixFFT_Mono() /****** Class FastFFTC/FixFFT_Mono ****************************************** * * NAME * FastFFTC::FixFFT_Mono -- Calcule une FFT * * SYNOPSIS * FastFFTC::FixFFT_Mono(pd_FFTIn, iFFTSize) * * void FastFFTC::FixFFT_Mono(REAL *, UWORD); * * FUNCTION * Calcule une FFT sur un groupe d'échantillons (fichier mono). * * INPUTS * pd_FFTIn - Tableau de sortie de la FFT. * iFFTSize - Taille de la fenêtre FFT (est égale à la moitié du nombre d' * échantillons à analyser). * * NOTES * Les échantillons sont placés dans un ordre particulier afin de limiter * le nombre de calculs (échantillons d'index pair comme partie réelle, * échantillons d'index impair comme partie imaginaire). * Cette solution permet de limiter la taille de la fenêtre FFT à la moitié * de sa taille normale. * Il suffit de réordonner les nombres complexes en sortie de FFT pour * obtenir le résultat recherché. * * BUGS * Fonction contrôlée conforme le 05/07/98 * * SEE ALSO * FastFFTC::FastFFT() * ***************************************************************************** * */ #if _USE_ASM_FCT_ == 0 void FastFFTC::FixFFT_Mono(T_FFTDATAS *pd_FFTIn, UWORD iFFTSize) { T_FFTDATAS *pd_FFTIn2 = &pd_FFTIn[iFFTSize * 2]; // &(*FFT_pdWinSmplShuffled)[0][iFFTSize][0]; T_FFTDATAS *pd_FFTOrg = pd_FFTIn; ULONG i; ULONG iFFTLen = iFFTSize >> 1; double wr, wi; #if _USE_SINCOS_TABLE_ == 0 double wpr, wpi; #else ULONG delta = 0; float *Cosin = &FFT_dSinCos[0]; float *Sinus = &FFT_dSinCos[SINCOSSIZE / 3]; #endif double h1r, h1i, h2r, h2i; double dOneHalf = 0.5; // Décoder les données en sortie de FFT pour obtenir la table finale. #if _USE_SINCOS_TABLE_ == 0 wpr = FFT_dSinCos[0][9]; // cos(PI/iFFTSize); wpi = FFT_dSinCos[1][9]; // sin(PI/iFFTSize); wr = wpr; // = cos(theta) wi = wpi; // = sin(theta) #endif pd_FFTIn += 2; for (i = 1; i <= iFFTLen; i++) { #if _USE_SINCOS_TABLE_ == 1 ++delta; wr = Cosin[delta]; wi = Sinus[delta]; #endif pd_FFTIn2 += (-2); h2r = dOneHalf * (pd_FFTIn[1] + pd_FFTIn2[1]); h2i = dOneHalf * (pd_FFTIn2[0] - pd_FFTIn[0]); h1r = h2r; h2r = wr * h1r - wi * h2i; h2i = wr * h2i + wi * h1r; h1r = dOneHalf * (pd_FFTIn[0] + pd_FFTIn2[0]); h1i = dOneHalf * (pd_FFTIn[1] - pd_FFTIn2[1]); *pd_FFTIn++ = h2r + h1r; // + wr * h2r - wi * h2i; *pd_FFTIn++ = h2i + h1i; // + wr * h2i + wi * h2r; pd_FFTIn2[0] = h1r - h2r; // - wr * h2r + wi * h2i; pd_FFTIn2[1] = h2i - h1i; // + wr * h2i + wi * h2r; #if _USE_SINCOS_TABLE_ == 0 h1r = wr; wr = h1r * wpr - wi * wpi; // = cos(x).cos(theta) - sin(x).sin(theta) = cos(x+theta) wi = wi * wpr + h1r * wpi; // = sin(x).cos(theta) + cos(x).sin(theta) = sin(x+theta) #endif } iFFTSize *= 2; pd_FFTOrg[iFFTSize] = pd_FFTOrg[0] - pd_FFTOrg[1]; pd_FFTOrg[0] += pd_FFTOrg[1]; ++pd_FFTOrg; pd_FFTOrg[0] = pd_FFTOrg[iFFTSize] = 0.0; } #endif // _USE_ASM_FCT_ /// /// FastFFTC::SetupWinSamples_1024() /****** Class FastFFTC/SetupWinSamples_1024 ********************************* * * NAME * FastFFTC::SetupWinSamples_1024 -- Préparer les données pour la FFT. * * SYNOPSIS * FastFFTC::SetupWinSamples_1024(pSamples, wNumChannels) * * void FastFFTC::SetupWinSamples_1024(WORD *, UWORD); * * FUNCTION * Initialise le tampon d'entrée de la FFT à partir d'un bloc contenant 1024 * échantillons. * Lorsque le fichier audio traité est stéréo, les deux canaux sont préparés * simultanément. * * INPUTS * pSamples - Pointeur vers les échantillons. * wNumChannels - Nombre de canaux. * * NOTES * Afin de tirer partie des symétries de la FFT, les deux canaux sont * préparés de sorte que l'un corresponde à la partie réelle, et l'autre * à la partie imaginaire. De fait, la FFT effectue son calcul sur les deux * canaux simultanément. * * Pour un fichier mono, les échantillons sont alternativement placés dans * la partie réelle (échantillons d'indices pairs) ou la partie imaginaire * (échantillons d'indices impairs) du tampon d'entrée. De cette façon, la * fonction FFT peut être fortement optimisée. * * BUGS * * SEE ALSO * FastFFTC::ComputeFastFFT_Stereo() * ***************************************************************************** * * Benchmarks (12 cycles par seconde) : * 21.06.98 : 0'25.42 * */ #if _USE_ASM_FCT_ == 0 void FastFFTC::SetupWinSamples_1024(WORD *pSamples, UWORD wNumChannels) { const T_FFTDATAS *pdHANNWindow = &(*FFT_pdHannWindowShuffled)[0]; T_FFTDATAS *pdWSmpls = &(*FFT_pdWinSmplShuffled)[0][0][0]; double f; WORD *pwSmpl; ULONG i, k; if (wNumChannels == 2) { for (i = 0; i < 1024 / 4; i++) { k = GetShuffleIndex_8(i) * 2; pwSmpl = &pSamples[k]; f = *pdHANNWindow++; *pdWSmpls++ = pwSmpl[0] * f; // Canal gauche. *pdWSmpls++ = pwSmpl[1] * f; // Canal droit. f = *pdHANNWindow++; *pdWSmpls++ = pwSmpl[512 * 2] * f; // Canal gauche. *pdWSmpls++ = pwSmpl[512 * 2 + 1] * f; // Canal droit. f = *pdHANNWindow++; *pdWSmpls++ = pwSmpl[256 * 2] * f; // Canal gauche. *pdWSmpls++ = pwSmpl[256 * 2 + 1] * f; // Canal droit. f = *pdHANNWindow++; *pdWSmpls++ = pwSmpl[768 * 2] * f; // Canal gauche. *pdWSmpls++ = pwSmpl[768 * 2 + 1] * f; // Canal droit. } } else { for (i = 0; i < 1024 / 4; i++) { k = GetShuffleIndex_8(i) * 2; pwSmpl = &pSamples[k]; *pdWSmpls++ = pwSmpl[0] * *pdHANNWindow++; *pdWSmpls++ = pwSmpl[1] * pdHANNWindow[511]; *pdWSmpls++ = pwSmpl[512] * *pdHANNWindow++; *pdWSmpls++ = pwSmpl[513] * pdHANNWindow[511]; } } } #endif // _USE_ASM_FCT_ ///