[Из песочницы] Преобразование Фурье. The Fast and the Furious28.12.2018 20:04

Зачастую при разработке алгоритмов мы упираемся в предел вычислительной сложности, который, казалось бы, преодолеть невозможно. Преобразование Фурье имеет сложность $O(n^2)$ , а быстрый вариант, предложенный около 1805 года Гаусом¹ (и переизобретенный в 1965 году Джеймсом Кули и Джоном Тьюки) $O(nlog(n))$ . В данной статье хочу вам показать, что можно получить результаты преобразования за линейное время $O(n)$ или даже достичь константной сложности $O(1)$ при определенных условиях, которые встречаются в реальных задачах.

Когда передо мной встала задача написания программы анализа передаточных функций звуковых систем в реальном времени, я, как и все, сперва обратился к быстрому преобразованию. Всё было хорошо, но при больших размерах временного окна загрузка CPU становилась неприлично большой и с этим надо было что-то делать. Было решено сделать паузу и изучить преобразование еще раз, а заодно поискать способы решения проблемы. Вернемся к оригинальному преобразованию Жозефа Фурье2:

$f(x) = \sum\limits_{-\infty}^{+\infty} c_ke^{2\pi ikx/T} \\ c_k = \frac{1}{T}\int\limits_0^Tf(x)e^{-2\pi ikx / T}dx$

Посмотрим внимательно, что же тут происходит. Каждое выходное значение в частотной области $c_k$ является суммой всех входных значений сигнала $f(x)$ умноженных на $e^{-2\pi ikx / T}$ . Чтобы произвести вычисления нам надо для каждого выходного значения пройтись по всем входным данным, т.е. выполнить те самые $n^2$ операций.

Избавляемся от n

Напомню, что изначально была задача анализа звуковых данных в реальном времени. Для этого выбранное временное окно (по сути буфер) размером N заполняется данными с частотой fd соответствующей частоте дискретизации. С периодом Т происходят преобразования входных данных из временного окна в частотное. Если посмотреть на реальные числа, то N варьируется от 214 (16 384) до 216 (65 536) сэмплов (значения остались в наследство от БПФ, где размер окна обязан быть степенью двойки). Время T = 80 мс (12,5 раз в секунду), что позволяет видеть достаточно удобно изменения и не слишком нагружать CPU и GPU. Частота дискретизации fd стандартная и составляет 48кГц. Давайте посчитаем, сколько же данных во временном окне изменяется между измерениями. За время Т в буфер поступает $48000 * 0,08 = 3840$ сэмплов. Таким образом в окне обновляется только от 5% до 23% данных. В худшем случае 95% (а в лучшем 73%, что тоже немало!) обрабатываемых сэмплов будут снова и снова попадать в преобразование несмотря на то, что уже были обработаны в прошлых итерациях.

Внимательный читатель в этот момент поднимет руку и скажет: «постойте, а как же коэффициент $e^{-2\pi ikx / T}$ ? Ведь при каждом новом преобразовании те же данные будут располагаться на новых позициях ряда, а как следствие иметь разные коэффициенты?» Каждому пятерку за внимательность, давайте вспоминать немаловажную деталь в преобразовании о которой часто забывают. При исследовании значений функции $f(t)$ на интервале от 0 до t, функция считается периодической, что позволяет безболезненно производить сдвиг функции влево или вправо во времени. Как следствие, мы имеем полное право не вставлять в конец новое значение и удалять из начала старое, а циклически заменять данные в буфере.

Для наглядности, можно записать в табличном виде как будет изменяться буфер:

t=0	f (0)	f (1)	f (2)	f (3)	f (4)	f (5)	f (6)	f (7)	f (8)	f (9)
t=1	f (10)	f (1)	f (2)	f (3)	f (4)	f (5)	f (6)	f (7)	f (8)	f (9)
t=2	f (10)	f (11)	f (2)	f (3)	f (4)	f (5)	f (6)	f (7)	f (8)	f (9)
t=3	f (10)	f (11)	f (12)	f (3)	f (4)	f (5)	f (6)	f (7)	f (8)	f (9)
t=4	f (10)	f (11)	f (12)	f (13)	f (4)	f (5)	f (6)	f (7)	f (8)	f (9)

Можно записать как изменяется преобразование во времени от t1 до t2:

$F_t = F_{t-1} + \Delta F \\ \Delta F : \Delta c_k = \frac{1}{T}\int\limits_{t_1}^{t_2}(f_t(x) - f_{t-1}(x))e^{-2\pi ikx / T}dx$

Значение $F_{t-1}(x)$ является результатом предыдущего преобразования, а сложность вычисления $\Delta f(x)$ не зависит от размера временного окна и, следовательно, является константным. В итоге сложность преобразования будет $O(n)$ *, т.к. всё что нам остается, это один раз пройтись по частотному окну и применить изменения для изменившихся за время Т сэмплов. Отдельно хочу обратить ваше внимание, что коэффициенты $e^{-2\pi ikx / T}$ могут быть посчитаны заранее, что дает дополнительный выигрыш в производительности, а в цикле останется только две операции: вычитание вещественных чисел и умножение вещественного числа на комплексное, на практике обе этих операции просты и дёшевы.

