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

Глава 6. Компьютерные библиотеки ДПФ
80 для двумерного комплексного массива размера N0×N1 имеет следующий прототип: fftw_plan fftw_plan_dft_2d(int N0, int N1,fftw_complex *in, fftw_complex *out, int sign, unsigned flags). где N0 – число строк, а N1 – число столбцов в двумерном массиве.
Предполагается, что массив комплексных чисел типа fftw_complex занимает в памяти непрерывную область и располагается в ней в порядке, характерном для языка Си – быстрее всего меняется последний индекс, медленнее всего – первый. Фактически обрабатываемые данные структурируются в виде трехмерного массива типа double in[N0][N1][2], где in[i][j][0] это реальная часть комплексного числа в i-ой строке и j-ом столбце двумерной матрицы размера N0×N1, а in[i][j][1] – мнимая часть.
Если декларировать массивы на этапе компиляции программы на языке Си, то они автоматически получаются в нужном формате, например,
{ fftw_complex data[N0][N1]; fftw_plan plan; plan = fftw_plan_dft_2d(N0, N1, &data[0][0], &data[0][0], FFTW_FORWARD,
FFTW_ESTIMATE);
}
В этом примере создается план двумерного комплексного преобразования массива размера N0×N1 с размещением выходных данных поверх (в том же массиве data) входных. Однако такое статическое выделение памяти является не очень удачным. Во-первых, могут возникнуть проблемы с переполнение стека для локальных автоматических массивов. А во-вторых, элементы массива могут быть не оптимальным образом выравнены для использования SIMD- инструкций. Вместо этого, память лучше выделять динамически с использованием fftw_malloc, как и в одномерном случае. Например, для комплексного массива размером 8 на 16, это можно сделать так: fftw_complex *in; in = (fftw_complex*) fftw_malloc(8*16*sizeof(fftw_complex));

Глава 6. Компьютерные библиотеки ДПФ
81 при этом обращение к (i, j)-му элементу исходного двумерного массива происходит как к элементу одномерного массива in[i*16+j].

§6.2. ДПФ в пакете Wolfram Mathematica
Дискретное преобразование Фурье можно найти во многих пакетах, предназначенных для научных и технических расчетов.
Рассмотрим в качестве примера пакет Wolfram Mathematica®, который используется на физическом факультете МГУ имени М.В. Ломоносова в рамках дисциплины «Компьютерная физика». Для выполнения дискретного преобразования Фурье в пакете имеются функции
Fourier[List] (прямое преобразование) и InverseFourier[List] – обратное преобразование. Аргументом обеих функций является список (List) комплексных чисел, над которыми происходит соответствующее преобразование. По умолчанию прямым преобразованием в пакете
Математика называют преобразование вида:

