来自 1986 的 DFT 改进算法（YUNSWJ 数值分析实现版）

import numpy as np

def fourier_linear_interp(x, dt, freqs, angular_freq=False):
    """
    Improved DFT-like Fourier transform based on piecewise-linear interpolation.

    Parameters
    ----------
    x : array_like, shape (N,)
        Time-domain samples x[n] of a real or complex signal.
        Assumed uniformly spaced with sampling interval dt.
        You should provide N >= 2 samples.
    dt : float
        Sampling interval in seconds.
    freqs : array_like
        Frequencies at which to evaluate the Fourier transform.
        If angular_freq=False (default), this is in Hz.
        If angular_freq=True, this is in rad/s.
    angular_freq : bool, optional
        If True, `freqs` is interpreted as angular frequency ω (rad/s).
        If False (default), `freqs` is in Hz and ω = 2π f.

    Returns
    -------
    X : ndarray, shape (len(freqs),)
        Approximation to the continuous-time Fourier transform
        X(ω) = ∫ x(t) e^{-j ω t} dt
        computed using piecewise-linear interpolation between samples.

    Notes
    -----
    - This is *not* the same as the standard DFT / FFT.
      Standard DFT corresponds roughly to a rectangular-rule
      numerical integral; this function corresponds to integrating
      a piecewise-linear interpolant exactly.
    - For ω = 0, the formula is handled using exact limits:
        A(0) = Δt
        B(0) = Δt^2 / 2
      which avoids division-by-zero and improves numerical stability.
    """
    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 < 2:
        raise ValueError("Need at least 2 samples for linear interpolation.")

    # Angular frequency array
    if angular_freq:
        omega = freqs
    else:
        omega = 2.0 * np.pi * freqs  # rad/s

    N = x.size
    # Segment endpoints: x_n on [t_n, t_{n+1}]
    x0 = x[:-1]          # shape (N-1,)
    x1 = x[1:]           # shape (N-1,)
    a = x0               # intercept
    b = (x1 - x0) / dt   # slope
    t_n = dt * np.arange(N - 1, dtype=np.float64)  # left endpoints

    # Prepare output
    X = np.zeros_like(omega, dtype=np.complex128)

    # Handle ω = 0 separately (exact limits)
    eps = 1e-12  # small tolerance
    mask_zero = np.abs(omega) < eps
    mask_nonzero = ~mask_zero

    # ---- ω = 0 branch ----
    if np.any(mask_zero):
        # ∫ x(t) dt over each linear segment:
        # x(t) = a + b (t - t_n), t ∈ [t_n, t_n + dt]
        # Integral is: a*dt + b*(dt^2 / 2)  (shift does not matter for area)
        seg_int = a * dt + b * (dt**2 / 2.0)  # shape (N-1,)
        X[mask_zero] = np.sum(seg_int)

    # ---- ω ≠ 0 branch ----
    if np.any(mask_nonzero):
        w = omega[mask_nonzero]            # shape (M,)
        # Common factors in closed-form integrals over [0, dt]
        # A(w) = ∫_0^dt e^{-j w t} dt
        A = (1.0 - np.exp(-1j * w * dt)) / (1j * w)  # shape (M,)
        # B(w) = ∫_0^dt t e^{-j w t} dt
        # from sympy: ((j*dt*w + 1)*e^{-j*w*dt} - 1)/w^2
        B = (((1j * dt * w + 1.0) * np.exp(-1j * w * dt)) - 1.0) / (w**2)

        # Broadcast to (M, N-1)
        A = A[:, None]   # (M,1)
        B = B[:, None]   # (M,1)

        a_col = a[None, :]  # (1, N-1)
        b_col = b[None, :]  # (1, N-1)
        t_col = t_n[None, :]  # (1, N-1)

        # e^{-j w t_n} factor for each segment and frequency
        exp_term = np.exp(-1j * w[:, None] * t_col)  # (M, N-1)

        # Segment contributions:
        # I_n(w) = e^{-j w t_n} [ a_n A(w) + b_n B(w) ]
        seg_contrib = exp_term * (a_col * A + b_col * B)  # (M, N-1)

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

    return X

