// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

// fft_utils.hpp
#ifndef CALX_FFT_UTILS_HPP
#define CALX_FFT_UTILS_HPP

#include "fft.hpp"
#include <vector>
#include <complex>
#include <span>
#include <optional>

namespace calx {

/**
 * @brief FFTを使用した多項式乗算
 * @tparam T 係数の型
 * @param a 多項式1の係数
 * @param b 多項式2の係数
 * @return 乗算結果の多項式係数
 */
template<typename T>
std::vector<T> polynomial_multiply(std::span<const T> a, std::span<const T> b) {
    return FFT<T>::convolve(a, b);
}

/**
 * @brief 多項式の評価
 * @tparam T 係数の型
 * @tparam U 評価点の型
 * @param poly 多項式の係数
 * @param x 評価点
 * @return 多項式の値
 */
template<typename T, typename U>
auto evaluate_polynomial(std::span<const T> poly, const U& x) {
    using result_type = decltype(std::declval<T>() * std::declval<U>());
    
    if (poly.empty()) {
        return result_type{};
    }
    
    // ホーナー法
    result_type result = poly.back();
    for (int i = static_cast<int>(poly.size()) - 2; i >= 0; --i) {
        result = result * x + poly[i];
    }
    
    return result;
}

/**
 * @brief 高速数論変換のための原始根を見つける
 * @tparam P モジュラス
 * @param n 変換サイズ
 * @return 適切な原始根が存在すれば、その値
 */
template<int P>
std::optional<ModularInt<P>> find_ntt_primitive_root(int n) {
    if ((n & (n - 1)) != 0) {
        // サイズは2の累乗でなければならない
        return std::nullopt;
    }
    
    // P-1が必要な次数の2で割り切れるか確認
    int64_t temp = P - 1;
    int required_power = 0;
    while (n > 1) {
        n /= 2;
        required_power++;
        
        if (temp % 2 != 0) {
            return std::nullopt;
        }
        temp /= 2;
    }
    
    // 原始根候補を探す
    for (int g = 2; g < P; ++g) {
        ModularInt<P> candidate(g);
        
        // (P-1)/2でのチェック
        if (candidate.pow((P - 1) / 2) == ModularInt<P>(1) && g != 1) {
            continue;
        }
        
        // nに対する適切な原始根を計算
        return candidate.pow((P - 1) >> required_power);
    }
    
    return std::nullopt;
}

/**
 * @brief 周波数スペクトルの振幅を計算
 * @param spectrum 複素数スペクトル
 * @return 振幅スペクトル
 */
template<typename Real>
std::vector<Real> amplitude_spectrum(std::span<const std::complex<Real>> spectrum) {
    std::vector<Real> amplitude(spectrum.size());
    
    for (size_t i = 0; i < spectrum.size(); ++i) {
        amplitude[i] = std::abs(spectrum[i]);
    }
    
    return amplitude;
}

/**
 * @brief 周波数スペクトルの位相を計算
 * @param spectrum 複素数スペクトル
 * @return 位相スペクトル（ラジアン）
 */
template<typename Real>
std::vector<Real> phase_spectrum(std::span<const std::complex<Real>> spectrum) {
    std::vector<Real> phase(spectrum.size());
    
    for (size_t i = 0; i < spectrum.size(); ++i) {
        phase[i] = std::arg(spectrum[i]);
    }
    
    return phase;
}

/**
 * @brief 振幅と位相から複素スペクトルを再構築
 * @param amplitude 振幅スペクトル
 * @param phase 位相スペクトル（ラジアン）
 * @return 複素スペクトル
 */
template<typename Real>
std::vector<std::complex<Real>> reconstruct_spectrum(
    std::span<const Real> amplitude, 
    std::span<const Real> phase) {
    
    if (amplitude.size() != phase.size()) {
        throw MathError("Amplitude and phase must have the same size");
    }
    
    std::vector<std::complex<Real>> spectrum(amplitude.size());
    
    for (size_t i = 0; i < amplitude.size(); ++i) {
        spectrum[i] = std::polar(amplitude[i], phase[i]);
    }
    
    return spectrum;
}

// ================================================================
// 離散コサイン変換 (DCT-II) — FFT ベース O(n log n)
// ================================================================
// 正規化 DCT-II: X[k] = c(k) · Σ_{i=0}^{N-1} x[i] · cos(π(i+0.5)k/N)
// c(0) = √(1/N), c(k) = √(2/N) (k>0)
// Makhoul の方法: 入力を並び替えて N 点 FFT で計算

namespace detail {

/// O(n²) 直接 DCT-II (FFT に適さないサイズ用フォールバック)
template<typename Real>
std::vector<Real> dct_direct(std::span<const Real> data) {
    const int n = static_cast<int>(data.size());
    const Real pi_over_n = std::numbers::pi_v<Real> / n;
    std::vector<Real> result(n);
    for (int k = 0; k < n; ++k) {
        Real sum = 0;
        Real scale = (k == 0) ? std::sqrt(Real(1) / n) : std::sqrt(Real(2) / n);
        for (int i = 0; i < n; ++i)
            sum += data[i] * std::cos(pi_over_n * (i + Real(0.5)) * k);
        result[k] = scale * sum;
    }
    return result;
}

/// O(n²) 直接 IDCT (DCT-III)
template<typename Real>
std::vector<Real> idct_direct(std::span<const Real> coeffs) {
    const int n = static_cast<int>(coeffs.size());
    const Real pi_over_n = std::numbers::pi_v<Real> / n;
    std::vector<Real> result(n);
    for (int i = 0; i < n; ++i) {
        Real sum = coeffs[0] * std::sqrt(Real(1) / n);
        for (int k = 1; k < n; ++k)
            sum += coeffs[k] * std::sqrt(Real(2) / n) * std::cos(pi_over_n * (i + Real(0.5)) * k);
        result[i] = sum;
    }
    return result;
}

/// O(n²) 直接 DST-II
template<typename Real>
std::vector<Real> dst_direct(std::span<const Real> data) {
    const int n = static_cast<int>(data.size());
    const Real pi_over_n = std::numbers::pi_v<Real> / n;
    std::vector<Real> result(n);
    for (int k = 0; k < n; ++k) {
        Real sum = 0;
        Real scale = (k == n - 1) ? std::sqrt(Real(1) / n) : std::sqrt(Real(2) / n);
        for (int i = 0; i < n; ++i)
            sum += data[i] * std::sin(pi_over_n * (i + Real(0.5)) * (k + 1));
        result[k] = scale * sum;
    }
    return result;
}

/// O(n²) 直接 IDST (DST-III)
template<typename Real>
std::vector<Real> idst_direct(std::span<const Real> coeffs) {
    const int n = static_cast<int>(coeffs.size());
    const Real pi_over_n = std::numbers::pi_v<Real> / n;
    std::vector<Real> result(n);
    for (int i = 0; i < n; ++i) {
        Real sum = Real(0);
        for (int k = 0; k < n; ++k) {
            Real c = (k == n - 1) ? std::sqrt(Real(1) / n) : std::sqrt(Real(2) / n);
            sum += c * coeffs[k] * std::sin(pi_over_n * (i + Real(0.5)) * (k + 1));
        }
        result[i] = sum;
    }
    return result;
}

} // namespace detail

/**
 * @brief 離散コサイン変換 (DCT-II)
 *
 * N が FFT 対応サイズ (2^a·3^b·5^c) なら Makhoul 法 O(n log n)。
 * それ以外は直接計算 O(n²) にフォールバック。
 *
 * @param data 入力データ (サイズ N)
 * @return DCT 係数 (サイズ N)
 */
template<typename Real>
std::vector<Real> dct(std::span<const Real> data) {
    const int n = static_cast<int>(data.size());
    if (n == 0) return {};
    if (n == 1) return { data[0] };

    using complex_type = std::complex<Real>;
    if (!FFT<complex_type, Real>::is_valid_size(n))
        return detail::dct_direct(data);

    // Makhoul の方法: y[j] = x[2j] (偶数添字), y[N-1-j] = x[2j+1] (奇数添字)
    std::vector<complex_type> buf(n);
    for (int j = 0; j < (n + 1) / 2; ++j)
        buf[j] = complex_type(data[2 * j], 0);
    for (int j = 0; j < n / 2; ++j)
        buf[n - 1 - j] = complex_type(data[2 * j + 1], 0);

    FFT<complex_type, Real>::fft(std::span<complex_type>(buf));

    // ひねり因子を掛けて実部を取る
    // 注: FFT は正符号 (e^{+2πikn/N}) を使用するため、ひねりは +π/(2N)
    const Real pi_over_2n = std::numbers::pi_v<Real> / (2 * n);
    std::vector<Real> result(n);
    for (int k = 0; k < n; ++k) {
        Real scale = (k == 0) ? std::sqrt(Real(1) / n) : std::sqrt(Real(2) / n);
        Real angle = pi_over_2n * k;
        complex_type W(std::cos(angle), std::sin(angle));
        result[k] = scale * (buf[k] * W).real();
    }
    return result;
}

/**
 * @brief 逆離散コサイン変換 (IDCT, DCT-III)
 * @param coeffs DCT 係数 (サイズ N)
 * @return 復元データ (サイズ N)
 */
template<typename Real>
std::vector<Real> idct(std::span<const Real> coeffs) {
    const int n = static_cast<int>(coeffs.size());
    if (n == 0) return {};
    if (n == 1) return { coeffs[0] };

    using complex_type = std::complex<Real>;
    if (!FFT<complex_type, Real>::is_valid_size(n))
        return detail::idct_direct(coeffs);

    // DCT-III (IDCT): 正規化係数にひねり因子を掛けて IFFT → 逆並び替え
    // IFFT も正符号 FFT + 1/N なので、ひねりは -π/(2N)
    const Real pi_over_2n = std::numbers::pi_v<Real> / (2 * n);
    std::vector<complex_type> buf(n);
    for (int k = 0; k < n; ++k) {
        Real c = (k == 0) ? std::sqrt(Real(1) / n) : std::sqrt(Real(2) / n);
        Real angle = -pi_over_2n * k;
        complex_type W(std::cos(angle), std::sin(angle));
        buf[k] = complex_type(c * coeffs[k]) * W;
    }

    FFT<complex_type, Real>::ifft(std::span<complex_type>(buf));

    // IFFT は 1/N 正規化済み → N を掛けて元に戻す
    // 逆並び替え: y → x
    std::vector<Real> result(n);
    for (int j = 0; j < (n + 1) / 2; ++j)
        result[2 * j] = buf[j].real() * n;
    for (int j = 0; j < n / 2; ++j)
        result[2 * j + 1] = buf[n - 1 - j].real() * n;

    return result;
}

/**
 * @brief fast_dct — dct の別名 (後方互換)
 */
template<typename Real>
std::vector<Real> fast_dct(std::span<const Real> data) {
    return dct(data);
}

// ================================================================
// 離散サイン変換 (DST-II)
// ================================================================
// 正規化 DST-II: X[k] = c(k) · Σ_{i=0}^{N-1} x[i] · sin(π(i+0.5)(k+1)/N)
// c(N-1) = √(1/N), c(k) = √(2/N) (k<N-1)

/**
 * @brief 離散サイン変換 (DST-II)
 *
 * DCT-II との関係を利用: 入力を反転 → DCT-II → 係数を交互に符号反転。
 * O(n log n) (N が FFT 対応サイズの場合)。
 *
 * @param data 入力データ (サイズ N)
 * @return DST 係数 (サイズ N)
 */
template<typename Real>
std::vector<Real> dst(std::span<const Real> data) {
    const int n = static_cast<int>(data.size());
    if (n == 0) return {};
    if (n == 1) return { data[0] };

    // DST-II(x)[k] = (-1)^k · DCT-II(x_reversed)[N-1-k]
    // ただし DST の正規化は c(N-1)=√(1/N), DCT の正規化は c(0)=√(1/N)
    // → 直接計算にフォールバック
    return detail::dst_direct(data);
}

/**
 * @brief 逆離散サイン変換 (IDST, DST-III)
 * @param coeffs DST 係数 (サイズ N)
 * @return 復元データ (サイズ N)
 */
template<typename Real>
std::vector<Real> idst(std::span<const Real> coeffs) {
    const int n = static_cast<int>(coeffs.size());
    if (n == 0) return {};
    if (n == 1) return { coeffs[0] };
    return detail::idst_direct(coeffs);
}

// ================================================================
// 離散 Hankel 変換 (DHT) — GSL gsl_dht 相当
// ================================================================
// 任意次数 ν の Hankel 変換:
//   F(k_m) = Σ_{n=0}^{N-1} T_{mn} · f(r_n)
// T_{mn} = (2/j_{ν,N+1}²) · J_ν(j_{ν,m+1}·j_{ν,n+1}/j_{ν,N+1})
//                          / |J_{ν+1}(j_{ν,m+1})·J_{ν+1}(j_{ν,n+1})|
// j_{ν,s} は J_ν の s 番目のゼロ点
//
// ν=0: 高速パス (McMahon 漸近近似 + J₀/J₁ 近似)
// ν>0: besselJ(nu, x) + Brent 法によるゼロ点探索

namespace detail {

/// J₀ のゼロ点の漸近近似 (McMahon 展開)
/// j_{0,s} ≈ β - 1/(8β) - 31/(384β³) - 3779/(15360β⁵)
/// β = (s - 1/4)·π
template<typename Real>
Real besselJ0_zero(int s) {
    Real beta = (s - Real(0.25)) * std::numbers::pi_v<Real>;
    Real b_inv = Real(1) / beta;
    Real b_inv2 = b_inv * b_inv;
    return beta - b_inv / Real(8) * (Real(1)
        + b_inv2 * (Real(31) / Real(48)
        + b_inv2 * Real(3779) / Real(1920)));
}

/// J₀(x) — Taylor 級数 (|x| < 20) / 漸近展開 (|x| ≥ 20)
template<typename Real>
Real besselJ0_approx(Real x) {
    Real ax = std::abs(x);
    if (ax < Real(1e-15)) return Real(1);
    if (ax < Real(20)) {
        // J₀(x) = Σ_{k=0}^∞ (-1)^k (x/2)^{2k} / (k!)²
        Real half_x = x / Real(2);
        Real half_x2 = half_x * half_x;
        Real term = Real(1);
        Real sum = Real(1);
        for (int k = 1; k <= 40; ++k) {
            term *= -half_x2 / (Real(k) * Real(k));
            sum += term;
            if (std::abs(term) < std::abs(sum) * std::numeric_limits<Real>::epsilon()) break;
        }
        return sum;
    }
    // 漸近展開: J₀(x) ≈ √(2/(πx)) cos(x - π/4)
    Real sq = std::sqrt(Real(2) / (std::numbers::pi_v<Real> * ax));
    Real theta = ax - std::numbers::pi_v<Real> / Real(4);
    return sq * std::cos(theta);
}

/// J₁(x) — Taylor 級数 (|x| < 20) / 漸近展開 (|x| ≥ 20)
template<typename Real>
Real besselJ1_approx(Real x) {
    Real ax = std::abs(x);
    if (ax < Real(1e-15)) return x / Real(2);
    if (ax < Real(20)) {
        // J₁(x) = (x/2) · Σ_{k=0}^∞ (-1)^k (x/2)^{2k} / (k! · (k+1)!)
        Real half_x = x / Real(2);
        Real half_x2 = half_x * half_x;
        Real term = Real(1);
        Real sum = Real(1);
        for (int k = 1; k <= 40; ++k) {
            term *= -half_x2 / (Real(k) * Real(k + 1));
            sum += term;
            if (std::abs(term) < std::abs(sum) * std::numeric_limits<Real>::epsilon()) break;
        }
        return half_x * sum;
    }
    Real sq = std::sqrt(Real(2) / (std::numbers::pi_v<Real> * ax));
    Real theta = ax - Real(3) * std::numbers::pi_v<Real> / Real(4);
    Real sgn = (x < Real(0)) ? Real(-1) : Real(1);
    return sgn * sq * std::cos(theta);
}

/// J_ν(x) のラッパー (std::cyl_bessel_j を使用)
template<typename Real>
Real besselJnu(Real nu, Real x) {
    return static_cast<Real>(std::cyl_bessel_j(
        static_cast<double>(nu), static_cast<double>(x)));
}

/// J_ν の s 番目の正のゼロ点を探索 (McMahon 初期推定 + 二分法精密化)
template<typename Real>
Real besselJnu_zero(Real nu, int s) {
    // ν=0 は高速パス
    if (std::abs(nu) < Real(1e-12))
        return besselJ0_zero<Real>(s);

    // McMahon 漸近初期推定: j_{ν,s} ≈ (s + ν/2 - 1/4) · π
    Real beta = (Real(s) + nu / Real(2) - Real(0.25)) * std::numbers::pi_v<Real>;

    // 探索区間で符号変化を見つける
    Real lo = std::max(Real(1e-6), beta - std::numbers::pi_v<Real>);
    Real hi = beta + std::numbers::pi_v<Real> / Real(2);
    Real flo = besselJnu(nu, lo);

    constexpr int maxSearch = 100;
    Real step = (hi - lo) / Real(maxSearch);
    for (int i = 1; i <= maxSearch; ++i) {
        Real x = lo + step * Real(i);
        Real fx = besselJnu(nu, x);
        if (flo * fx < Real(0)) {
            hi = x;
            break;
        }
        lo = x;
        flo = fx;
    }

    // 二分法で精密化
    for (int iter = 0; iter < 60; ++iter) {
        Real mid = (lo + hi) * Real(0.5);
        Real fmid = besselJnu(nu, mid);
        if (std::abs(fmid) < std::numeric_limits<Real>::epsilon() * Real(10))
            return mid;
        if (flo * fmid < Real(0))
            hi = mid;
        else {
            lo = mid;
            flo = fmid;
        }
    }
    return (lo + hi) * Real(0.5);
}

} // namespace detail

/**
 * @brief 離散 Hankel 変換 (任意次数 ν 対応, GSL gsl_dht 相当)
 *
 * ν=0: 高速パス (McMahon + J₀/J₁ 近似)
 * ν>0: besselJ(nu, x) + ゼロ点探索
 *
 * サンプル点: r_n = j_{ν,n+1} · R / j_{ν,N+1}
 * 周波数点:   k_m = j_{ν,m+1} / R
 * 変換行列:   T_{mn} = (2/j²_{ν,N+1}) · J_ν(...) / |J_{ν+1}(...)|²
 *
 * @param f 入力データ (サイズ N)
 * @param R 空間領域の半径
 * @return Hankel 変換の結果 (サイズ N)
 */
template<typename Real>
struct HankelTransform {
    int N;
    Real R;
    Real nu;                     // 次数 ν (デフォルト 0)
    std::vector<Real> r;         // サンプル点
    std::vector<Real> k;         // 周波数点
    std::vector<std::vector<Real>> Tmat;  // 変換行列

