Main menu

Алгоритмы вычисления Быстрых Преобразований Фурье (БПФ) в многомерных пространствах

Алгоритм, изменивший цифровую эпоху

Преобразование Фурье — это математическая основа обработки любых сигналов. Оно позволяет раскладывать сложные волновые формы на набор простых синусоид. Дискретное преобразование Фурье (ДПФ) выполняет эту задачу для цифровых (сэмплированных) данных. Однако прямое, «школьное» вычисление ДПФ по формуле имеет квадратичную алгоритмическую сложность O(N^2). Для аудиофайла длиной в одну секунду (44 100 отсчетов) потребуется около двух миллиардов операций умножения. Если бы мы пользовались этим прямым методом, современная потоковая передача видео, мобильная связь 4G/5G и магнитно-резонансная томография (МРТ) были бы просто невозможны из-за вычислительных ограничений.

Проблема была решена в 1965 году, когда Джеймс Кули и Джон Тьюки опубликовали алгоритм Быстрого Преобразования Фурье (БПФ, или FFT). (Справедливости ради, похожие идеи высказывал еще Карл Фридрих Гаусс в 1805 году, но они опередили свое время). Алгоритм БПФ радикально снижает количество операций с O(N^2) до O(N log N). Для массива в миллион точек ускорение составляет почти 50 000 раз! Журнал Computing in Science & Engineering заслуженно включил БПФ в десятку величайших алгоритмов XX века.

Разделяй и властвуй: принцип бабочки

Самая популярная версия алгоритма Кули-Тьюки базируется на принципе «разделяй и властвуй» (divide and conquer) по основанию 2 (Radix-2). Алгоритм требует, чтобы количество точек N в массиве было степенью двойки (например, 256, 512, 1024). Если это не так, массив дополняется нулями (zero-padding).

Алгоритм рекурсивно разбивает задачу вычисления одного ДПФ размера N на вычисление двух ДПФ размера N/2 (одно для четных индексов массива, другое — для нечетных). Затем каждый из этих подмассивов снова делится пополам, и так далее, пока размер массива не достигнет 1. Затем результаты начинают объединяться (синтезироваться) обратно. На этапе объединения используется базовая вычислительная структура, которая из-за графа информационных потоков получила название «операция бабочки» (Butterfly operation). Она комбинирует два комплексных числа, используя умножение на специальные комплексные экспоненты (поворачивающие множители). Перед выполнением БПФ элементы массива переставляются в специальном порядке, называемом битовой инверсией (Bit-reversal permutation), что позволяет выполнять все вычисления прямо на месте (in-place), не выделяя дополнительную оперативную память.

Многомерное БПФ: от звука к изображениям

Обычное одномерное БПФ применяется для обработки звука или одномерных радиосигналов. Но как анализировать цифровые фотографии (поиск краев, сжатие JPEG, фильтрация шума) или решать трехмерные дифференциальные уравнения Пуассона на сетке? Для этого используется многомерное БПФ.

Счастье для инженеров заключается в том, что многомерное дискретное преобразование Фурье математически сепарабельно (разделимо). Это означает, что для вычисления двумерного БПФ изображения (матрицы пикселей) не нужно писать новый сложный алгоритм. Достаточно сначала применить обычное одномерное БПФ к каждой строке матрицы независимо. Получив промежуточную матрицу комплексных чисел, мы затем применяем то же самое одномерное БПФ к каждому ее столбцу. Этот строчный-столбцовый подход (Row-Column algorithm) легко масштабируется на 3D и 4D пространства, идеально распараллеливается на тысячи ядер графических процессоров (с помощью библиотек вроде cuFFT от NVIDIA) и является бьющимся сердцем современной компьютерной томографии и молекулярной химии.

Подробнее

Вычисление определителей и обращение матриц: численные аспекты

Разрыв между теоретической алгеброй и вычислительной реальностью

В курсе высшей алгебры студенты учатся решать системы линейных алгебраических уравнений (СЛАУ) вида Ax = b с помощью правила Крамера (через определители) или путем явного вычисления обратной матрицы (x = A^-1 * b). Эти методы красивы в теории и отлично подходят для ручного решения простейших систем 3x3 на листке бумаги. Однако в реальной вычислительной математике при написании программного обеспечения инженеры избегают этих методов как огня. Почему же аналитические фавориты терпят крах при встрече с кремниевыми процессорами?

Все дело в вычислительной сложности. Вычисление определителя по классической формуле разложения по строке (формула Лейбница) требует факториального количества операций (O(N!)). Для матрицы размером всего 20x20 вычисление определителя таким способом потребует 20! операций. Если суперкомпьютер будет выполнять миллиард операций в секунду, ему потребуется несколько сотен тысяч лет, чтобы решить эту крошечную задачу. Использование правила Крамера для СЛАУ требует вычисления N+1 таких определителей, что делает этот метод абсолютно бесполезным для практических расчетов в физике и инженерии.

