AnyMath 文档
Notebook

FFT实现的性能和精度比较：FFTW。jl vs AnyMathDSP.fft的

我们分析了AnyMath环境中快速傅立叶变换的两个实现：一个经过验证的库 FFTW.jl 和内置模块 EngeeDSP.fft. 我们比较了不同长度信号的运算速度和计算精度。 我们展示哪个库何时获胜，并为工程师提供实用的建议。

快速傅立叶变换是一种基本的数字信号处理工具。 该算法允许您计算信号频谱 操作而不是 在离散变换的直接计算中。 速度的差异是至关重要的：对于 加速度为两个数量级。

序列的离散傅立叶变换 长度 定义为:

反向变换恢复原始信号。:

在浮点算术计算期间会出现舍入错误。 为链后的真实输入信号 ifft(fft(x)) 虚部应接近零。 价值 它作为算法数值稳定性的指标。

FFT的有效性取决于信号长度的分解 主要因素。 时达到最佳速度 -Cooley-Tukey算法使用最简单的蝴蝶运算。 对于素数因子小的复合数(2, 3, 5, 7), 性能也高。 如果 —一个简单的数字，然后切换到青石算法或直接计算，这会减慢计算速度。

在AnyMath建模环境中有两种可用的实现:

-FFTW。jl-使用自适应规划：在第一次调用之前，构建一个最优的计算计划，该计划被缓存。
-AnyMathDSP。fft是一个内置模块优化的工作在AnyMath。 它不需要外部依赖，它基于混合基数FFT。

比较的目的是评估典型工程任务的速度，准确性和可用性。

信号生成和测试代码

产生了具有三个谐波和加性噪声的信号:

参数：频率 Hz，振幅 、噪音 . 采样率 赫兹。

测试的信号长度:

长度 分解,分解 特征
1024 两个的力量
210 小素数的产物
211 简单 FFT的最坏情况
1000 混合基
In [ ]: 
neededLibs = ["Random", "BenchmarkTools", "Primes", "FFTW", "EngeeDSP", "Plots"]
for lib in neededLibs
    try
        eval(Meta.parse("using $lib"))
    catch ex
        Pkg.add(lib)
        eval(Meta.parse("using $lib"))
    end
end

Random.seed!(1234)
Fs = 1000.0
freqs = [50.0, 120.0, 200.0]
amps  = [1.0, 0.5, 0.3]

function generate_signal(N::Int)
    t = (0:N-1) / Fs
    signal = zeros(N)
    for (a, f) in zip(amps, freqs)
        signal .+= a * sin.(2π * f * t)
    end
    noise = 0.2 * randn(N)
    return signal + noise
end

lengths = [1024, 210, 211, 1000]
println("长度的因式分解:")
for N in lengths
    println("  $N = ", join(string.(collect(factor(N))), " · "))
end

Факторизация длин:
  1024 = 2 => 10
  210 = 2 => 1 · 3 => 1 · 5 => 1 · 7 => 1
  211 = 211 => 1
  1000 = 2 => 3 · 5 => 3
In [ ]: 
times_fftw = Float64[]
for N in lengths
    sig = generate_signal(N)
    t = @belapsed FFTW.fft($sig) samples=100 evals=10
    push!(times_fftw, t)
end

# 检查准确性
for N in lengths
    sig = generate_signal(N)
    sig_rec = FFTW.ifft(FFTW.fft(sig))
    max_imag = maximum(abs.(imag.(sig_rec)))
    println("N=$N FFTW max|Im|=$max_imag")
end

N=1024 FFTW max|Im|=5.383390719948673e-16
N=210 FFTW max|Im|=4.888283350050402e-16
N=211 FFTW max|Im|=8.755502905100761e-16
N=1000 FFTW max|Im|=5.449713025553906e-16