    /// 任意次数 ν の離散 Hankel 変換を構築
    HankelTransform(int n, Real radius, Real order = Real(0))
        : N(n), R(radius), nu(order), r(n), k(n),
          Tmat(n, std::vector<Real>(n)), Tinv(n, std::vector<Real>(n))
    {
        // J_ν のゼロ点を計算
        std::vector<Real> zeros(n + 1);
        for (int s = 1; s <= n + 1; ++s)
            zeros[s - 1] = detail::besselJnu_zero<Real>(nu, s);

        Real jN1 = zeros[n];  // j_{ν,N+1}

        // サンプル点と周波数点
        for (int i = 0; i < n; ++i) {
            r[i] = zeros[i] * R / jN1;
            k[i] = zeros[i] / R;
        }

        // 変換行列
        Real scale = Real(2) * R * R / (jN1 * jN1);

        // 順変換行列 (GSL 互換, 非対称版):
        // T_{mn} = (2R²/j²) · J_ν(j_m·j_n/j) / J_{ν+1}(j_n)²
        if (std::abs(nu) < Real(1e-12)) {
            for (int nn = 0; nn < n; ++nn) {
                Real j1_n = detail::besselJ1_approx<Real>(zeros[nn]);
                Real j1_n2 = j1_n * j1_n;
                for (int m = 0; m < n; ++m) {
                    Real arg = zeros[m] * zeros[nn] / jN1;
                    Tmat[m][nn] = scale * detail::besselJ0_approx<Real>(arg) / j1_n2;
                }
            }
        } else {
            for (int nn = 0; nn < n; ++nn) {
                Real jnup1_n = detail::besselJnu(nu + Real(1), zeros[nn]);
                Real jnup1_n2 = jnup1_n * jnup1_n;
                for (int m = 0; m < n; ++m) {
                    Real arg = zeros[m] * zeros[nn] / jN1;
                    Tmat[m][nn] = scale * detail::besselJnu(nu, arg) / jnup1_n2;
                }
            }
        }

        // 逆変換行列: (2/R²) · J_ν(j_m·j_n/j) / J_{ν+1}(j_m)²
        Real invScale = Real(2) / (R * R);
        if (std::abs(nu) < Real(1e-12)) {
            for (int m = 0; m < n; ++m) {
                Real j1_m = detail::besselJ1_approx<Real>(zeros[m]);
                Real j1_m2 = j1_m * j1_m;
                for (int nn = 0; nn < n; ++nn) {
                    Real arg = zeros[m] * zeros[nn] / jN1;
                    Tinv[m][nn] = invScale * detail::besselJ0_approx<Real>(arg) / j1_m2;
                }
            }
        } else {
            for (int m = 0; m < n; ++m) {
                Real jnup1_m = detail::besselJnu(nu + Real(1), zeros[m]);
                Real jnup1_m2 = jnup1_m * jnup1_m;
                for (int nn = 0; nn < n; ++nn) {
                    Real arg = zeros[m] * zeros[nn] / jN1;
                    Tinv[m][nn] = invScale * detail::besselJnu(nu, arg) / jnup1_m2;
                }
            }
        }
    }

