// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later
//
// example_doc_fft_samples.cpp
// Verification harness for samples in api/FFT.html (ja+en).

#define _USE_MATH_DEFINES
#include <math/fft/fft.hpp>
#include <math/fft/fft2d.hpp>
#include <math/fft/fft_utils.hpp>
#include <math/core/matrix.hpp>
#include <math/core/mp/Float.hpp>
#include <cmath>
#include <complex>
#include <iomanip>
#include <iostream>
#include <vector>

int main() {
    std::cout << std::setprecision(15);

    using C = sangi::Complex<double>;  // NOTE: sangi::FFT requires sangi::Complex
                                        // (NOT std::complex — see doc note)

    // ---- basic 8-point FFT round-trip ----------------------------
    {
        std::cout << "[basic FFT 8-pt]\n";
        std::vector<C> signal = {C(1), C(2), C(3), C(4), C(4), C(3), C(2), C(1)};
        std::vector<C> orig = signal;
        sangi::FFT<C>::fft(signal);
        std::cout << "  X[0].re=" << signal[0].re << "   (sum of samples = 20)\n";
        std::cout << "  X[1].re=" << signal[1].re << "  .im=" << signal[1].im << "\n";
        std::cout << "  X[4].re=" << signal[4].re << "   (Nyquist bin)\n";
        sangi::FFT<C>::ifft(signal);
        double err = 0;
        for (std::size_t i = 0; i < signal.size(); ++i) {
            err += std::abs(signal[i].re - orig[i].re) + std::abs(signal[i].im - orig[i].im);
        }
        std::cout << "  round-trip L1 error = " << err << "\n";
    }

    // ---- spectrum peak (440 Hz @ 44.1k, N=1024) ------------------
    {
        std::cout << "[spectrum analysis 440Hz]\n";
        constexpr int N = 1024;
        constexpr double fs = 44100.0;
        constexpr double f0 = 440.0;
        std::vector<double> signal(N);
        for (int i = 0; i < N; ++i)
            signal[i] = std::sin(2.0 * M_PI * f0 * i / fs);
        auto spectrum = sangi::RealFFT<double>::fft(signal);
        auto amp = sangi::amplitude_spectrum<double>(spectrum);
        int peak = 0;
        for (int k = 1; k < (int)amp.size(); ++k)
            if (amp[k] > amp[peak]) peak = k;
        double peakFreq = peak * fs / N;
        std::cout << "  peak bin = " << peak << "  peakFreq = " << peakFreq << " Hz\n";
        std::cout << "  (expected ~ 440 Hz; bin resolution = "
                  << (fs / N) << " Hz so ~ 430.66 Hz is the closest)\n";
    }

    // ---- convolution (polynomial multiplication) ----------------
    {
        std::cout << "[convolve]\n";
        std::vector<C> p = {C(1), C(2), C(3)};
        std::vector<C> q = {C(4), C(5)};
        auto r = sangi::FFT<C>::convolve(p, q);
        std::cout << "  conv = [";
        for (std::size_t i = 0; i < r.size(); ++i) {
            if (i) std::cout << ", ";
            std::cout << r[i].re;
        }
        std::cout << "]   (expected 4, 13, 22, 15)\n";
    }

    // ---- 2D FFT round-trip ---------------------------------------
    {
        std::cout << "[2D FFT 4x4]\n";
        using SC = sangi::Complex<double>;
        sangi::Matrix<SC> img(4, 4);
        for (int r = 0; r < 4; ++r)
            for (int c = 0; c < 4; ++c)
                img(r, c) = SC(double(r + c), 0.0);
        auto freq = sangi::fft2d(img);
        auto rec = sangi::ifft2d(freq);
        std::cout << "  Recovered (0,0).re = " << rec(0, 0).re
                  << "   (expected 0)\n";
        std::cout << "  Recovered (3,3).re = " << rec(3, 3).re
                  << "   (expected 6)\n";
    }

    // ---- DCT round-trip ------------------------------------------
    {
        std::cout << "[DCT 8-pt]\n";
        std::vector<double> data = {1, 2, 3, 4, 5, 6, 7, 8};
        auto coeffs = sangi::dct<double>(data);
        auto rec = sangi::idct<double>(coeffs);
        std::cout << "  coeffs[0..2] = " << coeffs[0] << ", "
                  << coeffs[1] << ", " << coeffs[2] << "\n";
        double err = 0;
        for (std::size_t i = 0; i < data.size(); ++i) err += std::abs(rec[i] - data[i]);
        std::cout << "  round-trip L1 error = " << err << "\n";
    }

    // ---- FFT with multi-precision Float (sangi::Float) ----------
    {
        std::cout << "[FFT<Complex<Float>, Float> 8-pt]\n";
        using sangi::Float;
        Float::setDefaultPrecision(170);  // ≈ 50 decimal digits
        using FC = sangi::Complex<Float>;
        std::vector<FC> signal = {FC(Float(1)), FC(Float(2)), FC(Float(3)), FC(Float(4)),
                                   FC(Float(4)), FC(Float(3)), FC(Float(2)), FC(Float(1))};
        std::vector<FC> orig = signal;
        sangi::FFT<FC, Float>::fft(signal);
        std::cout << "  X[0].re = " << signal[0].re.toDecimalString(15) << "   (sum=20)\n";
        std::cout << "  X[1].re = " << signal[1].re.toDecimalString(15) << "\n";
        std::cout << "  X[1].im = " << signal[1].im.toDecimalString(15) << "\n";
        std::cout << "  X[4].re = " << signal[4].re.toDecimalString(15) << "   (Nyquist)\n";
        sangi::FFT<FC, Float>::ifft(signal);
        Float err(Float(0));
        for (std::size_t i = 0; i < signal.size(); ++i) {
            Float dr = signal[i].re - orig[i].re;
            Float di = signal[i].im - orig[i].im;
            if (dr.isNegative()) dr = -dr;
            if (di.isNegative()) di = -di;
            err = err + dr + di;
        }
        std::cout << "  round-trip L1 err = " << err.toDecimalString(8)
                  << "   (multi-precision rounding)\n";
    }

    // ---- Hilbert / instantaneous amplitude ----------------------
    {
        std::cout << "[Hilbert envelope]\n";
        constexpr int N = 256;
        std::vector<double> signal(N);
        for (int i = 0; i < N; ++i) {
            double t = double(i) / N;
            signal[i] = (1.0 + 0.5 * std::cos(2 * M_PI * 5 * t))
                       * std::cos(2 * M_PI * 50 * t);
        }
        auto env = sangi::instantaneousAmplitude<double>(signal);
        std::cout << "  env[0]   = " << env[0] << "   (modulator at t=0: 1.5)\n";
        std::cout << "  env[N/2] = " << env[N / 2] << "\n";
    }

    return 0;
}