blob: 0e040ad706bb0ce21c1631d50a486281e4874a0a [file] [log] [blame]
 // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2009 Mark Borgerding mark a borgerding net // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. #include "main.h" #include template inline std::complex RandomCpx() { return std::complex((T)(rand() / (T)RAND_MAX - .5), (T)(rand() / (T)RAND_MAX - .5)); } using namespace std; using namespace Eigen; template inline complex promote(complex x) { return complex((long double)x.real(), (long double)x.imag()); } inline complex promote(float x) { return complex((long double)x); } inline complex promote(double x) { return complex((long double)x); } inline complex promote(long double x) { return complex((long double)x); } template long double fft_rmse(const VT1& fftbuf, const VT2& timebuf) { long double totalpower = 0; long double difpower = 0; long double pi = acos((long double)-1); for (size_t k0 = 0; k0 < (size_t)fftbuf.size(); ++k0) { complex acc = 0; long double phinc = (long double)(-2.) * k0 * pi / timebuf.size(); for (size_t k1 = 0; k1 < (size_t)timebuf.size(); ++k1) { acc += promote(timebuf[k1]) * exp(complex(0, k1 * phinc)); } totalpower += numext::abs2(acc); complex x = promote(fftbuf[k0]); complex dif = acc - x; difpower += numext::abs2(dif); // cerr << k0 << "\t" << acc << "\t" << x << "\t" << sqrt(numext::abs2(dif)) << endl; } // cerr << "rmse:" << sqrt(difpower/totalpower) << endl; return sqrt(difpower / totalpower); } template long double dif_rmse(const VT1 buf1, const VT2 buf2) { long double totalpower = 0; long double difpower = 0; size_t n = (min)(buf1.size(), buf2.size()); for (size_t k = 0; k < n; ++k) { totalpower += (long double)((numext::abs2(buf1[k]) + numext::abs2(buf2[k])) / 2); difpower += (long double)(numext::abs2(buf1[k] - buf2[k])); } return sqrt(difpower / totalpower); } enum { StdVectorContainer, EigenVectorContainer }; template struct VectorType; template struct VectorType { typedef vector type; }; template struct VectorType { typedef Matrix type; }; template void test_scalar_generic(int nfft) { typedef typename FFT::Complex Complex; typedef typename FFT::Scalar Scalar; typedef typename VectorType::type ScalarVector; typedef typename VectorType::type ComplexVector; FFT fft; ScalarVector tbuf(nfft); ComplexVector freqBuf; for (int k = 0; k < nfft; ++k) tbuf[k] = (T)(rand() / (double)RAND_MAX - .5); // make sure it DOESN'T give the right full spectrum answer // if we've asked for half-spectrum fft.SetFlag(fft.HalfSpectrum); fft.fwd(freqBuf, tbuf); VERIFY((size_t)freqBuf.size() == (size_t)((nfft >> 1) + 1)); VERIFY(T(fft_rmse(freqBuf, tbuf)) < test_precision()); // gross check fft.ClearFlag(fft.HalfSpectrum); fft.fwd(freqBuf, tbuf); VERIFY((size_t)freqBuf.size() == (size_t)nfft); VERIFY(T(fft_rmse(freqBuf, tbuf)) < test_precision()); // gross check if (nfft & 1) return; // odd FFTs get the wrong size inverse FFT ScalarVector tbuf2; fft.inv(tbuf2, freqBuf); VERIFY(T(dif_rmse(tbuf, tbuf2)) < test_precision()); // gross check // verify that the Unscaled flag takes effect ScalarVector tbuf3; fft.SetFlag(fft.Unscaled); fft.inv(tbuf3, freqBuf); for (int k = 0; k < nfft; ++k) tbuf3[k] *= T(1. / nfft); // for (size_t i=0;i<(size_t) tbuf.size();++i) // cout << "freqBuf=" << freqBuf[i] << " in2=" << tbuf3[i] << " - in=" << tbuf[i] << " => " << (tbuf3[i] - // tbuf[i] ) << endl; VERIFY(T(dif_rmse(tbuf, tbuf3)) < test_precision()); // gross check // verify that ClearFlag works fft.ClearFlag(fft.Unscaled); fft.inv(tbuf2, freqBuf); VERIFY(T(dif_rmse(tbuf, tbuf2)) < test_precision()); // gross check } template void test_scalar(int nfft) { test_scalar_generic(nfft); // test_scalar_generic(nfft); } template void test_complex_generic(int nfft) { typedef typename FFT::Complex Complex; typedef typename VectorType::type ComplexVector; FFT fft; ComplexVector inbuf(nfft); ComplexVector outbuf; ComplexVector buf3; for (int k = 0; k < nfft; ++k) inbuf[k] = Complex((T)(rand() / (double)RAND_MAX - .5), (T)(rand() / (double)RAND_MAX - .5)); fft.fwd(outbuf, inbuf); VERIFY(T(fft_rmse(outbuf, inbuf)) < test_precision()); // gross check fft.inv(buf3, outbuf); VERIFY(T(dif_rmse(inbuf, buf3)) < test_precision()); // gross check // verify that the Unscaled flag takes effect ComplexVector buf4; fft.SetFlag(fft.Unscaled); fft.inv(buf4, outbuf); for (int k = 0; k < nfft; ++k) buf4[k] *= T(1. / nfft); VERIFY(T(dif_rmse(inbuf, buf4)) < test_precision()); // gross check // verify that ClearFlag works fft.ClearFlag(fft.Unscaled); fft.inv(buf3, outbuf); VERIFY(T(dif_rmse(inbuf, buf3)) < test_precision()); // gross check } template void test_complex(int nfft) { test_complex_generic(nfft); test_complex_generic(nfft); } template void test_complex2d() { typedef typename Eigen::FFT::Complex Complex; FFT fft; Eigen::Matrix src, src2, dst, dst2; src = Eigen::Matrix::Random(); // src = Eigen::Matrix::Identity(); for (int k = 0; k < ncols; k++) { Eigen::Matrix tmpOut; fft.fwd(tmpOut, src.col(k)); dst2.col(k) = tmpOut; } for (int k = 0; k < nrows; k++) { Eigen::Matrix tmpOut; fft.fwd(tmpOut, dst2.row(k)); dst2.row(k) = tmpOut; } fft.fwd2(dst.data(), src.data(), ncols, nrows); fft.inv2(src2.data(), dst.data(), ncols, nrows); VERIFY((src - src2).norm() < test_precision()); VERIFY((dst - dst2).norm() < test_precision()); } inline void test_return_by_value(int len) { VectorXf in; VectorXf in1; in.setRandom(len); VectorXcf out1, out2; FFT fft; fft.SetFlag(fft.HalfSpectrum); fft.fwd(out1, in); out2 = fft.fwd(in); VERIFY((out1 - out2).norm() < test_precision()); in1 = fft.inv(out1); VERIFY((in1 - in).norm() < test_precision()); } EIGEN_DECLARE_TEST(FFTW) { CALL_SUBTEST(test_return_by_value(32)); CALL_SUBTEST(test_complex(32)); CALL_SUBTEST(test_complex(32)); CALL_SUBTEST(test_complex(256)); CALL_SUBTEST(test_complex(256)); CALL_SUBTEST(test_complex(3 * 8)); CALL_SUBTEST(test_complex(3 * 8)); CALL_SUBTEST(test_complex(5 * 32)); CALL_SUBTEST(test_complex(5 * 32)); CALL_SUBTEST(test_complex(2 * 3 * 4)); CALL_SUBTEST(test_complex(2 * 3 * 4)); CALL_SUBTEST(test_complex(2 * 3 * 4 * 5)); CALL_SUBTEST(test_complex(2 * 3 * 4 * 5)); CALL_SUBTEST(test_complex(2 * 3 * 4 * 5 * 7)); CALL_SUBTEST(test_complex(2 * 3 * 4 * 5 * 7)); CALL_SUBTEST(test_scalar(32)); CALL_SUBTEST(test_scalar(32)); CALL_SUBTEST(test_scalar(45)); CALL_SUBTEST(test_scalar(45)); CALL_SUBTEST(test_scalar(50)); CALL_SUBTEST(test_scalar(50)); CALL_SUBTEST(test_scalar(256)); CALL_SUBTEST(test_scalar(256)); CALL_SUBTEST(test_scalar(2 * 3 * 4 * 5 * 7)); CALL_SUBTEST(test_scalar(2 * 3 * 4 * 5 * 7)); #if defined EIGEN_HAS_FFTWL || defined EIGEN_POCKETFFT_DEFAULT CALL_SUBTEST(test_complex(32)); CALL_SUBTEST(test_complex(256)); CALL_SUBTEST(test_complex(3 * 8)); CALL_SUBTEST(test_complex(5 * 32)); CALL_SUBTEST(test_complex(2 * 3 * 4)); CALL_SUBTEST(test_complex(2 * 3 * 4 * 5)); CALL_SUBTEST(test_complex(2 * 3 * 4 * 5 * 7)); CALL_SUBTEST(test_scalar(32)); CALL_SUBTEST(test_scalar(45)); CALL_SUBTEST(test_scalar(50)); CALL_SUBTEST(test_scalar(256)); CALL_SUBTEST(test_scalar(2 * 3 * 4 * 5 * 7)); CALL_SUBTEST((test_complex2d())); CALL_SUBTEST((test_complex2d())); CALL_SUBTEST((test_complex2d())); CALL_SUBTEST((test_complex2d())); // fail to build since Eigen limit the stack allocation size,too big here. // CALL_SUBTEST( ( test_complex2d () ) ); #endif #if defined EIGEN_FFTW_DEFAULT || defined EIGEN_POCKETFFT_DEFAULT || defined EIGEN_MKL_DEFAULT CALL_SUBTEST((test_complex2d())); CALL_SUBTEST((test_complex2d())); CALL_SUBTEST((test_complex2d())); CALL_SUBTEST((test_complex2d())); #endif #if defined EIGEN_FFTW_DEFAULT || defined EIGEN_POCKETFFT_DEFAULT || defined EIGEN_MKL_DEFAULT CALL_SUBTEST((test_complex2d())); CALL_SUBTEST((test_complex2d())); CALL_SUBTEST((test_complex2d())); CALL_SUBTEST((test_complex2d())); #endif }