−
=
=
1 0
)
(
1
)
(
N
j
kj
N
W
j
u
N
k
U
,
(6.3) и обратным:

−
=
−
=
1 0
)
(
1
)
(
N
k
kj
N
W
k
U
N
j
u
,
(6.4) т.е. со знаком «+» в показателе экспоненты прямого преобразования, знаком
«-» при обратном преобразовании и симметричным коэффициентом
N
/
1
перед суммой в обоих преобразованиях. С помощью опции FourierParameters -> {a, b} можно изменить поведение функций и привести преобразования к более удобному для конкретного применения виду. В общем случае, с использованием этой опции прямое и обратное преобразования имеют, соответственно, вид:

−
=
π
+
+
=
1 0
2 2
/
)
1
(
)
(
1
)
(
N
j
N
ibjk
a
e
j
u
N
k
U
,
(6.5)

−
=
π
−
−
=
1 0
2 2
/
)
1
(
)
(
1
)
(
N
k
N
ibjk
a
e
k
U
N
j
u
(6.6)

Глава 6. Компьютерные библиотеки ДПФ
82
В частности, значение b = –1 меняет знак показателя экспоненты на противоположный и совместно со значением a = 1 приводит ДПФ к виду (5.1), принятому в настоящем пособии
6
В качестве примера рассмотрим
ДПФ уже ранее использованного ряда из N = 4 чисел: u(0) = 1, u(1) = 2, u(2) = 3, u(3) = 4:
Результатом прямого ДПФ является список комплексных амплитуд дискретных гармоник, который в данном примере сохраняется в переменной S1. Первый элемент этого списка – это амплитуда нулевой гармоники U(0), затем последовательно через запятую располагаются амплитуды первой - U(1), второй – U(2) и третьей (в терминах отрицательных частот – минус первой) – U(3) гармоник. Далее в программе вычисляется спектральная плотность мощности анализируемого ряда чисел и, в конце, выполняется обратное дискретное преобразование Фурье, которое, будучи применённое к результату прямого преобразования, даёт исходный ряд чисел.
§6.3. ДПФ на языке Python
7
В настоящее время одним из популярных языков программирования является Питон (Python) - интерпретируемый мультипарадигмальный язык, с помощью которого можно решать
6
В пакете Математика такой выбор параметров связывают с типовым применением ДПФ в задачах обработки сигналов.
7
Авторы выражают искреннюю признательность Е.В. Васильеву за предоставление материалов для этого параграфа.

Глава 6. Компьютерные библиотеки ДПФ
83 широкий круг задач. Если у вас установлен дистрибутив
Python(x,y), специально предназначенный для численных расчётов, анализа и визуализации данных, то ничего дополнительного не потребуется. Если вариант дистрибутива другой, то надо проверить наличие, и при необходимости установить пакеты NumPy
(Numeric Python) и SciPy
(Scientific
Python).
Эти пакеты предоставляют множество математических операций и алгоритмов на компилируемых языках, благодаря чему обеспечивается высокая скорость вычислений.
Установка пакета
NumPy в виртуальное окружение
Python осуществляется командой: pip install numpy.
1. ДПФ из модуля NumPy
Прямое и обратное дискретные преобразования Фурье из пакета NumPy соответствуют формулам (5.1). Рассмотрим опять в качестве входных данных одномерный массив из четырех чисел arr_in = [1, 2, 3, 4]:

Глава 6. Компьютерные библиотеки ДПФ
84
Здесь исходный массив arr_in создается с помощью команды array, в которую в явном виде передается тип complex64. Затем функция fft, имя которой переопределено на fft_numpy, выполняет прямое одномерное комплексное ДПФ, результат которого сохраняется в forward_fft_numpy_res.
Спектральная интенсивность гармоник вычисляется как сумма квадратов действительной и мнимой частей полученного преобразования и сохраняется в массиве spectral_density типа float64. Обратное ДПФ осуществляется функцией ifft
(с псевдонимом ifft_numpy).
Она возвращает массив forward_fft_numpy_res типа complex64.
Аналогично выглядит использование ДПФ из библиотеки
SciPy. Разница заключается в том, что модуль, который реализует
Фурье-преобразование, необходимо импортировать из scipy.fftpack.
2. Использование библиотеки FFTW
Рассмотренная выше библиотека
FFTW доступна для использования в Питоне через интерфейс («обертку»), которая называется PyFFTW
Устанавливается эта библиотека командой: pip install pyfftw.
Пользователям предлагается несколько вариантов интерфейсов.
Они отличаются степенью необходимой детализации при организации выполнения ДПФ. Минимальную детализацию (более простой вариант
– близкий по смыслу к numpy.fft/scipy.fftpack.fft) поддерживают модули pyfftw.interfaces.numpy_fft и
pyfftw.interfaces.scipy_fftpack.fft
Выполнение ДПФ с функциями из этих модулей осуществляется аналогично п. 1:
Увеличить эффективность
ДПФ возможно при предварительном создании планов преобразований, доступных в модуле pyfftw.FFTW
Некоторый компромисс между детальной настройкой и простым numpy-подобным интерфейсом достигается в модуле pyfftw.builders.

Глава 6. Компьютерные библиотеки ДПФ
85 3. Сравнение быстродействия ДПФ в Python
Заметим, что различные модули и библиотеки не только предлагают разный пользовательский интерфейс, но и характеризуются различной эффективностью выполнения преобразований. Общая закономерность проявляется в том, что чем более простой интерфейс и чем меньше предварительных настроек необходимо производить, тем больше время выполнения преобразования. Сравним быстродействие
ДПФ из следующих библиотек:
1) Библиотека NumPy.
2) Библиотека SciPy.
3) Библиотека PyFFTW:
- модуль с интерфейсом numpy_fft;
- модуль с интерфейсом scipy_fftpack;
- модуль pyfftw.builders;
- модуль pyfftw.FFTW с созданием планов Фурье-преобразований.

