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

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

(2)

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

(3)

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

(4)

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

, (5)

Только функции sin и cos c частотами k/T будут взаимно ортогональны, а это является непременным условием преобразования Фурье. Набор первых функций разложения в ряд Фурье приведен на рисунке 1. При этом длительность функций совпадает с длительностью анализа T .

Рисунок 1. Функции разложения в ряд Фурье

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

Рисунок 2. Спектр функции x (t ) при анализе на ограниченном интервале времени

В данном случае формула вычисления прямого преобразования Фурье (4) преобразуется к следующему виду:

(6)

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

(7)

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

(8)

Теперь рассмотрим как изменились свойства дискретного преобразования Фурье (ДПФ) по сравнению с прямым преобразованием Фурье на ограниченном интервале времени. Когда мы рассматривали дискретизацию аналогового сигнала, мы выяснили, что спектр входного сигнала должен быть ограничен по частоте. Это требование ограничивает количество дискретных составляющих спектра сигнала. Первоначально может показаться, что мы можем ограничить спектр сигнала частотой f д /2, что соответствует количеству частотных составляющих K = N /2 . Однако это не так. Несмотря на то, что спектр сигнала для действительных отсчетов сигнала для положительных частот и отрицательных частот симметричен относительно 0, отрицательные частоты могут потребоваться для некоторых алгоритмов работы со спектрами, например, для . Еще больше отличие получается при выполнении дискретного преобразования Фурье над комплексными отсчетами входного сигнала. В результате для полного описания спектра цифрового сигнала требуется N частотных отсчетов (k = 0, ..., N/2 ).

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

Рассмотрим непрерывный сигнал конечной длительности с числом степеней свободы, равным Для этого сигнала можно записать разложение в ряд Котельникова:

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

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

Рассмотрим спектр который определяется выражением

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

Очевидно, справедливо и обратное соотношение

Применяя теорему о запаздывании, можно записать

Подставляя (3.2) в (3.1), получим окончательное выражение для спектра

Чтобы перейти к дискретному преобразованию Фурье, значения спектра в выражении (3.3) нужно вычислять не для всех значений частоты, а для дискретных (выборочных):

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

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

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

Рассмотрим значение определяемое формулой (3.4) для где целое число:

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

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

Спектральную плотность сигнала запишем в виде ряда Котельникова

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

Интеграл в выражении аналогичен вычисленному ранее интегралу (3.2). Пользуясь этой аналогией, запишем

Подставляя (3.6) в (3.5), получим выражение для временной функции

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

где А принимает значения от 0 до

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

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

Подставим это выражение дискретного сигнала в формулу дискретного преобразования Фурье (3.4)

Для сравнения найдем спектральную плотность исходного треугольного импульса:

Легко видеть, что дискретный спектр (3.11) неточно описывает спектральную плотность треугольного импульса (3.12). Значения несколько отличаются от соответствующих значений спектра треугольного импульса (рис. 3.1, б).

Теперь подставим дискретные значения спектра (3.11) в выражение для обратного дискретного преобразования Фурье (3.8):

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

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

Рис. 3.1. Дискретное преобразование Фурье дискретизированного треугольного импульса

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

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

При этом обратное дискретное преобразование Фурье имеет вид

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

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

Фурье (3.15). Получившуюся при этом последовательность обозначим

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

Внутренняя сумма выражения (3,16) равна нулю, если и равна если Следовательно, при т. е. числовые последовательности совпадают друг с другом. Таким образом, при последовательном применении к любой числовой последовательности прямого и обратного дискретного преобразования Фурье получают в результате ту же последовательность.

Проиллюстрируем это положение простейшими примерами.

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

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

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

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

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

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

Введение

На лабораторном занятии были изучены возможности по дискретному тригонометрическому преобразованию (ДТП) со следующих точек зрения:

1. Проверили свойство обратимости заданного ДТП.

2. Исследовали линейность предложенного ДТП.

3. Изучили особенности повтора спектра у проверяемого ДТП.

4. Определили наличие симметричного отражения спектра у ДТП, а именно

4.1. наличие центральной симметрии,