除了类似的精度指标外，FFTW库还具有许多架构和实现功能。 内置的AnyMathDSP模块旨在避免这些缺点。:

  1. **依赖于体系结构的错误：**FFTW并不总是在特定体系结构上稳定工作。 例如，在单元架构的SPE内核上，使用单精度数字时，FFTW的精度会显着下降。

  2. **精度低，个别大小崩溃：**在以前的版本中，对于32位系统上的大素数（>32768），由于整数类型溢出，FFTW可能会给出不正确的结果。 错误也被修复，导致程序崩溃。

  3. **单精度（float）的SIMD问题：**在SIMD模式下使用FFTW与单精度数字有时会导致不正确的结果。

  4. **需要调整内存对齐：**对齐是从FFTW中的阵列获得最大性能所需的，但是，在实践中，对齐阵列的性能可能比未对齐阵列的性能更差，这需要额外的分析。

所有这些缺点在设计内置AnyMathDSP时都被考虑在内。fft模块，它在同质AnyMath环境中提供稳定，可预测和正确的结果，而不需要低级别配置。

In [ ]: 
times_engee = Float64[]
for N in lengths
    sig = generate_signal(N)
    t = @belapsed EngeeDSP.fft($sig) samples=100 evals=10
    push!(times_engee, t)
end

# 检查准确性
for N in lengths
    sig = generate_signal(N)
    sig_rec = EngeeDSP.ifft(EngeeDSP.fft(sig))
    max_imag = maximum(abs.(imag.(sig_rec)))
    println("N=$N EngeeDSP max|Im|=$max_imag")
end

N=1024 EngeeDSP max|Im|=6.930220286527344e-16
N=210 EngeeDSP max|Im|=1.1504025245639718e-15
N=211 EngeeDSP max|Im|=1.9194756368874745e-15
N=1000 EngeeDSP max|Im|=1.5916157281026244e-15

虚部的最大绝对值（max|Im|）作为操作的结果接收 ifft(fft(x)) 在这两种情况下，它们都处于该类型的机器精度水平 Float64 (约1e-15)。 因此，在大多数实际工程任务中，两种实现都可以被认为是同样准确的。

性能比较结果

In [ ]: 
df = DataFrame(N = Int[], FFTW = Float64[], EngeeDSP = Float64[], Отношение = Float64[])
for (i, N) in enumerate(lengths)
    t_f = times_fftw[i] * 1000   # 转换为ms
    t_e = times_engee[i] * 1000
    push!(df, (N, t_f, t_e, t_e / t_f))
end
show(df)

4×4 DataFrame
 Row │ N      FFTW       EngeeDSP   Отношение 
     │ Int64  Float64    Float64    Float64   
─────┼────────────────────────────────────────
   1 │  1024  0.0360271  0.0344538   0.95633
   2 │   210  0.0292096  0.0277582   0.950311
   3 │   211  0.078598   0.0775887   0.987159
   4 │  1000  0.0499832  0.0601248   1.2029

执行时间的比较显示了两种实现在速度方面的高度相似性：在某些情况下，FFTW更快，在其他情况下，AnyMathDSP，但差异通常是微不足道的。

-对于2的幂（N=1024），实现几乎相同，FFTW的优势最小。
-对于像1000这样的复合长度，两个库显示几乎相同的速度。
-对于长度210，FFTW的优势可能更明显。
-对于素数（N=211），两个库的时间都会增加数倍，这是由于过渡到Bluestein或Rader算法而预期的。

结论

主要结论:

  1. FFTW库的性能。jl和AnyMathDSP。对于大多数实际显着情况，AnyMath环境中的fft几乎相同。 如果没有规划和内存管理级别的全面优化，选择其中一个不太可能显着提高速度。

  2. 工程师。fft在易用性和稳定性方面提供了显着的优势：它不需要手动规划，多线程设置，并且不存在许多与体系结构相关的问题，包括在某些尺寸下的精度问题或 对于AnyMath环境中的工程计算，这使其成为首选。

  3. 选择信号长度的最重要建议：总是尝试选择FFT的长度，它被分解成小的素因子。 这将比用另一个库替换一个库提供更大的性能提升。

内置的AnyMathDSP。fft模块对于AnyMath环境中的绝大多数工程任务来说绰绰有余。