Для полноты картины осталось лишь обозначить начальное состояние, но тут всё просто:

$F_0(x) = 0$

*- конечно, конечная сложность всего преобразования так и останется $O(n^2)$ , но оно будет выполнено постепенно, за n итераций, пока обновляется буфер. $O(n)$ — это сложность обновления данных, но именно это нам и нужно (при использовании БПФ сложность каждого преобразования $O(nlog(n))$ ).

А что, если копнуть глубже. Или избавляемся от второго n

Сразу хочу оговориться, что дальнейшие шаги применимы только в том случае, если вы не планируете для результата производить обратное преобразование (с целью коррекции сигнала или получения импульсной характеристики). Для начала, хочу напомнить, что в результате преобразования мы получаем симметричный массив данных, что сходу позволяет уменьшить количество преобразований в два раза.

А теперь давайте проанализируем результирующий набор данных, учитывая условия задачи. Мы имеем набор комплексных чисел, каждое из которых описывает амплитуду и фазу колебаний на определенной частоте. Частоту можно определить по формуле: $f[j] = j \frac{fd}{N}$ для $j < \frac{N}{2}$ . Давайте оценим шаг частотного окна на наших данных: $\Delta f = \frac{fd}{N}$ Для N = 214: 2,93Гц (а для 216: 0,73Гц). Таким образом в диапазоне от 1кГц до 2кГц мы получим 341 результат. Попробуйте самостоятельно оценить сколько данных будет в диапазоне от 8кГц до 16кГц для N = 65536. Много, правда? Очень много! Нужно ли нам столько данных? Конечно в задачах отображения частотных характеристик звуковых систем, ответ: нет. А с другой стороны, для анализа в области низких частот, маленький шаг очень полезен. Не стоит забывать и о том, что впереди еще графика, которая должна эти объемы ( $\frac{N}{2}$ ) преобразовать к понятному человеку виду (усреднение, сплайн или сглаживание) и вывести их на экран. Да и на высоких частотах даже имея 4К экран и выводя график в полноэкранном режиме при логарифмической оси частот размер шага быстро окажется намного меньше 1 пикселя.

Опытным путем можно выяснить, что вполне достаточно иметь лишь 48 точек на октаву, а чтобы иметь данные чуть более сглаженными и усредненными, предлагаю остановиться на значении 96. В звуковом диапазоне частот от 20Гц до 20кГц нетрудно насчитать всего 10 октав: 20, 40, 80, 160, 320, 640, 1280, 2560, 5120, 10240, 20480, каждую из которых можно разделить на заданное количество поддиапазонов (не забывайте, что разбиение следует производить геометрически, а не арифметически), следовательно, более чем достаточно произвести преобразование только для 960 частот, чтобы получить результат, что в 16…65 раз меньше изначального варианта.

Таким образом, комбинируя оба подхода, получаем константную сложность выполнения алгоритма обновления данных $O(1)$ .

Мёд в квадрате и ложка дёгтя

Вот теперь можно смело заявить, что от сложности $O(n^2)$ мы пришли к сложности $O(1)$ использовав два простых лайфхака:

проанализировав задачу, заметили, что данные добавляются постепенно, а период полного обновления временного окна значительно выше периода преобразований и перешли к вычислению разницы преобразования Фурье.
ушли от арифметического шага в частотном окне к ограниченному лишь заданными значениями, что позволяет разительно снизить количество преобразований.

Но, конечно, жизнь была бы действительно сказкой, если бы не одно но. Применение этих двух подходов позволило действительно разгрузить CPU так, что догадаться о том, что он рассчитывает преобразование Фурье и выводит результаты на экран даже при $N=2^{20}$ было сложно. Но кара не заставила себя ждать, когда ваши сигналы в реальности не периодичны (а это обязательно для получения верных результатов преобразования) и подобрать подходящий размер окна не представляется возможным, возникает необходимость использования различных оконных функций, что уже не позволяет использовать полноценно первый шаг. Практика показала, что применение оконных функций критично при исследовании сигналов с частотой меньше, чем $0,1f_d$ . На больших частотах количество периодов, попавших во временное окно, значительно и ослабляет искажения, возникающие в следствии наличия разрыва первого рода (между f (0) и f (N-1)) в изначальной функции.

От второго шага я тоже в итоге отказался и вернулся обратно к БПФ, т.к. выигрыш в данной задаче уже был невелик.

В заключение

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

Увы, для меня в данной задаче, это осталось лишь небольшим математическим развлечением, но я надеюсь, что оно вдохновит вас заняться на праздниках исследованием других алгоритмов с точки зрения изменений входных данных во времени :)

Литература

Cooley–Tukey FFT algorithm
Е.А. Власова Ряды. Издательство МГТУ им Н.Э. Баумана. Москва. 2002

Изображение взято из манги Митио Сибуя «ЗАНИМАТЕЛЬНАЯ МАТЕМАТИКА. АНАЛИЗ ФУРЬЕ»