4.2. наличие осевой (вертикальной) симметрии.

5. Рассмотрели влияние фазовых сдвигов сигнала на результирующее ДТП.

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

7. Исследовали возможность фильтрации сигналов с помощью заданного ДТП.

8. Проверили экспериментально сохранение энергии исследуемым ДТП.

9. Обнаружили связь данного ДТП с дискретным преобразованием Фурье.

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

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

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

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

ДПФ служит для анализа периодических функций, и его можно получить исходя из теории рядов Фурье. Пусть x0(t) - непрерывная периодическая функция с периодом Р и частотой f0 = 1/P так что

Функцию x0(t) можно разложить в ряд Фурье:

где коэффициенты разложения Х0(n) заданы формулой

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

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

где n1 - целое значение n, задающее частоту f1.

На фиг. 1 показаны такой ограниченный спектр и колебание, которому он соответствует.

интервал дискретизации Т равен

так что число отсчетов на период будет

Фиг. 1. Периодическая функция x0(t) с ограниченной полосой частот и ее спектр X0(n).

1В результате дискретизации получаем периодическое, нормализованное относительно Т колебание вида

Это колебание определено на интервале, равном его периоду, т. е.

Поскольку x(t/T) – периодическая функция для расчета коэффициентов ряда Фурье используется соотношение (2)

(Замена Р на /V в делителе и пределах интегрирования соответствует переходу к нормализованной переменной.) Подставляя выражение (3), получаем

Известно, что

Окончательно с учетом того, что по определению

Соотношение, связывающее x(k) с Х(n), может быть получено непосредственно из формулы (1), если подставить t=kT и учесть, что при ограниченной ширине спектра функции x0(t) сумма содержит конечное число членов. Итак,

Следует заметить, что x(k) -периодическая функция, т. е.

и аналогично

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

Итак, при дискретизации периодической функции x0(t) соотношение (4) позволяет по выборкам x0(t) найти спектр Х(n), который на интервале 0 ≤ n ≤ N - 1 в точности равен спектру Х0(n) исходной периодической функции. Функция x(k) и ее спектр графически представлены на фиг. 2. Поскольку соотношение (5.4) получено на основании теоремы отсчетов, оно является точным и экономичным (при расчетах) эквивалентом исходного интегрального соотношения (2) и может быть использовано для вычисления коэффициентов разложения на ЭВМ. Соотношения (4) и (5) будем называть дискретным преобразованием Фурье (ДПФ) и обратным дискретным преобразованием Фурье (ОДПФ) соответственно. Заметим, что переменная n меняется здесь от нуля до N-1. Получаемый спектр можно интерпретировать следующим образом. Первые (N/2-1) точек Х(n) -соответствуют (N/2 - 1) спектральным линиям Х0(n) на положительных частотах, как показано на фиг. 5.3, а последние (N/2-1) точек Х(n) соответствуют (N/2-1) спектральным линиям на отрицательных частотах.

Пара преобразований, заданная соотношениями (4) и (5), встречается и в другом виде. Например, множитель 1 / N и знак минус у экспоненты могут быть записаны как в прямом, так и в обратном преобразовании, общий смысл при этом не меняется.

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

Фиг. 2. Дискретизированная периодическая функция x(k) и ее периодический спектр Х(n).

Фиг. 3. Соотношение между коэффициентами ряда Фурье и ДПФ.

Свойства ДПФ

Некоторые свойства ДПФ играют в практических вопросах обработки сигналов важную роль.

Линейность

Если xр(n) и ур(n) - периодические последовательности (с периодом в N отсчетов каждая), а Хр(k) и Yp(k) - их ДПФ, то дискретное преобразование Фурье последовательности хр(n) + + ур(n) равно Хр(k) + Yp(k). Это положение справедливо и для последовательностей конечной длины.

Сдвиг

Если последовательность хр(n) периодическая с периодом в N отсчетов, а ее ДПФ равно Хр(k), то ДПФ периодической последовательности вида хр(n-n0) будет равно.

Фиг. 4. К определению ДПФ сдвинутой последовательности.

