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

Быстрое преобразование Фурье (БПФ) - это алгоритм вычисления преобразования Фурье для дискретного  случая. В отличие от простейшего  алгоритма, который имеет сложность  порядка O(N²), БПФ имеет сложность всего лишь O(Nlog₂N). Алгоритм БПФ был впервые опубликован в 1965 году в статье Кули (Cooly) и Тьюки (Tukey).

Данное пособие  содержит исходный код работающей программы  для вычисления БПФ, подробное объяснение принципа ее работы и теоретическое  обоснование. Все это можно найти  и на других ресурсах, но трудно найти именно в таком комплекте: и программа, и объяснения, и теория, и на русском языке.

Если у вас  нет времени и желания разбираться  с теорией, то можете сразу скопировать  текст программы на C++. Здесь находится заголовочный файл fft.h и исходник fft.cpp для быстрого преобразования Фурье для числа отсчетов, равного степени двойки. Вызывать надо функцию fft. А здесь находится заголовочный файл и исходник для произвольного (!) числа отсчетов. Он чуть медленнее, но скорость там тоже порядка Nlog₂N. Вызывать надо функцию universal_fft.

Определение 1.

Дана конечная последовательность x₀, x₁, x₂,...,x_N-1 (в общем случае комплексных). Дискретное преобразование Фурье (ДПФ) заключается в поиске другой последовательности X₀, X₁, X₂,...,X_N-1 элементы которой вычисляются по формуле:

    (1).  

Определение 2.

Дана конечная последовательность X₀, X₁, X₂,...,X_N-1 (в общем случае комплексных). Обратное дискретное преобразование Фурье (ДПФ) заключается в поиске другой последовательности x₀, x₁, x₂,...,x_N-1 элементы которой вычисляются по формуле:

    (2).  

Основным свойством  этих преобразований (которое доказывается в соответствующих разделах математики) является тот факт, что из последовательности {x} получается (при прямом преобразовании) последовательность {X}, а если потом применить к {X} обратное преобразование, то снова получится исходная последовательность {x}.  

Определение 3.

Величина 

называется поворачивающим множителем.

Рассмотрим ряд  свойств поворачивающих множителей, которые нам понадобятся в  дальнейшем.

Верхняя цифра  в поворачивающем множителе не является индексом, это - степень. Поэтому, когда  она равна единице, мы не будем  ее писать:

Прямое преобразование Фурье можно выразить через поворачивающие множители. В результате формула (1) примет вид:

    (3).

Эти коэффициенты действительно оправдывают свое название. Нарисуем на комплексной  плоскости любое комплексное  число, в виде вектора, исходящего из начала координат. Представим это комплексное число в показательной форме: re^jφ, где r - модуль числа, а φ - аргумент. Модуль соответствует длине вектора, а аргумент - углу поворота:

Теперь возьмем  какой-нибудь поворачивающий множитель  . Его модуль равен единице, а фаза - 2π/N. Как известно, при умножении комплексных чисел, представленных в показательной форме, их модули перемножаются, а аргументы суммируются. Тогда умножение исходного числа на поворачивающий множитель не изменит длину вектора, но изменит его угол. То есть, произойдет поворот вектора на угол 2π/N (см. предыдущий рисунок).

Если теперь посмотреть на формулу (3), то станет ясен геометрический смысл преобразования Фурье: он состоит в том, чтобы  представить N комплексных чисел-векторов из набора {x}, каждое в виде суммы векторов из набора {X}, повернутых на углы, кратные 2π/N.  

Теорема 0.

Если комплексное  число представлено в виде e ^j2πN, где N - целое, то это число e ^j2πN = 1.

Доказательство:

По формуле  Эйлера, и ввиду периодичности синуса и косинуса:

e ^j2πN = cos(2πN) + j sin(2πN) = cos 0 + j sin 0 = 1 + j0 = 1  

Теорема 1.

Величина  периодична по k и по n с периодом N. То есть, для любых целых l и m выполняется равенство:

    (4).

Доказательство:

    (5)

Величина -h = -(nl+mk+mlN) - целая, так как все множители целые, и все слагаемые целые. Значит, мы можем применить Теорему 0:

Что и требовалось  доказать по (4).  

Теорема 2.

Для величины справедлива формула:

Доказательство:

Алгоритм быстрого преобразования Фурье (БПФ) - это оптимизированный по скорости способ вычисления ДПФ. Основная идея заключается в двух пунктах.

1. Необходимо разделить сумму (1) из N слагаемых на две суммы по N/2 слагаемых, и вычислить их по отдельности. Для вычисления каждой из подсумм, надо их тоже разделить на две и т.д.
2. Необходимо повторно использовать уже вычисленные слагаемые.

Применяют либо "прореживание по времени" (когда  в первую сумму попадают слагаемые  с четными номерами, а во вторую - с нечетными), либо "прореживание по частоте" (когда в первую сумму попадают первые N/2 слагаемых, а во вторую - остальные). Оба варианта равноценны. В силу специфики алгоритма приходится применять только N, являющиеся степенями 2. Рассмотрим случай прореживания по времени.  

Теорема 3.

Определим еще две последовательности: {x_[even]} и {x_[odd]} через последовательность {x} следующим образом:

x_[even]n = x_2n,  
x_[odd]n = x_2n+1,     (6)  
n = 0, 1,..., N/2-1

Пусть к этим последовательностям применены  ДПФ и получены результаты в виде двух новых последовательностей {X_[even]} и {X_[odd]} по N/2 элементов в каждой.

Утверждается, что  элементы последовательности {X} можно  выразить через элементы последовательностей {X_[even]} и {X_[odd]} по формуле:

    (7).

Доказательство:

