From 18c33edc2fdf6abc9f8f36ea67b256d4a885493a Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Thu, 8 May 2025 17:32:13 -0300 Subject: Revert "Work to debug log-fourier" This reverts commit 77e8c94e5b41287e97e36e34a97478400637d102. --- log-fourier.cpp | 28 ++++++++++++++++++++++------ 1 file changed, 22 insertions(+), 6 deletions(-) (limited to 'log-fourier.cpp') diff --git a/log-fourier.cpp b/log-fourier.cpp index bc2dd87..9d1f2cb 100644 --- a/log-fourier.cpp +++ b/log-fourier.cpp @@ -64,13 +64,15 @@ std::vector LogarithmicFourierTransform::fourier(const std::vector= (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::exp(II * s(n)) * std::pow(II * σ, II * s(n) - k) * Γ(k - II * s(n)); + â[(pad*N / 2 + n) % (pad*N)] *= std::pow(II * σ, II * s(n) - k) * Γ(k - II * s(n)); } FFTW_EXECUTE(â_to_a); for (unsigned n = 0; n < N; n++) { @@ -85,16 +87,30 @@ std::vector LogarithmicFourierTransform::inverse(const std::vector c(N); std::vector σ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)); + 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; + } } 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)); + 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)); } FFTW_EXECUTE(â_to_a); for (unsigned n = 0; n < N; n++) { - c[n] += std::exp(-k * τ(n)) * a[n].real() / (2 * M_PI * N); + c[n] += std::exp(-k * τ(n)) * a[(pad - 1)*N+n].real() / (Real)(pad*N) / (2 * M_PI); } } -- cgit v1.2.3-70-g09d2