Разработка программно-математического обеспечения корреляционного совмещения изображений с использованием быстрого преобразования Фурье

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

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

Реферат

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

Результаты будут использоваться в НИР проводимой на кафедре ЭВМ ФГБОУ ВПО «РГРТУ».

The abstract

The degree project contains development of a program complex, designed for the study and conduct of direct and reverse Fourier transform. The software provides the ability to make image geometric transformation, find the Fourier transformation of images

Results will be used in research and development.

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

Определение местоположения летательного аппарата (ЛА) — это одна из первостепенных задач современной авионики. Эксплуатация ЛА невозможна без быстродействующей, надежной навигационной системы. Одним из наиболее динамично развивающихся и перспективных направлений в данной области являются корреляционно-экстремальные навигационные системы (КЭНС).

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

Основой работы корреляционно-экстремальных навигационных систем является сравнение изображения совокупности ориентиров (текущего изображения) с эталонным изображением, полученным ранее. Разница в положении этих изображений в принятой системе координат позволяет формировать команды для удержания объекта управления на заданной траектории движения.

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

1. Технико-экономическое обоснование

Актуальность выбранной темы дипломного проекта обуславливается общими тенденциями широкого развития и внедрения корреляционно экстремальных навигационных систем в современных летательных аппаратах.

В настоящее время КЭНС применяются в качестве систем навигации (наведения) пилотируемых самолетов, дистанционно-пилотируемых летательных аппаратов, ракет. При работе КЭНС происходит сравнение ТИ подстилающей поверхности с ЭИ, хранящимся в бортовой базе данных. На этапе предполетной подготовки в бортовую базу данных помещаются ЭИ той местности, над которой будет происходить полет с учетом возможных отклонений от намеченного курса. Общая схема функционирования системы сравнения ТИ с ЭИ показана на рисунке 1.1.

Рисунок 1.1 — Общая схема функционирования системы сравнения ТИ с ЭИ

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

· оптическое

· радиолокационное

· радиотепловое

· магнитное

Для представления окружающей обстановки на борту ЛА наряду с другими системами могут присутствовать подсистемы технического зрения (СТЗ), такие как:

· бортовая радиолокационная станция

· телекамера

· тепловизионный датчик

· лазерный локатор

В автоматизированных системах снижаются и требования к «разрешению» распознающих алгоритмов. Системе информационной поддержки достаточно привлечь внимание оператора к определённому участку сцены, после чего распознавание точного типа объектов и принятие решения о необходимости тех или иных действий осуществит сам оператор. При такой постановке задачи нет необходимости поддерживать сверхподробную базу моделей возможных целей. База моделей может включать лишь общее описание крупных классов целей. В то же время уменьшение подробности выдаваемых оператору «подсказок» позволяет резко увеличить скорость обработки информации, что ведёт к высвобождению вычислительных ресурсов для решения других задач управления ЛА.

Задача автоматического или автоматизированного обнаружения ориентиров и навигации над местностью является базовой во всём комплексе задач машинного зрения в перспективных ЛА. Указанные задачи могут быть сформулированы следующим образом:

· обнаружение объектов и изменений в сцене наблюдения;

· высокоточные измерения элементов сцены;

· слежение за объектами;

· самоориентация и самопозиционирование ЛА;

· реконструкция наблюдаемых поверхностей и обнаружение трёхмерных структур;

· описание сцены и идентификация объектов.

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

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

К достоинствам ГИС относятся:

· Удобные инструменты визуализации данных

· Наиболее естественное отображение пространственной информации

· Мощные возможности пространственного моделирования

· Полноценная работа со стандартными СУБД и др.

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

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

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

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

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

Для решения данной задачи необходимо:

1. Разработать программно — математическое обеспечение, позволяющее произвести разложение изображения в Фурье-образ, используя быстрое преобразование Фурье.

2. Разработать программно — математическое обеспечение, позволяющее произвести обратное Фурье-преобразование, получить изображение из его Фурье-образа.

3. Произвести анализ результатов, полученных с помощью разработанного программно — математического обеспечения.

Основным назначением разрабатываемого программного стенда является реализация разложения изображения в Фурье-образ с помощью быстрого преобразования Фурье с проверкой правильности разработанного программного продукта на реальных образцах РЛИ.

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

1) Программный стенд должен обеспечить:

· Выбор и загрузку исходного ТИ.

