diff options
Diffstat (limited to 'log-fourier.cpp')
| -rw-r--r-- | log-fourier.cpp | 28 |
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); } } |