libicr: FFTW.cpp Source File

Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra. Eigen itself is part of the KDE project.
00003 //
00004 // Copyright (C) 2009 Mark Borgerding mark a borgerding net
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #include "main.h"
00026 #include <unsupported/Eigen/FFT>
00027 
00028 template <typename T> 
00029 std::complex<T> RandomCpx() { return std::complex<T>( (T)(rand()/(T)RAND_MAX - .5), (T)(rand()/(T)RAND_MAX - .5) ); }
00030 
00031 using namespace std;
00032 using namespace Eigen;
00033 
00034 float norm(float x) {return x*x;}
00035 double norm(double x) {return x*x;}
00036 long double norm(long double x) {return x*x;}
00037 
00038 template < typename T>
00039 complex<long double>  promote(complex<T> x) { return complex<long double>(x.real(),x.imag()); }
00040 
00041 complex<long double>  promote(float x) { return complex<long double>( x); }
00042 complex<long double>  promote(double x) { return complex<long double>( x); }
00043 complex<long double>  promote(long double x) { return complex<long double>( x); }
00044     
00045 
00046     template <typename VT1,typename VT2>
00047     long double fft_rmse( const VT1 & fftbuf,const VT2 & timebuf)
00048     {
00049         long double totalpower=0;
00050         long double difpower=0;
00051         long double pi = acos((long double)-1 );
00052         for (size_t k0=0;k0<(size_t)fftbuf.size();++k0) {
00053             complex<long double> acc = 0;
00054             long double phinc = -2.*k0* pi / timebuf.size();
00055             for (size_t k1=0;k1<(size_t)timebuf.size();++k1) {
00056                 acc +=  promote( timebuf[k1] ) * exp( complex<long double>(0,k1*phinc) );
00057             }
00058             totalpower += norm(acc);
00059             complex<long double> x = promote(fftbuf[k0]); 
00060             complex<long double> dif = acc - x;
00061             difpower += norm(dif);
00062             //cerr << k0 << "\t" << acc << "\t" <<  x << "\t" << sqrt(norm(dif)) << endl;
00063         }
00064         cerr << "rmse:" << sqrt(difpower/totalpower) << endl;
00065         return sqrt(difpower/totalpower);
00066     }
00067 
00068     template <typename VT1,typename VT2>
00069     long double dif_rmse( const VT1 buf1,const VT2 buf2)
00070     {
00071         long double totalpower=0;
00072         long double difpower=0;
00073         size_t n = (min)( buf1.size(),buf2.size() );
00074         for (size_t k=0;k<n;++k) {
00075             totalpower += (norm( buf1[k] ) + norm(buf2[k]) )/2.;
00076             difpower += norm(buf1[k] - buf2[k]);
00077         }
00078         return sqrt(difpower/totalpower);
00079     }
00080 
00081 enum { StdVectorContainer, EigenVectorContainer };
00082 
00083 template<int Container, typename Scalar> struct VectorType;
00084 
00085 template<typename Scalar> struct VectorType<StdVectorContainer,Scalar>
00086 {
00087   typedef vector<Scalar> type;
00088 };
00089 
00090 template<typename Scalar> struct VectorType<EigenVectorContainer,Scalar>
00091 {
00092   typedef Matrix<Scalar,Dynamic,1> type;
00093 };
00094 
00095 template <int Container, typename T>
00096 void test_scalar_generic(int nfft)
00097 {
00098     typedef typename FFT<T>::Complex Complex;
00099     typedef typename FFT<T>::Scalar Scalar;
00100     typedef typename VectorType<Container,Scalar>::type ScalarVector;
00101     typedef typename VectorType<Container,Complex>::type ComplexVector;
00102 
00103     FFT<T> fft;
00104     ScalarVector tbuf(nfft);
00105     ComplexVector freqBuf;
00106     for (int k=0;k<nfft;++k)
00107         tbuf[k]= (T)( rand()/(double)RAND_MAX - .5);
00108 
00109     // make sure it DOESN'T give the right full spectrum answer
00110     // if we've asked for half-spectrum
00111     fft.SetFlag(fft.HalfSpectrum );
00112     fft.fwd( freqBuf,tbuf);
00113     VERIFY((size_t)freqBuf.size() == (size_t)( (nfft>>1)+1) );
00114     VERIFY( fft_rmse(freqBuf,tbuf) < test_precision<T>()  );// gross check
00115 
00116     fft.ClearFlag(fft.HalfSpectrum );
00117     fft.fwd( freqBuf,tbuf);
00118     VERIFY( (size_t)freqBuf.size() == (size_t)nfft);
00119     VERIFY( fft_rmse(freqBuf,tbuf) < test_precision<T>()  );// gross check
00120 
00121     if (nfft&1)
00122         return; // odd FFTs get the wrong size inverse FFT
00123 
00124     ScalarVector tbuf2;
00125     fft.inv( tbuf2 , freqBuf);
00126     VERIFY( dif_rmse(tbuf,tbuf2) < test_precision<T>()  );// gross check
00127 
00128 
00129     // verify that the Unscaled flag takes effect
00130     ScalarVector tbuf3;
00131     fft.SetFlag(fft.Unscaled);
00132 
00133     fft.inv( tbuf3 , freqBuf);
00134 
00135     for (int k=0;k<nfft;++k)
00136         tbuf3[k] *= T(1./nfft);
00137 
00138 
00139     //for (size_t i=0;i<(size_t) tbuf.size();++i)
00140     //    cout << "freqBuf=" << freqBuf[i] << " in2=" << tbuf3[i] << " -  in=" << tbuf[i] << " => " << (tbuf3[i] - tbuf[i] ) <<  endl;
00141 
00142     VERIFY( dif_rmse(tbuf,tbuf3) < test_precision<T>()  );// gross check
00143 
00144     // verify that ClearFlag works
00145     fft.ClearFlag(fft.Unscaled);
00146     fft.inv( tbuf2 , freqBuf);
00147     VERIFY( dif_rmse(tbuf,tbuf2) < test_precision<T>()  );// gross check
00148 }
00149 
00150 template <typename T>
00151 void test_scalar(int nfft)
00152 {
00153   test_scalar_generic<StdVectorContainer,T>(nfft);
00154   //test_scalar_generic<EigenVectorContainer,T>(nfft);
00155 }
00156 
00157 
00158 template <int Container, typename T>
00159 void test_complex_generic(int nfft)
00160 {
00161     typedef typename FFT<T>::Complex Complex;
00162     typedef typename VectorType<Container,Complex>::type ComplexVector;
00163 
00164     FFT<T> fft;
00165 
00166     ComplexVector inbuf(nfft);
00167     ComplexVector outbuf;
00168     ComplexVector buf3;
00169     for (int k=0;k<nfft;++k)
00170         inbuf[k]= Complex( (T)(rand()/(double)RAND_MAX - .5), (T)(rand()/(double)RAND_MAX - .5) );
00171     fft.fwd( outbuf , inbuf);
00172 
00173     VERIFY( fft_rmse(outbuf,inbuf) < test_precision<T>()  );// gross check
00174     fft.inv( buf3 , outbuf);
00175 
00176     VERIFY( dif_rmse(inbuf,buf3) < test_precision<T>()  );// gross check
00177 
00178     // verify that the Unscaled flag takes effect
00179     ComplexVector buf4;
00180     fft.SetFlag(fft.Unscaled);
00181     fft.inv( buf4 , outbuf);
00182     for (int k=0;k<nfft;++k)
00183         buf4[k] *= T(1./nfft);
00184     VERIFY( dif_rmse(inbuf,buf4) < test_precision<T>()  );// gross check
00185 
00186     // verify that ClearFlag works
00187     fft.ClearFlag(fft.Unscaled);
00188     fft.inv( buf3 , outbuf);
00189     VERIFY( dif_rmse(inbuf,buf3) < test_precision<T>()  );// gross check
00190 }
00191 
00192 template <typename T>
00193 void test_complex(int nfft)
00194 {
00195   test_complex_generic<StdVectorContainer,T>(nfft);
00196   test_complex_generic<EigenVectorContainer,T>(nfft);
00197 }
00198 /*
00199 template <typename T,int nrows,int ncols>
00200 void test_complex2d()
00201 {
00202     typedef typename Eigen::FFT<T>::Complex Complex;
00203     FFT<T> fft;
00204     Eigen::Matrix<Complex,nrows,ncols> src,src2,dst,dst2;
00205 
00206     src = Eigen::Matrix<Complex,nrows,ncols>::Random();
00207     //src =  Eigen::Matrix<Complex,nrows,ncols>::Identity();
00208 
00209     for (int k=0;k<ncols;k++) {
00210         Eigen::Matrix<Complex,nrows,1> tmpOut;
00211         fft.fwd( tmpOut,src.col(k) );
00212         dst2.col(k) = tmpOut;
00213     }
00214 
00215     for (int k=0;k<nrows;k++) {
00216         Eigen::Matrix<Complex,1,ncols> tmpOut;
00217         fft.fwd( tmpOut,  dst2.row(k) );
00218         dst2.row(k) = tmpOut;
00219     }
00220 
00221     fft.fwd2(dst.data(),src.data(),ncols,nrows);
00222     fft.inv2(src2.data(),dst.data(),ncols,nrows);
00223     VERIFY( (src-src2).norm() < test_precision<T>() );
00224     VERIFY( (dst-dst2).norm() < test_precision<T>() );
00225 }
00226 */
00227 
00228 
00229 void test_return_by_value(int len)
00230 {
00231     VectorXf in;
00232     VectorXf in1;
00233     in.setRandom( len );
00234     VectorXcf out1,out2;
00235     FFT<float> fft;
00236 
00237     fft.SetFlag(fft.HalfSpectrum );
00238 
00239     fft.fwd(out1,in);
00240     out2 = fft.fwd(in);
00241     VERIFY( (out1-out2).norm() < test_precision<float>() );
00242     in1 = fft.inv(out1);
00243     VERIFY( (in1-in).norm() < test_precision<float>() );
00244 }
00245 
00246 void test_FFTW()
00247 {
00248   CALL_SUBTEST( test_return_by_value(32) );
00249   //CALL_SUBTEST( ( test_complex2d<float,4,8> () ) ); CALL_SUBTEST( ( test_complex2d<double,4,8> () ) );
00250   //CALL_SUBTEST( ( test_complex2d<long double,4,8> () ) );
00251   CALL_SUBTEST( test_complex<float>(32) ); CALL_SUBTEST( test_complex<double>(32) ); 
00252   CALL_SUBTEST( test_complex<float>(256) ); CALL_SUBTEST( test_complex<double>(256) ); 
00253   CALL_SUBTEST( test_complex<float>(3*8) ); CALL_SUBTEST( test_complex<double>(3*8) ); 
00254   CALL_SUBTEST( test_complex<float>(5*32) ); CALL_SUBTEST( test_complex<double>(5*32) ); 
00255   CALL_SUBTEST( test_complex<float>(2*3*4) ); CALL_SUBTEST( test_complex<double>(2*3*4) ); 
00256   CALL_SUBTEST( test_complex<float>(2*3*4*5) ); CALL_SUBTEST( test_complex<double>(2*3*4*5) ); 
00257   CALL_SUBTEST( test_complex<float>(2*3*4*5*7) ); CALL_SUBTEST( test_complex<double>(2*3*4*5*7) ); 
00258 
00259   CALL_SUBTEST( test_scalar<float>(32) ); CALL_SUBTEST( test_scalar<double>(32) ); 
00260   CALL_SUBTEST( test_scalar<float>(45) ); CALL_SUBTEST( test_scalar<double>(45) ); 
00261   CALL_SUBTEST( test_scalar<float>(50) ); CALL_SUBTEST( test_scalar<double>(50) ); 
00262   CALL_SUBTEST( test_scalar<float>(256) ); CALL_SUBTEST( test_scalar<double>(256) ); 
00263   CALL_SUBTEST( test_scalar<float>(2*3*4*5*7) ); CALL_SUBTEST( test_scalar<double>(2*3*4*5*7) ); 
00264   
00265   #ifdef EIGEN_HAS_FFTWL
00266   CALL_SUBTEST( test_complex<long double>(32) );
00267   CALL_SUBTEST( test_complex<long double>(256) );
00268   CALL_SUBTEST( test_complex<long double>(3*8) );
00269   CALL_SUBTEST( test_complex<long double>(5*32) );
00270   CALL_SUBTEST( test_complex<long double>(2*3*4) );
00271   CALL_SUBTEST( test_complex<long double>(2*3*4*5) );
00272   CALL_SUBTEST( test_complex<long double>(2*3*4*5*7) );
00273   
00274   CALL_SUBTEST( test_scalar<long double>(32) );
00275   CALL_SUBTEST( test_scalar<long double>(45) );
00276   CALL_SUBTEST( test_scalar<long double>(50) );
00277   CALL_SUBTEST( test_scalar<long double>(256) );
00278   CALL_SUBTEST( test_scalar<long double>(2*3*4*5*7) );
00279   #endif
00280 }