· Нахождение Фурье-образа загруженного исходного ТИ.

· Восстановление исходного изображения по полученному Фурье-образу.

· Визуальное отображение найденного Фурье-образа исходного ТИ в виде пары изображений.

В качестве исходного РЛИ должно использоваться растровое изображение в формате битовой карты (BMP) с 256 градациями яркости серого цвета.

Программный стенд должен быть разработан на языке C++ в среде программирования C++ Builder версии 6.0 для операционной системы Windows XP (Windows 7), иметь удобный пользовательский интерфейс.

2. Анализ проблем, возникающих при совмещении изображений в корреляционно-экстремальных навигационных системах

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

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

Перечислим эти факторы более подробно:

*шумовые эффекты — имеют десятки видов источников возникновения, к числу которых можно отнести несовершенство сенсоров приёмо-передающей аппаратуры, аппаратуры оцифровки изображений, трудные условия съёмки, недостаток освещения и ряд других;

*сложный текстурированный фон, на котором должно происходить обнаружение объектов;

*эффекты загораживания (заслонения) одних объектов другими, как правило, не определённой заранее формы, например — облако на космофотоснимке и т. п., загораживающие помехи;

*искажающие оптические эффекты в виде различных расфокусировок и дисторсий, ракурсные искажения и др.;

*эффекты резкой смены освещения, блики, тени, особенно в динамически меняющихся сценах;

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

*эффекты изменения среды между сенсорами и объектами наблюдения — задымления, атмосферные осадки, пыль, искусственные помехи и многие другие;

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

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

Локализация.

Важное отличие, присущее собственно проблеме обнаружения объектов на изображениях по сравнению с задачами распознавания или интерпретации заранее сегментированного образа, заключается в том, что обнаружение в практических задачах всегда связано с процедурой поиска объекта. Именно реализация процедуры поиска объекта связана с угрозой взрывообразного роста необходимого объёма вычислений.

Проиллюстрируем это на примере простой задачи поиска объекта путём сравнения текущего изображения сцены с растровым эталоном или шаблоном формы объекта. Если построить какой-либо функционал соответствия между объектом размером M? M и фрагментом M? M из изображения размера N? N, то простой перебор фрагментов требует количества вычислений не менее чем операций, что, например, при размере объекта 50?50, а изображения — 2000?2000 элементов составляет 10 миллиардов операций. Даже с учётом значительного увеличения мощности современных БЦВМ (бортовая центральная вычислительная машина), такие объёмы вычислений по-прежнему далеко выходят за пределы возможностей реализации бортовых систем реального времени, предназначенных для таких задач как навигация и наведение ЛА.

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

Это понятие означает, что необходимо не только обнаружить объект, но и точно указать в системе координат изображения (или сцены) его положение в каком-либо смысле. Несколько неясное толкование понятия «локализации», приведенное выше, связано с тем, что по сравнению со своей эталонной моделью объект на реальном изображении может быть заметно искажён геометрически, причём аналитическая модель искажения может отсутствовать. В этих случаях локализация объекта является нетривиальной задачей. В более простой ситуации, при аналитически известной с точностью до параметров геометрии искажений, под точной локализацией можно понимать знание о положении какой-либо характерной точки объекта и параметрах геометрии искажения (углы поворота, элементы проективного преобразования и т. п.). При этом встречающиеся случаи ошибок локализации целесообразно разделить на две группы — нормальные и аномальные ошибки.

Нормальная ошибка — это правильная локализация объекта с некоторой позиционной или параметрической неточностью, характеризуемой количественными оценками. Для объектов, характеризуемых габаритными размерами, большими 3?3…5?5 элементов изображения, позиционные нормальные ошибки могут быть значительно меньше размера элемента изображения, уменьшаясь с величиной объекта. В этом случае принято говорить о возможности субпиксельной локализации. Это особенно важно для задач стереообнаружения, так как при малых параллаксах 3D-объектов субпиксельная локализация самым существенным образом определяет точность их пространственного положения.

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

Вычислительная реализуемость.

Несмотря на отмеченный ранее колоссальный прогресс вычислительной техники и создание обширной специализированной процессорной базы для обработки изображений, для основной массы бортовых приложений реального времени характеристики вычислителей и их свойства всё ещё далеки от желаемых. Даже в случае реализации простейших алгоритмов оконной фильтрации изображения с минимальной апертурой 3?3 элемента объём вычислений составляет десятки операций на точку изображения. При обработке более высокого уровня необходимый объём вычислений колеблется в пределах от сотен до тысяч операций на пиксел.

