diff options
| -rw-r--r-- | fourier.cpp | 30 | ||||
| -rw-r--r-- | fourier.hpp | 5 |
2 files changed, 22 insertions, 13 deletions
diff --git a/fourier.cpp b/fourier.cpp index c2ed600..07f8fdc 100644 --- a/fourier.cpp +++ b/fourier.cpp @@ -24,26 +24,32 @@ Real ddf(Real λ, unsigned p, unsigned s, Real q) { return (1 - λ) * ddfP(p, q) + λ * ddfP(s, q); } -FourierTransform::FourierTransform(unsigned n, Real Δω, Real Δτ, unsigned flags) : a(2 * n), â(n + 1), Δω(Δω), Δτ(Δτ) { +FourierTransform::FourierTransform(unsigned n, Real Δω, Real Δτ, unsigned flags) : n(n), Δω(Δω), Δτ(Δτ) { + a = fftw_alloc_real(2 * n); + â = reinterpret_cast<Complex*>(fftw_alloc_complex(n + 1)); fftw_init_threads(); fftw_plan_with_nthreads(FFTW_THREADS); fftw_import_wisdom_from_filename("fftw.wisdom"); - plan_r2c = fftw_plan_dft_r2c_1d(2 * n, a.data(), reinterpret_cast<fftw_complex*>(â.data()), flags); - plan_c2r = fftw_plan_dft_c2r_1d(2 * n, reinterpret_cast<fftw_complex*>(â.data()), a.data(), flags); + plan_r2c = fftw_plan_dft_r2c_1d(2 * n, a, reinterpret_cast<fftw_complex*>(â), flags); + plan_c2r = fftw_plan_dft_c2r_1d(2 * n, reinterpret_cast<fftw_complex*>(â), a, flags); fftw_export_wisdom_to_filename("fftw.wisdom"); } FourierTransform::~FourierTransform() { fftw_destroy_plan(plan_r2c); fftw_destroy_plan(plan_c2r); + fftw_free(a); + fftw_free(â); fftw_cleanup(); } std::vector<Complex> FourierTransform::fourier(const std::vector<Real>& c) { - a = c; + for (unsigned i = 0; i < 2 * n; i++) { + a[i] = c[i]; + } fftw_execute(plan_r2c); - std::vector<Complex> ĉ(â.size()); - for (unsigned i = 0; i < â.size(); i++) { + std::vector<Complex> ĉ(n + 1); + for (unsigned i = 0; i < n + 1; i++) { ĉ[i] = â[i] * (Δτ * M_PI); } return ĉ; @@ -51,18 +57,20 @@ std::vector<Complex> FourierTransform::fourier(const std::vector<Real>& c) { std::vector<Complex> FourierTransform::fourier() { fftw_execute(plan_r2c); - std::vector<Complex> ĉ(â.size()); - for (unsigned i = 0; i < â.size(); i++) { + std::vector<Complex> ĉ(n+1); + for (unsigned i = 0; i < n+1; i++) { ĉ[i] = â[i] * (Δτ * M_PI); } return ĉ; } std::vector<Real> FourierTransform::inverse(const std::vector<Complex>& ĉ) { - â = ĉ; + for (unsigned i = 0; i < n + 1; i++) { + â[i] = ĉ[i]; + } fftw_execute(plan_c2r); - std::vector<Real> c(a.size()); - for (unsigned i = 0; i < a.size(); i++) { + std::vector<Real> c(2*n); + for (unsigned i = 0; i < 2*n; i++) { c[i] = a[i] * (Δω / (2 * M_PI)); } return c; diff --git a/fourier.hpp b/fourier.hpp index 1ebb7bc..9451f69 100644 --- a/fourier.hpp +++ b/fourier.hpp @@ -19,10 +19,11 @@ Real ddf(Real λ, unsigned p, unsigned s, Real q); class FourierTransform { private: - std::vector<Real> a; - std::vector<Complex> â; + Real* a; + Complex* â; fftw_plan plan_r2c; fftw_plan plan_c2r; + unsigned n; Real Δω; Real Δτ; public: |