Быстрое преобразование Фурье (БПФ)

Узнаём, как умный алгоритм превращает медленное ДПФ в молниеносный БПФ — без потери точности.

БПФ (быстрое преобразование Фурье, FFT) — алгоритм вычисления того же ДПФ за O(N log N) операций вместо O(N²), основанный на стратегии «разделяй и властвуй».

Прямое ДПФ корректно, но медленно: для миллиона отсчётов нужен триллион умножений. БПФ, переоткрытый Кули и Тьюки в 1965 году (идею знал ещё Гаусс), считает тот же результат за считанные миллионы операций. Без него не было бы ни цифрового звука в реальном времени, ни Wi-Fi, ни МРТ. Разберём идею и реализуем рекурсивный БПФ на чистом Python.

Идея: разделяй и властвуй

Ключевое наблюдение: ДПФ от N точек можно собрать из двух ДПФ по N/2 точек — отдельно по чётным и отдельно по нечётным индексам. Эти половинки комбинируются с помощью «поворачивающих множителей» exp(-2j*pi*k/N). А каждую половинку, в свою очередь, снова делят пополам — и так до отдельных точек. Деление пополам log2(N) раз и даёт логарифмический множитель.

import math, cmath

def fft(x):
    N = len(x)
    if N == 1:
        return [complex(x[0])]
    even = fft(x[0::2])      # ДПФ чётных индексов
    odd = fft(x[1::2])       # ДПФ нечётных индексов
    T = [cmath.exp(-2j * math.pi * k / N) * odd[k] for k in range(N // 2)]
    return ([even[k] + T[k] for k in range(N // 2)] +
            [even[k] - T[k] for k in range(N // 2)])

x = [0, 1, 2, 3, 4, 5, 6, 7]
X = fft(x)
for k, c in enumerate(X):
    re = round(c.real, 3) or 0.0
    im = round(c.imag, 3) or 0.0
    print(f"X[{k}] = {re:+.3f} {im:+.3f}j")

Вывод:

X[0] = +28.000 +0.000j
X[1] = -4.000 +9.657j
X[2] = -4.000 +4.000j
X[3] = -4.000 +1.657j
X[4] = -4.000 +0.000j
X[5] = -4.000 -1.657j
X[6] = -4.000 -4.000j
X[7] = -4.000 -9.657j

X[0] = 28 — сумма чисел от 0 до 7. Это настоящий БПФ, посчитавший спектр за O(N log N) операций.

БПФ даёт тот же результат, что и ДПФ

Важно убедиться: БПФ — это не приближение, а точно то же ДПФ, только быстрее. Сравним их выходы.

import math, cmath

def fft(x):
    N = len(x)
    if N == 1:
        return [complex(x[0])]
    even = fft(x[0::2]); odd = fft(x[1::2])
    T = [cmath.exp(-2j * math.pi * k / N) * odd[k] for k in range(N // 2)]
    return ([even[k] + T[k] for k in range(N // 2)] +
            [even[k] - T[k] for k in range(N // 2)])

def dft(x):
    N = len(x)
    return [sum(x[n] * cmath.exp(-2j * math.pi * k * n / N) for n in range(N))
            for k in range(N)]

x = [1, 0, 2, 0, 1, 0, 2, 0]
max_err = max(abs(a - b) for a, b in zip(fft(x), dft(x)))
print(f"Максимальное расхождение БПФ и ДПФ: {max_err:.1e}")

Вывод:

Максимальное расхождение БПФ и ДПФ: 4.8e-15

Расхождение — на уровне машинной точности (десятые доли от 10**-15), то есть это буквально один и тот же результат. БПФ ничего не теряет.

Насколько быстрее

import math

def dft_ops(N):
    return N * N
def fft_ops(N):
    return int(N * math.log2(N))

print(f"{'N':>9} | {'ДПФ ~':>15} | {'БПФ ~':>12} | ускорение")
for N in [8, 64, 1024, 1048576]:
    sp = dft_ops(N) // max(fft_ops(N), 1)
    print(f"{N:>9} | {dft_ops(N):>15} | {fft_ops(N):>12} | {sp}x")

Вывод:

        N |           ДПФ ~ |        БПФ ~ | ускорение
        8 |              64 |           24 | 2x
       64 |            4096 |          384 | 10x
     1024 |         1048576 |        10240 | 102x
  1048576 |   1099511627776 |     20971520 | 52428x

Для миллиона отсчётов БПФ быстрее в 50 тысяч раз. Именно поэтому спектр аудио считают в реальном времени.

Как работает под капотом

Магия в поворачивающих множителях. Заметим, что exp(-2j*pi*(k+N/2)/N) = -exp(-2j*pi*k/N) — половинный сдвиг просто меняет знак. Это позволяет посчитать множитель один раз и использовать его и со знаком «+», и со знаком «-» (схема «бабочка»). Деление на чётные/нечётные индексы плюс «бабочки» — и из N² операций остаётся N log N. Классический алгоритм Кули-Тьюки требует, чтобы N было степенью двойки (наша реализация — тоже). Если длина не степень двойки, сигнал дополняют нулями (zero-padding) до ближайшей степени или применяют более общие варианты БПФ. В реальных библиотеках (FFTW, numpy) используют сильно оптимизированные итеративные версии без рекурсии.

Частые ошибки

Подавать длину не степень двойки в этот БПФ. Рекурсия делит пополам; для произвольной длины нужен zero-padding или общий алгоритм.
Думать, что БПФ точнее ДПФ. Он не точнее и не «другой» — это тот же результат, посчитанный экономнее.
Игнорировать накладные расходы на малых N. На очень коротких сигналах прямое ДПФ может быть не медленнее из-за простоты.

Итог

БПФ считает то же ДПФ за O(N log N) вместо O(N²).
Идея — «разделяй и властвуй»: ДПФ из двух половинных ДПФ через поворачивающие множители.
Классический БПФ требует длину — степень двойки; иначе дополняют нулями.
Ускорение колоссально: для миллиона точек — в десятки тысяч раз.