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

Синтезатор играет арпеджио, а картинка показывает, какие частоты звучат в каждый момент: светлая полоса — громкая частота. Чтобы её нарисовать, браузер режет шесть секунд звука на 513 окон по 1024 отсчёта и каждое раскладывает на частоты. По определению это 1024 × 1024 умножений на окно, больше полумиллиарда на картинку. Быстрое преобразование Фурье делает то же самое за 5120 умножений на окно — в двести раз меньше, — и картинка пересчитывается быстрее, чем вы отпустите кнопку. Меняйте форму волны, фильтр, длину окна.

У синуса на картинке одна линия — одна частота. У «пилы» — лестница: основной тон и все кратные ему обертоны, у меандра — только нечётные, и лестница реже. FM-колокол звучит обертонами, которые не кратны основному тону, поэтому лестница у него кривая. Фильтр, который плавно открывается и закрывается, виден как волна яркости по верху. А длина окна — это выбор между временем и частотой: окно в 256 отсчётов видит каждую ноту отдельно, но размазывает частоты, окно в 4096 разрешает частоты до 5 Гц и размазывает ноты друг в друга.

Идея

Дискретное преобразование Фурье (ДПФ) N отсчётов x₀, …, x_N−1 — это N чисел

Посмотрите на эту сумму как на многочлен: P(z) = x₀ + x₁z + x₂z² + … + x_N−1z^N−1. Тогда X_k = P(ω^k) — значение многочлена в k-м корне N-й степени из единицы. ДПФ — это вычисление многочлена в N точках, и в лоб на это уходит N² умножений.

Разделим многочлен на чётные и нечётные степени: P(z) = E(z²) + z · O(z²), где у E и O вдвое меньше коэффициентов. Квадраты N корней из единицы — это всего N/2 разных чисел, корни степени N/2. Значит, E и O достаточно вычислить в N/2 точках — это два ДПФ половинной длины, — а затем собрать ответ: X_k = E_k + ω^kO_k и X_k+N/2 = E_k − ω^kO_k. Одна задача размера N превратилась в две размера N/2 и N/2 умножений.

Повторяя деление, пока длина не станет единицей, получаем log₂ N уровней по N/2 умножений: N/2 · log₂ N. Знак минус во второй формуле — потому что ω^N/2 = −1: точки z и −z дают одинаковые квадраты, и одно вычисление E и O обслуживает обе. Вывод с доказательством обратимости и с историей Гаусса, который придумал этот приём около 1805 года, — в главе «Ряды Фурье» «Царицы наук». Здесь посмотрим, как алгоритм работает на деле и что с ним делают.

Бабочки

Рекурсию удобно развернуть. Если на каждом уровне чётные отсчёты уходят налево, а нечётные направо, то в самом низу отсчёты окажутся в порядке, где номер записан двоичными цифрами задом наперёд: 0, 4, 2, 6, 1, 5, 3, 7. Переставим их так сразу — и дальше алгоритм идёт снизу вверх, уровень за уровнем, сливая соседние блоки операцией, которую за форму схемы называют бабочкой: из пары a, b получаются a + w·b и a − w·b.

Каждый кружок — комплексное число, стрелка показывает его направление. На входе четыре единицы и четыре нуля. На выходе X₀ = 4 — сумма всех отсчётов, а X₂, X₄, X₆ равны нулю: сигнал из блока единиц не содержит этих частот. Двенадцать бабочек — двенадцать умножений вместо шестидесяти четырёх, и всё вычисление идёт на месте, в одном массиве.

Код

Итеративный вариант: перестановка в обратном порядке битов, затем уровни бабочек длиной 2, 4, 8, … Обратное преобразование — то же самое с ω, заменённым на ω⁻¹, и делением на N. Программа раскладывает сигнал из четырёх синусоид и проверяет себя: на 128 отсчётах сравнивает ответ с ДПФ по определению, а обратным преобразованием возвращает исходный сигнал.

n = 1024: 5120 multiplications in the FFT, 1048576 by the definition
bin   55 =    440 Hz, amplitude 1.000
bin  110 =    880 Hz, amplitude 0.500
bin  165 =   1320 Hz, amplitude 0.250
bin  250 =   2000 Hz, amplitude 0.300

