Документация Engee
Notebook

Сравнение производительности и точности реализаций БПФ: FFTW.jl против EngeeDSP.fft

Разбираем две реализации быстрого преобразования Фурье в среде Engee: проверенную библиотеку FFTW.jl и встроенный модуль EngeeDSP.fft. Сравниваем скорость работы и точность расчётов для сигналов разной длины. Показываем, когда какая библиотека выигрывает, и даём практические рекомендации для инженеров.

Быстрое преобразование Фурье — базовый инструмент цифровой обработки сигналов. Алгоритм позволяет вычислить спектр сигнала за операций вместо при прямом расчёте дискретного преобразования. Разница в скорости критична: для ускорение составляет два порядка.

Дискретное преобразование Фурье для последовательности длины определяется как:

Обратное преобразование восстанавливает исходный сигнал:

При вычислениях в арифметике с плавающей точкой возникают ошибки округления. Для вещественного входного сигнала после цепочки ifft(fft(x)) мнимая часть должна быть близка к нулю. Величина служит индикатором численной устойчивости алгоритма.

Эффективность БПФ зависит от разложения длины сигнала на простые множители. Наилучшая скорость достигается при — алгоритм Кули–Тьюки использует простейшие операции «бабочка». Для составных чисел с малыми простыми множителями (2, 3, 5, 7) производительность также высока. Если — простое число, библиотеки переключаются на алгоритм Блустейна или прямое вычисление, что замедляет расчёт.

В среде моделирования Engee доступны две реализации:

  • FFTW.jl — Использует адаптивное планирование: перед первым вызовом строится оптимальный план вычислений, который кэшируется.
  • EngeeDSP.fft — встроенный модуль, оптимизированный для работы внутри Engee. Не требует внешних зависимостей, основан на смешанно-радиксном БПФ.

Цель сравнения: оценить скорость, точность и удобство использования в типовых инженерных задачах.

Код для генерации сигнала и тестирования

Генерировали сигнал с тремя гармониками и аддитивным шумом:

Параметры: частоты Гц, амплитуды , шум . Частота дискретизации Гц.

Длины сигналов для тестов:

Длина Разложение Характеристика
1024 степень двойки
210 произведение малых простых
211 простое наихудший случай для БПФ
1000 смешанное основание
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("Факторизация длин:")
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 имеет ряд архитектурных и реализационных особенностей. Встроенный в среду модуль EngeeDSP спроектирован таким образом, чтобы избежать этих недостатков:

  1. Архитектурно-зависимые ошибки: FFTW не всегда стабильно работает на специфичных архитектурах. Например, на ядрах SPE архитектуры Cell точность FFTW может снижаться в разы при работе с числами одинарной точности.

  2. Низкая точность и крэши для отдельных размеров: В прошлых версиях для больших простых чисел (>32768) в 32-битных системах FFTW могла давать неверные результаты из-за переполнения целочисленных типов. Также были зафиксированы баги, приводящие к краху программ.

  3. Проблемы SIMD для одинарной точности (float): Использование FFTW в режиме SIMD с числами одинарной точности иногда приводит к некорректным результатам.

  4. Необходимость настройки выравнивания памяти: Для максимальной производительности от массивов в FFTW требуется выравнивание, однако на практике производительность с выровненными массивами может быть хуже, чем с невыровненными, что требует дополнительного профилирования.

Все эти недостатки учтены при проектировании встроенного модуля EngeeDSP.fft, который обеспечивает стабильные, предсказуемые и корректные результаты в гомогенной среде Engee без необходимости низкоуровневой настройки.

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 [ ]: 
using DataFrames
df = DataFrame(N = Int[], FFTW = Float64[], EngeeDSP = Float64[], Отношение = Float64[])
for (i, N) in enumerate(lengths)
    t_f = times_fftw[i] * 1000   # перевод в мс
    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.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

Сравнение времени выполнения показывает высокую схожесть обеих реализаций по скорости: в одних случаях быстрее FFTW, в других — EngeeDSP, но разница в целом незначительна.

  • Для степени двойки (N=1024) реализации практически идентичны, с минимальным преимуществом FFTW.
  • Для составных длин, таких как 1000, обе библиотеки показывают почти одинаковую скорость.
  • Для длины 210 преимущество FFTW может быть заметнее.
  • Для простого числа (N=211) время возрастает в несколько раз у обеих библиотек, что ожидаемо из-за перехода на алгоритм Блустейна или Радера.

Вывод

Основные выводы:

  1. Производительность библиотек FFTW.jl и EngeeDSP.fft в среде Engee практически идентична для большинства практически значимых случаев. Выбор в пользу одной из них вряд ли даст ощутимый прирост скорости без комплексной оптимизации на уровне планирования и управления памятью.

  2. EngeeDSP.fft предоставляет существенные преимущества с точки зрения простоты использования и стабильности: он не требует ручного планирования, настройки многопоточности и не имеет многих архитектурно-зависимых проблем, включая проблемы с точностью на определённых размерах или с SIMD-оптимизациями. Для инженерных расчётов в среде Engee это делает его предпочтительным выбором.

  3. Важнейший совет по выбору длины сигнала: Всегда старайтесь выбирать длину БПФ, которая разлагается на малые простые множители. Это даст гораздо больший прирост производительности, чем замена одной библиотеки на другую.

Для подавляющего большинства инженерных задач в среде Engee встроенного модуля EngeeDSP.fft более чем достаточно.