Алгоритмы вычисления Быстрых Преобразований Фурье (БПФ) в многомерных пространствах

Алгоритм, изменивший цифровую эпоху

Преобразование Фурье — это математическая основа обработки любых сигналов. Оно позволяет раскладывать сложные волновые формы на набор простых синусоид. Дискретное преобразование Фурье (ДПФ) выполняет эту задачу для цифровых (сэмплированных) данных. Однако прямое, «школьное» вычисление ДПФ по формуле имеет квадратичную алгоритмическую сложность O(N^2). Для аудиофайла длиной в одну секунду (44 100 отсчетов) потребуется около двух миллиардов операций умножения. Если бы мы пользовались этим прямым методом, современная потоковая передача видео, мобильная связь 4G/5G и магнитно-резонансная томография (МРТ) были бы просто невозможны из-за вычислительных ограничений.

Проблема была решена в 1965 году, когда Джеймс Кули и Джон Тьюки опубликовали алгоритм Быстрого Преобразования Фурье (БПФ, или FFT). (Справедливости ради, похожие идеи высказывал еще Карл Фридрих Гаусс в 1805 году, но они опередили свое время). Алгоритм БПФ радикально снижает количество операций с O(N^2) до O(N log N). Для массива в миллион точек ускорение составляет почти 50 000 раз! Журнал Computing in Science & Engineering заслуженно включил БПФ в десятку величайших алгоритмов XX века.

Разделяй и властвуй: принцип бабочки

Самая популярная версия алгоритма Кули-Тьюки базируется на принципе «разделяй и властвуй» (divide and conquer) по основанию 2 (Radix-2). Алгоритм требует, чтобы количество точек N в массиве было степенью двойки (например, 256, 512, 1024). Если это не так, массив дополняется нулями (zero-padding).

Алгоритм рекурсивно разбивает задачу вычисления одного ДПФ размера N на вычисление двух ДПФ размера N/2 (одно для четных индексов массива, другое — для нечетных). Затем каждый из этих подмассивов снова делится пополам, и так далее, пока размер массива не достигнет 1. Затем результаты начинают объединяться (синтезироваться) обратно. На этапе объединения используется базовая вычислительная структура, которая из-за графа информационных потоков получила название «операция бабочки» (Butterfly operation). Она комбинирует два комплексных числа, используя умножение на специальные комплексные экспоненты (поворачивающие множители). Перед выполнением БПФ элементы массива переставляются в специальном порядке, называемом битовой инверсией (Bit-reversal permutation), что позволяет выполнять все вычисления прямо на месте (in-place), не выделяя дополнительную оперативную память.

Многомерное БПФ: от звука к изображениям

Обычное одномерное БПФ применяется для обработки звука или одномерных радиосигналов. Но как анализировать цифровые фотографии (поиск краев, сжатие JPEG, фильтрация шума) или решать трехмерные дифференциальные уравнения Пуассона на сетке? Для этого используется многомерное БПФ.

Счастье для инженеров заключается в том, что многомерное дискретное преобразование Фурье математически сепарабельно (разделимо). Это означает, что для вычисления двумерного БПФ изображения (матрицы пикселей) не нужно писать новый сложный алгоритм. Достаточно сначала применить обычное одномерное БПФ к каждой строке матрицы независимо. Получив промежуточную матрицу комплексных чисел, мы затем применяем то же самое одномерное БПФ к каждому ее столбцу. Этот строчный-столбцовый подход (Row-Column algorithm) легко масштабируется на 3D и 4D пространства, идеально распараллеливается на тысячи ядер графических процессоров (с помощью библиотек вроде cuFFT от NVIDIA) и является бьющимся сердцем современной компьютерной томографии и молекулярной химии.