Множители ω^k программа получает, раз за разом умножая на одно и то же комплексное число, — так короче, но ошибки округления копятся. В длинных преобразованиях множители считают заново синусом и косинусом или берут из таблицы, как в следующих двух листингах.

Теперь настоящие данные: средняя суточная температура Алматы за 85 лет, 31 046 дней по реанализу ERA5. Какие периоды в ней спрятаны? Дополним ряд нулями до 131 072 отсчётов — больше чем вчетверо, чтобы частотная сетка была частой, — и посмотрим на самые высокие пики спектра.

31046 days, mean 8.76 C, FFT of 131072 samples: 1114112 butterflies instead of 17179869184 products
peak at bin 359: period 365.2 days, amplitude 13.2 C
peak at bin 4: period 33500.9 days, amplitude 1.0 C
peak at bin 718: period 182.6 days, amplitude 0.8 C

Главный пик — 365,2 дня: БПФ само нашло длину года, ничего о календаре не зная. Амплитуда 13,2 °C — лето теплее среднего примерно на столько же, на сколько зима холоднее. Второй пик — самая длинная волна, которая помещается в 85 лет: так в спектре выглядят медленные изменения средней температуры за десятилетия. Третий — полгода, 182,6 дня, с амплитудой 0,8 °C: годовой ход не чистая синусоида, и вторая гармоника описывает его неровность. Счёт в лоб занял бы 17 миллиардов умножений, БПФ обходится миллионом бабочек.

Умножение многочленов и больших чисел

Перемножим (1 + 2x + 3x²)(4 + 5x + 6x²). Коэффициент при x^k произведения — сумма a_ib_j по всем i + j = k:

Такая операция над последовательностями коэффициентов называется свёрткой, и в лоб она стоит n² умножений. Подставьте x = 10 — получится 4 + 130 + 2800 + 27 000 + 180 000 = 209 934 = 321 · 654: число — это многочлен от десяти, цифры — его коэффициенты, младшая первой. Умножение чисел столбиком и есть свёртка цифр с переносами.

БПФ делает свёртку быстрой. Многочлен степени n однозначно задаётся своими значениями в n + 1 точках. Вычислим оба множителя в корнях из единицы — прямым БПФ, — перемножим значения попарно — это n умножений, — и восстановим коэффициенты произведения обратным БПФ. Три преобразования по N log N плюс N умножений вместо n².

9543 digits times 10142 digits: FFT of length 8192
the product has 19684 digits: 399553154987076…763468471600001
rounding error below 0.01: yes
equals the built-in multiplication: yes
multiplications: column method 10754961, FFT about 167936

Цифры сгруппированы по три: каждый коэффициент — число от 0 до 999. Больше брать опасно: коэффициенты произведения доходят до n · 999², и вычисления в числах с плавающей точкой должны отличать их от соседних целых. Программа проверяет, насколько найденные коэффициенты далеки от целых, — здесь меньше чем на 0,01. Для чисел ещё длиннее ошибка растёт, и тогда переходят к арифметике остатков: теоретико-числовое преобразование делает то же самое без округлений.

Ниже замер: столбик, наше БПФ и встроенное умножение — int в Python и math/big в Go.

Столбик растёт как квадрат, БПФ — чуть быстрее линейного, и наше БПФ обгоняет столбик на Python уже к тысяче цифр, а на Go — между тысячей и тремя тысячами. А встроенное умножение оказывается быстрее обоих: в CPython и в Go это алгоритм Карацубы, асимптотически он хуже БПФ — d^1,58, — но написан на C и ассемблере и выигрывает за счёт константы на сотнях тысяч цифр. Библиотеки переходят к преобразованию Фурье там, где числа ещё длиннее: модуль decimal в Python использует теоретико-числовое преобразование для коэффициентов длиннее 1024 машинных слов, а проект GIMPS возводит в квадрат числа из десятков миллионов цифр специальным вариантом БПФ.

Сколько это стоит

