Файл: В. П. Кандидов и др. Москва физический факультет мгу.pdf

Глава 5. Быстрое преобразование Фурье
69
Переставим отчеты в соответствии с двоичной инверсией, что соответствует следующей замене номеров последовательности
0 0
=
→
=
j
j
,
2 1
=
→
=
j
j
,
1 2
=
→
=
j
j
,
3 3
=
→
=
j
j
. Учтем далее, что для N = 4 в силу периодичности гармонических функций W справедливы равенства
kj
n
kj
W
W
−
+
−
=
)
4
(
, где
K
,
2
,
1
,
0
±
±
=
n
. Тогда прямое преобразование Фурье принимает вид
[
]
)
3
(
)
1
(
)
2
(
)
0
(
4 1
)
0
(
u
u
u
u
U
+
+
+
=
,
[
]
3 1
2
)
3
(
)
1
(
)
2
(
)
0
(
4 1
)
1
(
−
−
−
+
+
+
=
W
u
W
u
W
u
u
U
,
[
]
2 2
)
3
(
)
1
(
)
2
(
)
0
(
4 1
)
2
(
−
−
+
+
+
=
W
u
W
u
u
u
U
,
[
]
5 3
2
)
3
(
)
1
(
)
2
(
)
0
(
4 1
)
3
(
−
−
−
+
+
+
=
W
u
W
u
W
u
u
U
Вычислим далее спектры разреженных отсчетов. В соответствии с (5.6), где
1
,
0 1
=
j
и
1
,
0 2
=
k
, имеем
[
]
0
,
0
,
)
1
,
0
(
)
0
,
0
(
2 1
)
0
,
0
(
2 1
=
=
+
=
k
j
u
u
U
,
[
]
0
,
1
,
)
1
,
1
(
)
0
,
1
(
2 1
)
0
,
1
(
2 1
=
=
+
=
k
j
u
u
U
,
[
]
1
,
0
,
)
1
,
0
(
)
0
,
0
(
2 1
)
1
,
0
(
2 1
2
=
=
+
=
−
k
j
W
u
u
U
,
[
]
1
,
1
,
)
1
,
1
(
)
0
,
1
(
2 1
)
1
,
1
(
2 1
2
=
=
+
=
−
k
j
W
u
u
U
Затем по (5.7) находим все фурье-гармоники U(k) где k = 0, 1, 2, 3:
[
]
)
0
,
1
(
)
0
,
0
(
2 1
)
0
(
U
U
U
+
=
,
[
]
1
)
1
,
1
(
)
1
,
0
(
2 1
)
1
(
−
+
=
W
U
U
U
,
[
]
2
)
0
,
1
(
)
0
,
0
(
2 1
)
2
(
−
+
=
W
U
U
U
,
[
]
3
)
1
,
1
(
)
1
,
0
(
2 1
)
3
(
−
+
=
W
U
U
U

Глава 5. Быстрое преобразование Фурье
70
Однократное вычисление сумм повторяющихся слагаемых и хранение их готовом виде в памяти ЭВМ позволяет достичь экономии числа операций при использовании БПФ по сравнению с ДПФ, причем эта экономия тем существеннее, чем больше N. Наглядное представление о структуре вычислений при использовании алгоритма БПФ можно получить с помощью рис. 5.5 (N = 4) и рис. 5.6 (N = 8).
Рис. 5.5. Графическая схема прямого ДПФ (слева) и БПФ (справа) при N = 4
Рис. 5.6. Графическая схема БПФ при N = 8

Глава 6. Компьютерные библиотеки ДПФ
71
Глава
6. Компьютерные библиотеки ДПФ
§6.1. Свободно распространяемая библиотека FFTW
FFTW – это библиотека подпрограмм на языке Cи для вычисления дискретного преобразования Фурье. Она содержит подпрограммы как для одномерного преобразования, так и многомерного преобразования с произвольным размером обрабатываемых массивов как действительных, так и комплексных чисел. Библиотека является свободно распространяемым программным продуктом и поэтому доступна для самого широкого круга пользователей.
Библиотека разработана в
Массачусетском технологическом институте Маттео Фриго и Стивеном Джонсоном.
Абревиатура FFTW расшифровывается как Fastest Fourier Transform in
the West. Домашней страницей библиотеки является сайт www.fftw.org
, где она доступна для скачивания, включая написанное авторами подробное пособие [
http://www.fftw.org/fftw3.pdf
], часть которого была использована для написания данной главы. Цифра 3 в названии файлов соответствует главному номеру версии библиотеки.
Будучи написанной на языке Си, библиотека является переносимым продуктом и на большинстве компьютеров её использование не потребует каких-либо модификаций, сохраняя высокую производительность – одну из лучших среди универсальных программ для быстрого преобразования Фурье. Эта производительность во многом достигается за счет того, что конкретные алгоритмы вычислений не являются раз и навсегда заданными, а подбираются исходя из имеющихся вычислительных ресурсов и структуры обрабатываемых данных. В частности, могут быть использованы такие расширения системы команд для микропроцессоров x86, как SSE или
AVX и другие. Такой подход обуславливает общую двухэтапную процедуру применения программ библиотеки. На первом этапе вызывается программный модуль планировщик, который оптимизирует внутренние параметры алгоритмов для последующих расчётов. Его результатом является некоторая структура данных, так называемый
план, который используется для более эффективного выполнения БПФ

Глава 6. Компьютерные библиотеки ДПФ
72 на втором этапе. Для массивов фиксированной длины планировщик достаточно вызвать только один раз
1
Несмотря на то, что библиотека поддерживает ДПФ массивов произвольной длины, включая случай, когда длина является простым числом, следует помнить, что стандартный дистрибутив FFTW обеспечивает наиболее эффективное преобразование для массивов, длина которых является целой степенью 2, 3, 5 или 7. Далее, следуя руководству [
http://www.fftw.org/fftw3.pdf
], мы рассмотрим примеры использования библиотеки для выполнения одномерного и двумерного дискретных преобразований Фурье на языке Си. Они представляют собой, вероятно, наиболее востребованные задачи применения ДПФ и дают общее представление о том, как устроена и как может быть использована данная библиотеки. Конечно, эти примеры далеко не исчерпывают все возможности библиотеки FFTW. Например, она позволяет выполнять
ДПФ произвольной размерности, есть параллельные версии библиотеки как для систем с разделяемой памятью (используют OpenMP), так и с распределённой (используют
MPI).
1. Установка библиотеки FFTW на компьютере
Если вы используете Unix систему, такую как, например,
GNU/Linux, то процесс установки библиотеки FFTW будет для вас, наверняка, прост. Сама библиотека изначально создавалась для таких систем. Последнюю версию
2
исходного кода в архивированном виде можно скачать со странички сайта http://www.fftw.org/download.html
Архив снабжен конфигурационным файлом для утилиты make, а подробные инструкции по установке можно найти здесь http://www.fftw.org/fftw3_doc/Installation-on-Unix.html#Installation-on-Unix
Несколько менее очевидным может оказаться процесс установки для пользователей Microsoft Visual Studio на Windows системах. Для упрощения ситуации авторы библиотеки разместили для скачивания предкомпилированную версию FFTW
3
в виде DLL-
1
На тот случай, если время выполнения планировщика является критичным, библиотека предлагает быстрые (но, возможно, с неоптимальным результатом) планировщики, основанные на эвристике или результатах предыдущих планирований.
2
На момент написания пособия это была версия FFTW 3.3.8 3
На этот же момент Windows-версия имела номер FFTW 3.3.5

Глава 6. Компьютерные библиотеки ДПФ
73 библиотеки. Далее мы рассмотрим процесс установки и пример использования библиотеки FFTW в среде Visual Studio 12.
Скачать DLL-библиотеку FFTW и заголовочные файлы к ней можно на странице http://www.fftw.org/install/windows.html
. На этой странице размещены два архива с библиотеками, соответственно, для
32-х и 64-х–битных приложений. Скачиваем нужный вариант и разархивируем его в подходящей директории. В архиве содержатся библиотеки для вычислений с одинарной, двойной (стандартный вариант) и четверной (long double) точностью. Мы будем использовать стандартный вариант – файл libfftw3-3.dll. Наличие этой библиотеки будет важно на момент запуска исполняемого файла приложения, использующего FFTW. Во время сборки приложения («построения решения») нужна другая, статическая библиотека
*.lib
(«импортированная библиотека»), содержащая информацию о библиотечных функциях FFTW. Эту библиотеку нужно предварительно создать.
Для создания импортированной библиотеки следует воспользоваться утилитой lib из пакета Visual Studio. Для её запуска в командном режиме можно, например, использовать утилиту
«Командная строка разработчика для VS2012», которая входит в состав дистрибутива с Visual Studio 12 и по умолчанию может быть запущена из той же папки меню «Пуск», что и сама студия.
После запуска командной строки разработчика, в ее консольном окне с помощью команд выбора диска и смены директории (cd) следует выбрать текущей ту директорию, в которой находятся ранее разархивированные файлы библиотеки FFTW. В этой директории исполняется команда lib: lib /machine:x86 /def:libfftw3-3.def
Опция /machine:x86 этой команды является опцией по умолчанию и предназначена для создания 32-х битных приложений.
Если вы скачали архив FFTW для 64-х битных приложений, то в этой опции «86» надо заменить на «64». При успешном выполнении команды в директории с FFTW должен появиться файл импортированной библиотеки libfftw3-3.lib, которую нужно будет включить в проект приложения.
Предположим, что мы хотим использовать библиотеку FFTW в консольном приложении, создаваемом в Visual Studio 12. Добавить в проект импортированную библиотеку можно, например, используя следующие опции меню:
Проект → Свойства → Свойства
конфигурации → Компоновщик → Все параметры → Дополнительные

Глава 6. Компьютерные библиотеки ДПФ
74
зависимости, где в отдельном окне вводится имя библиотеки libfftw3-
3.lib. Кроме этого, в окно ввода опции Все параметры → Дополнитель-
ные каталоги библиотек следует добавить имя директории, содержа- щей эту библиотеку.
На этом процесс подключения библиотеки FFTW к проекту можно считать завершённым. Саму же динамическую библиотеку libfftw3-3.dll удобно просто скопировать в ту директорию, откуда будет запускаться приложение.
Если запуск будет происходить непосредственно из среды Visual Studio, то скопировать в текущую директорию проекта.
2. Одномерное ДПФ массива комплексных чисел
Дискретное преобразование Фурье в библиотеке FFTW выполняется в ненормализованном виде. Прямое преобразование имеет вид:

−
=
−
=
1 0
)
(
)
(
N
j
kj
N
W
j
u
k
U
,
(6.1) где






π
=
−
=
−
=
N
i
W
N
j
N
k
N
2
exp
,
1
,
1
,
0
,
1
,
1
,
0
K
K
В отличие от формулы (5.1) оно не содержит множителя
N
/
1
, при этом обратное преобразование алгоритмически отличается от прямого только знаком в показателе экспоненты:

−
=
=
1 0
)
(
)
(
N
k
kj
N
W
k
U
j
u
(6.2)
Таким образом, последовательное применение прямого и обратного преобразований приведёт в итоге к умножению входных данных на N.
Ниже приведён фрагмент типовой программы вычисления ДПФ одномерного массива комплексных чисел из N элементов на языке Си и комментарии к ней.
#include
{ fftw_complex *in, *out;

Глава 6. Компьютерные библиотеки ДПФ
75 fftw_plan plan; in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); plan = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD,
FFTW_ESTIMATE); fftw_execute(plan); /* repeat as needed */ fftw_destroy_plan(plan); fftw_free(in); fftw_free(out);
}
Эту программу нужно скомпоновать (слинковать) вместе с библиотекой
FFTW и стандартной математической библиотекой. На Unix системах для этого следует использовать параметры -lfftw3 -lm.
Использование библиотеки FFTW требует включение в программу заголовочного файла fftw3.h. Собственно вычисление ДПФ предполагает следующую последовательность действий, которые подробно описаны ниже:
1.
Выделение памяти для входного и выходного массивов.
2.
Вызов планировщика.
3.
Задание входного массива и вычисление ДПФ необходимое количество раз.
4.
Освобождение выделенных ресурсов.
1. Комплексное ДПФ производится над массивом, тип элементов которых fftw_complex или double[2]
4
. Т.е., фактически, обрабатываемые данные структурируются в виде двумерного массива типа double in[N][2], где in[i][0] это реальная часть i-го комплексного числа, а in[i][1] – мнимая часть. Память для входного массива выделяется с помощью вызова библиотечной функции fftw_malloc(), которая делает это аналогично стандартной Си-функции malloc(), но
4
Используется следующее определение для синонима типа: typedef double fftw_complex[2].

