Work to debug log-fourier

1 files changed, 6 insertions, 22 deletions

diff --git a/log-fourier.cpp b/log-fourier.cpp
index 9d1f2cb..bc2dd87 100644
--- a/log-fourier.cpp
+++ b/log-fourier.cpp
@@ -64,15 +64,13 @@ std::vector<Complex> LogarithmicFourierTransform::fourier(const std::vector<Real
     for (unsigned n = 0; n < pad*N; n++) {
       if (n < N) {
         a[n] = c[n] * exp((1 - k) * τ(n));
-      } else if (n >= (pad - 1) * N) {
-        a[n] = c[pad*N-n-1] * exp((1 - k) * τ(pad*N-n-1));
       } else {
         a[n] = 0;
       }
     }
     FFTW_EXECUTE(a_to_â);
     for (unsigned n = 0; n < pad*N; n++) {
-      â[(pad*N / 2 + n) % (pad*N)] *= std::pow(II * σ, II * s(n) - k) * Γ(k - II * s(n));
+      â[(pad*N / 2 + n) % (pad*N)] *= std::exp(II * s(n)) * std::pow(II * σ, II * s(n) - k) * Γ(k - II * s(n));
     }
     FFTW_EXECUTE(â_to_a);
     for (unsigned n = 0; n < N; n++) {
@@ -87,30 +85,16 @@ 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 < pad * N; n++) {
-      if (n < N) {
-        if (σ < 0) {
-          a[n] = std::conj(ĉ[n]) * std::exp((1 - k) * ω(n));
-        } else {
-          a[n] = ĉ[n] * std::exp((1 - k) * ω(n));
-        }
-      } else if (n >= (pad - 1) * N) {
-        if (σ < 0) {
-          a[n] = ĉ[pad*N-n-1] * std::exp((1 - k) * ω(pad*N-n-1));
-        } else {
-          a[n] = std::conj(ĉ[pad*N-n-1]) * std::exp((1 - k) * ω(pad*N-n-1));
-        }
-      } else {
-        a[n] = 0;
-      }
+    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 < pad*N; n++) {
-      â[(pad*N / 2 + n) % (pad*N)] *= std::pow(-II * σ, II * s(n) - k) * Γ(k - II * s(n));
+    for (unsigned n = 0; n < N; 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[(pad - 1)*N+n].real() / (Real)(pad*N) / (2 * M_PI);
+      c[n] += std::exp(-k * τ(n)) * a[n].real() / (2 * M_PI * N);
     }
   }