Преобразование Фурье самый подробный разбор


Оглавление (нажмите, чтобы открыть):

Реконструкция изображений \ Image Processing Toolbox

Реконструкция изображений

Преобразование Фурье

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

Далее рассмотрим следующие вопросы:

  1. Основные определения преобразования Фурье
  2. Дискретное преобразование Фурье, включая быстрое преобразование Фурье
  3. Применение Фурье-преобразования (некоторые примеры практического применения преобразования Фурье)

Основные определения преобразования Фурье

Если ƒ(m,n) представляет собой функцию двух дискретных пространственных переменных m и n, тогда двумерное преобразование Фурье функции ƒ(m,n) может быть представлено следующим выражением

Переменные представляют собой угловые частоты. Таким образом, представляет собой функцию ƒ(m,n) в частотной области. является комплекснозначной функцией с соответствующими частотами . Частоты находятся в пределах диапазона , . Отметим, что F (0,0) представляется в виде суммы всех переменных ƒ(m,n). По этой причине F (0,0) часто называют постоянной составляющей преобразования Фурье.

Обратное двумерное преобразование Фурье представляется выражением

Т.е. это выражение представляет ƒ(m,n) в виде суммы бесконечного числа сложных экспоненциальных функций (синусоид) с различными частотами. Амплитуда и фаза определяют вклад частот в представление .

При иллюстрации Фурье-преобразования допустим, что функция ƒ(m,n) равна 1 и представлена в виде прямоугольника. Для упрощения диаграммы, функция ƒ(m,n) будет представляться непрерывной функцией двух дискретных переменных m и n.

На рисунке внизу с использованием функции mesh визуализировано значения амплитуд, которые получены при Фурье-преобразовании прямоугольной функции, представленной на предыдущем рисунке. Визуализацию амплитуды еще называют визуализацией преобразований Фурье.

Амплитуда изображения прямоугольной функции

Пик функции находится в центре и отображает значение F (0,0), которое является суммой всех значений ƒ(m,n). Все остальные составляющие представляют собой распределение энергии по вертикальным и горизонтальным частотам.

Другой путь визуализации Фурье-преобразования заключается в отображении значений в виде изображения.

Логарифмическое представление Фурье-преобразования прямоугольной функции

Рассмотрим примеры преобразования Фурье функций различных простых форм.

Примеры Фурье преобразования функций различных простых форм

Дискретное косинусное преобразование

Дискретные косинусные преобразования представляют изображение в виде суммы синусоид с различной амплитудой и частотой. Функция dct2 в приложении Image Processing Toolbox реализует двумерные дискретные косинусные преобразования изображений. Одна из особенностей дискретного преобразования Фурье состоит в том, что некоторые локальные участки изображения можно охарактеризовать небольшим количеством коэффициентов дискретного преобразования Фурье. Это свойство очень часто используется при разработке методов сжатия изображений. Например, дискретное косинусное преобразование является основой международного стандарта, который используется в алгоритме сжатия изображений с потерями JPEG. Название формата “JPEG” состоит из первых букв названия рабочей группы, которая принимала участие в разработке этого стандарта (Joint Photographic Experts Group).

Двумерное дискретное косинусное преобразование матрицы A с размерами реализуется согласно следующему выражению

Значения Bpq называют коэффициентами дискретного косинусного преобразования матрицы A.

(Следует отметить, что индексы матрицы в MATLAB всегда начинаются с 1, а не с 0. Поэтому элементы матрицы, которые представлены в MATLAB как A(1,1) и B(1,1), будут соответствовать элементам A00 и B00 из приведенной выше формулы.)

Обратное дискретное косинусное преобразование реализуется согласно выражениям

Выражение обратного дискретного косинусного преобразования может интерпретироваться как представление матрицы A с размерами в виде суммы следующих функций

Эти функции называются основными (базовыми) функциями дискретного косинусного преобразования. Коэффициенты дискретного косинусного преобразования Bpq можно рассматривать как весовые при каждой базовой функции. Например, для матрицы с размером элементов существует 64 базовые функции, что продемонстрировано на изображении.

64 базовые функции, которые получены для матрицы с размерами элементов

Горизонтальные частоты увеличиваются слева направо, а вертикальные – сверху вниз.

Матрица дискретных косинусных преобразований

Приложение Image Processing Toolbox предлагает два различных пути реализации дискретных косинусных преобразований. Первый метод реализован в функции dct2. Функция dct2 использует быстрое преобразования Фурье для ускорения вычислений. Второй метод использует матрицу дискретных косинусных преобразований, которая возвращается функцией dctmtx. Матрица преобразований T формируется согласно следующего выражения

Для матрицы A с размерами представляет собой матрицу с размерами , где каждый столбец содержит одномерное дискретное косинусное преобразование A. Двумерное дискретное косинусное преобразование A вычисляется как B=T*A*T’. Обратное двумерное дискретное косинусное преобразование B вычисляется как T’*B*T.

Дискретные косинусные преобразования и сжатие изображений

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

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

На рисунке представлено два изображения – исходное и реконструированное. При реконструкции изображения использовалось только 15 % коэффициентов дискретных косинусных преобразований. Однако, следует отметить, что качество реконструированного изображения является довольно приемлемым. Для просмотра других свойств дискретного косинусного преобразования см. функцию dctdemo.

Преобразования Радона

Функция radon в приложении Image Processing Toolbox вычисляет матрицу проекций изображения вдоль заданных направлений. Проекция двумерной функции f(x,y) равна интегралу вдоль указанной линии. Функция Радона представляет собой вычисление проекций изображения на оси, которые задаются углами в градусах относительно горизонтали против часовой стрелки. На рисунке показана проекция некоторой фигуры под указанным углом

Параллельно-лучевая проекция с углом поворота theta

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

Горизонтальная и вертикальная проекции некоторой простой функции

Проекции могут вычисляться вдоль произвольного угла theta. Встроенная в приложение Image Processing Toolbox функция radon вычисляет проекции изображения вдоль определенных направлений. Проекция двумерной функции f(x,y) на ось x’ представляет собой линейный интеграл

Таким образом, оси x’ y’ задаются поворотом на угол против часовой стрелки.

На изображении внизу проиллюстрировано геометрию преобразования Радона.

Геометрия преобразования Радона

Визуализация преобразований Радона

При проведении преобразований Радона необходимо указать исходное изображение и вектор углов theta.

R представляет собой матрицу, в которой каждый столбец является преобразованием Радона для одного из углов, который содержится в векторе theta. Вектор xp содержит соответствующие координаты вдоль оси x. Центральный пиксель I определяется согласно выражению floor((size(I)+1)/2).

Рассмотрим, как в преобразованиях Радона вычисляются проекции. Рассмотрим проекции под углом 0° и 45°.

Преобразования Радона при 0°

Преобразования Радона при 45°

Преобразования Радона при большом числе углов часто отображается в виде изображения. В данном примере рассмотрено преобразования Радона для изображения в виде квадрата при диапазоне углов от 0° до 180° с дискретностью 1°.

Преобразования радона с использованием 180 проекций

Использование преобразований Радона при детектировании линий


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

    Вычисление бинарного изображения с использованием функции edge

Исходное изображение Изображение с выделенными границами

Вычисление преобразований Радона на основе бинарного изображения границ