Глава 6. Компьютерные библиотеки ДПФ
76 гарантирует правильное выравнивание элементов массива, если для ускорения вычислений используются SIMD-инструкции процессора.
Аналогично выделяется память и для массива выходных данных out. При желании можно ограничиться одним и тем же массивом для входных и выходных данных. В этом случае результат
ДПФ будет записан поверх входных данных (они не сохранятся).
2. Следующим шагом идёт вызов функции-планировщика, результатом работы которой будет объект (план), содержащий все необходимые данные для выполнения БПФ: fftw_plan_dft_1d(int N, fftw_complex *in, fftw_complex *out, int sign, unsigned flags).
Первый аргумент N – это число элементов входного (и выходного) массива.
Далее идут два указателя на входной и выходной массивы (могут быть одинаковыми).
Четвёртый аргумент, sign, может быть либо FFTW_FORWARD (–1) для прямого преобразования, либо FFTW_BACKWARD (+1) – для обратного преобразования.
Последний аргумент, flags, как правило, имеет значение либо
FFTW_MEASURE, либо FFTW_ESTIMATE. В случае FFTW_MEASURE планировщик выполняет тестовые вычисления, используя массивы in и out с целью подбора оптимальных параметров. Поэтому вызов планировщика должен быть выполнен до инициализации массива in, т.к. в противном случае после вызова планировщика исходные данные во входном массиве будут утеряны. В случае FFTW_ESTIMATE никакие преобразования не выполняются, подходящий план строится быстрее, но он может оказаться не оптимальным.
3. После построения плана он может быть использован многократно для выполнения ДПФ с указанными массивами in и out посредством вызова функции fftw_execute(plan): void fftw_execute(const fftw_plan plan).
Результат преобразования размещается в массиве out в порядке возрастания номера гармоники от 0 до N – 1. Те, кому более привычно оперировать положительными и отрицательными частотами (см. например, формулу (3.3)), могут считать, что в первой половине массива располагаются неотрицательные гармоники от 0 до N/2, а затем

