1 files changed, 19 insertions, 11 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;