A = (1.0 - np.exp(-1j * w * dt)) / (1j * w)

B = (((1j * dt * w + 1.0) * np.exp(-1j * w * dt)) - 1.0) / (w**2)

a_col = a[None, :]    # (1, N-1)
b_col = b[None, :]
t_col = t_n[None, :]
exp_term = np.exp(-1j * w[:, None] * t_col)  # (M, N-1), M 个频点
seg_contrib = exp_term * (a_col * A + b_col * B)  # (M, N-1)
X[mask_nonzero] = np.sum(seg_contrib, axis=1)

seg_int = a * dt + b * (dt**2 / 2.0)
X[mask_zero] = np.sum(seg_int)

import numpy as np

def fourier_linear_interp(x, dt, freqs, angular_freq=False):
    """
    Improved DFT-like Fourier transform based on piecewise-linear interpolation.

    Parameters
    ----------
    x : array_like, shape (N,)
        Time-domain samples x[n] of a real or complex signal.
        Assumed uniformly spaced with sampling interval dt.
        You should provide N >= 2 samples.
    dt : float
        Sampling interval in seconds.
    freqs : array_like
        Frequencies at which to evaluate the Fourier transform.
        If angular_freq=False (default), this is in Hz.
        If angular_freq=True, this is in rad/s.
    angular_freq : bool, optional
        If True, `freqs` is interpreted as angular frequency ω (rad/s).
        If False (default), `freqs` is in Hz and ω = 2π f.

    Returns
    -------
    X : ndarray, shape (len(freqs),)
        Approximation to the continuous-time Fourier transform
        X(ω) = ∫ x(t) e^{-j ω t} dt
        computed using piecewise-linear interpolation between samples.

    Notes
    -----
    - This is *not* the same as the standard DFT / FFT.
      Standard DFT corresponds roughly to a rectangular-rule
      numerical integral; this function corresponds to integrating
      a piecewise-linear interpolant exactly.
    - For ω = 0, the formula is handled using exact limits:
        A(0) = Δt
        B(0) = Δt^2 / 2
      which avoids division-by-zero and improves numerical stability.
    """
    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 < 2:
        raise ValueError("Need at least 2 samples for linear interpolation.")

    # Angular frequency array
    if angular_freq:
        omega = freqs
    else:
        omega = 2.0 * np.pi * freqs  # rad/s

    N = x.size
    # Segment endpoints: x_n on [t_n, t_{n+1}]
    x0 = x[:-1]          # shape (N-1,)
    x1 = x[1:]           # shape (N-1,)
    a = x0               # intercept
    b = (x1 - x0) / dt   # slope
    t_n = dt * np.arange(N - 1, dtype=np.float64)  # left endpoints

    # Prepare output
    X = np.zeros_like(omega, dtype=np.complex128)

    # Handle ω = 0 separately (exact limits)
    eps = 1e-12  # small tolerance
    mask_zero = np.abs(omega) < eps
    mask_nonzero = ~mask_zero

    # ---- ω = 0 branch ----
    if np.any(mask_zero):
        # ∫ x(t) dt over each linear segment:
        # x(t) = a + b (t - t_n), t ∈ [t_n, t_n + dt]
        # Integral is: a*dt + b*(dt^2 / 2)  (shift does not matter for area)
        seg_int = a * dt + b * (dt**2 / 2.0)  # shape (N-1,)
        X[mask_zero] = np.sum(seg_int)

    # ---- ω ≠ 0 branch ----
    if np.any(mask_nonzero):
        w = omega[mask_nonzero]            # shape (M,)
        # Common factors in closed-form integrals over [0, dt]
        # A(w) = ∫_0^dt e^{-j w t} dt
        A = (1.0 - np.exp(-1j * w * dt)) / (1j * w)  # shape (M,)
        # B(w) = ∫_0^dt t e^{-j w t} dt
        # from sympy: ((j*dt*w + 1)*e^{-j*w*dt} - 1)/w^2
        B = (((1j * dt * w + 1.0) * np.exp(-1j * w * dt)) - 1.0) / (w**2)

        # Broadcast to (M, N-1)
        A = A[:, None]   # (M,1)
        B = B[:, None]   # (M,1)

        a_col = a[None, :]  # (1, N-1)
        b_col = b[None, :]  # (1, N-1)
        t_col = t_n[None, :]  # (1, N-1)

        # e^{-j w t_n} factor for each segment and frequency
        exp_term = np.exp(-1j * w[:, None] * t_col)  # (M, N-1)

        # Segment contributions:
        # I_n(w) = e^{-j w t_n} [ a_n A(w) + b_n B(w) ]
        seg_contrib = exp_term * (a_col * A + b_col * B)  # (M, N-1)

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

    return X

