#include "fftservintf.h" namespace SOPHYA { //// VOIR GRAND BLABLA EXPLICATIF A LA FIN DU FICHIER /*! \class FFTServerInterface \ingroup NTools Defines the interface for FFT (Fast Fourier Transform) operations. Definitions : - Sampling period \b T - Sampling frequency \b fs=1/T - Total number of samples \b N - Frequency step in Fourier space \b =fs/N=1/(N*T) - Component frequencies - k=0 -> 0 - k=1 -> 1/(N*T) - k -> k/(N*T) - k=N/2 -> 1/(2*T) (Nyquist frequency) - k>N/2 -> k/(N*T) (or negative frequency -(N-k)/(N*T)) For a sampling period T=1, the computed Fourier components correspond to : \verbatim 0 1/N 2/N ... 1/2 1/2+1/N 1/2+2/N ... 1-2/N 1-1/N 0 1/N 2/N ... 1/2 ... -2/N -1/N \endverbatim For complex one-dimensional transforms: \f[ out(i) = F_{norm} \Sigma_{j} \ e^{-2 \pi \sqrt{-1} \ i \ j} \ {\rm (forward)} \f] \f[ out(i) = F_{norm} \Sigma_{j} \ e^{2 \pi \sqrt{-1} \ i \ j} \ {\rm (backward)} \f] i,j= 0..N-1 , where N is the input or the output array size. For complex multi-dimensional transforms: \f[ out(i1,i2,...,id) = F_{norm} \Sigma_{j1} \Sigma_{j2} ... \Sigma_{jd} \ e^{-2 \pi \sqrt{-1} \ i1 \ j1} ... e^{-2 \pi \sqrt{-1} \ id \ jd} \ {\rm (forward)} \f] \f[ out(i1,i2,...,id) = F_{norm} \Sigma_{j1} \Sigma_{j2} ... \Sigma_{jd} \ e^{2 \pi \sqrt{-1} \ i1 \ j1} ... e^{2 \pi \sqrt{-1} \ id \ jd} \ {\rm (backward)} \f] For real forward transforms, the input array is real, and the output array complex, with Fourier components up to k=N/2. For real backward transforms, the input array is complex and the output array is real. */ /* --Methode-- */ FFTServerInterface::FFTServerInterface(string info) { _info = info; _fgnorm = true; } /* --Methode-- */ FFTServerInterface::~FFTServerInterface() { } // ----------------- Transforme pour les double ------------------- /* --Methode-- */ //! Forward Fourier transform for double precision complex data /*! \param in : Input complex array \param out : Output complex array */ void FFTServerInterface::FFTForward(TArray< complex > &, TArray< complex > &) { throw NotAvailableOperation("FFTServer::FFTForward(TArray...) Unsupported operation !"); } /* --Methode-- */ //! Backward (inverse) Fourier transform for double precision complex data /*! \param in : Input complex array \param out : Output complex array */ void FFTServerInterface::FFTBackward(TArray< complex > &, TArray< complex > &) { throw NotAvailableOperation("FFTServer::FFTBackward(TArray...) Unsupported operation !"); } /* --Methode-- */ //! Forward Fourier transform for double precision real input data /*! \param in : Input real array \param out : Output complex array */ void FFTServerInterface::FFTForward(TArray< r_8 > &, TArray< complex > &) { throw NotAvailableOperation("FFTServer::FFTForward(TArray...) Unsupported operation !"); } /* --Methode-- */ //! Backward (inverse) Fourier transform for double precision real output data /*! \param in : Input complex array \param out : Output real array \param usoutsz : if true, use the output array size for computing the inverse FFT. In all cases, the input/output array sizes compatibility is checked. if usoutsz == false, the size of the real array is selected based on the the imaginary part of the input complex array at the nyquist frequency. size_out_real = 2*size_in_complex - ( 1 or 2) */ void FFTServerInterface::FFTBackward(TArray< complex > &, TArray< r_8 > &, bool) { throw NotAvailableOperation("FFTServer::FFTBackward(TArray...) Unsupported operation !"); } // ----------------- Transforme pour les float ------------------- /* --Methode-- */ //! Forward Fourier transform for complex data /*! \param in : Input complex array \param out : Output complex array */ void FFTServerInterface::FFTForward(TArray< complex > &, TArray< complex > &) { throw NotAvailableOperation("FFTServer::FFTForward(TArray r_4 ... ) Unsupported operation !"); } /* --Methode-- */ //! Backward (inverse) Fourier transform for complex data /*! \param in : Input complex array \param out : Output complex array */ void FFTServerInterface::FFTBackward(TArray< complex > &, TArray< complex > &) { throw NotAvailableOperation("FFTServer::FFTBackward(TArray r_4 ... ) Unsupported operation !"); } /* --Methode-- */ //! Forward Fourier transform for real input data /*! \param in : Input real array \param out : Output complex array */ void FFTServerInterface::FFTForward(TArray< r_4 > &, TArray< complex > &) { throw NotAvailableOperation("FFTServer::FFTForward(TArray r_4 ... ) Unsupported operation !"); } /* --Methode-- */ //! Backward (inverse) Fourier transform for real output data /*! \param in : Input complex array \param out : Output real array \param usoutsz : if true, use the output array size for computing the inverse FFT. In all cases, the input/output array sizes compatibility is checked. if usoutsz == false, the size of the real array is selected based on the the imaginary part of the input complex array at the nyquist frequency. size_out_real = 2*size_in_complex - ( 1 or 2) */ void FFTServerInterface::FFTBackward(TArray< complex > &, TArray< r_4 > &, bool) { throw NotAvailableOperation("FFTServer::FFTBackward(TArray r_4 ... ) Unsupported operation !"); } /* --Methode-- */ /*! \class FFTArrayChecker \ingroup NTools Service class for checking array size and resizing output arrays, to be used by FFTServer classes */ template FFTArrayChecker::FFTArrayChecker(string msg, bool checkpack, bool onedonly) { _msg = msg + " FFTArrayChecker::"; _checkpack = checkpack; _onedonly = onedonly; } /* --Methode-- */ template FFTArrayChecker::~FFTArrayChecker() { } template T FFTArrayChecker::ZeroThreshold() { return(0); } DECL_TEMP_SPEC /* equivalent a template <> , pour SGI-CC en particulier */ r_8 FFTArrayChecker< r_8 >::ZeroThreshold() { return(1.e-39); } DECL_TEMP_SPEC /* equivalent a template <> , pour SGI-CC en particulier */ r_4 FFTArrayChecker< r_4 >::ZeroThreshold() { return(1.e-19); } /* --Methode-- */ template int FFTArrayChecker::CheckResize(TArray< complex > const & in, TArray< complex > & out) { int k; string msg; if (in.Size() < 1) { msg = _msg + "CheckResize(complex in, complex out) - Unallocated input array !"; throw(SzMismatchError(msg)); } if (_checkpack) if ( !in.IsPacked() ) { msg = _msg + "CheckResize(complex in, complex out) - Not packed input array !"; throw(SzMismatchError(msg)); } int ndg1 = 0; for(k=0; k 1) ndg1++; if (_onedonly) if (ndg1 > 1) { msg = _msg + "CheckResize(complex in, complex out) - Only 1-D array accepted !"; throw(SzMismatchError(msg)); } out.ReSize(in); // sa_size_t sz[BASEARRAY_MAXNDIMS]; // for(k=0; k int FFTArrayChecker::CheckResize(TArray< T > const & in, TArray< complex > & out) { int k; string msg; if (in.Size() < 1) { msg = _msg + "CheckResize(real in, complex out) - Unallocated input array !"; throw(SzMismatchError(msg)); } if (_checkpack) if ( !in.IsPacked() ) { msg = _msg + "CheckResize(real in, complex out) - Not packed input array !"; throw(SzMismatchError(msg)); } int ndg1 = 0; for(k=0; k 1) ndg1++; if (_onedonly) if (ndg1 > 1) { msg = _msg + "CheckResize(real in, complex out) - Only 1-D array accepted !"; throw(SzMismatchError(msg)); } sa_size_t sz[BASEARRAY_MAXNDIMS]; // if (ndg1 > 1) { sz[0] = in.Size(0)/2+1; for(k=1; k int FFTArrayChecker::CheckResize(TArray< complex > const & in, TArray< T > & out, bool usoutsz) { int k; string msg; if (in.Size() < 1) { msg = _msg + "CheckResize(complex in, real out) - Unallocated input array !"; throw(SzMismatchError(msg)); } if (_checkpack) if ( !in.IsPacked() ) { msg = _msg + "CheckResize(complex in, real out) - Not packed input array !"; throw(SzMismatchError(msg)); } int ndg1 = 0; for(k=0; k 1) ndg1++; if (_onedonly) if (ndg1 > 1) { msg = _msg + "CheckResize(complex in, real out) - Only 1-D array accepted !"; throw(SzMismatchError(msg)); } if (usoutsz) { // We have to use output array size bool fgerr = false; if (ndg1 > 1) { if (in.Size(0) != out.Size(0)/2+1) fgerr = true; } else { if (in.Size(in.MaxSizeKA()) != out.Size(in.MaxSizeKA())/2+1) fgerr = true; } if (fgerr) { msg = _msg + "CheckResize(complex in, real out) - Incompatible in-out sizes !"; throw(SzMismatchError(msg)); } } else { // We have to resize the output array T thr = ZeroThreshold(); // Seuil pour tester Imag(Nyquist) == 0 sa_size_t sz[BASEARRAY_MAXNDIMS]; if (ndg1 > 1) { T imnyq = in(in.Size(0)-1,0,0).imag(); // Rz+cmc/Nov07 : // Calcul de la taille SizeX/Sz[0] paire/impaire en fonction de Imag(Nyquist) sz[0] = ((imnyq < -thr)||(imnyq > thr)) ? 2*in.Size(0)-1 : 2*in.Size(0)-2; if (sz[0] < 1) sz[0] = 1; for(k=1; k thr) ) ? 2*n-1 : 2*n-2; if (ncs < 1) ncs = 1; sz[in.MaxSizeKA()] = ncs; } out.ReSize(in.NbDimensions(), sz); } return(ndg1); } #ifdef __CXX_PRAGMA_TEMPLATES__ #pragma define_template FFTArrayChecker #pragma define_template FFTArrayChecker #endif #if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES) template class FFTArrayChecker; template class FFTArrayChecker; #endif } // FIN namespace SOPHYA /********************************************************************** Memo uniquement destine aux programmeurs: (cmv 03/05/04) -- cf programme de tests explicatif: cmvtfft.cc ===================================================================== ===================================================================== ============== Transformees de Fourier et setNormalize ============== ===================================================================== ===================================================================== - si setNormalize(true): invTF{TF{S}} = S - si setNormalize(false): invTF{TF{S}} = N * S ===================================================================== ===================================================================== ============== Transformees de Fourier de signaux REELS ============= ===================================================================== ===================================================================== ------- --- FFT d'un signal REEL S ayant un nombre pair d'elements N=2p ------- taille de la FFT: Nfft = N/2 + 1 = p + 1 abscisses de la fft: | 0 | 1/N | 2/N | ..... | p/N=1/2 | ^continu ^frequence de Nyquist ... Ex: N=6 -> Nfft = 6/3+1 = 4 le signal a N elements reels, la fft a Nfft elements complexes cad 2*Nfft reels = 2*(p+1) reels = 2p + 2 reels = N + 2 reels soit 2 reels en trop: ce sont les phases du continu et de la frequence de Nyquist relations: - si setNormalize(true) : fac = N setNormalize(false) : fac = 1/N sum(i=0,N-1){S(i)^2} = fac* [[ 2* sum(j=0,Nfft-1){|TF{S}(j)|^2} - |TF{S}(0)|^2 - |TF{S}(Nfft-1)|^2 ]] (On ne compte pas deux fois le continu et la freq de Nyquist) ------- --- FFT d'un signal REEL ayant un nombre impair d'elements N=2p+1 ------- taille de la FFT: Nfft = N/2 + 1 = p + 1 abscisses de la fft: | 0 | 1/N | 2/N | ..... | p/N | ^continu (la frequence de Nyquist n'y est pas) ... Ex: N=7 -> Nfft = 7/3+1 = 4 le signal a N elements reels, la fft a Nfft elements complexes cad 2*Nfft reels = 2*(p+1) reels = 2p + 2 reels = N + 1 reels soit 1 reel en trop: c'est la phase du continu relations: - si setNormalize(true) : fac = N setNormalize(false) : fac = 1/N sum(i=0,N-1){S(i)^2} = fac* [[ 2* sum(j=0,Nfft-1){|TF{S}(j)|^2} - |TF{S}(0)|^2 ]] (On ne compte pas deux fois le continu) ------------ --- FFT-BACK d'un signal F=TF{S} ayant un nombre d'elements Nfft ------------ Sback = invTF{TF{S}} Remarque: Nfft a la meme valeur pour N=2p et N=2p+1 donc Nfft conduit a 2 possibilites: { N = 2*(Nfft-1) signal back avec nombre pair d'elements { N = 2*Nfft-1 signal back avec nombre impair d'elements Pour savoir quel est la longueur N du signal TF^(-1){F} on regarde si F(Nfft-1) est reel ou complexe (la frequence de Nyquist d'un signal reel est reelle) - Si F(Nfft-1) reel cad Im{F(Nfft-1)}=0: N = 2*(Nfft-1) - Si F(Nfft-1) complexe cad Im{F(Nfft-1)}#0: N = 2*Nfft-1 Si setNormalize(true): invTF{TF{S}} = S setNormalize(false): invTF{TF{S}} = N * S ========================================================================= ========================================================================= ============== Transformees de Fourier de signaux COMPLEXES ============= ========================================================================= ========================================================================= ------- --- FFT d'un signal COMPLEXE S ayant un nombre d'elements N ------- taille de la FFT: Nfft = N abscisses de la fft: | 0 | 1/N | 2/N | ..... | (N-1)/N | ^continu Frequence de Nyquist: si N est pair: la frequence de Nyquist est l'absicce d'un des bins abscisses de TF{S}: Nfft = N = 2p | 0 | 1/N | 2/N | ... | (N/2)/N=p/N=0.5 | ... | (N-1)/N | ^frequence de Nyquist si N est impair: la frequence de Nyquist N'est PAS l'absicce d'un des bins abscisses de TF{S}: Nfft = N = 2p+1 | 0 | 1/N | 2/N | ... | (N/2)/N=p/N | ((N+1)/2)/N=(p+1)/N | ... | (N-1)/N | ... Ex: N = 2p =6 -> Nfft = 2p = 6 abscisses de TF{S}: | 0 | 1/6 | 2/6 | 3/6=0.5 | 4/6 | 5/6 | ... Ex: N = 2p+1 = 7 -> Nfft = 2p+1 = 7 abscisses de TF{S}: | 0 | 1/7 | 2/7 | 3/7 | 4/7 | 5/7 | 6/7 | relations: - si setNormalize(true) : fac = N setNormalize(false) : fac = 1/N sum(i=0,N-1){S(i)^2} = fac* [[ sum(j=0,Nfft-1){|TF{S}(j)|^2} ]] ------------ --- FFT-BACK d'un signal F=TF{S} ayant un nombre d'elements Nfft ------------ taille du signal: N = Nfft Si setNormalize(true): invTF{TF{S}} = S setNormalize(false): invTF{TF{S}} = N * S **********************************************************************/