Если размер анализируемого изображения составляет порядка 1000?1000 элементов, что не является чем-либо необычным для современных видео датчиков (можно вспомнить, что бытовые цифровые фотоаппараты давно превзошли отметку 5 Мпикс. в ПЗС-матрице (специализированная аналоговая интегральная микросхема, состоящая из светочувствительных фотодиодов, выполненная на основе кремния, использующая технологию ПЗС — приборов с зарядовой связью)), мы получим оценку количества потребных вычислений порядка нескольких гигабайтов операций на кадр.

Между тем, для приложений реального времени необходимо выполнять эти вычисления в темпе кадровой развертки (не менее 25 кадров в секунду), что приводит к оценке быстродействия около 50 Гфлопс (флопс — внесистемная единица, используемая для измерения производительности компьютеров, показывающая, сколько операций с плавающей точкой в секунду выполняет данная вычислительная система). Сами по себе эти оценки сегодня не являются запредельными для ЭВМ последнего поколения, однако следует учесть, что в случае создания систем управления перспективных ЛА массогабаритные характеристики конструируемых вычислительных устройств должны быть весьма ограничены. Таким образом, вычислительная реализуемость алгоритмов по-прежнему относится к числу наиболее важных факторов, учитываемых при их разработке.

Исходя из названных выше требованиям к алгоритмам сравнения изображений и накладываемых на них ограничений из-за не идеальности условий полета и возникающих помех, одной из важнейших задач, решаемых КЭНС, становится приведение ТИ и ЭИ к максимально сравнимому виду. Дополнительная обработка изображений перед их сравнением позволяет существенно снизить вероятность возникновения ошибок в определении местонахождения как самого летательного аппарата, так и объектов на сцене наблюдения.

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

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

Первая программная реализация алгоритма БПФ была осуществлена в начале 60-х годов XX века Джоном Кули в вычислительном центре IBM под руководством тески Джона Тьюки, а в 1965 году ими же была опубликована статья, посвященная алгоритму быстрого преобразования Фурье. С этого момента начинается настоящая БПФ-мания. Публикуются тысячи работ посвященных алгоритму БПФ, одна за одной выходят монографии, программисты соревнуются в эффективности реализации алгоритма. БПФ становится основным инструментом спектрального анализа сигналов.

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

Рассмотрим физический смысл дискретного преобразования Фурье. Пусть есть функция синуса x = sin (t).

Рисунок 3.1 — График функции x = sin (t)

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

Частота колебания обратна периоду: н = 1/T. Также говорят о круговой частоте, которая вычисляется по формуле: щ= 2рн = 2рT. Откуда: x = A sin (щt).

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

Рисунок 3.2 — График гармонического колебания x=Asin (2рt/T+ц) Очень похоже выглядит и выражение гармоники через косинус:

Рисунок 3.3 — График гармонического колебания x=Acos (2рt/T+ц) Принципиальной разницы в приведенных представлениях нет. Достаточно изменить фазу на р/2, чтобы перейти от синуса к косинусу и обратно. Далее будем подразумевать под гармоникой функцию косинуса:

x = A cos (2рt/T + ц) = A cos (2рнt + ц) = A cos (щt + ц) (3.1)

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

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

x = A cos ц cos (2рt/T) — A sin ц sin (2рt/T) (3.2)

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

x = Re cos (2рt/T) — Im sin (2рt / T) (3.3)

Re = A cos ц, Im = A sin ц

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

и (3.4)

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

(3.5)

Раскладывая каждое комплексное Xk на мнимую и действительную составляющие Xk = Rek + j Imk; разкладывая экспоненту по формуле Эйлера на синус и косинус действительного аргумента; перемножая полученные выражения; внеся 1/N под знак суммы и перегруппировав элементы в две суммы, получаем:

(3.6)

Рассмотрим следующую ситуацию. Пусть у нас есть звуковое или какое-либо иное колебание в виде функции x = f (t). Пусть это колебание задано в виде графика для отрезка времени [0, T]. Для обработки средствами вычислительной техники необходимо выполнить дискретизацию. Отрезок делится на N-1 частей и сохраняются значения функции x0, x1, x2,…, xN для N точек на границах отрезков t0 = 0, t1 = T/N, t2 = 2T/N,…, tn =nT/N,…, tN = T.

