diff options
Diffstat (limited to 'log-fourier.cpp')
| -rw-r--r-- | log-fourier.cpp | 14 |
1 files changed, 10 insertions, 4 deletions
diff --git a/log-fourier.cpp b/log-fourier.cpp index 1e2e7d9..7461a70 100644 --- a/log-fourier.cpp +++ b/log-fourier.cpp @@ -1,4 +1,5 @@ #include "log-fourier.hpp" +#include <complex> LogarithmicFourierTransform::LogarithmicFourierTransform(unsigned N, Real k, Real Δτ, unsigned pad) : N(N), pad(pad), k(k), Δτ(Δτ) { τₛ = -0.5 * N; @@ -40,7 +41,7 @@ Real LogarithmicFourierTransform::ν(unsigned n) const { return exp(ω(n)); } -Complex gamma(Complex z) { +Complex Γ(Complex z) { gsl_sf_result logΓ; gsl_sf_result argΓ; @@ -52,6 +53,7 @@ Complex gamma(Complex z) { std::vector<Complex> LogarithmicFourierTransform::fourier(const std::vector<Real>& c, bool symmetric) { std::vector<Complex> ĉ(N); std::vector<Real> σs = {1}; + /* c is either even or zero for negative arguments */ if (symmetric){ σs.push_back(-1); } @@ -65,7 +67,7 @@ std::vector<Complex> LogarithmicFourierTransform::fourier(const std::vector<Real } fftw_execute(a_to_â); for (unsigned n = 0; n < pad*N; n++) { - â[(pad*N / 2 + n) % (pad*N)] *= pow(1i * σ, 1i * s(n) - k) * gamma(k - 1i * s(n)); + â[(pad*N / 2 + n) % (pad*N)] *= pow(1i * σ, 1i * s(n) - k) * Γ(k - 1i * s(n)); } fftw_execute(â_to_a); for (unsigned n = 0; n < N; n++) { @@ -82,14 +84,18 @@ std::vector<Real> LogarithmicFourierTransform::inverse(const std::vector<Complex for (Real σ : σs) { for (unsigned n = 0; n < pad * N; n++) { if (n < N) { - a[n] = ĉ[n] * exp((1 - k) * ω(n)); + if (σ < 0) { + a[n] = std::conj(ĉ[n]) * exp((1 - k) * ω(n)); + } else { + a[n] = ĉ[n] * exp((1 - k) * ω(n)); + } } else { a[n] = 0; } } fftw_execute(a_to_â); for (unsigned n = 0; n < pad*N; n++) { - â[(pad*N / 2 + n) % (pad*N)] *= pow(-1i * σ, 1i * s(n) - k) * gamma(k - 1i * s(n)); + â[(pad*N / 2 + n) % (pad*N)] *= pow(-1i * σ, 1i * s(n) - k) * Γ(k - 1i * s(n)); } fftw_execute(â_to_a); for (unsigned n = 0; n < N; n++) { |