При анализе последовательностей конечной длины необходимо учитывать специфический характер временного сдвига последовательности. Так, на фиг. 4, а изображена конечная последовательность х (п) длиной в N отсчетов. Там же крестиками изображены отсчеты эквивалентной периодической последовательности хр(n), имеющей то же ДПФ, что и х(n). Чтобы найти ДПФ сдвинутой последовательности х(n - n0), причем n0 < N, следует рассмотреть сдвинутую периодическую последовательность Хр(n - n0) и в качестве эквивалентной сдвинутой конечной последовательности (имеющей ДПФ j принять отрезок последовательности хр(n - n0) в интервале 0 ≤ n ≤ N - 1. Таким образом, с точки зрения ДПФ последовательность х(n – n0) получается путем кругового сдвига элементов последовательности х(n) на n0 отсчетов

Свойства симметрии

Если периодическая последовательность хр(n) с периодом в./V отсчетов является действительной, то ее ДПФ Хр(k) удовлетворяет следующим условиям симметрии:

Аналогичные равенства справедливы и для конечной действительной последовательности х(n), имеющей N-точечное ДПФ X(k). Если ввести дополнительное условие симметрии последовательности хp(n), т. е. считать, что

то окажется, что Хр(k) может быть только действительной.

Поскольку чаще всего приходится иметь дело с действительными последовательностями, то, вычислив одно ДПФ, можно получить ДПФ двух последовательностей, используя свойства симметрии (6). Рассмотрим действительные периодические последовательности хр(n) и ур(n) с периодами в N отсчетов и N-точечными ДПФ Хр(k) и Yp(k) соответственно. Введем комплексную последовательность zp(n) вида

Ее ДПФ равно

Выделяя действительную и мнимую части равенства (10), получим

Действительные части Хр(k) и Yp(k) симметричны, а мнимые - антисимметричны, поэтому их легко разделить, используя операции сложения и вычитания:

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

Похожая информация.

Исследуем особенности спектрального представления дискретного сигнала, который задан на отрезке своими отсчётами
, взятыми соответственно в моменты времени
; полное число отсчётов
(- интервал дискретизации).

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

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

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

(5.1)

Где
– выборочные значения аналогового сигнала.

- дискретное преобразование Фурье (ДПФ) (5.4)

Основные свойства ДПФ

1. ДПФ- линейное преобразование т.е. сумме сигналов отвечает сумма их ДПФ

2. Число различных коэффициентов
, вычисляемых по формуле (5.4), равно числу N за период; при коэффициент

3. Коэффициент (постоянная составляющая) является средним значением всех отсчётов:

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

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

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

(5.5)

Очевидно, что (5.5) представляет собой формулу обратного дискретного преобразования Фурье (ОДПФ) .

11 Алгоритм быстрого преобразования Фурье. Число вычислительных операций. Сравнение дискретного и быстрого преобразования Фурье.

=0, 1, 2,…,( /2)-1 (5.7)

Учтём, что последовательности коэффициентов, относящихся к чётной и нечётной частям входного массива, являются периодическими с периодом N/2:

Кроме того, входящий в формулу (5.7) множитель при
можно преобразовать так:

Отсюда находим выражение для второй половины множества коэффициентов ДПФ

(5.8)

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

Число операций, необходимых для вычисления БПФ оценивается как
.

Выигрыш в скорости вычислений по сравнению с традиционным ДПФ достигает сотен и даже тысяч при достаточных длинах входных массивов.

12.. Z - преобразование и его свойства. Применение Z - преобразования.

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

Пусть
– числовая последовательность, конечная или бесконечная, содержащая отсчётные значения некоторого сигнала. Поставим ей в однозначное соответствие сумму ряда по отрицательным степеням комплексной переменнойZ:

(5.9)

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

Данное выражение имеет смысл при |Z|> .

Обратное Z- преобразование

Пусть x(z) – функция комплексной переменной Z. Замечательное свойство Z-преобразование состоит в том, что функция x(z) определяет всю бесконечную совокупность отсчётов (
).

Действительно, умножим обе части ряда (5.9) на множитель
:

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

(5.11)

Данное выражение носит название обратное Z-преобразование.

Важнейшие свойства Z -преобразования:

1. Линейность . Если
и
- некоторые дискретные сигналы, причём известны соответствующиеZ-преобразования x(z) и y(z), то сигналу
будет отвечать преобразованиепри любых постоянныхи. Доказательство проводится путём подстановки суммы в формулу (5.9).

2. Z -преобразование смещённого сигнала . Рассмотрим дискретный сигнал
, получающийся из дискретного сигнала
путём сдвига на одну позицию в сторону запаздывания, т.е. когда
. Непосредственно вычисляяZ-преобразование, получаем следующий результат:

Таким образом, символ
служит оператором единичной задержки (на один интервал дискретизации) вZ-области.

3. Z -преобразование свёртки . Пусть x(z) и y(z) – непрерывные сигналы, для которых определена свёртка:

(5.13)

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

Подобную дискретную свёртку называют линейной

Вычислим Z-преобразование дискретной свёртки:

(5.15)

Итак, свёртке двух дискретных сигналов отвечает произведение Z-преобразований.

Программирование

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

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

d(X,Y) = SUM (X - Y)^2

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

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

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

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

• Пусть даны два изображения X и Y – изображение и образец, размеров (N1,N2) и (M1,M2) соответственно и Ni > Mi

Например, найти:

в изображении

Корреляция как мера между изображениями

Согласно определению, корреляцией двух функций F и G называется величина:

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

m(X,Y) = SUM (X * Y) / (SQRT (SUM X ^2) * SQRT (SUM Y ^2))

или

m(X,Y) = / (SQRT () * SQRT ())

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

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

Свёртка двух функций

Согласно определению, свёрткой двух функций F и G называется функция FхG:

Пусть G’(t) = G(-t) и F’(t) = F(-t), тогда, очевидна справедливость равенств:

• FхF’(0) = SUM F(i)^2 – скалярное произведение вектора F на самого себя
• GхG’(0) = SUM G(j)^2– скалярное произведение вектора G на самого себя
• FхG’(0) = SUM F(i)*G(i) – скалярное произведение двух векторов F и G

Так же очевидно, что FхG’(t) равна корреляции получаемой в результате сдвига одного вектора, относительно другого на шаг t (это легко проверить явной подстановкой значений в формулу корреляции).

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

Преобразование Фурье функции f вещественной переменной является интегральным и задаётся следующей формулой:

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

Кроме того, существуют разнообразные обобщения данного понятия.

Многомерное преобразование Фурье

Преобразование Фурье функций, заданных на пространстве ℝ^n, определяется формулой:

Обратное преобразование в этом случае задается формулой:

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

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

Дискретное преобразование Фурье (в англоязычной литературе DFT, Discrete Fourier Transform) - это одно из преобразований Фурье, широко применяемых в алгоритмах цифровой обработки сигналов (его модификации применяются в сжатии звука в MP3, сжатии изображений в JPEG и др.), а также в других областях, связанных с анализом частот в дискретном (к примеру, оцифрованном аналоговом) сигнале. Дискретное преобразование Фурье требует в качестве входа дискретную функцию. Такие функции часто создаются путём дискретизации (выборки значений из непрерывных функций). Дискретные преобразования Фурье помогают решать дифференциальные уравнения в частных производных и выполнять такие операции, как свёртки. Дискретные преобразования Фурье также активно используются в статистике, при анализе временных рядов. Существуют многомерные дискретные преобразования Фурье.

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

Прямое преобразование:

Обратное преобразование:

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

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

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

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

Где

• FFT – операция прямого преобразования Фурье
• BFT – операция обратного преобразования Фурье

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

Фурье-преобразования для вычисления корреляции

Пусть (t) равна корреляции получаемой в результате сдвига одного вектора, относительно другого на шаг t
Тогда, как уже показано ранее, выполняется

(t) = FхG’(t) = BFT (FFT(F)*FFT(G’))

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

FFT(G’) = CONJUGATE (FFT(G))

Где CONJUGATE (FFT(G)) – матрица, составленная из сопряжённых элементов матрицы FFT(G)
Таким образом, получаем

(t) = BFT (FFT(F)*CONJUGATE (FFT(G)))

Фурье-преобразования для решения задачи

m(X,Y) (i,j) = (i,j) / (|X|(i,j)) * |Y|(i,j)),

