来自 1986 的 DFT 改进算法（高阶插值实现）

import numpy as np

def fourier_quad_interp(x, dt, freqs, angular_freq=False):
    """
    Piecewise-quadratic interpolation FT.

    每个区间 [t_n, t_{n+1}] 上把信号近似成二次多项式：
        x(t) ≈ a_n + b_n * τ + c_n * τ^2
        其中 τ = t - t_n, 0 <= τ <= dt

    系数的求法：
      - 内部段 (1..N-3): 用 x[n-1], x[n], x[n+1] 三点
         在局部坐标 τ = -dt, 0, dt 上拟合二次多项式
      - 第 0 段: 用 x[0], x[1], x[2] (τ = 0, dt, 2dt)
      - 最后一段: 用 x[N-3], x[N-2], x[N-1] (τ = -dt, 0, dt)

    然后对每段做解析积分，再叠加得到频谱。
    """
    x = np.asarray(x, dtype=np.complex128)
    freqs = np.asarray(freqs, dtype=np.float64)

    if x.ndim != 1:
        raise ValueError("x must be 1-D array of samples.")
    if x.size < 3:
        raise ValueError("Need at least 3 samples for quadratic interpolation.")

    if angular_freq:
        omega = freqs
    else:
        omega = 2.0 * np.pi * freqs

    N = x.size
    n_seg = N - 1

    a = np.zeros(n_seg, dtype=np.complex128)
    b = np.zeros(n_seg, dtype=np.complex128)
    c = np.zeros(n_seg, dtype=np.complex128)

    # --- 段 0：用 τ=0, dt, 2dt 上的三点 (x0,x1,x2) 拟合 ---
    x0 = x[0]
    x1 = x[1]
    x2 = x[2]
    a[0] = x0
    u = (x1 - a[0]) / dt
    v = (x2 - a[0]) / dt
    c[0] = (v - 2 * u) / (2 * dt)
    b[0] = u - c[0] * dt

    # --- 中间段：用 τ=-dt, 0, dt 上的三点 (x[n-1],x[n],x[n+1]) 拟合 ---
    for n in range(1, n_seg - 1):
        x_minus = x[n - 1]
        x0n = x[n]
        x_plus = x[n + 1]
        a[n] = x0n
        c[n] = (x_plus + x_minus - 2 * x0n) / (2 * dt**2)
        b[n] = (x_plus - x_minus) / (2 * dt)

    # --- 最后一段：同样用 (x[N-3], x[N-2], x[N-1])，τ=-dt,0,dt ---
    n = n_seg - 1
    x_minus = x[N - 3]
    x0n = x[N - 2]
    x_plus = x[N - 1]
    a[n] = x0n
    c[n] = (x_plus + x_minus - 2 * x0n) / (2 * dt**2)
    b[n] = (x_plus - x_minus) / (2 * dt)

    t_n = dt * np.arange(n_seg, dtype=np.float64)

    X = np.zeros_like(omega, dtype=np.complex128)

    eps = 1e-12
    mask_zero = np.abs(omega) < eps
    mask_nonzero = ~mask_zero

    # -------- ω ≈ 0 时用精确极限 --------
    if np.any(mask_zero):
        # ∫_0^dt (a + bτ + cτ^2)dτ = a dt + b dt^2/2 + c dt^3/3
        seg_int = a * dt + b * (dt**2 / 2.0) + c * (dt**3 / 3.0)
        X[mask_zero] = np.sum(seg_int)

    # -------- ω ≠ 0 时用闭式表达式 --------
    if np.any(mask_nonzero):
        w = omega[mask_nonzero]

        # A(w) = ∫_0^dt e^{-j w τ} dτ
        A = (1.0 - np.exp(-1j * w * dt)) / (1j * w)
        # B(w) = ∫_0^dt τ e^{-j w τ} dτ
        B = (((1j * dt * w + 1.0) * np.exp(-1j * w * dt)) - 1.0) / (w**2)
        # C(w) = ∫_0^dt τ^2 e^{-j w τ} dτ
        # 用符号积分推出来的结果：
        # C = [ (j dt^2 w^2 + 2 dt w - 2j) e^{-j w dt} + 2j ] / w^3
        C = (((1j * dt**2 * w**2 + 2 * dt * w - 2j) * np.exp(-1j * w * dt) + 2j)
             / (w**3))

        A = A[:, None]
        B = B[:, None]
        C = C[:, None]

        a_col = a[None, :]
        b_col = b[None, :]
        c_col = c[None, :]
        t_col = t_n[None, :]

        exp_term = np.exp(-1j * w[:, None] * t_col)
        seg_contrib = exp_term * (a_col * A + b_col * B + c_col * C)

        X[mask_nonzero] = np.sum(seg_contrib, axis=1)

    return X