Рисунок 3.4 — Дискретизация непрерывного сигнала В результате прямого дискретного преобразования Фурье были получены N значений для Xk:

(3.7)

Если применить обратное дискретное преобразование Фурье, то получится исходная последовательность {x}. Исходная последовательность состояла из действительных чисел, а последовательность {X} в общем случае комплексная.

Вернемся к рассмотрению формулы (3.6). В левой части находится действительное число xn, а справа — две суммы, одна из которых помножена на мнимую единицу j. Сами же суммы состоят из действительных слагаемых. Отсюда следует, что вторая сумма равна нулю, если исходная последовательность {x} была действительной. Отбрасывая её, получаем:

(3.8)

Поскольку при дискретизации были выбраны tn = nT/N и xn = f (tn), то можно выполнить замену: n = tnN/T. В результате получим:

(3.9)

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

x = A cos (2рt/T + ц) = A cos (2рнt + ц) = A cos (щt + ц) (3.1)

x = Re cos (2рt/T) — Im sin (2рt / T) (3.3)

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

(3.10)

Функцию

Gk (t) = Ak cos (2рtk/T + цk) (3.11)

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

Амплитуда, фаза, частота и период каждой из гармоник связаны с коэффициентами Xk формулами:

(3.12)

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

3.2 Использование двумерного дискретного преобразования Фурье для обработки изображений

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

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

(3.13)

где — функция, описывающая исходное прямоугольное изображение, состоящее из N строк и M столбцов;

— Фурье_образ изображения ;

— мнимая единица;

.

Зависимость (4.1) может быть представлена следующим образом:

(3.14)

где, .

Число является комплексным, поэтому с учетом равенства зависимости (3.13) и (3.14) соответственно примут вид:

(3.15)

(3.16)

Из выражения (3.15) можно выразить коэффициенты и, определяющие действительную и мнимую части числа :

(3.17)

. (3.18)

3.3 Повышение быстродействия преобразования Фурье на основе быстрого двумерного преобразования Фурье

Использование алгоритма дискретного преобразования Фурье является не очень практичным вследствие больших затрат времени на его реализацию. Поэтому на практике используют алгоритмы быстрого преобразования Фурье (БПФ). В БПФ обычно число данных ограничено и выражено степенью с основанием 2 (2, 4, 8, 16, 32, 64, 128, …). Но, несмотря на это ограничение, БПФ всё равно используется благодаря практичности и высокой скорости вычислений.

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

При вычислении ДПФ для значений необходимо умножить раз и сложить раз. Если невелико, как в приведенном примере, то объем вычислений тоже мал, но если, например, равно 1000, то число операций достигает 1 000 000, что крайне затрудняет аппаратную реализацию алгоритма.

Алгоритм вычисления БПФ называется методом вычисления «бабочкой». Так, для одномерной 4_точечной функции он выглядит следующим образом:

(3.19)

(3.20)

(3.21)

(3.22)

где — комплексный член дискретного ряда Фурье;

;

— число точек функции .

Вычисление совокупности зависимостей (3.19) — (3.21) осуществляется в 2 этапа. На 1_м этапе осуществляется нахождение следующих коэффициентов:

(3.23)

(3.24)

(3.25)

. (3.26)

На 2_м этапе на основе значений () вычисляются коэффициенты () 4_точечного БПФ:

.

Одним из главных пунктов в алгоритме быстрого преобразования Фурье является метод вычисления «бабочкой». Еще один важный момент заключается в последовательных разбиениях ряда значений сигнала на две группы и перестановке значений сигнала таким образом, чтобы в последующем перейти к методу вычисления «бабочкой». Способ перестановки значений функции называется техникой сортировки. Техника сортировки основана на перестановке разрядов. Ряд значений функции расстанавливается в порядке, ,, (в случае БПФ из 4_х членов) и в порядке, ,, ,, ,, (в случае БПФ из 8_и членов). Таким образом, значение индекса получается из исходного перестановкой старших и младших разрядов. Это правило является универсальным для любого числа членов ряда.

Например, двумерное 32×32_точечное БПФ вычисляется путем вычисления коэффициентов для каждой строки изображения за счет последующего их умножения на величину :

(3.27)

(3.28)