Начинаем от формулы (2), в которую подставляем равенства из (6):

    (8)

Теперь обратим  внимание на то, что:

    (9)

Подставляя (9) в (8) получаем:

    (10)

Сравним с формулами  для X_[even]k и X_[odd]k при k = 0,1,…,N/2-1:

    (11)

Подставляя (11) в (10) получим первую часть формулы (7):

Это будет верно  при k = 0,1,…,N/2-1.

Согласно теореме 1:

    (12)

Подставим (12) в (10), и заменим  по теореме 2:

    (13)

Для k = N/2,…,N-1 по формуле (2):

    (14)

Подставляем (14) в (13):

Эта формула  верна для k = N/2,…,N-1 и, соответственно, (k - N/2) = 0,1,…,N/2-1 и представляет собой вторую и последнюю часть утверждения (7), которое надо было доказать.  

Формула (7) позволяет  сократить число умножений вдвое (не считая умножений при вычислении X_[even]k и X _[odd]k), если вычислять X_k не последовательно от 0 до N - 1, а попарно: X₀ и X_N/2, X₁ и X_N/2+1,..., X_N/2-1 и X_N. Пары образуются по принципу: X_k и X_N/2+k.  

Теорема 4.

ДПФ можно вычислить  также по формуле:

    (15)

Доказательство:

Согласно второй части формулы (7), получим:

  

Это доказывает второе равенство в утверждении  теоремы, а первое уже доказано в  теореме 3.  

Также по этой теореме  видно, что отпадает необходимость  хранить вычисленные X_[even]k и X_[odd]k после использования при вычислении очередной пары и одно вычисление можно использовать для вычисления двух элементов последовательности {X}.

На этом шаге будет выполнено N/2 умножений комплексных  чисел. Если мы применим ту же схему  для вычисления последовательностей {X_[even]} и {X_[odd]}, то каждая из них потребует N/4 умножений, итого еще N/2. Продолжая далее в том же духе log₂N раз, дойдем до сумм, состоящих всего из одного слагаемого, так что общее количество умножений окажется равно (N/2)log₂N, что явно лучше, чем N² умножений по формуле (2).

Рассмотрим БПФ для разных N. Для ясности добавим еще один нижний индекс, который будет указывать общее количество элементов последовательности, к которой этот элемент принадлежит. То есть X_{R}k - это k-й элемент последовательности {X_{R}} из R элементов. X_{R}[even]k - это k-й элемент последовательности {X_{R}[even]} из R элементов, вычисленный по четным элементам последовательности {X_{2R}}. X_{R}[odd]k - это k-й элемент последовательности {X_{R}[odd]}, вычисленный по нечетным элементам последовательности {X_{2R}}.

В вырожденном случае, когда слагаемое всего одно (N = 1) формула (1) упрощается до:

,

Поскольку в  данном случае k может быть равно  только 0, то X_{1}0 = x_{1}0, то есть, ДПФ над одним числом дает это же самое число.

Для N = 2 по теореме 4 получим:

Для N = 4 по теореме 4 получим:

Отсюда видно, что если элементы исходной последовательности были действительными, то при увеличении N элементы ДПФ становятся комплексными.

Для N = 8 по теореме 4 получим:

Обратите внимание, что на предыдущем шаге мы использовали степени W₄, а на этом - степени W₈. Лишних вычислений можно было бы избежать, если учесть тот факт, что

Тогда формулы  для N=4 будут использовать те же W-множители, что и формулы для N=8:

Подведем итог:

В основе алгоритма  БПФ лежат следующие  формулы:

    (16)

Теперь рассмотрим конкретную реализацию БПФ. Пусть имеется N=2^T элементов последовательности x_{N} и надо получить последовательность X_{N}. Прежде всего, нам придется разделить x_{N} на две последовательности: четные и нечетные элементы. Затем точно так же поступить с каждой последовательностью. Этот итерационный процесс закончится, когда останутся последовательности длиной по 2 элемента. Пример процесса для N=16 показан ниже:

Итого выполняется (log₂N)-1 итераций.

Рассмотрим двоичное представление номеров элементов и занимаемых ими мест. Элемент с номером 0 (двоичное 0000) после всех перестановок занимает позицию 0 (0000), элемент 8 (1000) - позицию 1 (0001), элемент 4 (0100) - позицию 2 (0010), элемент 12 (1100) - позицию 3 (0011). И так далее. Нетрудно заметить связь между двоичным представлением позиции до перестановок и после всех перестановок: они зеркально симметричны. Двоичное представление конечной позиции получается из двоичного представления начальной позиции перестановкой битов в обратном порядке. И наоборот.

Этот факт не является случайностью для конкретного N=16, а является закономерностью. На первом шаге четные элементы с номером n переместились в позицию n/2, а  нечетные из позиции в позицию N/2+(n-1)/2. Где n=0,1,…,N-1. Таким образом, новая позиция вычисляется из старой позиции с помощью функции:

ror(n,N) = [n/2] + N{n/2}

Здесь как обычно [x] означает целую часть числа, а {x} - дробную.

В ассемблере эта  операция называется циклическим сдвигом вправо (ror), если N - это степень двойки. Название операции происходит из того факта, что берется двоичное представление числа n, затем все биты, кроме младшего (самого правого) перемещаются на 1 позицию вправо. А младший бит перемещается на освободившееся место самого старшего (самого левого) бита.

 
рис. 1

Дальнейшие разбиения  выполняются аналогично. На каждом следующем шаге количество последовательностей  удваивается, а число элементов  в каждой из них уменьшается вдвое. Операции ror подвергаются уже не все  биты, а только несколько младших (правых). Старшие же j-1 битов остаются нетронутыми (зафиксированными), где j - номер шага: