Changeset 2540 in Sophya for trunk/SophyaLib/NTools
- Timestamp:
- May 25, 2004, 9:16:58 AM (21 years ago)
- Location:
- trunk/SophyaLib/NTools
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/SophyaLib/NTools/fftservintf.cc
r2344 r2540 1 1 #include "fftservintf.h" 2 2 3 //// VOIR GRAND BLABLA EXPLICATIF A LA FIN DU FICHIER 3 4 4 5 /*! … … 156 157 } 157 158 158 159 // ----------------- Transformation de normalisation pour les energies -------------------160 161 /* --Methode-- */162 //! Compute the transform to be applied to "fftx=FFT(x)" so that "x" and "FFT(x)" have the same energy163 /*!164 \return The factor to be applied to the FFT energy such that we get the same energy as for the x.165 \verbatim166 fftx is computed by: FFTForward(x,fftx)167 Energy of x : Ex = sum{|x(i)|^2} i=0,x.Size()-1168 Energy of fftx : Efftx = sum{|fftx(i)^2|} i=0,fftx.Size()-1169 ( usually x.Size() != fftx.Size() )170 -------------------------------------------------------------------171 | TransfEnergyFFT return A and B such that : Ex = A * Efftx + B |172 | and Norm such that : Ex = Norm * Efftx |173 -------------------------------------------------------------------174 To normalize the fftx vector apply : "fftx *= sqrt(Norm)"175 \endverbatim176 */177 r_8 FFTServerInterface::TransfEnergyFFT178 (TVector< complex<r_8> > const& x, TVector< complex<r_8> > const& fftx, r_8& A, r_8& B)179 {180 B=0.;181 if(getNormalize()) A = x.Size(); else A = 1./x.Size();182 r_8 norm = A;183 return norm;184 }185 186 //! Compute the transform to be applied to "fftx=FFT(x)" so that "x" and "FFT(x)" have the same energy187 r_8 FFTServerInterface::TransfEnergyFFT188 (TVector< complex<r_4> > const & x, TVector< complex<r_4> > const & fftx, r_8& A, r_8& B)189 {190 B=0.;191 if(getNormalize()) A = x.Size(); else A = 1./x.Size();192 r_8 norm = A;193 return norm;194 }195 196 //! Compute the transform to be applied to "fftx=FFT(x)" so that "x" and "FFT(x)" have the same energy197 r_8 FFTServerInterface::TransfEnergyFFT198 (TVector< r_8 > const & x, TVector< complex<r_8> > const & fftx, r_8& A, r_8& B)199 {200 A= 2.;201 B= - abs(fftx(0)*fftx(0));202 if(x.Size()%2 == 0) B -= abs(fftx(fftx.Size()-1)*fftx(fftx.Size()-1));203 if(getNormalize()) {A *= x.Size(); B *= x.Size();}204 else {A /= x.Size(); B /= x.Size();}205 r_8 norm = 0.;206 for(int_4 i=0;i<fftx.Size();i++) norm += abs(fftx(i)*fftx(i));207 if(norm>0.) norm = (A*norm+B)/norm;208 return norm;209 }210 211 //! Compute the transform to be applied to "fftx=FFT(x)" so that "x" and "FFT(x)" have the same energy212 r_8 FFTServerInterface::TransfEnergyFFT213 (TVector< r_4 > const & x, TVector< complex<r_4> > const & fftx, r_8& A, r_8& B)214 {215 A= 2.;216 B= - abs(fftx(0)*fftx(0));217 if(x.Size()%2 == 0) B -= abs(fftx(fftx.Size()-1)*fftx(fftx.Size()-1));218 if(getNormalize()) {A *= x.Size(); B *= x.Size();}219 else {A /= x.Size(); B /= x.Size();}220 r_8 norm = 0.;221 for(int_4 i=0;i<fftx.Size();i++) norm += abs(fftx(i)*fftx(i));222 if(norm>0.) norm = (A*norm+B)/norm;223 return norm;224 }225 226 227 159 /* --Methode-- */ 228 160 /*! … … 408 340 template class FFTArrayChecker<r_8>; 409 341 #endif 342 343 344 345 346 347 /********************************************************************** 348 349 Memo uniquement destine aux programmeurs: (cmv 03/05/04) 350 -- cf programme de tests explicatif: cmvtfft.cc 351 352 ===================================================================== 353 ===================================================================== 354 ============== Transformees de Fourier et setNormalize ============== 355 ===================================================================== 356 ===================================================================== 357 358 - si setNormalize(true): invTF{TF{S}} = S 359 - si setNormalize(false): invTF{TF{S}} = N * S 360 361 ===================================================================== 362 ===================================================================== 363 ============== Transformees de Fourier de signaux REELS ============= 364 ===================================================================== 365 ===================================================================== 366 367 ------- 368 --- FFT d'un signal REEL S ayant un nombre pair d'elements N=2p 369 ------- 370 taille de la FFT: Nfft = N/2 + 1 = p + 1 371 abscisses de la fft: | 0 | 1/N | 2/N | ..... | p/N=1/2 | 372 ^continu ^frequence de Nyquist 373 374 ... Ex: N=6 -> Nfft = 6/3+1 = 4 375 376 le signal a N elements reels, la fft a Nfft elements complexes 377 cad 2*Nfft reels = 2*(p+1) reels = 2p + 2 reels = N + 2 reels 378 soit 2 reels en trop: 379 ce sont les phases du continu et de la frequence de Nyquist 380 381 relations: 382 - si setNormalize(true) : fac = N 383 setNormalize(false) : fac = 1/N 384 sum(i=0,N-1){S(i)^2} 385 = fac* [[ 2* sum(j=0,Nfft-1){|TF{S}(j)|^2} 386 - |TF{S}(0)|^2 - |TF{S}(Nfft-1)|^2 ]] 387 (On ne compte pas deux fois le continu et la freq de Nyquist) 388 389 390 ------- 391 --- FFT d'un signal REEL ayant un nombre impair d'elements N=2p+1 392 ------- 393 taille de la FFT: Nfft = N/2 + 1 = p + 1 394 abscisses de la fft: | 0 | 1/N | 2/N | ..... | p/N | 395 ^continu 396 (la frequence de Nyquist n'y est pas) 397 398 ... Ex: N=7 -> Nfft = 7/3+1 = 4 399 400 le signal a N elements reels, la fft a Nfft elements complexes 401 cad 2*Nfft reels = 2*(p+1) reels = 2p + 2 reels = N + 1 reels 402 soit 1 reel en trop: c'est la phase du continu 403 404 relations: 405 - si setNormalize(true) : fac = N 406 setNormalize(false) : fac = 1/N 407 sum(i=0,N-1){S(i)^2} 408 = fac* [[ 2* sum(j=0,Nfft-1){|TF{S}(j)|^2} 409 - |TF{S}(0)|^2 ]] 410 (On ne compte pas deux fois le continu) 411 412 413 ------------ 414 --- FFT-BACK d'un signal F=TF{S} ayant un nombre d'elements Nfft 415 ------------ 416 Sback = invTF{TF{S}} 417 418 Remarque: Nfft a la meme valeur pour N=2p et N=2p+1 419 donc Nfft conduit a 2 possibilites: 420 { N = 2*(Nfft-1) signal back avec nombre pair d'elements 421 { N = 2*Nfft-1 signal back avec nombre impair d'elements 422 423 Pour savoir quel est la longueur N du signal TF^(-1){F} on regarde 424 si F(Nfft-1) est reel ou complexe 425 (la frequence de Nyquist d'un signal reel est reelle) 426 427 - Si F(Nfft-1) reel cad Im{F(Nfft-1)}=0: N = 2*(Nfft-1) 428 - Si F(Nfft-1) complexe cad Im{F(Nfft-1)}#0: N = 2*Nfft-1 429 430 Si setNormalize(true): invTF{TF{S}} = S 431 setNormalize(false): invTF{TF{S}} = N * S 432 433 ========================================================================= 434 ========================================================================= 435 ============== Transformees de Fourier de signaux COMPLEXES ============= 436 ========================================================================= 437 ========================================================================= 438 439 ------- 440 --- FFT d'un signal COMPLEXE S ayant un nombre d'elements N 441 ------- 442 taille de la FFT: Nfft = N 443 abscisses de la fft: | 0 | 1/N | 2/N | ..... | (N-1)/N | 444 ^continu 445 446 Frequence de Nyquist: 447 si N est pair: la frequence de Nyquist est l'absicce d'un des bins 448 abscisses de TF{S}: Nfft = N = 2p 449 | 0 | 1/N | 2/N | ... | (N/2)/N=p/N=0.5 | ... | (N-1)/N | 450 ^frequence de Nyquist 451 si N est impair: la frequence de Nyquist N'est PAS l'absicce d'un des bins 452 abscisses de TF{S}: Nfft = N = 2p+1 453 | 0 | 1/N | 2/N | ... | (N/2)/N=p/N | ((N+1)/2)/N=(p+1)/N | ... | (N-1)/N | 454 455 ... Ex: N = 2p =6 -> Nfft = 2p = 6 456 abscisses de TF{S}: | 0 | 1/6 | 2/6 | 3/6=0.5 | 4/6 | 5/6 | 457 ... Ex: N = 2p+1 = 7 -> Nfft = 2p+1 = 7 458 abscisses de TF{S}: | 0 | 1/7 | 2/7 | 3/7 | 4/7 | 5/7 | 6/7 | 459 460 relations: 461 - si setNormalize(true) : fac = N 462 setNormalize(false) : fac = 1/N 463 sum(i=0,N-1){S(i)^2} = fac* [[ sum(j=0,Nfft-1){|TF{S}(j)|^2} ]] 464 465 ------------ 466 --- FFT-BACK d'un signal F=TF{S} ayant un nombre d'elements Nfft 467 ------------ 468 taille du signal: N = Nfft 469 470 Si setNormalize(true): invTF{TF{S}} = S 471 setNormalize(false): invTF{TF{S}} = N * S 472 473 **********************************************************************/ -
trunk/SophyaLib/NTools/fftservintf.h
r2526 r2540 46 46 virtual void FFTBackward(TArray< complex<r_8> > const & in, TArray< complex<r_8> > & out); 47 47 48 //! Compute the factor to be applied to "fftx=FFT(x)" so that "x" and "FFT(x)" have the same energy49 virtual r_8 TransfEnergyFFT(TVector<r_4> const & x, TVector< complex<r_4> > const& fftx,r_8& A,r_8& B);50 virtual r_8 TransfEnergyFFT(TVector<r_8> const & x, TVector< complex<r_8> > const& fftx,r_8& A,r_8& B);51 virtual r_8 TransfEnergyFFT(TVector< complex<r_4> > const & x, TVector< complex<r_4> > const& fftx,r_8& A,r_8& B);52 virtual r_8 TransfEnergyFFT(TVector< complex<r_8> > const & x, TVector< complex<r_8> > const& fftx,r_8& A,r_8& B);53 54 virtual inline r_8 TransfEnergyFFT(TVector<r_4> const & x, TVector< complex<r_4> > const& fftx)55 {r_8 A,B; return TransfEnergyFFT(x,fftx,A,B);}56 virtual inline r_8 TransfEnergyFFT(TVector<r_8> const & x, TVector< complex<r_8> > const& fftx)57 {r_8 A,B; return TransfEnergyFFT(x,fftx,A,B);}58 virtual inline r_8 TransfEnergyFFT(TVector< complex<r_4> > const & x, TVector< complex<r_4> > const& fftx)59 {r_8 A,B; return TransfEnergyFFT(x,fftx,A,B);}60 virtual inline r_8 TransfEnergyFFT(TVector< complex<r_8> > const & x, TVector< complex<r_8> > const& fftx)61 {r_8 A,B; return TransfEnergyFFT(x,fftx,A,B);}62 63 //! Compute the size of the FFT vector given the input vector64 virtual uint_4 SizeFFT(TVector<r_4> const & x) {return (uint_4) (x.Size()/2 + 1);}65 virtual uint_4 SizeFFT(TVector<r_8> const & x) {return (uint_4) (x.Size()/2 + 1);}66 virtual uint_4 SizeFFT(TVector< complex<r_4> > const & x) {return (uint_4) x.Size();}67 virtual uint_4 SizeFFT(TVector< complex<r_8> > const & x) {return (uint_4) x.Size();}68 69 48 protected: 70 49
Note:
See TracChangeset
for help on using the changeset viewer.