Быстрое преобразование фурье алгоритм c

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

Есть и несколько подходов для подгонки алгоритмов под машинную архитектуру, один из которых – управление данными(data driven). Каждый, кто сталкивался с ручным программированием данных, знает, что это дело это не простое. Однако в большинстве случаев можно спроектировать пре-компилятор, который существенно упростит задачу. В статье описана методика построения двух управляемых данными алгоритмов БПФ и способы достижения максимальной производительности, превосходящей теоретическую.

Немного математики

N временных отсчетов любого сигнала могут быть преобразованы в N фазо-частот при помощи Прямого Дискретного Преобразования Фурье.

Для того чтобы получить комплексный фазо-частотный спектр нужно и временные отсчеты перевести в комплексный вид с нулевой мнимой частью. Поворачивающийся множитель ω представляет собой “вращающийся вектор” на единичной окружности в комплексной плоскости.
Для упрощения вычислений используется Быстрое Преобразование Фурье, где N является степенью от 2. В таком случае комплексный сигнал X=(x, x_{1 …} x_n-1) можно поделить на два A=(x, x_{2 …} x_n-2) и B=(x₁, x_{3 …} x_n-1) уже из N/2 отсчетов.

Теперь можно вывести нужную нам формулу: X_k = A_k + ω_N k * B_k;

Здесь N фазо-частот нашего сигнала можно получить, применяя БПФ по очереди для N/2 четных отсчетов и для N/2 нечетных отсчетов. На следующем шаге A и B опять разбиваются пополам. Стоит заметить, что этот подход применим к Обратному Преобразованию Фурье, которое здесь не рассматривается.

БПФ на C++

Рекурсивная реализация[1] БПФ по приведенным выше формулам проста и приведена ниже:

Действительные части комплексных значений входного вектора состоят из входных параметров, количество которых усекается до ближайшей степени 2-ки. Первый уровень рекурсии разбивает сигнал In из n-значений на 2 сигнала A и B, каждый из которых состоит из n/2 значений. Следующий уровень разбивает данные на 4 сигнала по n/4 значения. Рекурсия прекращается, когда остается только сигнал из одного значения. Сигналы A и B преобразуются в At и Bt, а затем в возвращаемый вектор сигналов Out. В конце программа выводит реальные и мнимые части фаза-частотного сигнала из выходного вектора.

Уровни абстракции

Чем выше уровень абстракции языка приложения, тем меньше усилий требуется для программирования. Но тем ниже производительность. Можно писать более эффективный код, снижая уровень абстракции, но тратя больше усилий. Разберем это на основе элементарной операции БПФ, которая может быть представлена как C = a + w * b в комплексном виде. Эта формула в общем случае реализуется через умножения с накоплением, поэтому имеет аббревиатуру MAC (multiply–accumulate operation):

Логарифм по основанию 2 от размера сигнала определяет число рекурсий. Несколько стилей n2ln функции приведены ниже.

В последней реализации наихудший случай исполнения равен среднему случаю, который С-компилятор может реализовать при помощи 3-х сдвигов и 4-х условных переходов. Кроме того эта реализация может быть представлена в виде макроса для расчета значений констант во время компиляции.

БПФ на С

Классическая реализация рекурсивного БПФ на С похожа на С++ вариант

Относительно С++ реализации здесь сделано несколько алгоритмических изменений. Во-первых, комплексные расчеты производятся вручную. Во-вторых, для экономии памяти входной и выходной вектор из реализации C++ объединены в один буфер. Буфера выделяются явно, потому размер преобразования нужно передавать в рекурсивную функцию.

С реализация БПФ на основе управления данными

Пример рекурсивного БПФ на основе управления данными, сделанный из предыдущей реализации на С.

Здесь буфера для вычислений содержит не сами комплексные значения, а указатели на них. В Mac массив фактически последовательно записываются отложенные элементарные БНФ операции, для того чтобы их сделать потом. Другими словами, это байт-код элементарных БНФ операций.

Двумерный массив DataTrace используется для поддержки этих операций. После вызова рекурсивной процедуры пользователь должен вызвать calculate_macs для исполнения сгенерированного байт-кода. Эта процедура имеет всего один цикл, но может применяться многократно. Это максимальная производительность для теоретической сложности БПФ n*log2 (n). Но проблема тут в памяти — Mac и DataTrace так же имеют n*log2 (n) элементов. И это слишком много для недорогих встроенных решений.