3.4 Определение числа коэффициентов, необходимых для анализа изображения

Для анализа изображений по их Фурье_образам необходимо определить достаточное для анализа количество коэффициентов. Это позволит сократить время анализа. Ниже приводятся примеры описания простых изображений посредством коэффициентов ряда Фурье.

1) Описание одиночной точки на черном фоне Воспользовавшись зависимостями (3.17) и (3.18) можно получить следующие выражения:, (3.29)

(3.30)

(3.31)

(3.32)

. (3.33)

Из совокупности зависимостей (3.29) — (3.33) можно выразить координаты x и y, описывающие одиночную точку на черном фоне:

(3.34)

. (3.35)

Как видно из зависимостей (3.34) — (3.35), одиночная точка на черном фоне описывается 3_мя коэффициентами ряда Фурье. Из этого следует, что для уменьшения числа коэффициентов ряда Фурье, необходимых для анализа изображения, следует проводить контрастирование исходного изображения.

2) Описание линии по горизонтали, состоящей из двух точек, на черном фоне Такое изображение согласно (3.5) и (3.6) имеет следующее описание:

(3.36)

(3.37)

(3.38)

(3.39)

(3.40)

где и — координаты точек линии.

Из выражений (3.24) — (3.28) следует, что координаты искомого объекта могут быть найдены как:

(3.41)

. (3.42)

Как видно из совокупности (3.41) — (3.42) линия по горизонтали, состоящая из двух точек, описывается четырьмя коэффициентами ряда Фурье.

3.5 Восстановление изображений на основе Фурье-образов

Для проверки правильности нахождения Фурье_образа необходимо осуществить восстановление исходного изображения путем обратного преобразования Фурье, которое осуществляется согласно следующей зависимости:

(3.43)

3.6 Нахождение корреляционной функции радиолокационного и моделируемого изображений на основе их Фурье образов

Корреляционный анализ радиолокационного и моделируемого изображений целесообразно осуществлять на основе их Фурье_образов. Так корреляционная функция Фурье_образов радиолокационного и моделируемого изображений примет вид:

(3.44)

(3.45)

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

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

и — коэффициенты ряда Фурье для радиолокационного изображения;

и — коэффициенты ряда Фурье для моделируемого изображения;

— размер сравниваемых изображений.

Алгоритм нахождения корреляционной функции при совмещении радиолокационного и моделируемого изображений, основанный на вычислении зависимостей (3.44) и (3.45) представлен на рисунке 3.5.

Рисунок 3.5 — Алгоритм нахождения корреляционной функции при совмещении радиолокационного и моделируемого изображений

3.7 Алгоритмы быстрого преобразования Фурье FFT (fast Fourier transform). Принцип построения

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

(3.46)

ДПФ N отсчетам сигнала s (n), n=0.N-1(в общем случае комплексным) ставит в соответствие N комплексных отсчетов спектра S (k), k=0.N-1, причем для вычисления одного спектрального отсчета требуется N операций комплексного умножения и сложения. Таким образом вычислительная сложность алгоритма ДПФ составляет N2 комплексных умножений и сложений. При этом можно заметить что если одно ДПФ на N точек (отсчетов) заменить вычислением двух ДПФ по N/2 точек, то это приведет к уменьшению количества операций в 2 раза. Замена N-точечного ДПФ двумя N/2 точечными представлено на рисунке 3.6.

Рисунок 3.5 — Замена N-точечного ДПФ двумя N/2-точечными ДПФ

При этом каждое из N/2 — точечных ДПФ также можно вычислить путем замены N/2 — точечного ДПФ на два N/4-точечных. В этом случае количество вычислительных операций равно 4*(N2/16)= N2/4. Таким образом, можно продолжать разбиение исходной последовательности до тех пор, пока возможно деление последовательности на две. Для N=8 (L=3) такое разбиение представлено на рисунке 3.6.

Рисунок 3.6 — Разбиение и объединение последовательностей при N = 8

Каждое разбиение делит последовательность на две половинной длительности (красную и синюю), а каждое объединение «собирает «из двух последовательностей одну удвоенной длительности.

Алгоритмы БПФ, которые используют выборки длиной N=2L, называются «алгоритмами БПФ по основанию 2». Данные алгоритмы получили наибольшее распространение, из-за того что в машинной арифметике N=2L является «круглым» числом. Далее мы будем рассматривать только алгоритмы по основанию 2.

