Python实现DFT时浮点三角函数精度问题的解决方案咨询

一、降低DFT浮点误差的简便优化方案 #

1. 用复数指数运算替代拆分的三角函数计算 #

原生math.cos和math.sin分开计算会引入额外舍入误差，改用cmath.exp直接计算复数指数（exp(-1j*phi)等价于cos(phi) - i*sin(phi)），Python内部对复数运算有更优的精度优化，代码改动极小：

import math
import cmath

def optimized_DFT(values):
    N = len(values)
    out = []
    for n in range(N):
        complex_sum = 0j
        for i, value in enumerate(values):
            phi = 2 * math.pi * n * i / N
            complex_sum += value * cmath.exp(-1j * phi)
        complex_sum /= N
        real = complex_sum.real
        imag = complex_sum.imag
        freq = n
        amplitude = abs(complex_sum)
        phase = cmath.phase(complex_sum)
        out.append((freq, amplitude, phase, real, imag))
    return out

2. 利用实输入DFT的对称性减少计算 #

如果输入是实数值，DFT结果满足共轭对称性：X[n]和X[N-n]互为共轭复数。利用这一点只需计算前半部分频率点，剩下的直接通过共轭推导，既减少计算量，也降低了一半的误差累积：

def symmetric_DFT(values):
    N = len(values)
    out = []
    # 计算前半部分频率（包含0和N/2，当N为偶数时）
    for n in range(N//2 + 1):
        complex_sum = 0j
        for i, value in enumerate(values):
            phi = 2 * math.pi * n * i / N
            complex_sum += value * cmath.exp(-1j * phi)
        complex_sum /= N
        real = complex_sum.real
        imag = complex_sum.imag
        freq = n
        amplitude = abs(complex_sum)
        phase = cmath.phase(complex_sum)
        out.append((freq, amplitude, phase, real, imag))
        # 推导共轭点（跳过0和N/2，避免重复）
        if n != 0 and n != N - n:
            conj_real = real
            conj_imag = -imag
            conj_amplitude = amplitude
            conj_phase = -phase
            out.append((N - n, conj_amplitude, conj_phase, conj_real, conj_imag))
    return out

3. 用Kahan求和优化累加精度 #

处理大规模数据时，循环累加会累积浮点舍入误差，Kahan求和算法通过引入补偿项修正这类误差，适合高精度要求的场景：

import math
import cmath

def kahan_DFT(values):
    N = len(values)
    out = []
    for n in range(N):
        complex_sum = 0j
        # Kahan求和的补偿项，分别跟踪实部和虚部的误差
        c_real = 0.0
        c_imag = 0.0
        for i, value in enumerate(values):
            phi = 2 * math.pi * n * i / N
            term = value * cmath.exp(-1j * phi)
            # 实部的Kahan求和
            y_real = term.real - c_real
            t_real = complex_sum.real + y_real
            c_real = (t_real - complex_sum.real) - y_real
            complex_sum = complex(t_real, complex_sum.imag)
            # 虚部的Kahan求和
            y_imag = term.imag - c_imag
            t_imag = complex_sum.imag + y_imag
            c_imag = (t_imag - complex_sum.imag) - y_imag
            complex_sum = complex(complex_sum.real, t_imag)
        complex_sum /= N
        real = complex_sum.real
        imag = complex_sum.imag
        freq = n
        amplitude = abs(complex_sum)
        phase = cmath.phase(complex_sum)
        out.append((freq, amplitude, phase, real, imag))
    return out

二、最优选择 #

如果要平衡精度、代码复杂度和性能，优先级如下：

1. 优先使用cmath.exp替代拆分三角函数的版本，改动最小，精度提升明显，性能也不会下降。
2. 若输入是实数值，叠加对称性优化，进一步降低误差和计算量。
3. 只有处理超大规模数据且对精度要求极高时，才考虑Kahan求和；普通场景下双精度浮点的精度已经足够。

三、decimal模块的适用性判断 #

不推荐用decimal模块处理DFT，原因如下：

• 性能瓶颈：decimal是软件实现的高精度类型，运算速度比原生双精度慢几个数量级，O(N²)的DFT会变得极慢。
• 精度冗余：双精度浮点已经有15-17位有效数字，完全满足绝大多数信号处理场景的精度需求，decimal的超高精度属于过度设计。
• 实现复杂度：需要手动适配所有复数运算、累加逻辑，代码复杂度陡增，反而容易引入新的人为误差。

只有在极少数需要几十位有效数字的特殊科学计算场景，且能接受性能损失时，才考虑decimal模块。

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