Табличная реализация

Теперь пришло время разделить предыдущую реализацию БПФ на две программы. Первая будет генерировать С структуры байт-кода, а вторая их исполнять. При таком подходе генератор С структур фактически является пре-компилятором, в котором можно, не жалея ресурсов, реализовывать различные оптимизационные стратегии, например, для оптимизации памяти RAM. Ранее и DataTrace массив N*log2 (N) элементов, в программе ниже его аналог массив XY имеет всего N+2 элементов.

Это пример БПФ для 16 отсчетов. Один элемент tbl массива содержит комплексное значение поворачивающего множителя и три смещения для организации вычислений на ячейках массива XY. При этом сам код имеет всего один “for” цикл.

Метод Гусеницы

Следующий пример основан на группировке элементарных БПФ операций относительно поворачивающего множителя.

Это пример БПФ для 32 отсчетов. Количество элементарных БПФ операций — N. Массив XY организован по схеме свопа и его размер — 2*N. Это дает возможность условной инструкции XMAC оперировать с 3-мя банками памяти одновременно, хотя на практике компиляторами это игнорируется. Тем не менее, XMAC может быть теоретически реализован даже одной машинной инструкцией.

Самый известный БПФ алгоритм использует сложную логику перестановки адресов, которая в графическом виде напоминает крылья бабочки. А вот в данном примере новые адреса получаются простым инкрементом от старых, потому с названием все просто — это метод гусеницы. Стоит заметить, что длинный switch оператор так же напоминает гусеницу.

Дальнейшие оптимизации

Элементарная комплексная операция БПФ (c = a + ω * b) состоит из 4-х обычных операций:

В литературе по БПФ [2] для этих вычислений можно найти несколько стратегий оптимизаций, которые классифицированы ниже:

• Оптимизации, основанные на особенностях входа/выхода БПФ

• Мнимые части входного сигнала нулевые
• Выходной фаза-частотный сигнал дублируется, по факту из него нужно N/2+1 значений

Почти все эти оптимизации БПФ применимы к методу гусеницы.

Возможны также архитектурно-машинные оптимизации. Например, элементарная БПФ операция в коде выше реализована как макрос XMAC. Она может быть распараллелена при помощи SIMD инструкций для x86-процессоров Intel AVX:

Приведенный макрос поддерживает две операции с плавающей точкой одновременно при помощи 128-битных регистров. Но стоит взглянуть на гусеницу еще раз – из-за особенностей адресации ближайшие XMAC-и могут быть объединены вместе, например, для реализации при помощи 512-регистров(AVX3).

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

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

Дискретное преобразование Фурье (ДПФ) периодического дискретного сигнала x ( n ) с периодом N определяется как

(12.1),

где  — основная частота преобразования ( бин ДПФ).

Выражение (12.1) можно переписать в виде

(12.2),

где (12.3).

Коэффициент W_N называется поворачивающим множителем. Легко показать, что является периодической функцией с периодом N

(12.4).

Поэтому ДПФ X ( k ) также является периодической функцией по аргументу k с периодом N .

Дискретное преобразование Фурье может быть использовано и для представления сигнала x ( n ) конечной длины N , определенного при n =0,1,…, N -1 и равного нулю вне интервала [0, N -1]. Действительно, такой сигнал можно рассматривать как один период соответствующего периодического сигнала и сипользовать преобразования (12.2). Следует только считать, что вне интнрвала [0, N -1] X ( k ) и x ( n ) равны нулю.

Если сравнить ДПФ конечного дискретного сигнала со спектром этого же сигнала, определяемым выражением

(12.5),

то очевидно, что ДПФ представляет собой N отсчетов спектра, взятых на периоде с интервалом дискретизации по частоте, равным .

В случае, когда x ( n ) является комплексным, при прямом вычислении N -точечного ДПФ нужно выполнить для каждого значения k ( N -1) умножений и ( N -1) сложений комплексных чисел или 4 ( N -1) умножений и 2 ( N -1) сложений действительных чисел. Для всех N значений k =0,1,…, N -1 требуется ( N -1) 2 умножений и N ( N -1) сложений комплексных чисел. Для больших значений N (порядка нескольких сотен или тысяч) прямое вычисление ДПФ по выражению (12.2) требует выполнения весьма большого числа арифметических операций умножения и сложения, что затрудняет реализацию вычислений в реальном масштабе времени.

