Some insigned gained

1 files changed, 18 insertions, 21 deletions

diff --git a/log-fourier.cpp b/log-fourier.cpp
index bc2dd87..e15f62a 100644
--- a/log-fourier.cpp
+++ b/log-fourier.cpp
@@ -7,21 +7,12 @@
 LogarithmicFourierTransform::LogarithmicFourierTransform(unsigned N, Real k, Real Δτ, unsigned pad) : N(N), pad(pad), k(k), Δτ(Δτ) {
   τₛ = -0.5 * N;
   ωₛ = -0.5 * N;
-  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");
+  sₛ = -0.5 * N;
+  Δs = 0.1;
+  Δω = Δτ;
 }
 
 LogarithmicFourierTransform::~LogarithmicFourierTransform() {
-  FFTW_DESTROY_PLAN(a_to_â);
-  FFTW_DESTROY_PLAN(â_to_a);
-  FFTW_FREE(a);
-  FFTW_FREE(â);
-  FFTW_CLEANUP();
 }
 
 Real LogarithmicFourierTransform::τ(unsigned n) const {
@@ -33,7 +24,7 @@ Real LogarithmicFourierTransform::ω(unsigned n) const {
 }
 
 Real LogarithmicFourierTransform::s(unsigned n) const {
-  return (n + sₛ) * 2*M_PI / (pad * N * Δτ);
+  return (n + sₛ) * Δs;
 }
 
 Real LogarithmicFourierTransform::t(unsigned n) const {
@@ -84,20 +75,26 @@ 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) {
-    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++) {
-      â[(N / 2 + n) % N] *= std::pow(-II * σ, II * s(n) - k) * Γ(k - II * s(n));
+    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));
+      }
     }
-    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 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));
     }
   }
 
+  for (unsigned i = 0; i < N; i++) {
+    c[i] -= c[N-1];
+  }
+
   return c;
 }