получаем, что

• = XxY’
• |X|^2 = = XxX’xE’ = BFT (FFT(X) * CONJUGATE (FFT(X)) * CONJUGATE (FFT(E)))
• |Y|^2 = = YxY’xE’ = BFT (FFT(Y) * CONJUGATE (FFT(Y)) * CONJUGATE (FFT(E)))

Где

• |X|(i,j) – норма общей части изображения X при сдвиге (i,j)
• |Y|(i,j) – норма общей части изображения Y при сдвиге (i,j)

Упрощение формул для решения поставленной задачи

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

m(X,Y) (i,j) = (i,j) / |X|^2(i,j),

получаем, что

• = BFT (FFT(X) * CONJUGATE (FFT(Y)))
• = BFT (SQUAREMAGNITUDE(FFT(X)) * CONJUGATE (FFT(E)))

Где

• (i,j) – скалярное произведение двух изображений, получаемых при сдвиге (i,j) относительно друг друга изображений X и Y
• E – изображение размера равному минимальным размерам X и Y, и заполненное единичными значениями (то есть “кадр” в котором сравниваются X и Y)
• (i,j) – норма (сумма яркостей пикселей) общей части изображения X при сдвиге (i,j)
• FFT – операция прямого двухмерного дискретного преобразования Фурье
• BFT – операция обратного двухмерного дискретного преобразования Фурье
• CONJUGATE – операция вычисления матрицы из сопряжённых элементов
• SQUAREMAGNITUDE– операция вычисления матрицы квадратов амплитуд элементов

