Revert "Some insigned gained"

This reverts commit 79aa3df97669ddb8836752d55faa16de44d8e3ba.

1 files changed, 21 insertions, 18 deletions

diff --git a/log-fourier.cpp b/log-fourier.cpp
index e15f62a..bc2dd87 100644
--- a/log-fourier.cpp
+++ b/log-fourier.cpp
@@ -7,12 +7,21 @@
 LogarithmicFourierTransform::LogarithmicFourierTransform(unsigned N, Real k, Real Δτ, unsigned pad) : N(N), pad(pad), k(k), Δτ(Δτ) {
   τₛ = -0.5 * N;
   ωₛ = -0.5 * N;
-  sₛ = -0.5 * N;
-  Δs = 0.1;
-  Δω = Δτ;
+  sₛ = -0.5 * pad * N;
+  a = reinterpret_cast<Complex*>(FFTW_ALLOC_COMPLEX(pad*N));
+  â = reinterpret_cast<Complex*>(FFTW_ALLOC_COMPLEX(pad*N));
+  FFTW_IMPORT_WISDOM("fftw.wisdom");
+  a_to_â = FFTW_PLAN_DFT_1D(pad*N, reinterpret_cast<FFTW_COMPLEX*>(a), reinterpret_cast<FFTW_COMPLEX*>(â), FFTW_BACKWARD, 0);
+  â_to_a = FFTW_PLAN_DFT_1D(pad*N, reinterpret_cast<FFTW_COMPLEX*>(â), reinterpret_cast<FFTW_COMPLEX*>(a), FFTW_BACKWARD, 0);
+  FFTW_EXPORT_WISDOM("fftw.wisdom");
 }
 
 LogarithmicFourierTransform::~LogarithmicFourierTransform() {
+  FFTW_DESTROY_PLAN(a_to_â);
+  FFTW_DESTROY_PLAN(â_to_a);
+  FFTW_FREE(a);
+  FFTW_FREE(â);
+  FFTW_CLEANUP();
 }
 
 Real LogarithmicFourierTransform::τ(unsigned n) const {
@@ -24,7 +33,7 @@ Real LogarithmicFourierTransform::ω(unsigned n) const {
 }
 
 Real LogarithmicFourierTransform::s(unsigned n) const {
-  return (n + sₛ) * Δs;
+  return (n + sₛ) * 2*M_PI / (pad * N * Δτ);
 }
 
 Real LogarithmicFourierTransform::t(unsigned n) const {
@@ -75,24 +84,18 @@ std::vector<Complex> LogarithmicFourierTransform::fourier(const std::vector<Real
 std::vector<Real> LogarithmicFourierTransform::inverse(const std::vector<Complex>& ĉ) {
   std::vector<Real> c(N);
   std::vector<Real> σs = {1, -1};
-
   for (Real σ : σs) {
-    std::vector<Complex> F1(N);
-    for (unsigned l = 0; l < N; l++) {
-      for (unsigned m = 0; m < N; m++) {
-        F1[l] +=  (ĉ[m].real() + II * σ * ĉ[m].imag()) * std::exp((1 - k) * ω(m) + II * ω(m) * s(l));
-      }
+    for (unsigned n = 0; n < N; n++) {
+      a[n] = (ĉ[n].real() + II * σ * ĉ[n].imag()) * std::exp((1 - k) * ω(n));
     }
+    FFTW_EXECUTE(a_to_â);
     for (unsigned n = 0; n < N; n++) {
-      for (unsigned l = 0; l < N; l++) {
-        c[n] += (F1[l] * std::exp(II * s(l) * τ(n) + II * σ * (M_PI / 2) * (II * s(l) - k)) * Γ(k - II * s(l))).real();
-      }
-      c[n] *= (Δs / (2 * M_PI)) * (Δω / (2 * M_PI)) * exp(-k * τ(n));
+      â[(N / 2 + n) % N] *= std::pow(-II * σ, II * s(n) - k) * Γ(k - II * s(n));
+    }
+    FFTW_EXECUTE(â_to_a);
+    for (unsigned n = 0; n < N; n++) {
+      c[n] += std::exp(-k * τ(n)) * a[n].real() / (2 * M_PI * N);
     }
-  }
-
-  for (unsigned i = 0; i < N; i++) {
-    c[i] -= c[N-1];
   }
 
   return c;