1 files changed, 75 insertions, 2 deletions
diff --git a/fourier.cpp b/fourier.cpp
index bc4b633..cf45ad1 100644
--- a/fourier.cpp
+++ b/fourier.cpp
@@ -1,4 +1,5 @@
 #include "fourier.hpp"
+#include <boost/multiprecision/fwd.hpp>
 
 inline Real fP(unsigned p, Real q) {
   return 0.5 * pow(q, p);
@@ -27,8 +28,8 @@ Real ddf(Real λ, unsigned p, unsigned s, Real q) {
 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_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, reinterpret_cast<fftw_complex*>(â), flags);
   plan_c2r = fftw_plan_dft_c2r_1d(2 * n, reinterpret_cast<fftw_complex*>(â), a, flags);
@@ -99,6 +100,78 @@ void FourierTransform::writeToA(unsigned i, Real ai) {
   a[i] = ai;
 }
 
+LogarithmicFourierTransform::LogarithmicFourierTransform(unsigned N, Real k, Real Δω, Real Δτ, Real Δs) : N(N), k(k), Δτ(Δτ), Δω(Δω), Δs(Δs) {
+  τₛ = -0.5 * N;
+  ωₛ = -0.5 * N;
+  sₛ = -0.5 * N;
+  a = reinterpret_cast<Complex*>(fftw_alloc_complex(N));
+  â = reinterpret_cast<Complex*>(fftw_alloc_complex(N));
+  fftw_import_wisdom_from_filename("fftw.wisdom");
+  a_to_â = fftw_plan_dft_1d(N, reinterpret_cast<fftw_complex*>(a), reinterpret_cast<fftw_complex*>(â), 1, 0);
+  â_to_a = fftw_plan_dft_1d(N, reinterpret_cast<fftw_complex*>(â), reinterpret_cast<fftw_complex*>(a), 1, 0);
+  fftw_export_wisdom_to_filename("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 {
+  return Δτ * (n + τₛ);
+}
+
+Real LogarithmicFourierTransform::ω(unsigned n) const {
+  return Δω * (n + ωₛ);
+}
+
+Real LogarithmicFourierTransform::s(unsigned n) const {
+  return Δs * (n + sₛ);
+}
+
+Real LogarithmicFourierTransform::t(unsigned n) const {
+  return exp(τ(n));
+}
+
+Real LogarithmicFourierTransform::ν(unsigned n) const {
+  return exp(ω(n));
+}
+
+Complex gamma(Complex z) {
+  gsl_sf_result logΓ;
+  gsl_sf_result argΓ;
+
+  gsl_sf_lngamma_complex_e(z.real(), z.imag(), &logΓ, &argΓ);
+
+  return exp(logΓ.val + 1i * argΓ.val);
+}
+
+std::vector<Complex> LogarithmicFourierTransform::fourier(const std::vector<Real>& c, bool symmetric){ 
+  std::vector<Complex> ĉ(N);
+  std::vector<Real> σs = {1};
+  if (symmetric){
+    σs.push_back(-1);
+  }
+  for (Real σ : σs) {
+    for (unsigned n = 0; n < N; n++) {
+      a[n] = c[n] * exp((1 - k) * τ(n));
+    }
+    fftw_execute(a_to_â);
+    for (unsigned n = 0; n < N; n++) {
+      â[(n + N / 2) % N] *= pow(-1i * σ, 1i * (s(n)) - k) * gamma(k - 1i * (s(n)));
+    }
+    fftw_execute(â_to_a);
+    for (unsigned n = 0; n < N; n++) {
+      ĉ[n] += exp(-k * ω(n) * 4 * M_PI) * a[n] / pow(2 * M_PI, 1);
+    }
+  }
+
+  return ĉ;
+}
+
 std::string fourierFile(std::string prefix, unsigned p, unsigned s, Real λ, Real τ₀, Real y, unsigned log2n, Real τₘₐₓ) {
   return prefix + "_" + std::to_string(p) + "_" + std::to_string(s) + "_" + std::to_string(λ) + "_" + std::to_string(τ₀) + "_" + std::to_string(y) + "_" + std::to_string(log2n) + "_" + std::to_string(τₘₐₓ) + ".dat";
 }