278 lines
9.6 KiB
C++
278 lines
9.6 KiB
C++
// 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 <unsupported/Eigen/FFT>
|
|
|
|
template <typename T>
|
|
inline std::complex<T> RandomCpx() {
|
|
return std::complex<T>((T)(rand() / (T)RAND_MAX - .5), (T)(rand() / (T)RAND_MAX - .5));
|
|
}
|
|
|
|
using namespace std;
|
|
using namespace Eigen;
|
|
|
|
template <typename T>
|
|
inline complex<long double> promote(complex<T> x) {
|
|
return complex<long double>((long double)x.real(), (long double)x.imag());
|
|
}
|
|
|
|
inline complex<long double> promote(float x) { return complex<long double>((long double)x); }
|
|
inline complex<long double> promote(double x) { return complex<long double>((long double)x); }
|
|
inline complex<long double> promote(long double x) { return complex<long double>((long double)x); }
|
|
|
|
template <typename VT1, typename VT2>
|
|
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<long double> 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<long double>(0, k1 * phinc));
|
|
}
|
|
totalpower += numext::abs2(acc);
|
|
complex<long double> x = promote(fftbuf[k0]);
|
|
complex<long double> 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 <typename VT1, typename VT2>
|
|
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 <int Container, typename Scalar>
|
|
struct VectorType;
|
|
|
|
template <typename Scalar>
|
|
struct VectorType<StdVectorContainer, Scalar> {
|
|
typedef vector<Scalar> type;
|
|
};
|
|
|
|
template <typename Scalar>
|
|
struct VectorType<EigenVectorContainer, Scalar> {
|
|
typedef Matrix<Scalar, Dynamic, 1> type;
|
|
};
|
|
|
|
template <int Container, typename T>
|
|
void test_scalar_generic(int nfft) {
|
|
typedef typename FFT<T>::Complex Complex;
|
|
typedef typename FFT<T>::Scalar Scalar;
|
|
typedef typename VectorType<Container, Scalar>::type ScalarVector;
|
|
typedef typename VectorType<Container, Complex>::type ComplexVector;
|
|
|
|
FFT<T> 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<T>()); // 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<T>()); // 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<T>()); // 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<T>()); // gross check
|
|
|
|
// verify that ClearFlag works
|
|
fft.ClearFlag(fft.Unscaled);
|
|
fft.inv(tbuf2, freqBuf);
|
|
VERIFY(T(dif_rmse(tbuf, tbuf2)) < test_precision<T>()); // gross check
|
|
}
|
|
|
|
template <typename T>
|
|
void test_scalar(int nfft) {
|
|
test_scalar_generic<StdVectorContainer, T>(nfft);
|
|
// test_scalar_generic<EigenVectorContainer,T>(nfft);
|
|
}
|
|
|
|
template <int Container, typename T>
|
|
void test_complex_generic(int nfft) {
|
|
typedef typename FFT<T>::Complex Complex;
|
|
typedef typename VectorType<Container, Complex>::type ComplexVector;
|
|
|
|
FFT<T> 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<T>()); // gross check
|
|
fft.inv(buf3, outbuf);
|
|
|
|
VERIFY(T(dif_rmse(inbuf, buf3)) < test_precision<T>()); // 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<T>()); // gross check
|
|
|
|
// verify that ClearFlag works
|
|
fft.ClearFlag(fft.Unscaled);
|
|
fft.inv(buf3, outbuf);
|
|
VERIFY(T(dif_rmse(inbuf, buf3)) < test_precision<T>()); // gross check
|
|
}
|
|
|
|
template <typename T>
|
|
void test_complex(int nfft) {
|
|
test_complex_generic<StdVectorContainer, T>(nfft);
|
|
test_complex_generic<EigenVectorContainer, T>(nfft);
|
|
}
|
|
|
|
template <typename T, int nrows, int ncols>
|
|
void test_complex2d() {
|
|
typedef typename Eigen::FFT<T>::Complex Complex;
|
|
FFT<T> fft;
|
|
Eigen::Matrix<Complex, nrows, ncols> src, src2, dst, dst2;
|
|
|
|
src = Eigen::Matrix<Complex, nrows, ncols>::Random();
|
|
// src = Eigen::Matrix<Complex,nrows,ncols>::Identity();
|
|
|
|
for (int k = 0; k < ncols; k++) {
|
|
Eigen::Matrix<Complex, nrows, 1> tmpOut;
|
|
fft.fwd(tmpOut, src.col(k));
|
|
dst2.col(k) = tmpOut;
|
|
}
|
|
|
|
for (int k = 0; k < nrows; k++) {
|
|
Eigen::Matrix<Complex, 1, ncols> 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<T>());
|
|
VERIFY((dst - dst2).norm() < test_precision<T>());
|
|
}
|
|
|
|
inline void test_return_by_value(int len) {
|
|
VectorXf in;
|
|
VectorXf in1;
|
|
in.setRandom(len);
|
|
VectorXcf out1, out2;
|
|
FFT<float> fft;
|
|
|
|
fft.SetFlag(fft.HalfSpectrum);
|
|
|
|
fft.fwd(out1, in);
|
|
out2 = fft.fwd(in);
|
|
VERIFY((out1 - out2).norm() < test_precision<float>());
|
|
in1 = fft.inv(out1);
|
|
VERIFY((in1 - in).norm() < test_precision<float>());
|
|
}
|
|
|
|
EIGEN_DECLARE_TEST(FFTW) {
|
|
CALL_SUBTEST(test_return_by_value(32));
|
|
CALL_SUBTEST(test_complex<float>(32));
|
|
CALL_SUBTEST(test_complex<double>(32));
|
|
CALL_SUBTEST(test_complex<float>(256));
|
|
CALL_SUBTEST(test_complex<double>(256));
|
|
CALL_SUBTEST(test_complex<float>(3 * 8));
|
|
CALL_SUBTEST(test_complex<double>(3 * 8));
|
|
CALL_SUBTEST(test_complex<float>(5 * 32));
|
|
CALL_SUBTEST(test_complex<double>(5 * 32));
|
|
CALL_SUBTEST(test_complex<float>(2 * 3 * 4));
|
|
CALL_SUBTEST(test_complex<double>(2 * 3 * 4));
|
|
CALL_SUBTEST(test_complex<float>(2 * 3 * 4 * 5));
|
|
CALL_SUBTEST(test_complex<double>(2 * 3 * 4 * 5));
|
|
CALL_SUBTEST(test_complex<float>(2 * 3 * 4 * 5 * 7));
|
|
CALL_SUBTEST(test_complex<double>(2 * 3 * 4 * 5 * 7));
|
|
|
|
CALL_SUBTEST(test_scalar<float>(32));
|
|
CALL_SUBTEST(test_scalar<double>(32));
|
|
CALL_SUBTEST(test_scalar<float>(45));
|
|
CALL_SUBTEST(test_scalar<double>(45));
|
|
CALL_SUBTEST(test_scalar<float>(50));
|
|
CALL_SUBTEST(test_scalar<double>(50));
|
|
CALL_SUBTEST(test_scalar<float>(256));
|
|
CALL_SUBTEST(test_scalar<double>(256));
|
|
CALL_SUBTEST(test_scalar<float>(2 * 3 * 4 * 5 * 7));
|
|
CALL_SUBTEST(test_scalar<double>(2 * 3 * 4 * 5 * 7));
|
|
|
|
#if defined EIGEN_HAS_FFTWL || defined EIGEN_POCKETFFT_DEFAULT
|
|
CALL_SUBTEST(test_complex<long double>(32));
|
|
CALL_SUBTEST(test_complex<long double>(256));
|
|
CALL_SUBTEST(test_complex<long double>(3 * 8));
|
|
CALL_SUBTEST(test_complex<long double>(5 * 32));
|
|
CALL_SUBTEST(test_complex<long double>(2 * 3 * 4));
|
|
CALL_SUBTEST(test_complex<long double>(2 * 3 * 4 * 5));
|
|
CALL_SUBTEST(test_complex<long double>(2 * 3 * 4 * 5 * 7));
|
|
|
|
CALL_SUBTEST(test_scalar<long double>(32));
|
|
CALL_SUBTEST(test_scalar<long double>(45));
|
|
CALL_SUBTEST(test_scalar<long double>(50));
|
|
CALL_SUBTEST(test_scalar<long double>(256));
|
|
CALL_SUBTEST(test_scalar<long double>(2 * 3 * 4 * 5 * 7));
|
|
|
|
CALL_SUBTEST((test_complex2d<long double, 2 * 3 * 4, 2 * 3 * 4>()));
|
|
CALL_SUBTEST((test_complex2d<long double, 3 * 4 * 5, 3 * 4 * 5>()));
|
|
CALL_SUBTEST((test_complex2d<long double, 24, 60>()));
|
|
CALL_SUBTEST((test_complex2d<long double, 60, 24>()));
|
|
// fail to build since Eigen limit the stack allocation size,too big here.
|
|
// CALL_SUBTEST( ( test_complex2d<long double, 256, 256> () ) );
|
|
#endif
|
|
#if defined EIGEN_FFTW_DEFAULT || defined EIGEN_POCKETFFT_DEFAULT || defined EIGEN_MKL_DEFAULT
|
|
CALL_SUBTEST((test_complex2d<float, 24, 24>()));
|
|
CALL_SUBTEST((test_complex2d<float, 60, 60>()));
|
|
CALL_SUBTEST((test_complex2d<float, 24, 60>()));
|
|
CALL_SUBTEST((test_complex2d<float, 60, 24>()));
|
|
#endif
|
|
#if defined EIGEN_FFTW_DEFAULT || defined EIGEN_POCKETFFT_DEFAULT || defined EIGEN_MKL_DEFAULT
|
|
CALL_SUBTEST((test_complex2d<double, 24, 24>()));
|
|
CALL_SUBTEST((test_complex2d<double, 60, 60>()));
|
|
CALL_SUBTEST((test_complex2d<double, 24, 60>()));
|
|
CALL_SUBTEST((test_complex2d<double, 60, 24>()));
|
|
#endif
|
|
}
|