Быстрым преобразованием Фурье (БПФ) называют набор алгоритмов, реализация которых приводит к существенному уменьшению вычислительной сложности ДПФ. Основная идея БПФ состоит в том, чтобы разбить исходный N -отсчетный сигнал x ( n ) на два более коротких сигнала, ДПФ которых могут быть скомбинированы таким образом, чтобы получить ДПФ исходного N -отсчетного сигнала.

Так, если исходный N -отсчетный сигнал разбить на два N /2-отсчетных сигнала, то для вычисления ДПФ каждого из них потребуется около ( N /2) 2 комплексных умножений. Тогда для вычисления искомого N -отсчетного ДПФ потребуется порядка 2 ( N /2) 2 = N 2 /2 комплексных умножений , т.е. вдвое меньше по сравнению с прямым вычислением. Операцию разбиения можно повторить, вычисляя вместо ( N /2) -отсчетного ДПФ два ( N /4) -отсчетных ДПФ и сокращая тем сасым объем вычислений еще в два раза. Выигрыш в два раза является приблизительным, поскольку не учитывается, каким образом из ДПФ меньшего размера образуется искомое N -отсчетное ДПФ.

Существует большое количество алгоритмов БПФ. Однако все они являются частными случаями единого алгоритма, базирующегося на задаче разбиения одного массива чисел на два. Тот факт, что это можно сделать более чем одним способом, определяет многообразие алгоритмов БПФ. Расмотрим два из них.

Первый алгоритм называется алгоритмом Б ПФ с пр ореживанием по времени.

Пусть задан N -отсчетный дискретный сигнал x ( n ). Примем, что N равно степени двойки. Если это не так, то всегда можно легко доплнить заданный сигнал нулевыми отсчетами до количества отсчетов, равного ближайшей степени двойки.

Разобъем исходный сигнал x ( n ) на два N /2-отсчетных сигнала x ₁ ( n ) и x ₂ ( n ), составленных соответственно из четных и нечетных отсчетов исходного сигнала x ( n )

(12.6).

N -точечное ДПФ сигнала x ( n ) можно записать как

(12.7).

С учетом того, что

(12.8)

(12.9)

(12.10),

где X ₁ ( k ) и X ₂ ( k ) – N /2-отсчетные ДПФ сигналов x ₁ ( n ) и x ₂ ( n ) соответственно.

Таким образом N -точечное ДПФ X ( k ) может быть разложено на два N /2-точечных ДПФ, результаты которых объединяются согласно (12.10).

Если бы N /2-точечные ДПФ вычислялись прямым способом, то для вычисления N -точечного ДПФ потребовалось бы ( N 2 /2+ N ) комплексных умножений. При больших N (когда N 2 /2>> N ) это позволяет сократить время вычислений на 50%.

Поскольку X ( k ) определено при , а X ₁( k ) и X ₂( k ) определены при , необходимо доопределить формулу (12.10) для . Учитывая, что X ₁( k ) и X ₂( k ) – периодические функции с периодом N /2, можно записать

(12.11),

поскольку .

Поэтому окончательно для N -точечного ДПФ можно записать

(12.12).

На рис.12.1 представлена последовательность операций при выполнении восьмиточечного ДПФ с использованием двух четырехточечных ДПФ.

Внимание! В данной версии класса используется прореживание по времени. Сигнал на входе БПФ формируется следующим образом: четные отсчеты подаются на действительную часть входа, нечетные – на мнимую, что позволяет значительно ускорить вычисления без потери качества. Если Вы используете класс, полученный с этого сайта до 23.04.2018, имейте в виду, что текст этой страницы скорректирован в соответствии с изменениями класса.

См. также примечание внизу, связанное с обновлением от 07.08.2017.

Класс FFT предназначен для встраивания в исходный код задач, реализованных на языке C ++, в которых необходимо проводить расчет спектра (АЧХ) сигнала по методу БПФ (быстрого преобразования Фурье).

Класс содержит следующие основные переменные и функции, используемые программистом при встраивании кода:

Переменные на входе БПФ:

N – количество точек спектра (спектр строится на основе 2* N измерений сигнала).

T – степень числа 2, определяющая количество точек спектра N (то есть N = 2 ^ T ).

Xre [ T ][], X im [ T ][] – массивы измерений сигнала на входе (действительная и мнимая части; в действительную часть подаются четные отсчеты, в мнимую — нечетные).

t m[ ] – массив значений моментов времени измерений сигнала.