Глава 6. Компьютерные библиотеки ДПФ
77 идут отрицательные гармоники с номерами
5
от –(N/2 – 1) до –1. Если in ≠ out, то входной массив остаётся неизменным.
Если есть необходимость выполнить обратное или такое же преобразование, но с другим массивом, то надо заново создать план с помощью функции fftw_plan_dft_1d. При возможности этот новый план будет создан на основе предыдущего.
4. Если данный план больше не будет использоваться, то выделенную для него память следует освободить с помощью вызова функции fftw_destroy_plan(plan). Если память для массивов была получена с помощью fftw_malloc(), то для ее освобождения необходимо использовать функцию fftw_free().
В качестве примера использования библиотеки
FFTW рассмотрим комплексное ДПФ следующего ряда из N = 4 чисел: u(0) = 1,
u(1) = 2,
u(2) = 3,
u(3) = 4.
Приведенная ниже программа последовательно делает прямое преобразование Фурье, вычисляет спектральную плотность мощности (с точностью до некоторого множителя) и делает обратное преобразование:
#include
"fftw3.h"
#include

#define
N
4 int main()
{ fftw_complex
*in, *out; fftw_plan plan, plan_inverse; int i; in = (
fftw_complex
*) fftw_malloc(
sizeof
(
fftw_complex
) *
N
); out = (
fftw_complex
*) fftw_malloc(
sizeof
(
fftw_complex
)
*
N
); plan = fftw_plan_dft_1d(
N
, in, out,
FFTW_FORWARD
,
FFTW_MEASURE
); plan_inverse = fftw_plan_dft_1d(
N
, out, in,
FFTW_BACKWARD
,
FFTW_ESTIMATE
); for
(i=0; i <
N
; i++)
5
В силу периодичности преобразования гармоника с номером –k (соответствует частоте –k/N) совпадает с гармоникой с номером N – k (частота (N – k)/N).

Глава 6. Компьютерные библиотеки ДПФ
78
{ in[i][0]= i+1.;
// {1.,2.,3.,4.}
in[i][1]=0.;
} fftw_execute(plan);
// DFT printf(
"\nForward DFT:\n"
); for
(i=0; i <
N
; i++)
{ printf(
"%f %fi\n"
,out[i][0],out[i][1]);
// DFT
{1.,2.,3.,4.}
} printf(
"\nPower spectral density\n"
); for
(i=0; i <
N
; i++)
{ printf(
"%f\n"
,out[i][0]*out[i][0]+out[i][1]*out[i][1]);
// PSD
} fftw_execute(plan_inverse);
// Inverse DFT printf(
"\nBackward(Forward DFT):\n"
); for
(i=0; i <
N
; i++)
{ printf(
"%f %fi\n"
,in[i][0],in[i][1]);
//
Inverse(Direct DFT)
} fftw_destroy_plan(plan); fftw_destroy_plan(plan_inverse); fftw_free(in); fftw_free(out); return
0;
}
В этой программе директива
#include
"fftw3.h" предполагает, что заголовочный файл находится в текущей директории проекта. Если это не так, и файл не размещен в системных директориях с заголовочными файлами, то к нему следует указать полный путь. В нашей программе используются два плана, один для прямого