Очевидно что делить последовательности на две можно по-разному, однако от этого зависит сможем ли мы при объединении получить неискаженный спектр сигнала и чего с точки зрения вычислительных затрат это будет нам стоить. Можно сказать, что эффективность алгоритма БПФ полностью зависит от способа разбиения и объединения последовательности, поскольку если не учитывать операции на разбиение-объединение, то для расчета спектра требуется N/2 раз посчитать ДПФ на 2 точки, в результате общее количество вычислительных операций составит 2*N, то есть количество операций линейно зависит от величины выборки.

Мы рассмотрим два способа разбиения — объединения: прореживание по времени и прореживание по частоте.

Перед рассмотрением способов разбиения — объединения требуется рассмотреть обратное дискретное преобразование Фурье (ОДПФ):

(3.47)

которое ставит в соответствие N комплексным отсчетам спектра S (k), k=0.N-1 N комплексных значений сигнала s (n), k=0.N-1.

Имея алгоритм вычисления БПФ, логично использовать его и для обратного преобразования. Для этого обратим внимание на то, что:

(3.48)

Другими словами комплексные экспоненты в выражении для прямого и обратного ДПФ являются комплексно-сопряженными (черта сверху означает комплексное сопряжение). Теперь рассмотрим два комплексных числа x=a+jb и y=c+jd.

Рассмотрим произведение x на комплексно-сопряженное y:

(3.49)

Теперь рассмотрим произведение комплексно-сопряженного на y:

(3.50)

Сравнивая (3.49) и (3.50) можно сделать вывод:. (3.51)

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

(3.52)

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

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

3.8 БПФ по основанию 2 с прореживанием по времени

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

Выражение для ДПФ имеет вид:

(3.53)

Процедура разбиения-объединения представлена на рисунке 3.8.

Рисунок 3.8 — Разбиение и объединение последовательностей при N = 8

Вначале комплексную экспоненту в выражении (3.53) обозначим как:

(3.54)

Тогда выражение (3.53) принимает вид:

(3.55)

Прореживание по времени заключается в разбиении исходной последовательности отсчетов s (n), n=0.N-1 на две последовательности длительности N/2 s0(n) и s1(n), n=0.N-1, таких что s0(n)=s (2n), а s1(n)=s (2n+1), n=0.N/2−1. Другими словами, последовательность s0(n) содержит отсчеты последовательности s (n) с четными индексами, а s1(n) — с нечетными. Прореживание по времени для N = 8 наглядно представлено на рисунке 3.9.

Рисунок 3.9 — Прореживание по времени для N=8

Рассмотрим ДПФ сигнала прореженного по времени:

(3.56)

Если рассмотреть только первую половину спектра S (k), k=0.N-1, а также учесть что

(3.57)

тогда (3.8.4) можно записать:

(3.58)

где S0(k) и S1(k) — точечные ДПФ последовательностей «четной» s0(n) и «нечетной» s1(n), n=0.N/2−1

(3.59)

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

Процедура объединения. Граф «Бабочка»

Теперь рассмотрим вторую половину спектра S (k+N/2), k=0.N/2−1:

(3.60)

Рассмотрим подробнее множитель.

(3.61)

Учтем что,

(3.62)

тогда выражение (3.60) справедливо для любого целого n, тогда (3.59) можно записать :.

(3.63)

Рассмотрим теперь поворотный коэффициент в (3.8.8):

. (3.64)

Тогда выражение (3.60) с учетом (3.63) и (3.64) принимает вид:

(3.65)

Таким образом окончательно можно записать:

(3.66)

Выражение (3.66) представляет собой алгоритм объединения при прореживании по времени. Данную процедуру объединения можно представить в виде графа (рисунок 3.10), который называется «Бабочка».

Рисунок 3.10 — Процедура объединения на основе графа «Бабочка «

Представим в виде графа алгоритм БПФ с прореживанием по времени основанный на разбиении — объединении при (рисунок 3.11).

Рисунок 3.11 — Граф алгоритма БПФ с прореживанием по времени при N=8