Время работы удовлетворяет T(N) = 2T(N/2) + O(N): две задачи вдвое меньше и линейная склейка. По основной теореме это O(N log N), а точнее N/2 · log₂ N умножений и N log₂ N сложений. Память — O(N): итеративный вариант работает на месте.

Умножение многочленов степени n — три преобразования длины N ≥ 2n + 1: O(n log n). Для больших чисел из d цифр добавляются переносы, тоже линейные, и при группах фиксированного размера выходит O(d log d) — пока хватает точности чисел с плавающей точкой.

Где на БПФ спотыкаются

• Длина — степень двойки. Простейший вариант требует N = 2^m. Сигнал дополняют нулями — как с погодой выше, — а для произвольной длины есть разложения по другим множителям и алгоритм Блюстейна; библиотеки вроде numpy.fft принимают любую длину.
• Свёртка круговая. Произведение значений в корнях из единицы даёт свёртку по кругу: хвост произведения заворачивается в начало. Чтобы получить обычную, оба множителя дополняют нулями до длины не меньше суммы длин минус один.
• Утечка спектра. Окно из целого числа периодов даёт острые пики, как в первом листинге. В жизни окно режет волну где попало, и обрыв на краях разбрызгивает энергию по соседним частотам. Переключите в первой сцене окно на прямоугольное — линии расплывутся. Окно Ханна плавно гасит края и почти убирает брызги ценой чуть более широких пиков.
• Время против частоты. Окно длины N различает частоты с шагом (частота дискретизации)/N: чем точнее частота, тем длиннее окно и тем хуже видно, когда именно она звучала. Это не недостаток алгоритма, а свойство самих волн.
• Округление. В числах с плавающей точкой ошибка растёт с длиной и с величиной коэффициентов; для точных свёрток целых чисел — теоретико-числовое преобразование.

Где это работает

• Связь. Wi-Fi стандартов 802.11a/g использует в канале 20 МГц 52 поднесущие — 48 для данных и 4 опорные: передатчик собирает символ обратным 64-точечным БПФ, приёмник разбирает прямым.
• Звук и изображения. Спектроанализаторы, тюнеры, эквалайзеры, свёрточная реверберация, шумоподавление. Сжатие JPEG и MP3 держится на родственных косинусных преобразованиях — о них в главе «Ряды Фурье».
• Научные вычисления. numpy.fft работает на библиотеке pocketfft; решатели уравнений в частных производных, кристаллография, обработка данных телескопов.
• Длинная арифметика. Поиск простых чисел Мерсенна: тест Люка — Лемера возводит в квадрат числа из десятков миллионов цифр, и программы GIMPS делают это вариантом БПФ, который предложили Ричард Крэндалл и Барри Фейгин в 1994 году. Модуль decimal в Python умножает длинные коэффициенты теоретико-числовым преобразованием.
• Олимпиады. Свёртки счётчиков: сколькими способами набрать каждую сумму, совпадения строк при всех сдвигах — обе задачи ниже.

История

Гаусс придумал этот приём около 1805 года, вычисляя орбиты астероидов, но при жизни его не опубликовал — подробнее в «Царице наук». Второе рождение алгоритма связано с сейсмологией. По воспоминаниям участников, в 1963 году физик Ричард Гарвин из IBM слушал на совещании научных советников доклады об обнаружении подземных ядерных испытаний и заметил, что у статистика Джона Тьюки есть способ считать преобразование Фурье гораздо быстрее. Анализ сейсмических сигналов упирался именно в это преобразование. Гарвин нашёл в IBM Джеймса Кули, который написал программу; в 1964 году она заработала на машине IBM 7094, а в апреле 1965 года вышла пятистраничная статья Кули и Тьюки «An Algorithm for the Machine Calculation of Complex Fourier Series». В 1971 году Арнольд Шёнхаге и Фолькер Штрассен построили на преобразовании Фурье умножение целых чисел за O(n log n log log n); о том, как эту оценку довели до n log n, рассказывает глава «Разделяй и властвуй» базового курса.

Задачи

В первой свёртка счётчиков отвечает сразу на все вопросы про суммы. Во второй нужно сравнить две строки при всех сдвигах — и для этого придётся превратить совпадение букв в произведение чисел.