Hann 窗：主瓣变宽，但旁瓣瞬间掉到 以下，然后非常快地衰减→ 典型的“主瓣变宽、副瓣下降”的窗函数 trade-off，在两种算法中都一致。
同一窗形里，新算法 vs FFT：
对于矩形窗：两条矩形窗曲线整体很接近，但线性插值的细节更平滑；对于 Hann 窗：两条 Hann 曲线几乎“扣在一起”，说明当窗函数已经把时间域信号变得相对光滑时，线性插值带来的改进主要是数值精度，不再根本改变泄漏形态。
窗的作用是“改变核函数本身”（决定泄漏结构）
线性插值 FT 是在既定核函数下做“更精确积分”两者是正交的；可以 “先用窗 + 再用改进积分算法”，效果叠加。
数值算法实现
我先放使用论文核心实现的函数(其实数值算法才是我的真爱）：
import numpy as np
deffourier_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)
# ---- ω = 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)
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)
# ---- ω = 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)
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.")
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)