if __name__ == "__main__":
    import matplotlib.pyplot as plt

    # 构造一个简单信号：非整数周期的正弦 + 线性趋势
    fs = 1000.0   # 采样率 Hz
    dt = 1.0 / fs
    T = 1.0       # 观测长度 1 s
    t = np.arange(0, T, dt)
    f0 = 57.3     # 非整数 bin 频率
    x = 0.8 * np.sin(2*np.pi*f0*t) + 0.1*t

    # 频率轴：和 FFT 一样取 0..fs/2
    N = len(t)
    freqs = np.linspace(0, fs/2, 2048, endpoint=False)

    # 新算法（线性插值积分近似傅里叶变换）
    X_lin = fourier_linear_interp(x, dt, freqs, angular_freq=False)

    # 标准 FFT（只看正频）
    X_fft = np.fft.rfft(x) * dt               # 简单把和式当积分
    f_fft = np.fft.rfftfreq(N, dt)

    # 画幅度谱比较
    plt.figure()
    plt.semilogy(freqs, np.abs(X_lin), label="Linear-interp FT")
    plt.semilogy(f_fft, np.abs(X_fft), "--", label="Rectangular DFT/FFT")
    plt.xlim(0, 200)
    plt.xlabel("Frequency (Hz)")
    plt.ylabel("|X(f)|")
    plt.legend()
    plt.grid(True)
    plt.show()

import numpy as np
import matplotlib.pyplot as plt

def fourier_linear_interp(x, dt, freqs, angular_freq=False):
    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 < 2:
        raise ValueError("Need at least 2 samples for linear interpolation.")

    if angular_freq:
        omega = freqs
    else:
        omega = 2.0 * np.pi * freqs  # rad/s

    N = x.size
    x0 = x[:-1]
    x1 = x[1:]
    a = x0
    b = (x1 - x0) / dt
    t_n = dt * np.arange(N - 1, dtype=np.float64)

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

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

    if np.any(mask_zero):
        seg_int = a * dt + b * (dt**2 / 2.0)
        X[mask_zero] = np.sum(seg_int)

    if np.any(mask_nonzero):
        w = omega[mask_nonzero]
        A = (1.0 - np.exp(-1j * w * dt)) / (1j * w)
        B = (((1j * dt * w + 1.0) * np.exp(-1j * w * dt)) - 1.0) / (w**2)

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

        a_col = a[None, :]
        b_col = b[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)

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

    return X

# Common signal parameters
fs = 1000.0
dt = 1.0 / fs
N = 256
t = np.arange(N) * dt
f0 = 57.3  # non-bin-centered frequency
x_base = np.sin(2 * np.pi * f0 * t)

# === Figure 1: Main lobe comparison ===
freqs_dense = np.linspace(f0 - 20, f0 + 20, 2000)
X_lin_dense = fourier_linear_interp(x_base, dt, freqs_dense, angular_freq=False)

X_fft = np.fft.rfft(x_base) * dt
f_fft = np.fft.rfftfreq(N, dt)