Алгоритм поиска фрагмента в полном изображении

• Пусть даны два изображения X и Y – изображение и образец, размеров (N1,N2) и (M1,M2) соответственно и Ni > Mi
• Требуется найти координаты образца Y в полном изображении X и вычислить оценочную величину - меру близости.

1. Расширить изображение Y до размера (N1,N2), дополнив его нулями
2. Сформировать изображение E из единиц размера (M1,M2) и расширить до размера (N1,N2), дополнив его нулями
3. Вычислить = BFT (FFT(X) * CONJUGATE (FFT(Y)))
4. Вычислить = BFT (SQUAREMAGNITUDE(FFT(X)) * CONJUGATE (FFT(E)))
5. Вычислить M = (f + )/(f + )
6. В матрице M найти элемент с максимальным значением – координаты этого элемента и являются искомой позицией образца в полном изображении, а значение равно оценке меры сравнения.

Примечание:
При использовании дискретного преобразования Фурье, матрица M содержит также элементы от циклического сдвига изображений между собой. Поэтому, если не требуется анализировать циклический сдвиг кадров, то поиск максимального элемента в матрице M нужно ограничить областью (0,0)-(N1-M1, N2-M2).

Примеры реализации

Реализованные алгоритмы являются частью библиотеки с открытым исходным кодом FFTTools. Интернет-адрес: github.com/dprotopopov/FFTTools

Используемое программное обеспечение

• Microsoft Visual Studio 2013 C# - среда и язык программирования
• EmguCV/OpenCV – C++ библиотека структур и алгоритмов для обработки изображений
• FFTWSharp/FFTW – C++ библиотека реализующая алгоритмы быстрого дискретного преобразования Фурье

///