No – номер крайней пары значений (чет ./ нечет.) в накопленном массиве измерений сигнала. Изменяется от 0 до N – 1; после появления N -й пары значений сигнала сохраняется значение No = N – 1.

CFreq – минимальная частота регистрации сигнала на входе, при которой производится расчет спектра. В некоторых случаях, из-за сбоев или задержек в каналах измерений, частота регистрации падает. Проводить расчет спектра в этих случаях не имеет смысла.

Otsechka – признак отказа от расчета спектра при частоте регистрации, меньшей CFreq .

Flt – простейший сглаживающий фильтр, предназначенный для сглаживания значений моментов времени. Такое сглаживание имеет смысл применять в том случае, если известно, что регистрация сигнала производится с постоянной частотой, но, по причине задержек в канале, частота регистрации «гуляет». Фильтр имеет передаточную функцию 1 / ( Tp + 1). Чтобы отключить фильтр, достаточно при расчете переменной dT вместо dT = Flt.Go ( tm [N - 1] — tm [0] + dt ); написать dT = tm [N - 1] — tm [0] + dt ;

fH , fT – параметры сглаживающего фильтра.

Переменные на выходе БПФ:

A[] – набор амплитуд спектра.

f req [ ] – набор частот спектра.

Все остальные переменные не используются программистом при работе с классом.

Prepare ( ) – функция расчета поворачивающих множителей; выполняется один раз перед началом работы с поступающим сигналом. Набор поворачивающих множителей зависит только от заданного числа T , поэтому данная функция запускается один раз перед расчетами спектра по поступающему сигналу.

PutVal ( ) – функция подготовки данных в массивах tm и Xre и в других переменных перед очередным расчетом спектра.

Spectrum ( ) – расчет спектра.

Встраивание класса FFT на конкретном примере

Проект testfft создан в среде Borland Builder 6 и предназначен для опробования класса FFT : пользователь может задавать амплитуды и частоты 3-х гармоник и наблюдать на экране спектр получившегося сигнала при разных количествах точек спектра и разных частотах регистрации. Параметр T , задаваемый с формы, определяет количество точек сигнала, используемых для построения спектра, то есть T спектра будет на 1 меньше, чем T сигнала.

Исходный код проекта состоит из 4-х unit ’ ов :

umain – пользовательский интерфейс; здесь производится ввод параметров, создание нити ( thread ), объекта fft класса FFT , их инициализация и запуск, прием и вывод на форму поступающей от этих объектов информации;

thr – нить, формирующая и передающая информацию объекту fft ;

fft – здесь находится класс FFT ;

flt – здесь находится сглаживающий фильтр.

В обработчике нажатия кнопки «СТАРТ» производится создание и инициализация объектов Thr (нити) и fft (БПФ), запуск функции Prepare () (расчет поворачивающих множителей) и запуск нити функцией Resume ().

Обработчик закрытия формы останавливает нить и удаляет созданные объекты.

Таймер перерисовывает график, отображая на нем текущие результаты расчета спектра.

Функция Look () выполняется с частотой, зависящей от частоты вашего процессора. В ней выполняется расчет гармоник и значений времени (по заданной частоте дискретизации samplingFR ) и передача этих данных алгоритму БПФ для очередного расчета спектра.

Примечание (обновление класса)

В практических приложениях часто возникает необходимость рассчитывать АЧХ, изменяющуюся в режиме онлайн. Расчет можно проводить с разным шагом. Если, например, получив данные сигнала X , …, X_{N –1} , мы получаем на выходе первый набор амплитуд и частот, а получив следующее наблюдение X_N , мы рассчитываем по X ₁ , …, X_N следующий набор – то, в таком случае, мы считаем АЧХ с шагом 1.

В случае достаточно большой частоты регистрации считать АЧХ на каждом шаге нерационально, т.к. это занимает слишком много времени. Поэтому функция класса FFT PutVal ( void ) заменена на PutVal ( bool c ) . Если c = true , функция сообщает алгоритму БПФ вновь поступившие данные и рассчитывает АЧХ. Если же c = false , происходит передача данных, но расчет не производится.

В класс TThr нашего примера добавлена переменная step , получающая значение с формы (поле Шаг расчета). Как можно увидеть в исходниках примера, функция PutVal ( c ) запускается с переменной c = true с частотой, заданной step , что ускоряет расчет и позволяет выводить на экран изменяющийся спектр сигнала в реальном времени или еще быстрее, в зависимости от задачи.