plt.figure()
plt.plot(freqs_dense, np.abs(X_lin_dense), label="Linear-interp FT")
plt.plot(f_fft, np.abs(X_fft), linestyle="--", marker="", label="Rectangular FFT")
plt.xlim(f0 - 10, f0 + 10)
plt.xlabel("Frequency (Hz)")
plt.ylabel("Amplitude")
plt.title("Main lobe comparison (non-bin-centered sine)")
plt.grid(True)
plt.legend()
plt.show()

# === Figure 2: Leakage / sidelobe comparison (wider view, log scale) ===
freqs_wide = np.linspace(0, fs/2, 4000)
X_lin_wide = fourier_linear_interp(x_base, dt, freqs_wide, angular_freq=False)

plt.figure()
plt.semilogy(freqs_wide, np.abs(X_lin_wide), label="Linear-interp FT")
plt.semilogy(f_fft, np.abs(X_fft), linestyle="--", marker="", label="Rectangular FFT")
plt.xlim(0, 300)
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude (log scale)")
plt.title("Leakage / sidelobes comparison")
plt.grid(True)
plt.legend()
plt.show()

# === Figure 3: Boundary + interpolation error vs high-res numeric FT ===
# Build a slightly more complicated signal with trend to emphasize boundary
x_trend = 0.8 * np.sin(2 * np.pi * f0 * t) + 0.2 * (t / t[-1])

# High-res numeric FT as "truth"
oversample = 50
t_fine = np.arange(0, t[-1] + dt/oversample, dt/oversample)
x_fine = 0.8 * np.sin(2 * np.pi * f0 * t_fine) + 0.2 * (t_fine / t[-1])

freqs_err = np.linspace(0, 200, 300)
omega_err = 2 * np.pi * freqs_err

X_true = np.zeros_like(freqs_err, dtype=np.complex128)
for i, w in enumerate(omega_err):
    X_true[i] = np.trapz(x_fine * np.exp(-1j * w * t_fine), t_fine)

X_fft_trend = np.fft.rfft(x_trend) * dt
f_fft_trend = np.fft.rfftfreq(N, dt)
X_lin_trend = fourier_linear_interp(x_trend, dt, freqs_err, angular_freq=False)

# Interpolate FFT result onto freqs_err grid for fair comparison
X_fft_interp = np.interp(freqs_err, f_fft_trend, np.abs(X_fft_trend))

err_fft = np.abs(X_fft_interp - np.abs(X_true))
err_lin = np.abs(np.abs(X_lin_trend) - np.abs(X_true))

plt.figure()
plt.semilogy(freqs_err, err_fft + 1e-18, label="Rectangular FFT error")
plt.semilogy(freqs_err, err_lin + 1e-18, label="Linear-interp FT error")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Absolute error vs high-res numeric FT")
plt.title("Boundary + interpolation error comparison")
plt.grid(True)
plt.legend()
plt.show()

# === Figure 4: Window effect (Hann) on both methods ===
hann = 0.5 - 0.5 * np.cos(2 * np.pi * np.arange(N) / (N - 1))
x_win = x_base * hann

X_fft_win = np.fft.rfft(x_win) * dt
X_lin_win = fourier_linear_interp(x_win, dt, freqs_wide, angular_freq=False)

plt.figure()
plt.semilogy(freqs_wide, np.abs(X_lin_wide), label="Rectangular window - Linear-interp")
plt.semilogy(f_fft, np.abs(X_fft), linestyle="--", label="Rectangular window - FFT")
plt.semilogy(freqs_wide, np.abs(X_lin_win), label="Hann window - Linear-interp")
plt.semilogy(f_fft, np.abs(X_fft_win), linestyle="--", label="Hann window - FFT")
plt.xlim(0, 300)
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude (log scale)")
plt.title("Window effect on leakage (Rectangular vs Hann)")
plt.grid(True)
plt.legend()
plt.show()