На первом этапе отсчеты входного сигнала переставляются местами и исходная последовательность делится на «четную» и «нечетную последовательности» (обозначены красными и синими стрелками). Потом «четная» и «нечетная» последовательности в свою очередь делятся на «четную» и «нечетную» последовательности. При N=2L такое деление можно делать (L-1) раз. В нашем случае L=3. Данная процедура называется двоично-инверсной перестановкой, так можно выполнить перенумерацию отсчетов переписав номер отсчета в двоичной системе в обратном направлении. Например s (4) имеет индекс в десятичной системе счисления 410=1002, если же 1002 переписать справа налево то получим 0012, то есть s (4) после разбиения на четные нечетные перед первой операцией «Бабочка» встанет на место s (1), которая в свою очередь встанет на место s (4). По аналогичному правилу поменяются местами все отсчеты, при этом некоторые останутся на месте, в частности s (2), так как если 210=0102 переписать справа налево то все равно останется 0102, аналогично s (0), s (5) и s (7). Очень важно понять, что данный метод перенумерации должен применяться при записи числа в двоичной системе состоящей из L разрядов.

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

После двоично-инверсной перестановки получаем четыре 2-точечных ДПФ:

(3.67)

На основе четырех 2-точечных ДПФ формируются два 4-точечных ДПФ:

(3.68)

И на последнем уровне формируется полный спектр входного сигнала.

Поворотные коэффициенты Рассмотрим подробнее поворотные коэффициенты .

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

На втором уровне имеем два поворотных коэффициента: ;

(3.69)

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

На третьем уровне имеем уже 4 поворотных коэффициента:

(3.70)

Графически поворотные коэффициенты можно представить как векторы на комплексной плоскости (рисунок 3.12):

Рисунок 3.12 — Поворотные коэффициенты алгоритма БПФ с прореживанием по времени при N=8

Можно заметить что на всех уровнях объединения количество поворотных коэффициентов удваивается, причем все поворотные коэффициенты предыдущего уровня объединения присутствуют и на следующем уровне. Таким образом для того чтобы перейти на следующий уровень необходимо между поворотными коэффициентами текущего уровня вставить поворотные коэффициенты следующего. Графически для перехода на следующий уровень при N =16 необходимо дополнить рисунок 3.12 как это показано на рисунке 3.13. Серые вектора показывают поворотные коэффициенты, которые будут присутствовать на последнем уровне при N=16, которых нет при N=8 .

Рисунок 3.13 — Поворотные коэффициенты алгоритма БПФ с прореживанием по времени при N=16

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

Вычислительная эффективность алгоритма БПФ с прореживанием по времени Алгоритм с прореживанием по времени на каждом уровне требует N комплексных умножений и сложений. При N=2L количество уровней разложения — объединения равно L, таким образом общее количество операций умножения и сложения равно LN.

Рассмотрим во сколько раз алгоритм БПФ с прореживанием по времени эффективнее ДПФ. Для этого рассмотрим коэффициент ускорения K отношение количества комплексных умножение и сложений при использовании ДПФ и БПФ, при этом учтем что L=log2(N)

раз. (3.71)

В таблице ниже приведено требуемое количество операций Mоп для алгоритма ДПФ при различном N=2L и при использовании БПФ с прореживанием по времени.

Таблица 3.1. Эффективность БПФ с прореживанием по времени

Из таблицы хорошо видно, что использование БПФ приводит к существенному уменьшению требуемого количества вычислительных операций. Так например при Т=1024 БПФ требует в 100 раз меньше операций чем ДПФ, а при Т=4096 в 341 раз. Первое опубликование данных результатов произвело эффект разорвавшийся бомбы, всколыхнув всю научную общественность. При этом очень важно, что выигрыш по производительности тем больше, чем больше размер выборки N. Так например при N=216 выигрыш составит216/16=212=4096 раз.

Таким образом, для получения эффективного алгоритма БПФ по основанию 2 с прореживанием по времени, при N=2L необходимо:

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

Используя поворотные коэффициенты сделать N/2 операций «Бабочка» согласно выражению (3.68) для получения первого объединения Повторить операцию «Бабочка» согласно выражению (3.68) для объединения на следующий уровень еще L-1 раз, также используя поворотные коэффициенты.

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

3.9 Алгоритм БПФ с прореживанием по частоте

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

(3.72)

В алгоритме БПФ с прореживанием по времени производилось разделение исходного сигнала в соответствии с двоично-инверсной перестановкой. В алгоритме с прореживанием по частоте наоборот исходный сигнал s (n) делится пополам, т. е. s0(n)=s (n) и s1(n)=s (n+N/2), n=0.N/2−1. Тогда выражение (3.9.1) можно переписать:

