使用fftpack.fft做离散傅里叶变换时多数频率极小非零的原因是什么

首先明确结论：你观测到的极小非零频谱值确实属于数值舍入误差，常规场景下不需要关注。

误差来源说明 #

• 信号生成阶段：np.sin是基于数值近似实现的三角函数计算，双精度浮点下本身就存在1e-16量级的计算误差，无法输出完全理想的正弦信号采样值。
• FFT计算阶段：FFT的蝶形运算包含大量浮点乘加操作，每一步运算都会引入极小的舍入误差，双精度场景下整体误差量级通常在1e-12~1e-15区间，因此不会出现理论上的零值（对应对数坐标下的负无穷），只会呈现一个极低的噪声底。
  补充：你的采样参数设置中，两个正弦波的频率都是频率分辨率delta_f=0.01Hz的整数倍，不存在频谱泄漏导致的旁瓣干扰，因此可以排除频谱泄漏的影响，确认这些小值就是舍入误差。

是否需要关注 #

• 常规频谱分析、频率/幅值提取场景下，这些误差的幅值比你的主信号（1Hz、22Hz分量）小10个数量级以上，完全不会对计算结果产生可观测的影响，不需要做额外处理。
• 仅当你的应用需要极高的动态范围（比如要检测幅值比主信号小10个数量级以上的超微弱信号）时，才需要针对性优化：比如改用更高精度的浮点运算类型、优化计算流程减少误差积累。

你使用的测试代码整理如下：

import numpy as np
from scipy import fftpack
from matplotlib import pyplot as plt

def signal_samples(t):
    A1 = 4  #amplitude
    f1 = 1  #frequency [Hz]
    p1 = 0  #phase [rad]
    
    A2 = 2  #amplitude
    f2 = 22 #frequency [Hz]
    p2 = 0  #phase [rad]
    
    # sum of two sine waves        
    signal = (A1 * np.sin(2*np.pi*f1*t + p1) + A2 * np.sin(2*np.pi*f2*t + p2)) 
    
    return(signal)

B = 30                      #max frequency
f_s = 2*B                   #sampling frequency twice as high
delta_f = 0.01              #frequency resolution in spectrum
N = int(f_s / delta_f)      #number of samples
T = N / f_s 

t = np.linspace(0,T,N,endpoint=False)
s = signal_samples(t)

fig,axs = plt.subplots(1,2,figsize=(8,3),sharey=True)
axs[0].plot(t,s)
axs[0].set_xlabel("time (s)")
axs[0].set_ylabel("signal")
axs[1].plot(t,s)
axs[1].set_xlim(0,5)
axs[1].set_xlabel("time (s)")
plt.show()

# Perform Discrete Fourier Transform
F = fftpack.fft(s)               # FFT
f = fftpack.fftfreq(N,1.0/f_s)   # Get frequency bins to match FFT  

mask = np.where(f>=0) # spectrum is symmetric around 0; only take positive frequencies

fig, axs = plt.subplots(3,1,figsize=(8,6))
axs[0].plot(f[mask],np.log(abs(F[mask])), label="real")
axs[0].plot(B,0,'r*', markersize=10)
axs[0].set_ylabel("log |F|", fontsize=14)

axs[1].plot(f[mask],np.abs(F[mask])/N, label="real")
axs[1].set_ylabel("|F| / N", fontsize=14)
axs[2].plot(f[mask],np.abs(F[mask])/N, label="real")

axs[2].set_xlim(21,23)
axs[2].set_xlabel("Frequency (Hz)", fontsize=14)
axs[2].set_ylabel("|F| / N", fontsize=14)
plt.show()

内容的提问来源于stack exchange，提问作者Antoine