From c244e44630e611e038d1f18992bea400f7d7b6cb Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Sat, 5 Apr 2025 21:13:36 -0300
Subject: Switched to fftw allocators to ensure correct alignment for SIMD
 instructions

---
 fourier.cpp | 30 +++++++++++++++++++-----------
 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:
-- 
cgit v1.2.3-70-g09d2