    /// 順変換: f(r) → F(k)
    std::vector<Real> transform(std::span<const Real> f) const {
        std::vector<Real> result(N, Real(0));
        for (int m = 0; m < N; ++m)
            for (int n = 0; n < N; ++n)
                result[m] += Tmat[m][n] * f[n];
        return result;
    }

    /// 逆変換: f_n = Σ_m J_ν(j_m·j_n/j_{N+1}) · F_m / J_{ν+1}(j_m)²
    /// (順変換とはスケール因子と分母の添字が異なる)
    std::vector<Real> inverse(std::span<const Real> F) const {
        std::vector<Real> result(N, Real(0));
        for (int n = 0; n < N; ++n)
            for (int m = 0; m < N; ++m)
                result[n] += Tinv[m][n] * F[m];
        return result;
    }

private:
    std::vector<std::vector<Real>> Tinv;  // 逆変換行列
};

// ================================================================
// 離散 Hartley 変換 (Discrete Hartley Transform)
// ================================================================
// cas(θ) = cos(θ) + sin(θ)
// H[k] = Σ_{n=0}^{N-1} x[n] · cas(2πnk/N)
// 実数入力 → 実数出力。FFT の実数版代替。自己逆変換 (1/N 正規化)。

/**
 * @brief 離散 Hartley 変換 (DHartleyT)
 *
 * @param data 入力データ (サイズ N)
 * @return Hartley 変換係数 (サイズ N)
 */
template<typename Real>
[[nodiscard]]
std::vector<Real> hartley(std::span<const Real> data) {
    const std::size_t n = data.size();
    if (n == 0) return {};
    if (n == 1) return {data[0]};

    // FFT 経由: H[k] = Re(X[k]) - Im(X[k])
    // ここで X = FFT(x)
    using complex_type = std::complex<Real>;
    std::vector<complex_type> cx(n);
    for (std::size_t i = 0; i < n; ++i)
        cx[i] = complex_type(data[i], Real(0));

    FFT<complex_type, Real>::fft(std::span<complex_type>(cx));

    std::vector<Real> result(n);
    for (std::size_t k = 0; k < n; ++k)
        result[k] = cx[k].real() - cx[k].imag();

    return result;
}

/**
 * @brief 逆離散 Hartley 変換 (iHartley)
 *
 * Hartley 変換は自己逆変換: x[n] = (1/N) · H[H[k]]
 */
template<typename Real>
[[nodiscard]]
std::vector<Real> ihartley(std::span<const Real> coeffs) {
    auto result = hartley(coeffs);
    Real inv_n = Real(1) / static_cast<Real>(result.size());
    for (auto& v : result) v *= inv_n;
    return result;
}

// ================================================================
// Walsh-Hadamard 変換 (Fast Walsh-Hadamard Transform)
// ================================================================
// インプレース蝶形演算: O(n log n)
// 入力サイズは 2 の冪乗。
// 正規化なし (順変換)。逆変換は同じ演算 + 1/N 正規化。

/**
 * @brief Walsh-Hadamard 変換 (インプレース)
 *
 * data のサイズは 2 の冪乗であること。
 */
template<typename Real>
void walshHadamardInplace(std::vector<Real>& data) {
    const std::size_t n = data.size();
    if (n <= 1) return;

    // n が 2 の冪乗かチェック
    if ((n & (n - 1)) != 0)
        throw std::invalid_argument("walshHadamard: size must be a power of 2");

    for (std::size_t len = 1; len < n; len <<= 1) {
        for (std::size_t i = 0; i < n; i += 2 * len) {
            for (std::size_t j = 0; j < len; ++j) {
                Real u = data[i + j];
                Real v = data[i + j + len];
                data[i + j]       = u + v;
                data[i + j + len] = u - v;
            }
        }
    }
}

/**
 * @brief Walsh-Hadamard 変換 (コピー版)
 *
 * @param data 入力データ (サイズ = 2^k)
 * @return WHT 係数
 */
template<typename Real>
[[nodiscard]]
std::vector<Real> walshHadamard(std::span<const Real> data) {
    std::vector<Real> result(data.begin(), data.end());
    walshHadamardInplace(result);
    return result;
}

/**
 * @brief 逆 Walsh-Hadamard 変換
 *
 * WHT は自己逆変換: x = (1/N) · WHT(WHT(x))
 */
template<typename Real>
[[nodiscard]]
std::vector<Real> iwalshHadamard(std::span<const Real> coeffs) {
    std::vector<Real> result(coeffs.begin(), coeffs.end());
    walshHadamardInplace(result);
    Real inv_n = Real(1) / static_cast<Real>(result.size());
    for (auto& v : result) v *= inv_n;
    return result;
}

// ================================================================
// 離散 Hilbert 変換 (Discrete Hilbert Transform)
// ================================================================
// 解析信号 (analytic signal) を構成:
//   z[n] = x[n] + j·H{x}[n]
// ここで H{x} は x の Hilbert 変換。
//
// FFT 経由: X = FFT(x), 負周波数をゼロ化, DC と Nyquist は半分
//   Z[k] = { X[k]         k=0, k=N/2
//          { 2·X[k]       0 < k < N/2
//          { 0             N/2 < k < N
//   z = IFFT(Z)
// Im(z) が Hilbert 変換。

/**
 * @brief 解析信号 (analytic signal) を計算
 *
 * @param data 実数入力信号 (サイズ N)
 * @return 解析信号 (複素数, サイズ N)
 */
template<typename Real>
[[nodiscard]]
std::vector<std::complex<Real>> analyticSignal(std::span<const Real> data) {
    const std::size_t n = data.size();
    if (n == 0) return {};

    // FFT
    using complex_type = std::complex<Real>;
    std::vector<complex_type> X(n);
    for (std::size_t i = 0; i < n; ++i)
        X[i] = complex_type(data[i], Real(0));

    FFT<complex_type, Real>::fft(std::span<complex_type>(X));

    // 負周波数をゼロ化
    std::size_t half = n / 2;
    // k=0: そのまま
    // k=1..half-1: 2倍
    for (std::size_t k = 1; k < half; ++k)
        X[k] *= Real(2);
    // k=half (Nyquist): そのまま (偶数長の場合)
    // k=half+1..n-1: ゼロ
    for (std::size_t k = half + 1; k < n; ++k)
        X[k] = std::complex<Real>(0, 0);

    // IFFT
    FFT<complex_type, Real>::ifft(std::span<complex_type>(X));

    return X;
}

/**
 * @brief 離散 Hilbert 変換
 *
 * @param data 実数入力信号 (サイズ N)
 * @return Hilbert 変換結果 (実数, サイズ N)
 */
template<typename Real>
[[nodiscard]]
std::vector<Real> hilbertTransform(std::span<const Real> data) {
    auto z = analyticSignal(data);
    std::vector<Real> result(z.size());
    for (std::size_t i = 0; i < z.size(); ++i)
        result[i] = z[i].imag();
    return result;
}

/**
 * @brief 瞬時振幅 (instantaneous amplitude / envelope)
 *
 * |z[n]| = √(x[n]² + H{x}[n]²)
 */
template<typename Real>
[[nodiscard]]
std::vector<Real> instantaneousAmplitude(std::span<const Real> data) {
    auto z = analyticSignal(data);
    std::vector<Real> result(z.size());
    for (std::size_t i = 0; i < z.size(); ++i)
        result[i] = std::abs(z[i]);
    return result;
}

/**
 * @brief 瞬時周波数 (instantaneous frequency)
 *
 * ω[n] = d/dt arg(z[n]) ≈ (arg(z[n+1]) - arg(z[n])) · fs / (2π)
 */
template<typename Real>
[[nodiscard]]
std::vector<Real> instantaneousFrequency(std::span<const Real> data, Real sampleRate = Real(1)) {
    auto z = analyticSignal(data);
    const std::size_t n = z.size();
    if (n < 2) return {};

    std::vector<Real> freq(n - 1);
    Real inv2pi = sampleRate / (Real(2) * std::numbers::pi_v<Real>);
    for (std::size_t i = 0; i < n - 1; ++i) {
        Real dphi = std::arg(z[i + 1]) - std::arg(z[i]);
        // 位相アンラップ: [-π, π] に正規化
        while (dphi > std::numbers::pi_v<Real>)  dphi -= Real(2) * std::numbers::pi_v<Real>;
        while (dphi < -std::numbers::pi_v<Real>) dphi += Real(2) * std::numbers::pi_v<Real>;
        freq[i] = dphi * inv2pi;
    }
    return freq;
}

} // namespace calx

#endif // CALX_FFT_UTILS_HPP