Глава 6. Компьютерные библиотеки ДПФ
86
На рис. 6.1 приведена зависимость среднего времени последовательного выполнения прямого и обратного ДПФ для каждого из рассмотренных случаев, нормированное на время выполнения ДПФ в библиотеке NumPy, для одномерного (рис. 1а) и двумерного (рис. 1б) массивов.
(а) (б)
Рис. 6.1. Зависимость среднего времени последовательного выполнения прямого и обратного ДПФ для разных модулей и библиотек в Python, нормированное на время выполнения ДПФ в библиотеке NumPy, для одномерного массива типа complex64 размером 1024 точки (усреднение по 50000 реализациям) (а) и двумерного массива типа complex64 размером 1024 x 1024 (усреднение по 500 реализациям) (б).
Видно, что в случае одномерного массива сильно проигрывают в скорости модули библиотеки PyFFTW, имитирующие интерфейсы
NumPy и SciPy, а также модуль pyfftw.builders. Это связано с тем, что для одномерного массива время предварительной внутренней настройки параметров, которая всегда происходит в PyFFTW, сопоставимо со временем всего преобразования. В случае двумерного массива относительный вклад времени настройки (построения плана) мал, поэтому выполнение преобразований в указанных модулях происходит быстрее, чем в тех, интерфейс которых они повторяют. Заметим, что
ДПФ из SciPy немного более эффективен, чем из NumPy, а скорость выполнения преобразований в PyFFTW с предварительным созданием планов, является наибольшей
8 8
Построение планов не входит в измеряемое время.

Литература
87
Литература
1.
А.А. Харкевич. Спектры и анализ, 4-е изд. М.: ЛКИ, 2007. 236 с.
2.
Е. Скучик. Основы акустики, т.1, М.: Мир, 1976, 520 с.
3.
Ж. Макс. Методы и техника обработки сигналов при физических измерениях, т.1, М.: Мир, 1983, 312 с.
4.
Д. Рабинер, Б. Гоулд. Теория и применение цифровой обработки сигналов, М.: Мир, 1978, 848 с.
5.
Р. Отнес, Л. Эноксон. Прикладной анализ временных рядов, М.:
Мир, 1982, 428 с.
6.
Р. Брейсуэлл. Преобразование Фурье //В мире науки. Scientific
American. №8, с.48, 1989.

Учебно-методическое издание
КАНДИДОВ Валерий Петрович
ЧЕСНОКОВ Сергей Сергеевич
ШЛЕНОВ Святослав Александрович
ДИСКРЕТНОЕ
ПРЕОБРАЗОВАНИЕ ФУРЬЕ
Оригинал-макет и компьютерные иллюстрации: Чесноков С.С.
Подписано в печать 9.12 2019 г. Формат 60х90 1/16
Бумага офсетная. Печать офсетная. Объем 5,5 п.л.
Тираж 100 экз. Заказ № 192.
Физический факультет МГУ
119991, ГСП-1, Москва, Ленинские горы, МГУ им. М.В. Ломоносова
Напечатано с готового оригинал-макета
Отдел оперативной печати физического факультета
МГУ имени М.В. Ломоносова
119991, ГСП-1, Москва
Ленинские горы, д. 1, стр. 2
Тираж 100 экз. Заказ № 192