Глава 6. Компьютерные библиотеки ДПФ
79 преобразования (plan) и другой для обратного (plan_inverse). При этом для прямого преобразования входным является файл in, а выходным (с результатом преобразования) – файл out. Обратное преобразование выполняется с массивом данных, полученных в итоге прямого преобразования, поэтому эти файлы меняются местами. Так же эта программа может служить иллюстрацией использования в языке Си комплексного типа данных fftw_complex из библиотеки FFTW.
Результат выполнения программы в консольном окне следующий:
Forward DFT:
10.000000 0.000000i
-2.000000 2.000000i
-2.000000 0.000000i
-2.000000 -2.000000i
Power spectral density
100.000000 8.000000 4.000000 8.000000
Backward(Forward DFT):
4.000000 0.000000i
8.000000 0.000000i
12.000000 0.000000i
16.000000 0.000000i
Обратите внимание, что после обратного преобразования получается массив, числа в котором в N раз больше исходных данных.
Для приведения результата обратного преобразования к виду (5.2) необходимо выполнить деление элементов на N.
3. Комплексное ДПФ двумерных массивов
Двумерное дискретное преобразование Фурье выполняется аналогично одномерному. Также требуется выделить память для fftw_complex массивов (предпочтительно используя fftw_malloc), создать план, выполнить необходимое число раз преобразование с созданным планом и освободить выделенные ресурсы. Планировщик