Визуализация применения преобразований Радона к изображению границ

  • Поиск локальных пиков на матрице преобразований Радона. Месторасположение пиков соответствует расположению линий на исходном изображении
  • Наибольший пик в матрице R соответствует =1° и x´= -80. Из центра исходного изображения проводится линия под углом на расстояние x’. Перпендикулярно к этой линии проводится прямая, которая соответствует прямой на исходном изображении. Кроме того, на изображении присутствуют и другие линии, которые представлены в матрице R соответствующими пиками.

    Геометрия преобразования Радона при детектировании прямых линий

    Свойства преобразования

    DSPL-2.0 — свободная библиотека алгоритмов цифровой обработки сигналов

    Распространяется под лицензией LGPL v3

    Страница проекта на GitHub.

    Содержание

    Свойство линейности

    Пусть даны непериодические сигналы и , а также их спектральные плотности и соответственно. Везде далее мы будем предполагать, что и — абсолютно интегрируемые сигналы, тогда преобразование Фурье сигнала равно

    Следствием является свойство умножения на константу :

    Свойство временного сдвига

    Рассмотрим сигнал как результат временного сдвига исходного сигнала на произвольную величину . Тогда преобразование Фурье сигнала имеет вид:

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

    Преобразование Фурье свертки сигналов

    Пусть сигнал представляет собой свертку сигналов и :

    Поменяем порядок интегрирования, и используем свойство (4) временного сдвига:

    Таким образом, спектральная плотность свертки двух сигналов равна произведению их спектральных плотностей.

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

    Преобразование Фурье произведения сигналов

    Пусть сигнал представляет собой произведение сигналов и . Преобразование Фурье сигнала равно:

    Подставим в (8) вместо сигнала обратное преобразование Фурье его спектральной плотности :

    Поменяем в (9) операции интегрирования и получим:

    Масштабирование и инверсия во времени

    Пусть сигнал представляет собой масштабированный во времени сигнал , — вещественная константа, отличная от нуля.

    Тогда преобразование Фурье сигнала равно:

    Преобразование Фурье производной исходного сигнала

    Пусть сигнал представляет собой непрерывный на всей числовой оси абсолютно интегрируемый сигнал, чья спектральная плотность равна . Тогда сигнал также является абсолютно интегрируемым, и его преобразование Фурье равно:

    Используем правило интегрирования по частям [1] [2, стр. 330]:

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

    Как и в случае с периодическими сигналами, наличие множителя приводит к тому, что с ростом частоты затухает слабее чем спектральная плотность исходного сигнала . Поэтому изначально мы наложили ограничение на исходный сигнал: он должен быть непрерывным, тогда его спектральная плотность будет затухать быстрее чем , и умножение на не приведет к росту с увеличением частоты, т.е. обеспечит сходимость (17).

    Свойство интегрирования исходного сигнала

    Пусть теперь представляет собой сигнал с нулевой постоянной составляющей. Спектральная плотность сигнала равна нулю при .

    Обратим внимание, что при , сигнал является абсолютно интегрируемым.

    Рассмотрим спектральную плотность сигнала . Для этого заметим, что сигнал ничто иное, как производная сигнала . Тогда используя свойство преобразования Фурье производной сигнала (19) можно записать:

    При , спектральная плотность рассчитывается без особого труда. Однако на частоте получаем неопределенность вида , раскрытие которой по правилу Лопиталя [2, стр. 257] по аналогии со свойством ряда Фурье приводит к окончательному выражению вида:

    Анализируя (22) можно заключить, что интегрирование сигнала устраняет разрывы и приводит к более быстрому затуханию спектральной плотности, ввиду наличия дополнительного множителя .

    Преобразование Фурье комплексно-сопряженного сигнала

    Пусть исходный сигнал представляет собой комплексный абсолютно интегрируемый сигнал, где — синфазная и — квадратурная компоненты. Рассмотрим преобразование Фурье , используя формулу Эйлера представления комплексных экспонент:

    Рассмотрим теперь преобразование Фурье комплексно-сопряженного сигнала :

    Таким образом, спектральная плотность комплексно-сопряженного сигнала равна инверсной по частоте комплексно-сопряженной спектральной плотности исходного сигнала.

    Важным следствием (27) является свойство симметрии спектральной плотности вещественного сигнала.

    Пусть имеется вещественный сигнал , чья спектральная плотность равна . Поскольку сигнал вещественный, то комплексное сопряжение его не меняет, т.е. . Перейдя в частотную область, с учетом (27) получаем равенство:

    Для вещественного сигнала выражение (24) с учетом можно представить:

    Тогда определяет реальную часть спектральной плотности и является четной функцией частоты, а — мнимая часть спектральной плотности является нечетной функцией частоты.


    Амплитудно- и фазочастотная характеристики сигнала

    По аналогии с понятиями амплитудного и фазового спектра периодического сигнала можно ввести понятия амплитудно-частотной (АЧХ) и фазочастотной (ФЧХ) характеристики сигнала:

    Тогда, в случае вещественного сигнала, с учетом свойств симметрии спектральной плотности (28), можно заключить, что АЧХ является четной функцией частоты , а ФЧХ — нечетной: .

    Свойство симметрии наглядно показано на рисунке 1 для одностороннего экспоненциального импульса.

    Двойственность преобразования Фурье

    Пусть сигнал имеет спектральную плотность . Рассмотрим что произойдет, если мы возьмем преобразование Фурье от спектральной плотности :

    Обратите внимание, интегрирование идет по переменной , хотя выражение (32) представляет собой прямое преобразование Фурье. Тогда можно преобразовать:

    Можно сделать вывод, что преобразование Фурье от спектральной плотности снова возвращает сигнал, инверсный во времени, умноженный на . Это свойство носит название двойственности (дуальности) преобразования Фурье.

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

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

    Вводя замену переменной получаем , . Пределы интегрирования остаются неизменными, и выражение (34) принимает вид:

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

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

    Мы говорили, что спектральная плотность вещественного сигнала является симметричной относительно нулевой частоты: . Тогда мы можем переформулировать это свойство и в другую сторону: если сигнал (в общем случае комплексный) обладает свойством симметрии во временной области: , то его спектральная плотность чисто вещественна.

    Для доказательства данного утверждения представим сигнал в виде . Тогда . Если выполняется условие симметрии , то:

    Тогда с учетом выражения (23) и (24), мнимая часть спектральной плотности равна:

    Заметим, что (38) представляет собой интегралы в бесконечных симметричных пределах от произведения четной и нечетной функции (четность и нечетность мы установили выше). Тогда можно заключить, что оба интеграла (38) равны нулю, и мнимая часть спектральной плотности , симметричного во времени сигнала , также равна нулю согласно (37). Что и требовалось доказать.

    Убывание спектральной плотности сигнала по частоте

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

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

    Прежде всего, нам потребуется обратиться к лемме Римана-Лебега.

    Лемма (Римана-Лебега). Преобразование Фурье абсолютно-интегрируемой функции является ограниченной функцией частоты , при этом стремится к нулю при .

    Доказательство данной леммы приведено в [3, стр. 83–84]. Лемма Римана-Лебега доказывает качественное свойство убывания спектральной плотности абсолютно-интегрируемого сигнала с ростом частоты , но не дает количественной оценки скорости убывания .

    Пусть исходный сигнал является непрерывной, абсолютно-интегрируемой функцией времени, которая может быть дифференцируема раз, причем все первых производных также будут представлять собой абсолютно-интегрируемые функции. Тогда производную порядка сигнала можно обозначить как . Если все первых производных являются абсолютно-интегрируемыми, то мы можем перейти в частотную область и согласно свойству преобразования Фурье, спектральная плотность равна:

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

    Наличие в сигнале разрыва первого рода (например скачка в прямоугольном импульсе) приводит к убыванию спектральной плотности со скоростью с ростом частоты.

    Выводы

    В данном разделе мы рассмотрели некоторые свойства преобразования Фурье.

    В следующем разделе мы рассмотрим спектральные плотности некоторых распространенных сигналов.

    Смотри также

    Примечания

    [1] Преобразование Фурье представляет собой интеграл в бесконечных пределах. Применение правила интегрирования по частям возможно [1, стр. 374] при соблюдении условия сходимости выражения (17), которое обеспечивается при непрерывном .

    Физический смысл БПФ

    Для чего нужно быстрое преобразование Фурье или вообще дискретное преобразование Фурье (ДПФ)? Давайте попробуем разобраться.

    Пусть у нас есть функция синуса .

    Максимальная амплитуда этого колебания равна . Если умножить его на некоторый коэффициент , то получим тот же график, растянутый по вертикали в раз: .

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

    Частота колебания обратна периоду: . Также говорят о круговой частоте, которая вычисляется по формуле: . Откуда: .

    И, наконец, есть фаза, обозначаемая как φ. Она определяет сдвиг графика колебания влево. В результате сочетания всех этих параметров получается гармоническое колебание или просто гармоника:

    Очень похоже выглядит и выражение гармоники через косинус:

    Большой разницы нет. Достаточно изменить фазу на , чтобы перейти от синуса к косинусу и обратно. Далее будем подразумевать под гармоникой функцию косинуса:

    В природе и технике колебания, описываемые подобной функцией чрезвычайно распространены. Например, маятник, струна, водные и звуковые волны и прочее, и прочее.

    Преобразуем (18) по формуле косинуса суммы:

    Выделим в (19) элементы, независимые от , и обозначим их как и :

    По величинам и можно однозначно восстановить амплитуду и фазу исходной гармоники:

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

    В результате прямого дискретного преобразования Фурье были получены значений для :

    Теперь возьмем обратное преобразование Фурье:

    Выполним над этой формулой следующие действия: разложим каждое комплексное на мнимую и действительную составляющие ; разложим экспоненту по формуле Эйлера на синус и косинус действительного аргумента; перемножим; внесем под знак суммы и перегруппируем элементы в две суммы:

    Это была цепочка равенств, которая начиналась с действительного числа . В конце получилось две суммы, одна из которых помножена на мнимую единицу . Сами же суммы состоят из действительных слагаемых. Отсюда следует, что вторая сумма должна быть равна нулю. Отбросим ее и получим:

    Поскольку при дискретизации мы брали и , то можем выполнить замену: . Следовательно, в синусе и косинусе вместо можно написать . В результате получим:

    Сопоставим эту формулу с формулой (20) для гармоники:

    Слагаемые суммы (26) аналогичны формуле (20), а формула (20) описывает гармоническое колебание. Значит сумма (26) представляет собой сумму из гармонических колебаний разной частоты, фазы и амплитуды.

    Выше объяснялось, каким образом формула вида (20) может быть преобразована в формулу вида (18):


    Выполним такое же преобразование для слагаемых суммы (26), преобразуем их из вида (20) в вид (18). Получим:

    Далее будем функцию

    называть k-й гармоникой.

    Для вычисления и надо использовать формулу (21). Теперь выпишем в одном месте все формулы, которые связывают амплитуду, фазу, частоту и период каждой из гармоник с коэффициентами :

    Функция — это аргумент комплексного числа. В языке C++ для ее вычиления удобно использовать функцию .

    Итак. Физический смысл дискретного преобразования Фурье состоит в том, чтобы представить некоторый дискретный сигнал в виде суммы гармоник. Параметры каждой гармоники вычисляются прямым преобразованием, а сумма гармоник — обратным.

    Преобразование Фурье: самый подробный разбор

    Задачи распознавания сцен и изображений традиционно входят в число классических проблем искусственного интеллекта [3]. Причем, если решение задачи распознавания изображений обычно не вызывает серьезных трудностей [1, 2], то в распознавании сцен успехов достигнуто значительно меньше [3]. Способность определять, имеется ли на изображении тот или иной объект, важна для практики, в частности для систем видеонаблюдения. Обнаружив, что на изображении имеются люди, можно приступить к детекции и распознаванию лиц, а в случае наличия автомобилей, перейти на детекцию и распознавание автомобильных номерных знаков. Таким образом, нет необходимости постоянно отрабатывать алгоритмы для обработки объектов других видов. Когда наблюдение за территорией ведется круглосуточно, информация о наличии, либо отсутствии объекта на изображении позволит сократить время обработки изображений и избавит от лишней нагрузки на процессор, что, в свою очередь, позволит быстрее получить нужные данные. Например, можно вести учет въехавших и выехавших автомобилей на охраняемую территорию, учет отгруженных товаров и даже определять – где и в какое время находился определенный сотрудник.

    Проблема заключается в том, что не существует универсальных методов определения наличия объекта на изображении. Можно лишь использовать специальные алгоритмы для поиска конкретных объектов, но и в данном случае выбор алгоритмов слишком велик и среди них нет универсального. Например, для поиска лиц существует много алгоритмов, такие как: метод главных компонент (Principal component analysis, PCA) [7], линейный дискриминантный анализ (Linear Discriminant Analysis, LDA) [7], скрытые Марковские модели [7], гибкое сопоставление (Elastic matching, EM) [7], метод опорных векторов [7] и активная модель внешнего вида (Active appearance model, AAM) [7]. И это не полный список методов поиска лиц на изображении. Более обобщенным методом поиска объектов можно считать подход «оптический поток» (Оptical flow) [8], но он годится только для движущихся объектов – в общем случае определяют направление движения объекта и уже работают с примерной формой обнаруженного объекта.

    Преобразование Фурье было выбрано вместо вейвлет анализа, потому что дискретное вейвлет преобразование не сохраняет постоянство при сдвиге [6, с. 62]. У преобразования Фурье, согласно свойству сдвига во времени, спектр объекта не меняется, а меняется только фаза, поэтому спектр всегда получается один, где бы на изображении объект ни находился.

    Из готовых реализаций преобразования Фурье стоит упомянуть функцию fft в пакете прикладных программ MATLAB [9] и библиотеку FFTW, написанную на С [10]. Все существующие решения не пригодны для полного исследования преобразования Фурье, и их интеграция с нашими проектами стоила бы неоправданное количество времени и сил.

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

    Выбор типа преобразования Фурье

    Из анализа литературы по вопросу преобразования Фурье следует, что многие исследователи не имеют единого мнения о предмете. Если рассматривать разных авторов, то одни и те же обозначения в формулах означают разные параметры. Из приведенных ниже формул (1) [4, с. 6], (2) [5, с. 14] и (3) [6, с. 35, с. 38] будем использовать формулу (3), чтобы избежать путаницы в обозначениях.

    где – сигнал в области времени, – преобразование Фурье, – время, – частота, – мнимая единица.

    где функции и составляют пару преобразования Фурье, – преобразование Фурье, – сигнал в области времени, – частота, – время, – мнимая единица.

    где – преобразование Фурье, – угловая частота, измеряется в рад/с (где частота – обратная величина от времени , ), – мнимая единица, непрерывный сигнал (непрерывно изменяющийся с течением времени), .

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

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

    1. Не периодичный непрерывный. К таким сигналам относится, например экспонента и кривая Гаусса. Эти сигналы устремляются в бесконечность в обе стороны и не представляют собой совокупности повторяющихся периодов. Этот тип преобразования называют преобразованием Фурье (Fourier Transform).

    2. Периодичный непрерывный. К таким сигналам можно отнести синусоиды, квадратные волны и любые волны с повторяющимся периодом от минус бесконечности до плюс бесконечности. Этот тип преобразования называют ряды Фурье (Fourier Series).

    3. Не периодичный дискретный. Эти сигналы только определены дискретными точками между плюс бесконечностью и минус бесконечностью и не имеют повторяющегося периода. Этот тип преобразования называют дискретным преобразованием Фурье во времени (Discrete Time Fourier Transform).

    4. Периодичный дискретный. Это дискретные периодические сигналы от минус бесконечности до плюс бесконечности. Этот тип преобразования называют дискретными рядами Фурье (Discrete Fourier Series) или дискретным преобразованием Фурье (Discrete Fourier Transform).

    Заметим, что в каждой категории преобразования Фурье сигналы изменяются от минус бесконечности до плюс бесконечности. На практике приходиться работать с дискретными сигналами конечной длины. Так как волны синуса и косинуса определены от минус бесконечности до плюс бесконечности, то невозможно, используя группу бесконечно длинных сигналов, создать что-то конечной длины. Чтобы обойти эту проблему, нужно вообразить, что наш сигнал имеет бесконечное число точек слева и справа от наших реальных данных. Если все эти воображаемые точки имеют значение ноль, то сигнал выглядит как дискретный не периодический, тогда применяется дискретное преобразование Фурье во времени. Если наши воображаемые точки будут копиями исходных N точек сигнала, то сигнал будет выглядеть дискретным периодическим с периодом в N точек. В этом случае применяется дискретное преобразование Фурье (ДПФ).

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

    Каждая категория преобразования Фурье разбивается на два вида: в вещественных числах и в комплексных числах. Категория в комплексных числах более сложна для восприятия, но именно на ней основано быстрое преобразование Фурье (БПФ), которое повсеместно и применяется для обработки цифровых сигналов.

    Преобразование Фурье для получения спектра изображения

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

    Мастер Йода рекомендует:  Программирование для начинающих – куда сделать первый шаг

    Опишем некоторые свойства преобразования Фурье, которые нам понадобятся для интерпретации спектра изображений.

    1. Линейность. Если складывать или вычитать два сигнала в области времени, то их спектры также складываются или вычитаются. Это означает, что можно раскладывать изображения по спектру. Например, имеется в наличие кровавый отпечаток на ткани. Очень сложно отделить отпечаток от ткани, анализируя целое изображение. Однако если с помощью преобразования Фурье перевести изображение в область частот, то будут получены как сильные компоненты, представляющие текстуру ткани, так и слабые разбросанные компоненты, представляющие отпечаток. Если подавить частотные компоненты ткани и провести обратное преобразование Фурье, то ткань исчезнет с изображения и останется только отпечаток пальца, который можно легко рассмотреть [6, с. 57].

    2. Сдвиг во времени. Это означает, что при смещении сигнала в области времени, его спектр не изменяется, меняется только фаза.

    3. Поворот. Это означает, что при повороте сигнала в области времени, сигнал в области частот поворачивается точно также.

    4. Масштабируемость. Это означает, что при расширении сигнала в области времени, его спектр становится уже и наоборот.

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

    Поясним, как некоторые свойства преобразования Фурье влияют на изображения.

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

    Согласно свойству Фурье, «сдвиг во времени» при смене положения объекта на изображении, в области частот меняется только фаза, таким образом, по амплитуде можно понять, что это именно этот объект, пусть даже находящийся в другом месте на изображении. Если же на изображении имеется множество разных объектов, то их спектры сложатся, и будет невозможно узнать, присутствовал ли на изображении искомый объект. Разумеется, согласно свойству линейности, можно из спектра изображения вычесть спектр искомого объекта (рисунок 1). Тогда, если объект на изображении присутствовал, то он исчезнет с изображения. Если на изображении искомого предмета не было, то результат будет зависеть от того, как много частот занимает вычитаемый объект. А именно, либо будут небольшие помехи на изображении, либо изображение станет неразборчивым (рисунок 2).

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

    Рисунок 2. Вычитание спектра прямоугольника из изображения, где прямоугольник не присутствовал

    Все приведенные выше операции проводились с помощью созданного в процессе исследования комплекса программ.

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

    Рецензенты:

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

    Ясницкий Леонид Нахимович, доктор технических наук, профессор, заведующий кафедрой прикладной информатики, Пермский государственный гуманитарно-педагогический универ-ситет, г. Пермь.

    ДИСКРЕТНОЕ И БЫСТРОЕ ПРЕОБРАЗОВАНИЕ ФУРЬЕ

    Среди прямых можно выделить алгоритмы с непосредственным вычислением выражения (13.5) и быстрые прямые алгоритмы (БПА), которые позволяют для некоторых частных случаев (четных значений, периодический последовательностей и др.) уменьшить требуемые вычислительные затраты [3].

    Большое распространение получили алгоритмы на основе обобщенного дискретного преобразования Фурье, которые объединяют широкую совокупность алгоритмов, использующих: дискретное преобразование Фурье; теоретико-числовое преобразование [3,9]; полиноминальную алгебру [1]; специальные методы; комбинацию перечисленных алгоритмов.

    В цифровой обработке сигналов обобщённое дискретное преобразование Фурье имеет важное значение не только в связи с вычислением цифровых свёрток, но и в ряде других приложений, причём наиболее полно его аппарат разработан для метода дискретного преобразования Фурье (ДПФ). Последний представляет собой специальную версию общего Фурье-преобразования, позволяющего представить сигналы в виде суперпозиции гармонических функций и подробно рассмотренного во многих литературных источниках, например [6,11,13]. ДПФ, описываемое парой соотношений (прямым и обратным):

    где выполняет блочный анализ временного сигнала.

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

    Если временная функция содержит только действительные значения, что является случаем, часто встречающимся на практике, то спектр сигнала является сопряжённой линейчатой функцией, т.е. компоненты положительной и отрицательной областей частот имеют одинаковые амплитуды, но противоположные фазы. Иначе говоря, только N/2 отсчётов в частотной области являются независимыми, а именно те, что расположены на отрезке

    В общем случае при вычислении ДПФ от комплексной входной последовательности необходимо производить вычисление всех N спектральных компонент, так как в данном случае спектр не обладает свойством комплексно-сопряженной симметрии.

    Ещё одна особенность ДПБ связана с представлением анализируемого сигнала как одного из интервалов бесконечной периодической функции. Если начальное и конечное значения сигнала на периоде различны, то в моменты времени 0, T, 2T,… происходит так называемый скачок функции, который может привести к существенному искажению спектра. Для устранения негативного влияния на спектр данной разрывности и получения лучшей избирательности необходимо использовать сглаживание с помощью одного из видов весовой функции (окна) [6,13]. Подробное описание различных окон, наиболее часто применяемых в гармоническом анализе методом ДПФ, представлено, например, в работе [16].

    Возможность использования ДПФ при вычислении ЦС вытекает из фундаментального свойства преобразования Фурье, гласящего, что Фурье-образ произведения двух временных сигналов эквивалентен свёртке преобразований от каждого из сомножителей и, наоборот, свертка сигналов во временной области соответствует умножению их Фурье-отображений в частотной области [4, 5]. Поэтому возможно вычисление циклической (периодической) свёртки косвенным путём – посредством выполнения обратного преобразования Фурье от произведения ДПФ исходных периодических временных последовательностей. Такое решение имеет существенное практическое значение в том случае, когда преобразования могут быть выполнены быстрее, чем сама свёртка.

    Эффективными процедурами вычисления ДПФ являются алгоритмы быстрого преобразования Фурье (БПФ), начало разработки которых было положено исследованиями Кули и Тьюки [19]. Общий принцип построения алгоритмов БПФ заключается в разложении процесса нахождения ДПФ на вычислении преобразований последовательностей меньшего размера. Действительно, как видно из выражения (13.6), суть ДПФ заключается в получении суммы произведений вида


    Весовая функция является периодической с интервалом N. Для заданного k при изменении n от нуля до N-1 произведение (13.8) также меняется периодически, причём число этих периодов равно k. Более того, даже в пределах одного периода могут образоваться комплексно-сопряжённые пары. Используя данные свойства, можно существенно уменьшить необходимое для вычислений число умножений. Для этого входная последовательность разбивается на две вспомогательные последовательности — и (рис. 13.3, а) – так, что в попадают отсчёты только с чётными номерами, а в — с нечётными, т.е.

    Тогда выражение (13.6) записывается

    где и — ДПФ последовательностей и соответственно.

    Вычислительные затраты на реализацию N – точечного ДПФ в соответствии с выражением (13.9) составляют по операций на вычисление N/2 – точечных ДПФ и и дополнительно N операций на их соединение, тогда как прямое вычисление по (13.6) требует операций, что почти в два раза больше.

    В выражении (13.9) индекс k изменяется от 0 до N-1. Однако и имеют период N/2 и вычисляются только в диапазоне от 0 до N/2-1. Используя свойство периодичности функции

    можно найти искомое соотношение и для Окончательно получаем [4]:

    Это соотношение между , и изображено на рис 13.3,б, который показывает, как 8-точечное ДПФ сведено к двум 4-точечным.

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

    Данный алгоритм получил название БПФ с прореживанием по времени. Базовая операция его, называемая «бабочкой» (баттерфляй), состоит в получении выходных чисел P и Q из выходных A и B в соответствии с выражением

    Граф базовой операции изображён на рис. 13.3,в. Очевидно, что значения, получаемые в результате выполнения базовой операции, могут быть хранимы в запоминающем устройстве на месте исходных операндов, что существенно сокращает требуемый объём запоминающих устройств. Алгоритм БПФ, в котором реализуется такая переадресация данных, получил название алгоритма с замещением.

    Второй распространённой формой алгоритма БПФ является метод прореживания по частоте, суть которого заключается в выполнении ДПФ от вспомогательных подпоследовательностей и половинной длины, получаемых из исходной в соответствии с выражением:

    Базовая операция БПФ с прореживанием по частоте имеет вид:

    и может быть представлена графом, изображённым на рис. 13.3,г.

    Легко заметить сходство между алгоритмами БПФ с прореживанием по времени и по частоте. В обоих случаях при выполнении ДПФ требуется равное количество операций, и в том, и в другом случае вычисления могут быть проведены с замещением, иногда совпадают даже конфигурации графов преобразования, хотя значения весовых коэффициентов различны. Ещё одно общее свойство обоих видов алгоритма заключается в необходимости перестановки следования входных либо выходных отсчётов преобразования, выполняемое обычно с помощью двоичной инверсии их номеров. Потребность в такой перестановке вытекает из самой структуры алгоритмов. Действительно, при прореживании по времени для получения выходной последовательности в естественном порядке на выходе устройства необходимо выполнить перестановку данных таким образом, чтобы первым наступал нулевой отсчёт, вторым – (N-1)/2-й и т.д. В алгоритме прореживания по частоте аналогичная перестановка производится с выходными отсчётами.

    Важной особенностью алгоритма БПФ является возможность использования его и для вычисления обратного дискретного преобразования Фурье (ОДПФ) [4,5], которое определяется выражением (13.7). Выполнив комплексное сопряжение последнего и разделив его на N, получим

    Правая часть формулы (13.14) представляет собой ДПФ последовательности и может быть вычислена с использованием алгоритма БПФ. Тогда искомую последовательность можно получить, взяв комплексно-сопряжённое с (13.14) выражение и умножив его на N:

    На практике здесь обычно используется перестановка всех, кроме , членов последовательности таким образом, что меняется местами с ; — с и т.д.

    Существуют различные алгоритмы БПФ [3], однако все они могут быть получены путём последовательного применения одной операции – представления одномерного массива двумерным [5]. При этом если число членов входной последовательности преобразования N является простым, т.е. не разложимым на произведения меньших целых чисел, то одномерный сигнал невозможно представить в виде двумерного массива, поэтому для такого сигнала не существует алгоритма БПФ. Для выхода из этого положения в большинстве практических задач исходную последовательность удлиняют путём дополнения её нулями до желаемого числа членов, являющегося составной величиной, после чего становится возможным выполнение преобразования с фиксированным или смешанным основанием. Кроме того, в последнее время разработан также ряд специальных быстрых алгоритмов вычисления ДПФ на основе вспомогательных преобразований, с использованием цифровых фильтров, на базе представления ДПФ в виде ЦС [3], позволяющих выполнять преобразования для различных размерностей массивов входных данных, в том числе и для простых N.

    С вычислением выражений вида свёртки, как указано выше, связана также широкая область цифровой обработки сигналов, например, цифровая фильтрация, под которой подразумевается выделение с использованием цифровых методов полезного сигнала на фоне флюктуирующих шумов в определённом диапазоне частотного спектра [10].

    Математически работа цифрового фильтра ЦФ во временной области описывается разностным уравнением (уравнением в конечных разностях)

    где и — n-е отсчёты входного и выходного сигналов фильтра; коэффициенты и представляют собой константы ЦФ с постоянными параметрами, или отсчёты решётчатых функций, зависящих от n или других показателей (ЦФ с переменными параметрами).

    Цифровые фильтры принято делить на два класса: нерекурсивные и рекурсивные [18]. Если в выражении (13.16) все коэффициенты равны нулю, то фильтр, реализующий этот алгоритм, называется нерекурсивным и алгоритм его работы записывается как

    Если в (13.16) хотя бы один из коэффициентов не равен нулю, то такой фильтр называется рекурсивным. Очевидно, что нерекурсивный цифровой фильтр (НРЦФ) представляет собой устройство без обратной связи, а рекурсивный ЦФ (РЦФ) – устройство с обратной связью. В отдельных случаях выделяют также в виде отдельных классов полурекурсивные, адаптивные и прогнозирующие ЦФ [10].

    Возможна и другая классификация цифровых фильтров: ЦФ с конечными (КИХ-фильтры) и с бесконечными (БИХ-фильтры) импульсными характеристиками [5,13]. При этом любой реальный нерекурсивный ЦФ является КИХ-системой, а рекурсивный фильтр представляет собой, как правило, БИХ-систему. Однако возможны варианты НРЦФ, являющиеся системами с импульсными характеристиками конечной длины [18].

    Структурная реализация ЦФ может иметь различный вид. Весьма широкое распространение получила каскадная форма представления фильтров, элементами которой являются биквадратные звенья [4]. Функционирование последних описывается рекуррентными соотношениями второго порядка

    Структурная схема возможной реализации биквадратного блока представлена на рис 13.3,д.

    Широкому распространению биквадратного звена в цифровой фильтрации служит не только факт относительно более низкой (в данном случае) чувствительности характеристик ЦФ к погрешности задания коэффициентов (связанных с физическими и стоимостными ограничениями на разрядность), но и гибкость в построении различных по сложности структур фильтров, возможность широкого использования эффективных методов распараллеливания, временного разделения и поточной обработки.

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

    Не нашли то, что искали? Воспользуйтесь поиском:

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

    Преобразование Фурье – преобразование, сопоставляющее функции некой вещественной переменной. Данная операция выполняется каждый раз, когда мы воспринимаем различные звуки. Ухо производит автоматическое «вычисление», выполнить которое наше сознание способно только после изучения соответствующего раздела высшей математики. Орган слуха у человека строит преобразование, в результате которого звук (колебательное движение условных частиц в упругой среде, которые распространяются в волновом виде в твердой, жидкой или газообразной среде) предоставляется в виде спектра последовательно идущих значений уровня громкости тонов разной высоты. После этого мозг превращает данную информацию в привычный всем звук.

    Математическое преобразование Фурье

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

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

    Историческая справка

    Первым человеком, применившим данный метод, стал французский математик Жан Батист Фурье. Преобразование, названное впоследствии его именем, изначально использовалось для описания механизма теплопроводности. Фурье всю свою сознательную жизнь занимался изучением свойств тепла. Он внес огромный вклад в математическую теорию определения корней алгебраических уравнений. Фурье являлся профессором анализа в Политехнической школе, секретарем Института египтологии, состоял на императорской службе, на которой отличился во время строительства дороги на Турин (под его руководством было осушено более 80 тысяч квадратных километров малярийных болот). Однако вся эта активная деятельность не помешала ученому заниматься математическим анализом. В 1802 году им было выведено уравнение, которое описывает распространение тепла в твердых телах. В 1807 году ученый открыл метод решения данного уравнения, которое и получило название «преобразование Фурье».

    Анализ теплопроводности

    Ученый применил математический метод для описания механизма теплопроводности. Удобным примером, в котором не возникает трудностей с вычислением, является распространение тепловой энергии по железному кольцу, погруженному одной частью в огонь. Для проведения опытов Фурье накалял докрасна часть этого кольца и закапывал его в мелкий песок. После этого проводил замеры температуры на противоположной его части. Первоначально распределение тепла является нерегулярным: часть кольца — холодная, а другая — горячая, между данными зонами можно наблюдать резкий градиент температуры. Однако в процессе распространения тепла по всей поверхности металла она становится более равномерной. Так, вскоре данный процесс приобретает вид синусоиды. Сначала график плавно нарастает и так же плавно убывает, точно по законам изменения функции косинуса или синуса. Волна постепенно выравнивается и в результате температура становится одинаковой на всей поверхности кольца.

    Автор данного метода предположил, что начальное нерегулярное распределение вполне можно разложить на ряд элементарных синусоид. Каждая из них будет иметь свою фазу (первоначальное положение) и свой температурный максимум. При этом каждая такая компонента изменяется от минимума к максимуму и обратно на полном обороте вокруг кольца целое число раз. Составляющая, имеющая один период, была названа основной гармоникой, а значение с двумя и более периодами – второй и так далее. Так, математическая функция, которая описывает температурный максимум, фазу или позицию называет преобразованием Фурье от функции распределения. Ученый свел единую составляющую, которая трудно поддается математическому описанию, к удобному в обращении инструменту – рядам косинуса и синуса, в сумме дающим исходное распределение.

    Суть анализа

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

    Вызов современникам

    Алгоритм преобразования Фурье стал вызовом теоретическим основам математики того времени. В начале девятнадцатого века большинство выдающихся ученых, в том числе и Лагранж, Лаплас, Пуассон, Лежандр и Био, не приняли его утверждение о том, что начальное распределение температуры раскладывается на составляющие в виде основной гармоники и более высокочастотные. Однако академия наук не могла проигнорировать результаты, полученные математиком, и удостоила его премии за теорию законов теплопроводности, а также проведение сравнения ее с физическими экспериментами. В подходе Фурье главное возражение вызывал тот факт, что разрывная функция представлена суммой нескольких синусоидальных функций, которые являются непрерывными. Ведь они описывают разрывающиеся прямые и кривые линии. Современники ученого никогда не сталкивались с подобной ситуацией, когда разрывные функции описывались комбинацией непрерывных, таких как квадратичная, линейная, синусоида либо экспонента. В том случае, если математик был прав в своих утверждениях, то сумма бесконечного ряда тригонометрической функции должна сводиться к точной ступенчатой. В то время подобное утверждение казалось абсурдным. Однако, несмотря на сомнения, некоторые исследователи (например Клод Навье, Софи Жермен) расширили сферу исследований и вывели их за пределы анализа распределения тепловой энергии. А математики тем временем продолжали мучиться вопросом о том, может ли сумма нескольких синусоидальных функций сводиться к точному представлению разрывной.

    200-летняя история

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

    Преобразование Фурье сегодня

    В наши дни изучение данного процесса в основном сводится к нахождению эффективных методов перехода от функции к ее преобразованному виду и обратно. Такое решение называется прямое и обратное преобразование Фурье. Что это значит? Для того чтобы определить интеграл и произвести прямое преобразование Фурье, можно воспользоваться математическими методами, а можно и аналитическими. Несмотря на то что при их использовании на практике возникают определенные трудности, большинство интегралов уже найдены и внесены в математические справочники. С помощью численных методов можно рассчитывать выражения, форма которых основывается на экспериментальных данных, либо функции, интегралы которых в таблицах отсутствуют и их сложно представить в аналитической форме.

    До появления вычислительной техники расчеты таких преобразований были весьма утомительными, они требовали ручного выполнения большого количества арифметических операций, которые зависели от числа точек, описывающих волновую функцию. Для облегчения расчетов сегодня существуют специальные программы, позволившие реализовать новые аналитические методы. Так, в 1965 году Джеймс Кули и Джон Тьюки создали программное обеспечение, получившее известность как «быстрое преобразование Фурье». Оно позволяет экономить время проведения расчетов за счет уменьшения числа умножений при анализе кривой. Метод «быстрое преобразование Фурье» основан на делении кривой на большое число равномерных выборочных значений. Соответственно количество умножений снижается вдвое при таком же снижении количества точек.

    Применение преобразования Фурье

    Данный процесс используется в различных областях науки: в теории чисел, физике, обработке сигналов, комбинаторике, теории вероятности, криптографии, статистике, океанологии, оптике, акустике, геометрии и других. Богатые возможности его применения основаны на ряде полезных особенностей, которые получили название «свойства преобразования Фурье». Рассмотрим их.

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

    2. Преобразование является обратимым. Причем обратный результат имеет практически аналогичную форму, как и при прямом решении.

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

    4. Согласно теореме «свертки», данный процесс превращает сложную операцию в элементарное умножение.


    5. Дискретное преобразование Фурье может быть быстро рассчитано на компьютере с использованием «быстрого» метода.

    Разновидности преобразования Фурье

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

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

    3. Дискретное преобразование Фурье. Этот метод используется в компьютерной технике для проведения научных расчетов и для цифровой обработки сигналов. Для проведения данного вида расчетов требуется иметь функции, определяющие на дискретном множестве отдельные точки, периодические или ограниченные области вместо непрерывных интегралов Фурье. Преобразование сигнала в таком случае представлено как сумма синусоид. При этом использование «быстрого» метода позволяет применять дискретные решения для любых практических задач.

    4. Оконное преобразование Фурье является обобщенным видом классического метода. В отличие от стандартного решения, когда используется спектр сигнала, который взят в полном диапазоне существования данной переменной, здесь особый интерес представляет всего лишь локальное распределение частоты при условии сохранения изначальной переменной (время).

    5. Двумерное преобразование Фурье. Данный метод используется для работы с двумерными массивами данных. В таком случае сначала преобразование производится в одном направлении, а затем — в другом.

    Заключение

    Сегодня метод Фурье прочно закрепился в различных областях науки. Например, в 1962 году была открыта форма двойной ДНК-спирали с использованием анализа Фурье в сочетании с дифракцией рентгеновских лучей. Последние фокусировались на кристаллах волокон ДНК, в результате изображение, которое получалось при дифракции излучения, фиксировались на пленке. Данная картинка дала информацию о значении амплитуды при использовании преобразования Фурье к данной кристаллической структуре. Данные о фазе получили путем сопоставления дифракционной карты ДНК с картами, которые получены при анализе подобных химических структур. В результате биологи восстановили кристаллическую структуру — исходную функцию.

    Мастер Йода рекомендует:  Объектно-ориентированное программирование с примерами на C# — всё по этой теме для программистов

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

    Обработка изображений. Часть 1.

    Эта статья написана для людей, которые хорошо знают что такое БПФ, знают как с его помощью отфильтровать сигнал, вычислить функцию взаимной корреляции в одномерном случае, но испытывают определённые трудности в переносе этих алгоритмов в 2D.
    Цикл статей будет посвящен обработке изображений с помощью двумерного быстрого преобразования Фурье.
    Будут рассмотрены вопросы корреляции, фильтрации и элементарных преобразованиях изображений (масштабирование, поворот).

    Прямое двумерное преобразование Фурье.
    Двумерное преобразование очень просто получить делая обычное одномерное FFT сначала по строкам исходной матрицы, а потом по столбцам.
    Программа прямого и обратного 2D FFT написана но Фортране. Использует кошерную библиотеку IMSL для вызова одномерного прямого/обратного одномерного комплексного FFT.
    Претензии любителей Си, Дельфи, PHP и 1С-бухгалтерия — не принимаются!
    Специально для любителей Дельфи отмечу, что применение Фортрановской библиотеки IMSL может дать прирост производительности до 15%
    N,M — размерность входной матрицы изображения
    U — входная матрица (комплексная)
    UFFT — преобразование Фурье входной матрицы (комплексная матрица)

    Обратное двумерное преобразование Фурье:

    Программа настолько элементарна, что не было смысла искать аналог в учебниках или на просторах Интернета.
    Некоторое преимущество (в моём понимании): программа на Фортране легко переводится на любой язык. Надо только помнить, что индексы в Фортране начинаются с 1. 🙂
    Естественно, в собственной реализации можно использовать любую, нравящуюся Вам, программу одномерного комплексного FFT.
    Для простоты будем использовать информацию из цветных (8 бит на цвет) *.bmp — файлов.
    Но никто не мешает Вам получить изображение с оптической мыши, как описано в здесь.
    Или использовать для «тупых» приложений, типа «Аналитическое продолжение карты силы тяжести на земной поверхности в верхнее полупространство на высоту 16 км с помощью двумерного преобразования Фурье».
    Мы на первом этапе, по совету классиков, будем тренироваться на кошках колонках. 🙂

    Возьмём обычный *.bmp — файл. Структура файла описана здесь.
    Использовать файлы будем цветные 24 бита (1 байт на цвет). При желании можно превратить их в градации серого:
    Y= 0.0721*B+ 0.7154*G+0.2125*R
    Сделаем двумерное преобразование Фурье над каждым цветом:

    UR,UG,UB — комплексные матрицы R,G,B. В действительной части- информация о цвете из *.bmp — файла. Мнимая часть =0.
    Нужно каждое значение выходных матриц комплексного спектра UFFTR, UFFTG, UFFTB разделить на Width*Height
    Вот так выглядит модуль двумерного спектра:

    Если сделать обратное двумерное БПФ

    в действительной части выходных массивов UR,UG,UB будет информация о соответствующем цвете. Округлив это значение до байта можно создать новый *.bmp — файл полностью совпадающий с исходным.
    Вот теперь самое интересное 🙂
    Делая элементарные преобразования комплексных матриц UFFTR, UFFTG, UFFTB полученных после двумерного преобразования Фурье легко масштабировать, вращать и коррелировать изображения.
    Сделаем элементарное преобразование:

    После обратного двумерного преобразования Фурье массивов UFFT1R, UFFT1G, UFFT1B получим картинку:

    Картинка отзеркалилась и повернулась на 90° против часовой стрелки.
    Как можно видеть цвет отдельных пикселей исказился. Дело в том, что после обратного 2D FFT цвета (после округления в 1 байт) не всегда будут лежать в диапазоне 0-255. Могут появиться значения -1 или 256. Но это легко компенсировать 🙂
    Поворот на произвольный угол:

    Масштабирование.
    Очень полезная статья о масштабировании здесь.
    К сожалению там не упоминается напрямую метод с использование 2D FFT.
    Пример уменьшения картинки в 2 раза методом обрезания высокочастотных составляющих двумерного спектра:

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

    Рассмотрим такую распространённую операцию как умножение двух целых чисел. Квадратичный алгоритм — умножения в столбик — все знают со школы. Долгое время предполагалось, что ничего быстрее придумать нельзя.

    Первым эту гипотезу опроверг Анатолий Карацуба. Его алгоритм сводит умножение двух \(n\) -значных чисел к трём умножениям \(\frac<2>\) -значных чисел, что даёт оценку времени работы

    \[ T(n)=3T\left(\dfrac n 2\right)+O(n)=O\left(n^<\log_2 3>\right)\approx O(n^<1.58>) \]

    Чтобы перейти к алгоритму с лучшей оценкой, нам нужно сначала установить несколько фактов о многочленах.

    Умножение многочленов

    Обратим внимание на то, что любое число можно представить многочленом:

    \[ \begin A(x) &= a_0 + a_1\cdot x + a_2 \cdot x^2 + \dots + a_n \cdot x^n \\ &= a_0 + a_1\cdot 2 + a_2 \cdot 2^2 + \dots + a_n \cdot 2^n \end \]

    Основание x при этом может быть выбрано произвольно.

    Чтобы перемножить два числа, мы можем перемножить соответствующие им многочлены, а затем произвести каррирование: пройтись от нижних разрядов получившегося многочлена и «сдвинуть» переполнившиеся разряды:

    Прямая формула для произведения многочленов имеет вид

    \[ \left(\sum_^n a_i x^i\right)\cdot\left(\sum_^m b_j x^j\right)=\sum_^x^k\sum_a_i b_j \]

    Её подсчёт требует \(O(n^2)\) операций, что нас не устраивает. Подойдём к этой задаче с другой стороны.

    Интерполяция

    Теорема. Пусть есть набор различных точек \(x_0, x_1, \dots, x_\) . Многочлен степени \(n\) однозначно задаётся своими значениями в этих точках. (Коэффициентов у этого многочлена столько же, сколько и точек — прим. К. О.)

    Доказательство будет конструктивным — можно явным образом задать многочлен, который принимает заданные значения \(y_0, y_1, \ldots, y_n\) в этих точках:

    Корректность. Проверим, что в точке \(x_i\) значение действительно будет равно \(y\) :

    Для \(i\) -го слагаемого внутреннее произведение будет равно единице, если вместо \(x\) подставить \(x_i\) : в этом случае просто перемножается \((n-1)\) единица. Эта единица помножится на \(y_i\) и войдёт в сумму.

    Для всех остальных слагаемых произведение занулится: один из множетелей будет равен \((x_i — x_i)\) .

    Уникальность. Предположим, есть два подходящих многочлена степени \(n\) — \(A(x)\) и \(B(x)\) . Рассмотрим их разность. В точках \(x_i\) значение получившегося многочлена \(A(x) — B(x)\) будет равняться нулю. Если так, то точки \(x_i\) должны являться его корнями, и тогда разность можно записать так:

    \[ A(x) — B(x) = \alpha \prod_^n (x-x_i) \]

    для какого-то числа \(\alpha\) . Тут мы получаем противоречие: если раскрыть это произведение, то получится многочлен степени \(n+1\) , который нельзя получить разностью двух многочленов степени \(n\) .

    Этот многочлен называется интерполяционным многочленом Лагранжа, а сама задача проведения многочлена через точки — интерполяцией.

    Примечание. На практике интерполяцию решают методом Гаусса: её можно свести к решению линейного уравнения \(aX = y\) , где \(X\) это матрица следующего вида:

    \[ \begin 1 & x_0 & x_0^2 & \ldots & x_0^n \\ 1 & x_1 & x_1^2 & \ldots & x_1^n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n & x_n^2 & \ldots & x_n^n \\ \end \]

    Важный факт: многочлен можно однозначно задать не только своими коэффициентами, но также корнями и значениями хотя бы в \((n+1)\) -ой точке.

    Умножение через интерполяцию

    Что происходит со значениями многочлена-произведения \(A(x) B(x)\) в конкретной точке \(x_i\) ? Оно просто становится равным \(A(x_i) B(x_i)\) .

    Основная идея алгоритма: если мы знаем значения в каких-то различных \(n + m\) точках для обоих многочленов \(A\) и \(B\) , то, попарно перемножив их, мы за \(O(n + m)\) операций можем получить значения в тех же точках для многочлена \(A(x) B(x)\) — а с их помощью можно интерполяцией получить исходный многочлен и решить задачу.

    Если притвориться, что evaluate и interpolate работают за линейное время, то умножение тоже будет работать за линейное время.

    К сожалению, непосредственное вычисление значений требует \(O(n^2)\) операций, а интерполяция — как методом Гаусса, так и через символьное вычисление многочлена Лагранжа — и того больше, \(O(n^3)\) .

    Но что, если бы мы могли вычислять значения в точках и делать интерполяцию быстрее?

    Комплексные числа


    Определение. Комплексные числа — это числа вида \(a + bi\) , где \(a\) и \(b\) это обычные вещественные числа, а \(i\) это так называемая мнимая единица: это число, для которого выполняется равенство \(i^2 = -1\) .

    Комплексные числа ввели в алгебре, чтобы работать с корнями из отрицательных чисел: \(i\) в каком-то смысле равно \(\sqrt<-1>\) . Так же, как и отрицательные числа, они как бы «не существуют» в реальном мире, а только в сознании математиков.

    С комплексными числами можно работать почти так же, как с действительными. Они даже удобнее: все квадратные корни всегда извлекаются, все корни многочленов всегда находятся.

    Комплексная плоскость

    Комплексные числа удобно изображать на плоскости в виде вектора \((a, b)\) и считать через них всякую геометрию.

    Модулем комплексного числа называется действительное число \(r = \sqrt\) . Геометрически, это длина вектора \((a, b)\) .

    Аргументом комплексного числа называется действительное число \(\phi \in (-\pi, \pi]\) , для которого выполнено \(\tan \phi = \frac\) . Геометрически, это значение угла между \((a, 0)\) и \((a, b)\) . Для нуля — вектора \((0, 0)\) — аргумент не определён.

    Таким образом комплексное число можно представить в полярных координатах:

    \[ a+bi = r ( \cos \phi + i \sin \phi ) \]

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

    Упражнение. Докажите это.

    Формула эйлера

    Определим число Эйлера \(e\) как число со следующим свойством:

    \[ e^ = \cos \phi + i \sin \phi \]

    Просто введём такую нотацию для выражения \(\cos \phi + i \sin \phi\) . Не надо думать, почему это так.

    Геометрически, все такие точки живут на единичном круге:

    Такая нотация удобна, потому что можно обращаться с \(e^\) как с обычной экспонентой. Пусть мы, например, хотим перемножить два числа на единичном круге с аргументами \(a\) и \(b\) . Тогда это можно записать так:

    \[ (\cos a + i \sin a) \cdot (\cos b + i \sin b) = e^ \]

    Упражнение. Проверьте это: раскройте скобки и проделайте немного алгебры.

    Корни из единицы

    У комплексных чисел есть много других замечательных свойств, но нам для алгоритма на самом деле потребуется только следующее:

    Утверждение. Для любого натурального \(n\) есть ровно \(n\) комплексных «корней из единицы», то есть чисел \(w_k\) , для которых выполнено:

    А именно, это будут числа вида:

    На комплексной плоскости эти числа располагаются на единичном круге на равном расстоянии друг от друга:

    Все 9 комплексных корней степени 9 из единицы

    Первый корень \(w_1\) (точнее второй — единицу считаем нулевым корнем) называют образующим корнем степени \(n\) из единицы. Возведение его в нулевую, первую, вторую и так далее степени порождает последовательность нужных корней единицы, при этом на \(n\) -ном элементе последовательность зацикливается:

    Будем обозначать \(w_1\) как просто \(w\) .

    Упражнение. Докажите, что других корней быть не может.

    Дискретное преобразование Фурье

    Дискретным преобразованием Фурье называется вычисление значений многочлена в комплексных корнях из единицы:

    Обратным дискретным преобразованием Фурье называется, как можно догадаться, обратная операция — интерполяция коэффициентов \(x_i\) по значениям \(X_i\) .

    Почему эта формула верна? При вычислении ПФ мы практически применяем матрицу к вектору:

    \[ \begin w^0 & w^0 & w^0 & w^0 & \dots & w^0 \\ w^0 & w^1 & w^2 & w^3 & \dots & w^ <-1>\\ w^0 & w^2 & w^4 & w^6 & \dots & w^ <-2>\\ w^0 & w^3 & w^6 & w^9 & \dots & w^ <-3>\\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ w^0 & w^ <-1>& w^ <-2>& w^ <-3>& \dots & w^1 \end\begin a_0 \\ a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_ \end = \begin y_0 \\ y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_ \end \]

    То есть преобразование Фурье — это просто линейная операция над вектором: \(W a = y\) . Значит, обратное преобразование можно записать так: \(a = W^<-1>y\) .

    Как будет выглядеть эта \(W^<-1>\) ? Автор не будет пытаться изображать логичный способ рассуждений о её получении и сразу её приведёт:

    \[ W^ <-1>= \dfrac 1 n\begin w^0 & w^0 & w^0 & w^0 & \dots & w^0 \\ w^0 & w^ <-1>& w^ <-2>& w^ <-3>& \dots & w^ <1>\\ w^0 & w^ <-2>& w^ <-4>& w^ <-6>& \dots & w^ <2>\\ w^0 & w^ <-3>& w^ <-6>& w^ <-9>& \dots & w^ <3>\\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ w^0 & w^ <1>& w^ <2>& w^ <3>& \dots & w^ <-1>\end \]

    Проверим, что при перемножении \(W\) и \(W^<-1>\) действительно получается единичная матрица:

    Значение \(i\) -того диагонального элемента будет равно \(\frac<1> \sum_k w^ w^ <-ki>= \frac<1> n = 1\) .

    Значение любого недиагонального ( \(i \neq j\) ) элемента \((i, j)\) будет равно \(\frac<1> \sum_k w^ w^ <-jk>= \frac<1> \sum_k w^k w^ = \frac> \sum_k w^k = 0\) , потому что все комплексные корни суммируются в ноль, то есть \(\sum w^k = 0\) (см. картинку — там всё симметрично).

    Внимательный читатель заметит симметричность форм \(W\) и \(W^<-1>\) , а также формул для прямого и обратного преобразования. На самом деле, эта симметрия нам сильно упростит жизнь: для обратного преобразования Фурье можно использовать тот же алгоритм, только вместо \(w^k\) использовать \(w^<-k>\) , а в конце результат поделить на \(n\) .

    Зачем это надо?

    Напомним, что мы изначально хотели перемножать многочлены следующим алгоритмом:

    Посчитаем значения в \(n+m\) каких-нибудь точках обоих многочленов

    Перемножим эти значения за \(O(n+m)\) .

    Интерполяцией получим многочлен-произведение.

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

    Соответствующий алгоритм называется быстрым преобразованием Фурье (англ. fast Fourier transform). Он использует парадигму «разделяй-и-властвуй» и работает за \(O(n \log n)\) .

    Схема Кули-Тьюки

    Обычно, алгоритмы «разделяй-и-властвуй» делят задачу на две половины: на первые \(\frac<2>\) элементов и вторые \(\frac<2>\) элементов. Здесь мы поступим по-другому: поделим все элементы на чётные и нечётные.

    Представим многочлен в виде \(P(x)=A(x^2)+xB(x^2)\) , где \(A(x)\) состоит из коэффициентов при чётных степенях \(x\) , а \(B(x)\) — из коэффициентов при нечётных.

    Пусть \(n = 2k\) . Тогда заметим, что

    Зная это, исходную формулу для значения многочлена в точке \(w^t\) можно записать так:

    Ключевое замечание: корней вида \(w^<2t>\) в два раза меньше, потому что \(w^n = w^0\) , и можно сказать, что.

    У нас по сути в два раза меньше корней (но они так же равномерно распределены на единичной окружности) и в два раза меньше коэффициентов — мы только что успешно уменьшили нашу задачу в два раз.

    Сам алгоритм заключается в следующем: рекурсивно посчитаем БПФ для многочленов \(A\) и \(B\) и объединим ответы с помощью формулы выше. При этом в рекурсии нам нужно считать значения на корнях степени не \(n\) , а \(k = \frac<2>\) , то есть на всех «чётных» корнях степени \(n\) (вида \(w^<2t>\) ).


    Заметим, что если \(w\) это образующий корень степени \(n = 2k\) из единицы, то \(w^2\) будет образующим корнем степени \(k\) , то есть в рекурсию мы можем просто передать другое значение образующего корня.

    Таким образом, мы свели преобразование размера \(n\) к двум преобразованиям размера \(\dfrac n 2\) . Следовательно, общее время вычислений составит

    \[ T(n)=2T\left(\dfrac n 2\right)+O(n)=O(n\log n) \]

    Заметим, что предположение о делимости \(n\) на \(2\) имело существенную роль. Значит, \(n\) должно быть чётным на каждом уровне, кроме последнего, из чего следует, что \(n\) должно быть степенью двойки.

    Реализация

    Приведём код, вычисляющий БПФ по схеме Кули-Тьюки:

    Как обсуждалось выше, обратное преобразование Фурье удобно выразить через прямое:

    Теперь мы умеем перемножать два многочлена за \(O(n \log n)\) :

    Примечание. Приведённый выше код, являясь корректным и имея асимптотику \(O(n\log n)\) , едва ли пригоден для использования на реальных контестах. Он имеет большую константу и далеко не так численно устойчивый, чем оптимальные варианты написания быстрого преобразования Фурье. Мы его приводим, потому что он относительно простой.

    Читателю рекомендуется самостоятельно задуматься о том, как можно улучшить время работы и точность вычислений. Из наиболее важных недостатков:

    внутри преобразования не должно происходить выделений памяти

    работать желательно с указателями, а не векторами

    корни из единицы должны быть посчитаны наперёд

    Следует избавиться от операций взятия остатка по модулю

    Вместо вычисления преобразования с \(w^<-1>\) можно вычислить преобразование с \(w\) , а затем развернуть элементы массива со второго по последний.

    Здесь приведена одна из условно пригодных реализаций.

    Но главная проблема в численной стабильности — мы нарушили первое правило действительных чисел. Однако, от неё можно избавиться.

    Number-theoretic transform

    Нам от комплексных чисел на самом деле нужно было только одно свойство: что у единицы есть \(n\) «корней». На самом деле, помимо комплексных чисел, есть и другие алгебраические объекты, обладающие таким свойством — например, элементы кольца вычетов по модулю.

    Нам нужно просто найти такую пару \(m\) и \(g\) (играющее роль \(w_n^1\) ), такую что \(g\) является образующим элементом, то есть \(g^n \equiv 1 \pmod m\) и при для всех остальных \(k все степени \(g^k\) различны по модулю \(m\) . В качестве \(m\) на практике часто специально берут «удобные» модули, например

    \[ m = 998244353 = 7 \cdot 17 \cdot 2^ <23>+ 1 \]

    Это число простое, и при этом является ровно на единицу больше числа, делящегося на большую степень двойки. При \(n=2^23\) подходящим \(g\) является число \(31\) . Заметим, что, как и для комплексных чисел, если для некоторого \(n=2^k\) первообразный корень — \(g\) , то для \(n=2^\) первообразным корнем будет \(g^2 (mod m)\) . Таким образом, для \(m=998244353\) и \(n=2^k\) первообразный корень будет равен \(g=31 \cdot 2^ <23-k>(mod m)\) .

    Реализация практически не отличается.

    С недавнего времени некоторые проблемсеттеры начали использовать именно этот модулю вместо стандартного \(10^9+7\) , чтобы намекнуть (или сбить с толку), что задача на FFT.

    Применения

    Даны две бинарные строки \(a\) и \(b\) . Нужно найти такой циклический сдвиг строки \(b\) , что количество совпадающих соответствующих символов с \(a\) станет максимально.

    Сперва научимся для каждого циклического сдвига \(i\) второй строки считать количество совпадающих единиц \(c_i\) . Это можно сделать за \(O(n^2)\) множеством разных способов, мы рассмотрим следующий: рассмотрим каждую единицу во втором числе, пусть она стоит на \(j\) -й позиции; для каждого \(l\) от \(0\) до \(n-1\) , если \(a_l\) равно 1, то прибавим один к \(c_\) (при этом \(i-j\) берётся по модулю \(n\) ). Такой алгоритм верный, потому что по сути мы перебираем пары единиц, которые могут совпадать, и прибавляем +1 к количеству совпадающих единиц для соответствующего циклического сдвига. И тут мы можем заметить очень важную вещь: если перемножить числа, соответствующие \(a\) и \(b\) , в столбик и не переносить разряды при сложении, то мы получим как раз массив \(c\) (с одним нюансом: его длина может быть больше \(n\) , тогда нам нужно для всех \(i \geq n\) прибавить \(c_i\) к \(c_\) )! А перемножать длинные числа мы уже научились: это легко сделать при помощи БПФ. Таким образом, мы научились искать число совпадающих единиц; заметим, что мы можем инвертировать биты в строках и применить эквивалентный алгоритм, получив в итоге количества совпадающих нулей. Сложим соответствующие элементы в двух массивах и найдём индекс максимального. Также очень часто в задачах на FFT требуется не явно перемножить два полинома, а посчитать свёртку двух векторов. Прямой свёрткой векторов \(a\) длины \(n\) и \(b\) длины \(m\) называется вектор \(s\) длины \(n+m-1\) такой, что \(s_k = \Sigma_^ a_i \cdot b_ (\forall k \in [0;n+m-2])\) (при этом считается, что несуществующие элементы равны нулю). Круговой (циклической) свёрткой векторов \(a\) и \(b\) длины \(n\) называется вектор \(s\) длины \(n\) такой, что \(s_k = \Sigma_^ a_i \cdot b_ (\forall k \in [0; n-1])\) (при этом \(\) берётся по модулю \(n\) ). Оказывается, что линейную свёртку можно считать через круговую: для этого дополним нулями оба вектора до одинаковой длины \(n+m-1\) . Это очень легко доказать: если для некоторого \(k\) \(i \geq k+1\) , то либо \(a_i\) , либо \(b_\) будут равны нулю. Если расписать выражение для прямого преобразования Фурье круговой свёртки и перенести множетили, то можно получить, что круговая свёртка равна вектору произведений многочленов с коэффициентами \(a\) и \(b\) в точках \(0,1,\dots n-1\) . Возможно, когда-нибудь я это распишу.

    PARALLEL.RU — Информационно-аналитический центр по параллельным вычислениям

    Вводные замечания

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

    Постановка задачи

    Пусть f = (f, . ,fN-1) T — исходный вектор, который нужно умножить на матрицу Фурье F c элементами

    для получения вектора a= (a, . aN-1) T по формуле a = Ff. Тогда для элемента искомого вектора умножение даст формулу:

    Как известно, в общем случае умножение матрицы нам вектор требует O(N 2 ) арифметических операций. Следует рассмотреть вопрос о сокращении этого числа.

    1 Вообще говоря, в формулах ДПФ для физической корректности часто присутствует и нормировочный множитель, которого у нас нет. Здесь мы будем использовать терминологию, принятую в § 17 справочника В.В.Воеводина и Ю.А.Кузнецова «Матрицы и вычисления» (М., Наука, 1984г.). Теоретические оценки, приводимые ниже, взяты тоже оттуда.

    Структура вычислений классического ДПФ

    Рассмотрим структуру вычислений классического ДПФ, предположив, что суммирование происходит по порядку от малых значений индексов к большим. Казалось бы, что сложного в банальном умножении матрицы на вектор? Однако, в отличие от стандартного умножения, где каждый элемент матрицы используется только один раз, в ДПФ матрица — специального вида. Это приводит к тому, что одни и те же элементы матрицы используются в совершенно разных местах схемы вычислений. На рисунке представлена структура ДПФ с N=6.

    Чёрными линиями обозначены собственно ветви вычислений, при этом передачи осуществляются слева направо, в крайнем справа столбце результатом будет соответствующий элемент вычисляемого вектора. В одном столбце расположены те операции, которые используют в качестве одного из данных один и тот же элемент исходного вектора. А вот цвет вершины (все они соответствуют операции a+bc) соответствует тому, какой из элементов матрицы Фурье умножается на элемент исходного вектора. Серым цветом окрашены вершины, где из матрицы Фурье в качестве множителя берётся нулевая степень комплексного числа, т.е. просто единица (и где, следовательно, умножение можно пропустить). Как видим, её как раз больше всего в структуре вычислений. Далее степени одного и того же комплексного множителя растут от 1 до 5 в вершинах соответственно зелёного, голубого, синего (элемент там умножается на -1 и потому можно, как и в серых вершинах, отказаться от умножения, заменив сложение вычитанием), фиолетового и жёлтого цветов. Как видно, структура использования этих степеней довольно «хаотична», хотя на самом деле в этом «хаосе» есть строго определённый порядок.

    В случае простого N структура будет отличаться тем, что «серых» вершин, кроме как на левой и верхней границах соответствующего рисунка, больше не будет вообще.

    В случае же N=2 у нас будут только 3 сложения и одно вычитание, так что использование двойки в качестве базового сомножителя БПФ несёт несомненные выгоды.

    Структура классического ДПФ в целом настолько проста, что её можно записать с помощью такой простой фортран-подпрограммы:

    (программа записана в нотации Фортрана 77).

    Случай составного N

    Рассмотрим случай N = p1 p2, причём сомножители не равны 1. Тогда, представим индексы элементов векторов и матрицы q, j в виде

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

    и то, что мнимая экспонента периодична с периодом 2π, получим, что

    При вычислении по представленному алгоритму требуется всего

    арифметических операций. В случае, когда сомножители равны друг другу, общее количество операций составит

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

    Сведение преобразования Фурье к последовательности преобразований меньшей размерности

    Приведённые выше формулы показывают интересную вещь. Оказывается, если записать исходный вектор f по столбцам прямоугольной таблицы размера p2 на p1, то его преобразование Фурье в вектор a может быть выполнено в следующие три действия (на рисунках-примерах мы будем считать, что исходная размерность равна 15):

      Вычисление преобразований Фурье порядка p1всех строк таблицы (красным на рисунке показаны группы элементов массива, над которыми выполняются ДПФ размерности p1=3).

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


    Вычисление преобразований Фурье порядка p2 всех столбцов таблицы (красным на рисунке показаны группы элементов массива, над которыми выполняются ДПФ размерности p2=5).

    Искомый вектор a будет в получившейся таблице записан по строкам. Если же мы хотим, чтобы он тоже был записан по столбцам, как и исходный, то нам надо будет добавить четвёртое действие — транспонирование таблицы.

    Указанная особенность — основными действиями являются те же преобразования Фурье, но меньшей размерности — позволяет применить приём БПФ и к ним, в том случае, если сомножители, образующие составное N, сами окажутся составными. Разработано много программ, которые доводят разложение до простых делителей числа N. Их детальный разбор можно прочитать, например в известной книге Блейхута. Самые известные из них, однако, получаются в следующем частном случае.

    Размерность — степень двойки

    Пусть N = 2 k . Тогда вышеизложенную схему сведения преобразования Фурье к преобразованиям меньшей размерности можно выполнить особенно эффективно. Этому способствуют следующие причины:

    • При выполнении преобразования размерности 2 вычисляемые экспоненты равны ±1, что дополнительно экономит на операциях умножения, оставляя только сложения и вычитания.
    • Поворачивающие множители первого шага дают нам полный набор множителей второго и последующих шагов, что также даёт нам экономию (хоть и меньшего порядка, чем главные преобразования) в вычислениях как тригонометрических функций, так и комплексных произведений. О влиянии этого на ошибки округления в ДПФ будет сказано отдельно. При этом не равны 1 только множители второго столбца.

    Общая структура одного шага БПФ (без учёта транспонирований, но с учётом умножений на поворачивающие множители) для N = 2 k даёт нам рисунок такого вида:

    Красными линиями обозначена передача начальных данных (массива из 2+1 пар действительных чисел, каждая пара — действительная и мнимая часть соответствующих комплексных чисел). Далее мы вверху видим простую запись результатов, а внизу — структуру с умножением на поворачивающий множитель, не равный 1. Зелёными стрелками обозначена запись данных в массив (и здесь мы видим в середине структуру, заменившую транспонирование матрицы 2х2). Квадратами же обозначены элементы массива (красными — начальные данные шага БПФ, зелёными — вычисленное подпрограммой их ДПФ с умножением на множители, розовыми — поворачивающий множитель).

    Мастер Йода рекомендует:  Фотогалерея для своего сайта – все способы реализации

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

    На рисунке — общая структурная схема БПФ для N=16. Зелёные вершины соответствуют выполняемым макрооперациям, описанным выше. Передача данных показана чёрными и синими стрелками. Ярусы параллельной формы макрографа БПФ расположены по вертикалям. Если всё делать таким образом, при выполнении БПФ для N = 2 k известно, что оно может быть реализовано за N log2N / 2 операций комплексного умножения и N log2N операций комплексного сложения.

    Здесь мы специально не приводим оценки через действительные умножения и сложения, поскольку для самих комплексных умножений есть разные реализации — как обычная (через 4 умножения и 2 сложения), так и «быстрая» (через 3 умножения и 3 сложения).

    Как видно из рисунка, структура пересылок нелинейная, она экспоненциально зависит от номера шага БПФ. Однако при правильной организации сочетания транспонирований структура пересылок будет неизменной от шага к шагу БПФ. Платой за это будет, однако, то, что половина дуг будет принадлежать к классу длинных, т.е. передача данных не будет локальной, а также то, что предварительно данные надо будет расположить в так называемом битинверсивном порядке (то есть в таком, где двоичное представление индекса массива идёт наоборот). Эта подзадача (начальной пересортировки) занимает примерно N пересылок данных, но сложна логически, поэтому её решение мы здесь приводить не будем.

    Для реализации на ПЛИСах есть пример фрагмента программы на языке Colamo, использующей БПФ (метод Кули-Тьюки для произвольной степени двойки) прямоугольного массива размера 1024х1024 по строкам, разработанный в НИИ МВС при ТРТУ. Вот её текст (мы будем периодически вставлять свои пояснения):

    var ReX_in, ImX_in, ReX_out, ImX_out: Array real [1024:Stream, 1024:Stream] Mem; // Входной и выходной массив строки, столбцы
    var ReX_in_com, ImX_in_com: Array real [1024:Stream, 1024:Stream, 10: Stream] Com; // Входной и выходной массив строки, столбцы, итерации
    var invert: Array integer [10: Stream,512: Stream] mem; // массив с адресами в битинверсивном порядке

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

    var N,i,j,k,m,p : Number;
    var ReImW : Array real [512:Stream] Mem; // массив коэффициентов БПФ
    var ReA_reg, ImA_reg, ReB_reg, ImB_reg : real Reg; // массив коэффициентов БПФ
    var ReA_com, ImA_com, ReB_com, ImB_com : real com; // массив коэффициентов БПФ

    Здесь видно, что и вычисление тригонометрии тоже оставлено на предварительную работу хост-компьютера.

    Define N=1024;
    ReadFile (inv, FF_Int_64, «inv.txt»); // чтение из файла массива с адресами в битинверсивном порядке
    ReadFile (ReImW, FF_Int_64, «ReImW.txt»); // чтение из файла массива коэффициентов
    ReadFile (AB, FF_Int_64, «In_com.txt»); // чтение из файла массива коэффициентов

    Начало описания структуры основной операции (так называемая «бабочка» БПФ):

    let But_FFT(in: ReA, ImA, ReB, ImB, ReW, ImW; // входные данные
    out: ReA_out, ImA_out, ReB_out, ImB_out); // выходные данные
    var ReA, ImA, ReB, ImB, ReA_out, ImA_out, ReB_out, ImB_out, ReW, ImW: real com;
    var Re_C, Im_C: real com; // Массив коммут. переменных соотв. структуре БПФ
    var mul_1, mul_2, mul_3, mul_4: real com; // Массив коммут. переменных соотв. структуре БПФ

    ReA_out := ReA + ReB;
    ImA_out := ImA + ImB;
    Re_C := ReA — ReB;
    Im_C := ImA — ImB;
    Mul_1 := Re_C*ReW;
    Mul_2 := Im_C*ImW;
    Mul_3 := Im_C*ReW;
    Mul_4 := Re_C*ImW;
    ReB_out := Mul_1 — Mul_2;
    ImB_out := Mul_3 + Mul_4;
    EndLet;

    Видно, что ключевая структура But_FFT реализует именно структуру одного шага БПФ, которую мы видели выше. Далее следует БПФ строк изображения:

    CADR Gor_FFT; // Кадр вычисления спектра строк излбражения
    for m := 0 to 2 do
    begin
    if m = 0 then // передача входных данных com — переменным
    begin
    for i := 0 to 1024 — 1 do // цикл по строкам
    begin
    for j := 0 to 1024 — 1 do // цикл по строкам
    begin
    ReX_in_com[i, j, 0] := ReX_in[i, j];
    ImX_in_com[i, j, 0] := ImX_in[i, j];
    end;
    end;
    end

    Здесь видно, что данные и результаты вычислений будут разнесены по итерациям в «фиктивном» массиве (который на самом деле не является массивом, поскольку имеет коммутационный тип, а только связывает операции друг с другом). Далее следует собственно основное гнездо циклов БПФ по строкам:

    else if m = 1 then // «рабочая» итерация
    begin
    for i := 0 to 1024 — 1 do // цикл по строкам
    begin
    for j := 0 to 10 — 1 do // цикл по итерациям
    begin
    for k := 0 to 512 — 1 do // цикл по Базовым операциям
    begin
    But_FFT(ReX_in_com[i,k*2,j], ImX_in_com[i,k*2,j], ReX_in_com[i,k*2+1,j], ImX_in_com[i,k*2+1,j], // входные данные
    ReImW[Invert[j,k]], ReImW[Invert[j,k]+1], // коэффициенты
    ReA_com, ImA_com, ReB_com, ImB_com); // выходные данные
    ReA_reg := ReA_com; // передача выходного потока через регистры для того что бы избежать «кольца»
    ImA_reg := ImA_com;
    ReB_reg := ReB_com;
    ImB_reg := ImB_com;
    for p := 0 to 1 do begin // цикл передачи выходных параметров
    if p = 0 then begin
    ReX_in_com[i,k,j+1] := ReA_reg;
    ImX_in_com[i,k,j+1] := ImA_reg;
    end
    else if p = 1 then
    begin
    ReX_in_com[i,k+512,j+1] := ReB_reg;
    ImX_in_com[i,k+512,j+1] := ImB_reg;
    end;
    end;
    end;
    end;
    end;
    end

    Рассмотрим основное тело программы. Видно, что внутренний цикл (по k) реализует всё необходимое множество базовых «бабочек» на каждом шаге БПФ, выполняя при этом перетасовку данных для следующего шага. Структура индексов массивов, как видно, линейна, но имеет включённые константы, сопоставимые с половиной этой размерности массива. Охватывающий его цикл по j реализует шаги БПФ, а внешний цикл по i — перебор всех строк.Собственно, БПФ практически уже выполнено.

    else if m = 2 then // итерация передачи данных в выходное КРП
    begin
    for i := 0 to 1024 — 1 do // цикл по строкам
    begin
    for j := 0 to 1024 — 1 do // цикл по столбцам
    begin
    ReX_out[i, j] := ReX_in_com[i, j, 10];
    ImX_out[i, j] := ImX_in_com[i, j, 10];
    end;
    end;
    end;
    end;
    ENDCADR;

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

    Всё, как видим, программа не очень сложна. Отметим ещё раз важную вещь: несмотря на наличие в библиотеке MyLibNewElement тригонометрических функций, все массивы коэффициентов для этой программы должны быть предвычислены заранее. Связано это с двумя факторами. Во-первых, размер массива изображения фиксирован и повторять одну и туже работу нет нужды, в то время как поток изображений, обрабатываемых на ПЛИСах, может идти довольно большим потоком. Во-вторых, использование тригонометрических функций даже в случае одиночных вычислений нарушило бы равномерность загрузки элементов ПЛИС и, как следствие, привело бы к простою значительной части ПЛИС. Кроме того, отметим, что использование в вычислениях массива, казалось бы, десятикратно более объёмного, чем исходный, вовсе не увеличивает память, используемую в ПЛИС-модулях, поскольку коммутационный тип массива исключает его хранение там, где оно не необходимо.

    Размерность — чётная степень двойки

    Описанные выше схемы и программа БПФ для произвольной степени двойки являются той реализацией метода Кули-Тьюки, который наиболее распространён среди реализаций БПФ. Однако этот алгоритм, несмотря на его преимущества (например, однотипность основных операций на каждом шаге), не является самым быстрым из реализаций БПФ на степенях двойки. В случае, если k (показатель степени двойки) чётно, то оказывается, что БПФ методом Кули-Тьюки для N = 2 k может быть реализовано даже за 3N log2 N / 8 операций комплексного умножения и N log2 N операций комплексного сложения. За счёт чего же чётная степень даёт нам дополнительный выигрыш?

    Рассмотрим преобразование Фурье порядка 4. Применим к нему разложение на 2 БПФ порядка 2 и увидим, что все поворачивающие множители равны ±1 либо ±i, что даёт нам возможность не умножать соответствующие элементы таблицы, а просто брать их с разными знаками и комбинировать их действительные и мнимые части. Подпрограмму прямого БПФ порядка 4 набора 4 комплексных чисел можно на Фортране записать, например, так:

    subroutine FT4D(A1,A2,A3,A4)
    real A1(2),A2(2),A3(2),A4(2),C1(2),C2(2),C3(2),C4(2)
    C1(1) = A1(1)+A3(1)
    C1(2) = A1(2)+A3(2)
    C2(1) = A2(1)+A4(1)
    C2(2) = A2(2)+A4(2)
    C3(1) = A1(1)-A3(1)
    C3(2) = A1(2)-A3(2)
    C4(1) = A2(1)-A4(1)
    C4(2) = A2(2)-A4(2)
    A1(1) = C1(1)+C2(1)
    A1(2) = C1(2)+C2(2)
    A3(1) = C1(1)-C2(1)
    A3(2) = C1(2)-C2(2)
    A2(1) = C3(1)-C4(2)
    A2(2) = C3(2)+C4(1)
    A4(1) = C3(1)+C4(2)
    A4(2) = C3(2)-C4(1)
    return
    end

    А в обратном БПФ порядка 4 будут несколько переставлены знаки и индексы:

    subroutine FT4R(A1,A2,A3,A4)
    real A1(2),A2(2),A3(2),A4(2),C1(2),C2(2),C3(2),C4(2)
    C1(1) = A1(1)+A3(1)
    C1(2) = A1(2)+A3(2)
    C2(1) = A2(1)+A4(1)
    C2(2) = A2(2)+A4(2)
    C3(1) = A1(1)-A3(1)
    C3(2) = A1(2)-A3(2)
    C4(1) = A2(1)-A4(1)
    C4(2) = A2(2)-A4(2)
    A1(1) = C1(1)+C2(1)
    A1(2) = C1(2)+C2(2)
    A3(1) = C1(1)-C2(1)
    A3(2) = C1(2)-C2(2)
    A4(1) = C3(1)-C4(2)
    A4(2) = C3(2)+C4(1)
    A2(1) = C3(1)+C4(2)
    A2(2) = C3(2)-C4(1)
    return
    end

    Но его структура в принципе — та же. На рисунке

    показана структура прямого преобразования Фурье размерности 4, т.е. одно из элементарных преобразований БПФ порядка чётной степени двойки. Красными линиями обозначена передача начальных данных (массива из 4 пар действительных чисел, каждая пара — действительная и мнимая часть соответствующих комплексных чисел). Далее мы вверху видим структуру, похожую на шаг БПФ любой степени двойки, а внизу — «запутанную» структуру похожего вида. Запутанность вызвана тем, что умножение на поворачивающие множители было заменено расстановкой знаков с взаимными заменами действительных и мнимых частей. Зелёными стрелками обозначена запись данных в массив (и здесь мы видим в середине структуру, заменившую транспонирование матрицы 2х2). Квадратами же обозначены элементы массива (красными — начальные данные БПФ, зелёными — вычисленное подпрограммой прямое ДПФ).

    Теперь вспомним про умножение на поворачивающие множители, которое завершает шаг БПФ порядка, разложенного по четвёркам, и увидим, что из 4 полученных после Фурье порядка 4 комплексных чисел умножению подвергаются 3. То есть из всего массива — 3/4 2 , в отличие от половины 2 в БПФ, разложенном по двойкам. Однако за счёт того, что самих шагов у нас вдвое меньше, мы и получаем выигрыш в множителе при N log2 N в количестве комплексных умножений.

    Если мы попытаемся применить такой же приём для уменьшения доли комплексных умножений, например, взяв в качестве базового ДПФ порядка 8 (степень двойки, таким образом, должна делиться на 3), то увидим, что на поворачивающие множители умножаются 7 чисел из 8, полученных после ДПФ порядка 8, и за счёт того, что этапов втрое меньше — и получим, казалось бы, коэффициент 7/24, то есть ещё меньше. Однако в базовом ДПФ порядка 8, как оказывается, нельзя уже полностью избавиться от комплексных умножений, и при их учёте мы снова получаем те же, что для чётных степеней, 3/8 в коэффициенте при N log2 N. Правда, при ближайшем рассмотрении оказывается, что часть из умножений является умножением на числа специального вида и может быть выполнена быстрее, но это усложняет общую схему. Однако, если усложнение не пугает, то можно довести оптимизацию до размерности 16.

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

    2 На самом деле — не 3/4 и не половина, а чуть меньше.

    Размерности, не являющиеся степенями двойки

    Реальные задачи таковы, что размерность ДПФ, которое подлежит «убыстрению», может оказаться вовсе не степенью двойки. Для тех специалистов, кто не хочет ломать голову над выбором алгоритма, в настоящее время распространён метод решения «от степени», при котором реальную задачу просто «подгоняют» под степень двойки (бывает даже, что это делают с реальными АЦП в устройствах). Однако такой метод во многих задачах неприемлем по разным причинам. Увеличение размера до ближайшей степени двойки может существенно увеличить размер задачи. Кроме этого, ряд задач традиционно ориентирован на размерности, заведомо не являющиеся степенями двойки (хоть и содержащие их в качестве делителей). Например, в геофизических климатических задачах одной из традиционных сеток является градусная (или кратные градусной). При этом получается, что размерность по долготе (с учётом периодичности) будет кратна 360, то есть содержать одну пятёрку и две тройки. Поэтому существуют такие алгоритмы БПФ, как методы Гуда-Томаса, Винограда и др. Надо сказать, что они, кроме прочего, являются и более быстрыми, чем алгоритм Кули-Тьюки для произвольной степени двойки. Поэтому во времена, когда ещё не было массовой компьютеризации и задачами БПФ занимались только специалисты, в программах БПФ царило многообразие. Ныне же большинство использует вовсе не самые быстрые БПФ, даже не подозревая о том, что вобще можно работать не со степенями двойки. Причём это связано не только с низким уровнем знаний о БПФ.

    Особенности реализации БПФ на параллельных вычислительных системах

    В случае, когда БПФ адаптируется для реализации на параллельной вычислительной системе (мы имеем в виду случай, когда нужно распараллеливать само преобразование Фурье, а не те алгоритмы, в которых преобразований Фурье разных векторов много и они сами могут выполняться параллельно или конвейерно), логическая сложность алгоритма, и, особенно, сложность деления данных по частям алгоритма, становится существенным тормозом для эффективной реализации. Кроме того, не секрет, что большинство многопроцессорных вычислительных систем с раздельной памятью (кластерная архитектура) часто предоставляет максимальное число узлов, являющееся степенью двойки. Системы, не имеющие таких ограничений, вроде ПЛИСов, имеют другие ограничения. Например, при программировании на COLAMO усложнение логической структуры при равной вычислительной мощности будет фактором, снижающим производительность базовых модулей , поскольку часть процессорных элементов при сложной логике программы будет либо простаивать, либо работать, но вхолостую (без использования результата). На кластерной архитектуре также, в случае разбиения схемы БПФ на блоки по процессорам, возникнут сложности, если разбивать придётся части схемы, реализующие ПФ меньшей размерности ( разных делителей N). Всё это также приводит к тому, что реализации для степеней двойки вытесняют остальные алгоритмы БПФ.

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

    Более правильным шагом по преодолению ограниченности ресурсов либо минимизации пересылок, видимо, является использование такой схемы, которая заложена сама в БПФ, при её конструировании. Идея, в общем-то, настолько очевидная, что, когда у программиста есть полное разложение N на простые множители N=2 k , она не всем приходит в голову: это использовать не полное разложение на простые множители, а разложение типа N=ML, где M, L — довольно большие числа, тоже являющиеся степенями двойки. Тогда мы видим, что все поворачивающие множители БПФ порядка N будут использованы только однажды — между БПФ порядка M и БПФ порядка L. Сами же эти БПФ порядков M, L можно, в зависимости от ресурсов, реализовать как БПФ степени двойки, либо тоже разложить на меньшие множители, и так — до той степени, пока ресурсов хватит, либо пересылок между компонентами МВС станет достаточно мало.

    Ошибки округления и преобразование Фурье

    Как нетрудно видеть по формулам, при выполнении преобразований Фурье большой размерности вероятно вычисление тригонометрических функций от больших аргументов. К сожалению, среди программ-реализаций преобразования Фурье автором было обнаружено сравнительно мало реализаций, выполняющих предварительное приведение аргумента с использованием целочисленной арифметики, а действительное приведение (используемое стандартными тригонометрическими функциями) даёт большую относительную ошибку. В самом деле, если порядок N и (следовательно) аргумента — 10 7 , то его приведение к диапазону первого периода тригонометрических функций даст ошибку, на 32-разрядой арифметике (относительная точность которой 10 -7 ) уже сравнимую с единицей, что делает преобразование неточным, а программу — ограниченно годной для небольших N.

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

    Заключение

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


    К сожалению, неумелое программирование, встречающееся даже в кодах фирменной поддержки БПФ на различных ускорителях, приводит к тому, что преимущества БПФ с порядком вектора, равным степени двойки, в известной степени снижаются из-за того, что тригонометрические функции вычисляются на каждом шаге и потому неизбежно увеличивают коэффициенты количества операций при главном члене (N log2 N). Между тем, предварительное предвычисление всей нужной для ДПФ тригонометрии даёт возможность выполнить тригонометрию за O(N) операций, причём даже это предвычисление может быть в основном занято не вызовами тригонометрических функций, а комплексными умножениями, как хотя бы в примере программы, использующей формулы приведения для разбивания всего диапазона на 16 частей.

    Также вызывает сожаление другое — редкое 3 использование разработчиками для современных параллельных вычислительных систем преимуществ чётной (а также кратной тройке и четвёрке) степени двойки. Схемы при этом получаются сложнее, но они не дают нам «холостых» ветвей и разбиений разной размерности. Поэтому разработчикам, желающим оптимизировать БПФ для степеней двойки, советуем обращать внимание на эту возможность.

    3 Редкое — не то слово. Нам такие реализации просто не встречались.

    Практические советы и замечания к выполнению БПФ

    Как выполнить преобразование на практике

    Берем сигнал длительностью , разбиваем его на частей, получаем мгновенные значения сигнала в точках:

    Если есть возможность выбрать , то следует выбрать степень двойки, так как соответствующий алгоритм немного быстрее. Если нет возможности выбрать, то не страшно — все равно есть алгоритм, не такой быстрый, но тоже имеющий сложность .

    Далее, создаем массив комплексных чисел, длиной :

    Теперь заполняем его. Допустим исходный сигнал у нас в массиве , и это действительная величина. Тогда просто копируем в реальную часть элементов в массиве , а мнимую часть обнуляем:

    Теперь осталось применить алгоритм БПФ. Если — т.е. степень двойки, то используем соответствующий алгоритм отсюда:

    Если же мы не можем гарантировать, что — степень двойки, тогда используем универсальный алгоритм отсюда:

    После исполнения этой функции массив изменится и будет теперь содержать результат прямого преобразования Фурье.

    Для обратного преобразования надо вызвать ту же функцию с последним аргументом . Если нужна действительная функция , то надо будет еще скопировать массив обратно, отбрасывая мнимую часть:

    Наглядное представление БПФ

    Полученный результат представляет собой последовательность комплексных чисел в форме пар: (реальная часть, мнимая часть). Но для понимания физической структуры сигнала желательно преобразовать их в амплитуды, частоты и фазы гармоник.

    Для преобразования можно использовать формулы:

    Самый первый элемент с индексом показывает постоянную составляющую сигнала, то есть, насколько он в среднем «приподнят» над осью времени. Остальные элементы представляют гармоники вида , где — это частота, которая для гармоники с индексом будет равна:

    Часто бывает дано не суммарное время сигнала , а частота дискретизации: сколько отсчетов приходится на одну секунду. Обозначим ее . Тогда частота для гармоники с индексом будет равна:

    В сумме постоянная составляющая и все гармоники должны дать исходный сигнал.

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

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

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

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

    В результате алгоритм получения наглядного представления будет таким:

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

    Преобразование туда и обратно

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

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

    В частности, если вы исправляете зеркальный эффект, вы (после прямого преобразование) обнуляете элементы с индексами выше и удваиваете элементы с индексами от до . Если теперь выполнить обратное преобразование, то реальная часть сигнала останется прежней, но мнимая часть сигнала станет ненулевой. Ее можно проигнорировать.

    Эффект размазывания

    Данный эффект подробно расматривался ранее, а здесь приведен среди других эффектов. Если частота исходной гармоники не равна в точности дной из частот , то вместо одной линии спектра получится смазанный пик — тем более смазанный, чем дальше отстоит исходная частота от ближайшей частоты .

    Здесь зеленым показана амплитуда спектральных гармоник, фиолетовым — фаза (как видите, в районе нужной частоты она делает резкий скачок). Исходной гармоникой здесь было колебание с амплитудой 30000, фазой 0 и частотой 440,2 Гц. Время дискретизации было равно 1 секунде при разбиении на 44100 отрезков. Ближайшие частоты вида равны 440 и 441 Гц, для них никакого размазывания не было бы.

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

    Эффекты модуляции

    А вот эта картинка —

    — получена совсем другим способом. Здесь была гармоника 440 Гц, но ее амплитуда в течение секунды линейно затухала от 30000 до 0. Картинки очень похожи, так что пойди разбери — то ли причина размазывания в том, что гармоника не попала точно в величину , то ли в том, что было затухание во времени.

    Различие там есть, просто оно не видно на картинке. При непопадании на фаза максимальной гармоники приблизительно совпадает с фазой одной из двух соседник гармоник и сильно отличается от фазы другой. А при модуляции фазы соседних гармоник или совпадают, или обе сильно отличаются от фазы максимальной. Чтобы показать это, на следующей паре картинок фаза и амплитуда отмечены не вертикальными рисками, а крестиками:

    На левой картинке — размазывание от несовпадения частот, на правой — от затухания. Так что в целом для различения этих случаев можно применять критерий: если фаза максимальной гармоники где-то посредине между фазами соседних гармоник, то это — модуляция, а если фаза максимальной гармоники гораздо близка к одной из «соседок», то это — несовпадение частот.

    Вот еще картинка:

    — здесь амплитуда линейно затухала от 30000 до 15000.

    — здесь амплитуда нелинейно нарастала от 0 до 30000.

    — эта скромная картинка получена с добавлением пульсации: амплитуда синусоидально увеличивается-уменьшается 10 раз за секунду. Две дополнительные частоты по бокам от основной обеспечили пульсацию основной частоты.

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

    — наконец, вот этот хаос возникает в результате «вибрато» — то есть, периодического изменения не амплитуды, а частоты. В данном случае на периоде 0.8 секунд произошло 8 периодических изменений частоты на четверть тона.

    НОВОСТИ ФОРУМА
    Рыцари теории эфира
    01.10.2020 — 05:20: ВОСПИТАНИЕ, ПРОСВЕЩЕНИЕ, ОБРАЗОВАНИЕ — Upbringing, Inlightening, Education ->
    [center][Youtube]69vJGqDENq4[/Youtube][/center]
    [center]14:36[/center]
    Osievskii Global News
    29 сент. Отправлено 05:20, 01.10.2020 г.’ target=_top>Просвещение от Вячеслава Осиевского — Карим_Хайдаров.
    30.09.2020 — 12:51: ВОСПИТАНИЕ, ПРОСВЕЩЕНИЕ, ОБРАЗОВАНИЕ — Upbringing, Inlightening, Education ->
    [center][Ok]376309070[/Ok][/center]
    [center]11:03[/center] Отправлено 12:51, 30.09.2020 г.’ target=_top>Просвещение от Дэйвида Дюка — Карим_Хайдаров.
    30.09.2020 — 11:53: ВОСПИТАНИЕ, ПРОСВЕЩЕНИЕ, ОБРАЗОВАНИЕ — Upbringing, Inlightening, Education ->
    [center][Youtube]VVQv1EzDTtY[/Youtube][/center]
    [center]10:43[/center]

    интервью Раввина Борода https://cursorinfo.co.il/all-news/rav.
    мой телеграмм https://t.me/peshekhonovandrei
    мой твиттер https://twitter.com/Andrey54708595
    мой инстаграм https://www.instagram.com/andreipeshekhonow/

    [b]Мой комментарий:
    Андрей спрашивает: Краснодарская синагога — это что, военный объект?
    — Да, военный, потому что имеет разрешение от Росатома на манипуляции с радиоактивными веществами, а также иными веществами, опасными в отношении массового поражения. Именно это было выявлено группой краснодарцев во главе с Мариной Мелиховой.

    [center][Youtube]CLegyQkMkyw[/Youtube][/center]
    [center]10:22 [/center]

    Доминико Риккарди: Россию ждёт страшное будущее (хотелки ЦРУ):
    https://tainy.net/22686-predskazaniya-dominika-rikardi-o-budushhem-rossii-sdelannye-v-2000-godu.html

    Завещание Алена Даллеса / Разработка ЦРУ (запрещено к ознакомлению Роскомнадзором = Жид-над-рус-надзором)
    https://av-inf.blogspot.com/2013/12/dalles.html

    [center][b]Сон разума народа России [/center]

    [center][Youtube]CLegyQkMkyw[/Youtube][/center]
    [center]10:22 [/center]

    Доминико Риккарди: Россию ждёт страшное будущее (хотелки ЦРУ):
    https://tainy.net/22686-predskazaniya-dominika-rikardi-o-budushhem-rossii-sdelannye-v-2000-godu.html

    Завещание Алена Даллеса / Разработка ЦРУ (запрещено к ознакомлению Роскомнадзором = Жид-над-рус-надзором)
    https://av-inf.blogspot.com/2013/12/dalles.html

    [center][b]Сон разума народа России [/center]

    Добавить комментарий