/// Matrix of values private Matrix Catch(Image image) { const double f = 1.0; int length = image.Data.Length; int n0 = image.Data.GetLength(0); int n1 = image.Data.GetLength(1); int n2 = image.Data.GetLength(2); Debug.Assert(n2 == 1); // Allocate FFTW structures var input = new fftw_complexarray(length); var output = new fftw_complexarray(length); fftw_plan forward = fftw_plan.dft_3d(n0, n1, n2, input, output, fftw_direction.Forward, fftw_flags.Estimate); fftw_plan backward = fftw_plan.dft_3d(n0, n1, n2, input, output, fftw_direction.Backward, fftw_flags.Estimate); var matrix = new Matrix(n0, n1); double[,] patternData = _patternImage.Data; double[,] imageData = image.Data; double[,] data = matrix.Data; var doubles = new double; // Calculate Divisor Copy(patternData, data); Buffer.BlockCopy(data, 0, doubles, 0, length*sizeof (double)); input.SetData(doubles.Select(x => new Complex(x, 0)).ToArray()); forward.Execute(); Complex complex = output.GetData_Complex(); Buffer.BlockCopy(imageData, 0, doubles, 0, length*sizeof (double)); input.SetData(doubles.Select(x => new Complex(x, 0)).ToArray()); forward.Execute(); input.SetData(output.GetData_Complex().Zip(complex, (x, y) => doubles1 = output.GetData_Complex().Select(x => x.Magnitude); if (_fastMode) { // Fast Result Buffer.BlockCopy(doubles1.ToArray(), 0, data, 0, length*sizeof (double)); return matrix; } // Calculate Divider (aka Power) input.SetData(doubles.Select(x => new Complex(x*x, 0)).ToArray()); forward.Execute(); complex = output.GetData_Complex(); CopyAndReplace(_patternImage.Data, data); Buffer.BlockCopy(data, 0, doubles, 0, length*sizeof (double)); input.SetData(doubles.Select(x => new Complex(x, 0)).ToArray()); forward.Execute(); input.SetData(complex.Zip(output.GetData_Complex(), (x, y) => x*Complex.Conjugate(y)).ToArray()); backward.Execute(); IEnumerable doubles2 = output.GetData_Complex().Select(x => x.Magnitude); // Result Buffer.BlockCopy(doubles1.Zip(doubles2, (x, y) => (f + x*x)/(f + y)).ToArray(), 0, data, 0, length*sizeof (double)); return matrix; }
///

/// Copy 3D array to 2D array (sizes can be different) /// Flip copied data /// Reduce last dimension ///

/// Input array /// Output array private static void Copy(double[,] input, double[,] output) { int n0 = output.GetLength(0); int n1 = output.GetLength(1); int m0 = Math.Min(n0, input.GetLength(0)); int m1 = Math.Min(n1, input.GetLength(1)); int m2 = input.GetLength(2); for (int i = 0; i < m0; i++) for (int j = 0; j < m1; j++) output = input; for (int k = 1; k < m2; k++) for (int i = 0; i < m0; i++) for (int j = 0; j < m1; j++) output += input; } ///

/// Copy 3D array to 2D array (sizes can be different) /// Replace items copied by value /// Flip copied data /// Reduce last dimension ///

/// Input array /// Output array /// Value to replace copied data private static void CopyAndReplace(double[,] input, double[,] output, double value = 1.0) { int n0 = output.GetLength(0); int n1 = output.GetLength(1); int m0 = Math.Min(n0, input.GetLength(0)); int m1 = Math.Min(n1, input.GetLength(1)); int m2 = input.GetLength(2); for (int i = 0; i < m0; i++) for (int j = 0; j < m1; j++) output = value; } ///

/// Find a maximum element in the matrix ///

/// Matrix of values /// Index of maximum element /// Index of maximum element /// Value of maximum element public void Max(Matrix matrix, out int x, out int y, out double value) { double[,] data = matrix.Data; int n0 = data.GetLength(0); int n1 = data.GetLength(1); value = data; x = y = 0; for (int i = 0; i < n0; i++) { for (int j = 0; j < n1; j++) { if (data < value) continue; value = data; x = j; y = i; } } } ///

/// Catch pattern bitmap with the Fastest Fourier Transform ///

/// Array of values public Matrix Catch(Bitmap bitmap) { using (var image = new Image(bitmap)) return Catch(image); }

Попался, который кусался

Литература

1. А.Л. Дмитриев. Оптические методы обработки информации. Учебное пособие. СПб. СПюГУИТМО 2005. 46 с.
2. А.А.Акаев, С.А.Майоров «Оптические методы обработки информации» М.:1988
3. Дж.Гудмен «Введение в Фурье-оптику» М.: Мир 1970