diff options
| -rw-r--r-- | log-fourier.cpp | 18 | ||||
| -rw-r--r-- | log-fourier.hpp | 4 |
2 files changed, 15 insertions, 7 deletions
diff --git a/log-fourier.cpp b/log-fourier.cpp index e6a84fe..cb57229 100644 --- a/log-fourier.cpp +++ b/log-fourier.cpp @@ -13,7 +13,7 @@ Complex Γ(Complex z) { return std::exp((Real)logΓ.val + II * (Real)argΓ.val); } -LogarithmicFourierTransform::LogarithmicFourierTransform(unsigned N, Real k, Real Δτ, unsigned pad, Real shift) : N(N), pad(pad), k(k), Δτ(Δτ), ts(N), νs(N), Γs(pad * N), shift(shift) { +LogarithmicFourierTransform::LogarithmicFourierTransform(unsigned N, Real k, Real Δτ, unsigned pad, Real shift) : N(N), pad(pad), k(k), Δτ(Δτ), ts(N), νs(N), Γs(pad * N), exp1kω(N), exp1kτ(N), expkω(N), expkτ(N), shift(shift) { τₛ = -0.5 * N; ωₛ = -0.5 * N; sₛ = -0.5 * pad * N; @@ -26,6 +26,10 @@ LogarithmicFourierTransform::LogarithmicFourierTransform(unsigned N, Real k, Rea for (unsigned n = 0; n < N; n++) { ts[n] = std::exp(τ(n)) / shift; νs[n] = std::exp(ω(n)) * shift; + exp1kτ[n] = std::exp((1 - k) * τ(n)); + exp1kω[n] = std::exp((1 - k) * ω(n)); + expkτ[n] = std::exp(-k * τ(n)); + expkω[n] = std::exp(-k * ω(n)); } for (unsigned n = 0; n < pad * N; n++) { Γs[n] = Γ(k - II * s(n)); @@ -70,9 +74,9 @@ std::vector<Complex> LogarithmicFourierTransform::fourier(const std::vector<Real for (Real σ : σs) { for (unsigned n = 0; n < pad*N; n++) { if (n < N) { - a[n] = c[n] * std::exp((1 - k) * τ(n)); + a[n] = c[n] * exp1kτ[n]; } else if (n >= (pad - 1) * N) { - a[n] = c[pad*N-n-1] * std::exp((1 - k) * τ(pad*N-n-1)); + a[n] = c[pad*N-n-1] * exp1kτ[pad*N-n-1]; } else { a[n] = 0; } @@ -83,7 +87,7 @@ std::vector<Complex> LogarithmicFourierTransform::fourier(const std::vector<Real } FFTW_EXECUTE(â_to_a); for (unsigned n = 0; n < N; n++) { - ĉ[n] += std::exp(-k * ω(n)) * a[(pad - 1)*N+n] / (Real)(pad*N); + ĉ[n] += expkω[n] * a[(pad - 1)*N+n] / (Real)(pad*N); } } @@ -102,9 +106,9 @@ 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].real() + II * σ * ĉ[n].imag()) * std::exp((1 - k) * ω(n)); + a[n] = (ĉ[n].real() + II * σ * ĉ[n].imag()) * exp1kω[n]; } else if (n >= (pad - 1) * N) { - a[n] = (ĉ[pad*N-n-1].real() - II * σ * ĉ[pad*N-n-1].imag()) * std::exp((1 - k) * ω(pad*N-n-1)); + a[n] = (ĉ[pad*N-n-1].real() - II * σ * ĉ[pad*N-n-1].imag()) * exp1kω[pad*N-n-1]; } else { a[n] = 0; } @@ -115,7 +119,7 @@ std::vector<Real> LogarithmicFourierTransform::inverse(const std::vector<Complex } 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] += expkτ[n] * a[(pad - 1)*N+n].real() / (Real)(pad*N) / (2 * M_PI); } } diff --git a/log-fourier.hpp b/log-fourier.hpp index 730639c..3abc379 100644 --- a/log-fourier.hpp +++ b/log-fourier.hpp @@ -23,6 +23,10 @@ private: std::vector<Real> ts; std::vector<Real> νs; std::vector<Complex> Γs; + std::vector<Real> exp1kω; + std::vector<Real> exp1kτ; + std::vector<Real> expkω; + std::vector<Real> expkτ; public: Real shift; LogarithmicFourierTransform(unsigned N, Real k, Real Δτ, unsigned pad = 4, Real shift = 0.5); |