(3.73)

Учитывая, что:

(3.74)

(3.75)

Рассмотрим теперь четные отсчеты спектра S (2k), K=0.N/2

(3.76)

Учитывая, что

(3.77)

Тогда

(3.78)

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

Рассмотрим теперь нечетные отсчеты спектра S (2k+1), k=0.N/2

(3.79)

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

(3.80)

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

Граф бабочка для алгоритма с прореживанием по частоте представлен на рисунке 3.14:

Рисунок 3.14 — Граф бабочка для алгоритма БПФ с прореживанием по частоте Поворотные коэффициенты в алгоритме с прореживанием по частоте полностью совпадают с поворотными коэффициентами алгоритма БПФ с прореживанием по времени.

Представим в виде графа алгоритм БПФ с прореживанием по частоте основанный на разбиении — объединении при N=8 (рисунок 3.15).

Рисунок 3.15 — Граф алгоритма БПФ с прореживанием по частоте для N=8

На первом этапе исходный сигнал делится на 2 половины (красные и синие стрелочки). Далее вычисляются

(3.81)

Тогда если выполнить ДПФ S0(n), то получим четные отсчеты спектра в соответствии с (3.9.9), а если ДПФ S1(n) — то нечетные отсчеты спектра. Таким образом одно ДПФ длительности N=8 заменили двумя ДПФ длительности N/2=4. Для вычисления каждой из ДПФ половинной длительности снова применим прореживание по частоте. В результате получим:

(3.82)

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

Сравнение алгоритмов БПФ по основанию 2 с прореживанием по времени и частоте Можно произвести сравнение алгоритма БПФ с прореживанием по частоте с алгоритмом БПФ с прореживанием по времени:

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

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

3. В связи с вышесказанным, вычислительная эффективность обоих алгоритмов практически идентична.

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

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

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

2. Выигрыши, получаемые альтернативными алгоритмами БПФ, незначительные по сравнению с алгоритмами по основанию два и нивелируются экспоненциально растущими вычислительными мощностями процессоров.

3. Алгоритмы по основанию два прекрасно «распараллеливаются» при использовании жесткой логики.

Все это переводит альтернативные алгоритмы БПФ в разряд «экзотики», и на практике алгоритмы по основанию два являются оптимальным выбором.

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

4. Разработка программно-математического обеспечения

4.1 Разработка алгоритмов

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

· ввод исходного изображения для Фурье-преобразования;

· Фурье-преобразование изображения;

· выполнение геометрических искажений изображения;

4.1.1 Нахождение Фурье-образа изображения

В данной работе был применен алгоритм Фурье-преобразования, состоящий из следующих этапов:

1) Предварительная перестановка элементов;

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

На первом шаге четные элементы с номером n переместились в позицию n/2, а нечетные из позиции в позицию N/2+(n-1)/2. Где n=0,1,…, N-1. Таким образом, новая позиция вычисляется из старой позиции с помощью функции:

ror (n, N) = [n/2] + N{n/2}

Здесь [x] означает целую часть числа, а {x} - дробную.

В ассемблере эта операция называется циклическим сдвигом вправо (ror), если N — это степень двойки. Название операции происходит из того факта, что берется двоичное представление числа n, затем все биты, кроме младшего (самого правого) перемещаются на 1 позицию вправо. А младший бит перемещается на освободившееся место самого старшего (самого левого) бита.

Рисунок 4.1 — Схема циклического сдвига вправо Дальнейшие разбиения выполняются аналогично. На каждом следующем шаге количество последовательностей удваивается, а число элементов в каждой из них уменьшается вдвое. Операции ror подвергаются уже не все биты, а только несколько младших (правых). Старшие же j-1 битов остаются нетронутыми (зафиксированными), где j — номер шага:

Рисунок 4.2 — Схема разбиения на шаге j

На рисунке 4.3 проиллюстрирован второй этап вычисления ДПФ. Линиями сверху вниз показано использование элементов для вычисления значений новых элментов. Очень удобно то, что два элемента на определенных позициях в массиве дают два элемента на тех же местах. Только принадлежать они будут уже другим, более длинным массивам, размещенным на месте прежних, более коротких. Это позволяет обойтись одним массивом данных для исходных данных, результата и хранения промежуточных результатов для всех T итераций.