Engee documentation
Notebook

Comparison of performance and accuracy of FFT implementations: FFTW.jl vs EngeeDSP.fft

We analyze two implementations of the fast Fourier transform in the Engee environment: a proven library FFTW.jl and the built-in module EngeeDSP.fft. We compare the speed of operation and the accuracy of calculations for signals of different lengths. We show when which library wins, and give practical recommendations for engineers.

Fast Fourier transform is a basic digital signal processing tool. The algorithm allows you to calculate the signal spectrum for operations instead of in the direct calculation of the discrete transformation. The difference in speed is critical: for the acceleration is two orders of magnitude.

Discrete Fourier transform for a sequence lengths defined as:

The reverse transformation restores the original signal.:

Rounding errors occur during calculations in floating-point arithmetic. For the real input signal after the chain ifft(fft(x)) the imaginary part should be close to zero. Value It serves as an indicator of the numerical stability of the algorithm.

The effectiveness of the FFT depends on the decomposition of the signal length by prime factors. The best speed is achieved when — the Cooley–Tukey algorithm uses the simplest butterfly operations. For composite numbers with small prime factors (2, 3, 5, 7), the performance is also high. If — a simple number, libraries switch to the Bluestone algorithm or direct calculation, which slows down the calculation.

There are two implementations available in the Engee modeling environment:

  • FFTW.jl — Uses adaptive planning: before the first call, an optimal calculation plan is built, which is cached.
  • EngeeDSP.fft is a built—in module optimized for working inside Engee. It does not require external dependencies, it is based on a mixed-radix FFT.

The purpose of the comparison is to evaluate the speed, accuracy, and usability of typical engineering tasks.

Code for signal generation and testing

A signal with three harmonics and additive noise was generated:

Parameters: frequencies Hz, amplitudes , noise . Sampling rate Hz.

Signal lengths for tests:

Length Decomposition Characteristic
1024 The power of two
210 the product of small primes
211 simple the worst case scenario for the FFT
1000 mixed base
In [ ]: 
using Random, BenchmarkTools, Primes, FFTW, EngeeDSP, Plots

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("Factorization of lengths:")
for N in lengths
    println("  $N = ", join(string.(collect(factor(N))), " · "))
end

Factorization of lengths:
  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

# Checking accuracy
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

In addition to similar accuracy indicators, the FFTW library has a number of architectural and implementation features. The built-in EngeeDSP module is designed to avoid these disadvantages.:

  1. Architecture-dependent errors: FFTW does not always work stably on specific architectures. For example, on SPE cores of the Cell architecture, the accuracy of FFTW can decrease significantly when working with single-precision numbers.

  2. Low accuracy and crashes for individual sizes: In previous versions, for large primes (>32768) on 32-bit systems, FFTW could give incorrect results due to integer type overflow. Bugs were also fixed that lead to the crash of programs.

  3. SIMD problems for single precision (float): Using FFTW in SIMD mode with single precision numbers sometimes leads to incorrect results.

  4. The need to adjust the memory alignment: Alignment is required for maximum performance from arrays in FFTW, however, in practice, performance with aligned arrays may be worse than with unaligned arrays, which requires additional profiling.

All these disadvantages are taken into account when designing the built-in EngeeDSP.fft module, which provides stable, predictable and correct results in a homogeneous Engee environment without the need for low-level configuration.

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

# Checking accuracy
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

Maximum absolute values of the imaginary part (max|Im|) received as a result of the operation ifft(fft(x)) in both cases , they are at the level of machine precision for the type Float64 (about 1e-15). Therefore, both implementations can be considered equally accurate in most practical engineering tasks.

Performance comparison results

In [ ]: 
using DataFrames
df = DataFrame(N = Int[], FFTW = Float64[], EngeeDSP = Float64[], Relation = Float64[])
for (i, N) in enumerate(lengths)
    t_f = times_fftw[i] * 1000   # conversion to ms
    t_e = times_engee[i] * 1000
    push!(df, (N, t_f, t_e, t_e / t_f))
end
show(df)

4×4 DataFrame
 [1m Row │ [1m N [1m FFTW [1m EngeeDSP [1m Ratio
     │ Int64  Float64    Float64    Float64
─────┼────────────────────────────────────────
   1 │  1024  0.022323   0.0231451   1.03683
   2 │   210  0.0202879  0.0250997   1.23718
   3 │   211  0.0644022  0.0657494   1.02092
   4 │  1000  0.0312717  0.0311902   0.997394

A comparison of the execution time shows the high similarity of both implementations in terms of speed: in some cases, FFTW is faster, in others, EngeeDSP, but the difference is generally insignificant.

  • For the power of two (N=1024), the implementations are almost identical, with a minimal advantage of FFTW.
  • For composite lengths like 1000, both libraries show almost the same speed.
  • For length 210, the advantage of FFTW may be more noticeable.
  • For a prime number (N=211), the time increases several times for both libraries, which is expected due to the transition to the Bluestein or Rader algorithm.

Conclusion

The main conclusions:

  1. Performance of FFTW libraries.jl and EngeeDSP.The fft in the Engee environment is almost identical for most practically significant cases. Choosing one of them is unlikely to give a significant increase in speed without comprehensive optimization at the planning and memory management levels.

  2. EngeeDSP.fft provides significant advantages in terms of ease of use and stability: it does not require manual planning, multithreading settings, and does not have many architecture-dependent problems, including problems with precision at certain sizes or with SIMD optimizations. For engineering calculations in the Engee environment, this makes it the preferred choice.

  3. The most important advice on choosing the signal length: Always try to choose the length of the FFT, which is decomposed into small prime factors. This will give a much greater performance boost than replacing one library with another.

The built-in EngeeDSP.fft module is more than enough for the vast majority of engineering tasks in the Engee environment.