Почему обращать матрицы явно — плохая идея?

Вторая распространенная ошибка новичков в программировании — явное вычисление обратной матрицы A^-1 для решения системы уравнений. Во-первых, вычисление обратной матрицы методом Гаусса-Жордана требует в 3 раза больше арифметических операций, чем решение самой системы методом исключения Гаусса (или с помощью LU-разложения). Мы тратим драгоценное время процессора впустую.

Во-вторых, обращение разреженных матриц приводит к катастрофическому заполнению памяти. Как мы помним из предыдущих статей, в методе конечных элементов матрица A на 99% состоит из нулей. Однако ее обратная матрица A^-1 почти всегда является полностью плотной (все элементы не равны нулю). Если исходная матрица из миллиона уравнений в формате CSR занимала несколько десятков мегабайт, то обратная матрица займет терабайты оперативной памяти. Именно поэтому золотое правило вычислительной линейной алгебры гласит: «Никогда не вычисляйте обратную матрицу явно, если вам нужно только решить систему уравнений». Вычисляйте LU-разложение или используйте итерационные методы.

Обусловленность и устойчивость обратной матрицы

Еще одной колоссальной проблемой является накопление ошибок округления. Каждая математическая матрица характеризуется числом обусловленности (Condition Number) — грубо говоря, это мера того, насколько матрица близка к вырожденной (к матрице с нулевым определителем).

Если число обусловленности велико (матрица плохо обусловлена), то микроскопические ошибки в исходных данных или ошибки округления с плавающей запятой внутри процессора при вычислении обратной матрицы будут многократно усилены (иногда в миллионы раз). Полученная компьютерная «обратная матрица» при умножении на исходную даст результат, крайне далекий от единичной матрицы. Анализ собственных и сингулярных чисел (SVD), а также использование алгоритмов предобуславливания (Preconditioning) стали обязательными инструментами для стабилизации процесса решения СЛАУ и борьбы с физической некорректностью инженерных моделей.

Подробнее

Численное моделирование в электродинамике: метод FDTD (сетка Йи)

Динамика электромагнитных полей в реальном времени

Система уравнений Максвелла является фундаментом всей современной классической электродинамики, описывая генерацию и распространение электромагнитных волн. Для расчета сложных антенн смартфонов, микроволновых печей, оптоволоконных кабелей и радарных систем-невидимок инженерам необходимо решать эти уравнения с высочайшей точностью. Одним из самых мощных и интуитивно понятных численных алгоритмов для этой задачи является Метод конечных разностей во временной области (Finite-Difference Time-Domain, FDTD).

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

Магия смещенной сетки Кейна Йи (Yee Grid)

Гениальность метода FDTD, предложенного Кейном Йи в 1966 году, заключается в особой топологии пространственной дискретизации. Если в обычном методе сеток все неизвестные величины (например, температура и давление) вычисляются в одних и тех же узлах, то в сетке Йи компоненты электрического (E) и магнитного (H) полей разнесены в пространстве.

Представьте себе кубическую ячейку. Векторы электрического поля располагаются строго на серединах ребер этого куба, а векторы магнитного поля — строго в центрах граней. Такое шахматное расположение идеально соответствует физической сути роторных уравнений Максвелла (закону Фарадея и закону Ампера-Максвелла): циркуляция электрического поля по замкнутому контуру равна изменению магнитного потока сквозь этот контур. Кроме того, вычисления разнесены не только в пространстве, но и во времени (метод «чехарды», leapfrog). Электрическое поле обновляется на целых шагах по времени, а магнитное — на полуцелых. Это обеспечивает высокую, вторую степень точности аппроксимации как по времени, так и по пространству, а также гарантирует безусловную физическую консервативность алгоритма (отсутствие фиктивных электрических зарядов).

Проблема открытого пространства: поглощающие условия (PML)

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

Для решения этой фундаментальной проблемы Жан-Пьер Беренжер в 1994 году изобрел Идеально Согласованные Слои (Perfectly Matched Layer, PML). Это специальный математический фиктивный материал, который располагается по краям расчетной сетки. Волна любой поляризации, падающая на PML под любым углом, проникает в него без отражения на границе раздела, а внутри слоя ее амплитуда экспоненциально затухает до нуля. Изобретение PML сделало метод FDTD индустриальным стандартом в вычислительной электродинамике, оптике и фотонике.

Подробнее

ГДЗ по математике

1 класс 2 класс 3 класс
4 класс 5 класс 6 класс

Книги

Фильмы

Старый сайт

Соц. сети