Обработка сигналов и изображений в пространственно частотной фильтрации

Размещено на http://www.allbest.ru/

Обработка сигналов и изображений в пространственно-частотной фильтрации



1. Многомерные сигналы и системы

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

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

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

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

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

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

s(x,y) = sin(x²+y²), -< x,y < (1.1)

График функции (в пределах одного периода) приведен на рис. 1.1.

Двумерный дискретный сигнал (цифровой массив) - это функция, определенная на совокупности пар числовых значений координат с определенным шагом дискретизации ?x и ?y. В общем случае, при различной физической размерности аргументов x и y, значения ?x и ?y не равны друг другу:



Элемент последовательности s_n,m представляет собой отсчет двумерной функции s в координатной точке (x=n?x,y=m?y), где значения x и y - независимые переменные (аргументы) функции. Для числовых массивов значения шага дискретизации по аргументам также могут приниматься равными 1 (независимо от размерности) и использоваться аргументация s(n,m) s_n,m. Результаты геофизических съемок какого-либо одного геофизического параметра по поверхности земли относятся к двумерным функциям: дискретным - если это отсчеты в отдельных точках по определенной координатной сети (x,y), или смешанным - если это непрерывная регистрация данных по профилям (например - мощности экспозиционной дозы гамма излучения горных пород при аэросъемке). Но в настоящее время геофизические съемки относятся даже не к двумерным, а к многомерным функциям, так как регистрируется, как правило, сразу несколько физических параметров геологических сред. Так, например, при спектрометрической съемке естественной радиоактивности горных пород регистрируется содержание в горных породах урана, тория и калия, в гравиразведке - трехкоординатный вектор силы тяжести, и т.п. Если на какой-либо площади проведена съемка нескольких видов геофизики, то их результаты также могут рассматриваться в совокупности, как многомерная функция физических параметров данной геологической среды.

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

Отметим некоторые двумерные последовательности (функции, сигналы), имеющие специальные названия.

Двумерный единичный импульс ?(n?x,m?y) = ?_n,m или единичный отсчет:

?_n,m = 1, при n = m = 0.

= 0, при n0, m0.

?_n,m = ?_n ?_m,

где ?_n, ?_m - одномерные единичные импульсы (импульсы Кронекера) по координатам n и m.

Произвольное расположение двумерного единичного импульса по координатам n₁, m₁ соответственно записывается в виде: ?((n-n₁)?x,(m-m₁)?y) = ?_n-n1,m-m1. Попутно напомним, что математическая запись импульса Кронекера обозначает не единичный отсчет, а функцию, определяющую место положения единичного отсчета и нулевые значения по остальным координатам (аргументам).

Двумерный линейный импульс представляет собой последовательность единичных отсчетов по одной координате: s(n,m) = ?(n) или s(n,m) = ?(m).

На рис. 1.3 приведены два двумерных линейных импульса, первый - по координате m = 0: s(n,m) = ?(m), и второй импульс по координате n = 2: s(n,m) = ?(n-2).

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



Двумерная единичная ступенька u(n,m), представленная на рис. 1.4, определяется выражением:

u(n,m) = 1, при n0 и m0,

= 0, в остальных случаях.

u(n,m) = u(n) u(m),

где u(j) представляют собой единичные ступеньки соответственно по координатам n и m: u(j)=1 при j0, u(j)=0 при j<0. Двумерная единичная ступенька отлична от нуля в одном квадранте (n,m)- плоскости.

Экспоненциальная последовательность: s(n,m) = aⁿb^m, -< n,m <, где а и b в общем случае комплексные числа. При а = exp(j?₁), b = exp(j?₂), |а|=1, |b|=1:

s(n,m) = exp(jn?₁+jm?₂) = cos(n?₁+m?₂)+jsin(n?₁+m?₂).

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

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

s(n,m) = s(n) s(m).

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



s(n,m) =s_in(n) s_im(m),

где N- число ненулевых строк или столбцов массива. В крайнем случае, для этого достаточно выразить s(n,m) в виде суммы отдельных строк:

s(n,m) = s(n,i) ?(m-i). (1.2)

Конечные последовательности. Важным классом сигналов являются последовательности конечной протяженности, для которых сигнал равен нулю вне определенной области, называемой опорной областью сигнала. На рис.1.5 условно представлена двумерная последовательность конечной протяженности, значения которой отличны от нулевых только внутри ограниченной прямоугольной области -3n2, -2m. Опорная область сигнала может быть произвольной формы и выходить за пределы сигнала, частично включая нулевые отсчеты. Отсчеты за пределами опорной области считаются равными нулю.

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

s(n,m+M) = s(n,m),

s(n+N,m) = s(n,m), (1.3)

обладает периодичностью в двух направлениях, по n и по m. Значения М и N называют интервалами периодичности сигнала соответственно по координатам m и n (горизонтальными и вертикальными интервалами периодичности). Прямоугольная форма области периода (пример на рис.1.6) наиболее удобна при обработке данных, но не является единственно возможной.

Для двумерных последовательностей условия (1.3) могут рассматриваться как частный случай общих условий периодичности:

s(n+N₁, m+M₁) = s(n,m), (1.4)

s(n+N₂, m+M₂) = s(n,m),

D = N₁M₂ - N₂M₁ 0.

Упорядоченные пары (N₁,M₁) и (N₂,M₂) представляют собой смещения от отсчетов одного периода к соответствующим отсчетам других периодов и могут рассматриваться как векторы N и M, которые образуют области периодов в форме параллелограмма.

Линейная независимость векторов обеспечивается при ненулевом определителе D, а количество отсчетов в пределах периода равно |D|. Пример периодической последовательности с векторами (4,4) и (3,-5) приведен на рис. 1.7.

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

s() = s(+), i = 1,2,3, ... ,P.



Столбцы векторов N_i образуют матрицу периодичности N размером P х P. Векторы периодичности матрицы линейно независимы при наличии у матрицы ненулевого определителя. Абсолютное значение определителя равно числу отсчетов в периоде. Последовательность s() прямоугольно периодична для случаев диагональной матрицы N. Если функция s() периодична с матрицей периодичности N, то для любого целочисленного вектора Р имеет место s(+) = s(), и матрица PN также будет матрицей периодичности для s(). Отсюда следует, что любая многомерная периодическая последовательность имеет не единственную матрицу периодичности.

2. Двумерные системы

Системы осуществляют преобразование сигналов. Формализованная система - это оператор (операция) отображения входного сигнала на выходной: z(x,y) = Т[s(x,y)].

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

z(n,m) = c s(n,m),

z(n,m) = s(n-N,m-M),

z(n,m) = s(n,m)+u(n,m).

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

s(n,m) =s(i,j) ?(n-i,m-j). (2.1)



Обобщением скалярного умножения является пространственное маскирование:

z(n,m) = c_n,m s(n,m). (2.2)

Правая часть равенства (2.2) представляет собой поэлементное произведение входного сигнала на совокупность чисел с_n,m.

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

z_n,m = (s_n,m)².

Линейные системы. Система считается линейной при выполнении двух условий:

1. Пропорциональное изменение входного сигнала вызывает пропорциональное изменение выходного сигнала.

2. Суммарный сигнал двух входных последовательностей дает суммарный сигнал двух соответствующих выходных последовательностей.

Другими словами, если оператор Т[s(x,y)] описывает линейную систему и имеет место z(x,y) = Т[s(x,y)], q(x,y) = Т[u(x,y)], то Т[as(x,y)+bu(x,y)] = az(x,y)+bq(x,y). Линейные системы подчиняются принципу суперпозиции сигналов.

В выражении (2.1) значения s(i,j) можно рассматривать как скалярные множители для соответствующих единичных импульсов. Применяя оператор преобразования Т[.] к левой и правой части (2.1), получаем:



Т[s(n,m)] = у(n,m) = s(i,j) T[??n-i,m-j)],

z(n,m) = s(i,j) h_ij(n,m), (2.3)

где h_ij(n,m) - отклик системы в точке (n,m) на единичный импульс в точке (i,j). Если импульсный отклик h_ij(n,m) определен для всех точек (i,j), то отклик системы на произвольный многомерный сигнал, как и для одномерных систем, находится с помощью суперпозиции.

Инвариантность к сдвигу. Система инвариантна к сдвигу, если сдвиг входной последовательности приводит к такому же сдвигу выходной последовательности:

Т[s(n-N,m-M)] = z(n-N,m-M).

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

В дальнейшем ограничимся рассмотрением только систем, широко распространенных при решении практических задач - линейных и инвариантных к сдвигу (ЛИС-системы).

Импульсный отклик на произвольно расположенный входной импульс, как следует из выражения (2.3), описывается выражением:

h_ij(n,m) = T[?(n-n_i,m-m_j)].

Для частного случая i = j = 0 имеем:

h_o(n,m) = T[?(n,m)].

Используя принцип инвариантности к сдвигу, получим:

h_ij(n,m) = h_o(n-i,m-j) = h(n-i,m-j), (2.4)

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

Двумерная свертка. Подставляя (2.4) в выражение (2.3), получаем:

z(n,m) = ?_i ?_j s(i,j) h(n-i,m-j). (.2.5)

Двумерная дискретная свертка (2.5), является аналогом одномерной дискретной свертки. При замене переменных n-i = k, m-j = l, получим:

z(n,m) = ?_k ?_l h(k,l) s(n-k,m-l), (2.5')

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

Для упрощения символьного аппарата двумерную свертку обозначают индексом (_**):

z(n,m) = h(k,l) _** s(n-k,m-l).



При обобщении этого выражения на многомерные системы, в векторной форме:

z()= h() _** s(-).

Разделимые системы. Если импульсный отклик системы может быть разделен:

h(k,l) = h(k) h(l), (.2.6)

то выражение (2.5') принимает вид:

z(n,m) = ?_k h(k) ?_l h(l) s(n-k,m-l), (2.7)

z(n,m) = ?_k h(k) g(n-k,m), g(n-k,m) = ?_l h(l) s(n,m-l).

Массив g(n,m) вычисляется одномерной сверткой столбцов массива s(n,m) при n = const (сечения массива по координатам n) с откликом h(l), с последующим вычислением выходного массива z(n,m) одномерной сверткой строк g(n,m) при m = const с откликом h(k). Результат не изменится, если сначала выполнять свертку по строкам, а затем по столбцам. Система с откликом вида (2.6) называется разделимой. Отметим, что в разделимой системе входной и выходной сигнал не обязаны быть разделимыми.

Аналогичные разделимые системы могут существовать и в многомерном варианте.

Устойчивость систем. Интерес для практики представляют только устойчивые системы, обеспечивающие определенный конечный результат системной операции на конечные входные сигналы. Необходимым и достаточным условием устойчивости системы является абсолютная суммируемость ее импульсного отклика: ?_k ?_l |h(k,l)| < .

Специальные двумерные системы. На практике используются также системы с несколькими входами и/или выходами.

Допустим, система имеет i-входы и j-выходы, линейна и инвариантна к сдвигу по переменной t. Если на i-вход системы поступает одномерный единичный импульс ?_i(t) при нулевых сигналах на остальных входах, то j-выходные сигналы будут импульсным откликом системы h_ij(t). При известном полном ансамбле значений h_ij для всех i-входов, для произвольной комбинации входных сигналов s_i(t) сигнал на j-выходе будет определяться выражением:

z_j(t) = ?_i ?_k h_ij(k) s_i(t-k). (2.8)

3. Частотные характеристики сигналов и систем

Частотный отклик системы. Допустим, что двумерная ЛИС-система имеет импульсный отклик h(k?x,l?y). Подадим на вход системы сигнал вида комплексной синусоиды:

s(n,m) = exp(jn?x?_x+jm?y?_y),

где ?_x и ?_y - значения частоты сигнала соответственно по координатам x и y. Принимая ?x = 1, ?y = 1 и выполняя двумерную свертку (2.5), получаем:

z(n,m) =h(k,l) exp[j?_x(n-k)+j?_y(m-l)] =

= exp(jn?_x+jm?_y) h(k,l) exp(-jk?_x-jl?_y) = H(?_x,?_y) exp(jn?_x+jm?_y).

H(?_x,?_y) = h(k,l) exp(-jk?_x-jl?_y). (3.1)

Таким образом, выходной сигнал представляет собой комплексную синусоиду с теми же значениями частоты, что и у входного сигнала, с изменением амплитуды и фазы за счет комплексного множителя H(?_x,?_y), который носит название частотного отклика (частотной характеристики) системы. Для дискретных сигналов частотный отклик периодичен с периодом 2? по обеим частотным переменным:

H(?_x+2?k,?_y+2?l) = H(?_x,?_y).

Пример расчета частотного отклика системы.

Определить частотную характеристику системы с импульсным откликом:

h(0,0) = 0.25, h(0,1) = 0.125, h(1,0) = 0.125, h(1,1) = 0.0625.

Частотный отклик:

H(?_x,?_y) =h(n,m)exp(-jn?_x-jm?_y) = 0.25+0.125[exp(-j?_x)+exp(j?_x)+exp(-?_y)+exp(j?_y)] + 0.0625[exp(-j?_x-j?_y)+exp(-j?_x+j?_y)+exp(j?_x-j?_y)+exp(j?_x+j?_y)] = 0.25(1+cos ?_x)(1+cos ?_y).

Система является примером двумерного фильтра нижних частот. Частотный отклик системы на плоскости (?_x,?_y), приведенный на рис. 3.1, имеет осевую симметрию с коэффициентом передачи 1 в центре (?_x=0, ?_y=0) со спадом до нуля при ?_x=? и ?_y=?.

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

h(k,l)= q(k)g(l) ??Q(?_x)G(?_y)= H(?_x,?_y)

Q(?_x) = ?_k q(k) exp(-jk?_x).

G(?_y) = ?_l g(l) exp(-jl?_y).

Импульсный отклик системы. Выражение (3.1) описывает разложение функции Н(?_x,?_y) в двумерный рад Фурье с коэффициентами разложения в виде отсчетов импульсного отклика h(k,l), т.е. прямое преобразование Фурье. Очевидно, что обратным преобразованием Фурье с интегрированием в пределах одного периода из частотного отклика H(?_x,?_y) можно получить импульсный отклик системы:

h(k,l) = H(?_x,?_y) exp(jk?_x+jl?_y) d?_xd?_y. (3.2)



Круговой низкочастотный фильтр (справа - сечения по координате m)

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

1. Фурье-преобразования сигналов.

S(?_x,?_y) = ?_n ?_m s(n,m) exp(-jn?_x-jm?_y). (3.3)

s(n,m) =S(?_x,?_y) exp(jn?_x+jm?_y) d?_xd?_y. (3.4)



2. Теорема о свертке.

z(n,m) = h(n,m) _** s(n,m) ? H(?_x,?_y) S(?_x,?_y) = Z(?_x,?_y).

z(n,m) = c(n,m) s(n,m) ? C(?_x,?_y) _** S(?_x,?_y) = Z(?_x,?_y).

3. Основные свойства Фурье-преобразования.

1) Линейность (в том числе для любых комплексных чисел a и b):

аs(n,m)+bz(n,m) ? aS(?_x,?_y)+bZ(?_x,?_y).

2) Пространственный сдвиг:

s(n-N,m-M) ? S(?_x,?_y) exp(-jN?_x-jM?_y).

3) Дифференцирование:

dS(?_x,?_y)/d?_x ? -jn s(n,m),

dS(?_x,?_y)/d?_y ? -jm s(n,m),

d²S(?_x,?_y)/(d?_x d?_y) ? -nm s(n,m).

4) Комплексное сопряжение:

х*(n,m) ? S*(-?_x,-?_y).

Вещественная и мнимая части Фурье-образов последовательностей s(n,m):

S(?_x,?_y) = S*(-?_x,-?_y).

Re [S(?_x,?_y)] = Re [S(-?_x,-?_y)].

Im [S(?_x,?_y)] = -Im [S(-?_x,-?_y)].

5) Теорема Парсеваля:

?_n ?_m s(n,m) s*(n,m) = S(?_x,?_y) S*(?_x,?_y) d?_x d?_y.

В частности, при s(n,m) = s(n,m):

?_n ?_m |s(n,m)|² = |S(?_x,?_y)|² d?_x d?_y,

где левая часть уравнения представляет собой полную энергию дискретного сигнала s(n,m), a функция |S(?_n,?_m)|² - спектральную плотность энергии сигнала.

4. Дискретизация двумерных сигналов

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

s(n,m) = s_a(n?x,m?y)

где ?x и ?y - горизонтальный и вертикальный интервалы дискретизации двумерного непрерывного сигнала s_a(x,y) с непрерывными координатами x и y. Ниже значения ?x и ?y, как и в одномерном случае, принимаются равными 1.

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

Интегральные преобразования Фурье аналоговых сигналов в непрерывной шкале частот ?_x и ?_y:

S_a(?_x,?_y) =s_a(x,y) exp(-j?_xx-j?_yy) dxdy. (4.1)

s_a(x,y) =S_a(?_x,?_y) exp(j?_xx+j?_yy₎ d?_xd?_y. (4.2)

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

S(k,l) =s(n,m) exp(-jn??k/N-jm2?l/M), (4.3)

S(k,l) =exp(-jn2?k/N) s(n,m) exp(-jm2?l/M), (4.3')

s(n,m) =S(k,l) exp(-jn2?k/N-jm2?l/M). (4.4)

s(n,m) =exp(-jn2?k/N) S(k,l) exp(-jm2?l/M). (4.4')

Выражения (4.3') и (4.4') показывают, что двумерное ДПФ по прямоугольному растру дискретизации данных может вычисляться с помощью одномерных последовательных ДПФ. Вторые суммы выражений являются одномерными ДПФ сечений функций s(n,m) и S(k,l) по линиям n и k соответственно, а первые - одномерными ДПФ вычисленных функций в сечениях по m и l. Другими словами, исходные матрицы значений s(n,m) и S(k,l) пересчитываются сначала в промежуточные матрицы с ДПФ по строкам (или по столбцам), а промежуточные - в окончательные с ДПФ по столбцам (или соответственно по строкам).

Интерполяционный ряд восстановления двумерного сигнала. Если непрерывный сигнал s_a(x,y) является сигналом с ограниченным спектром, а периоды дискретизации выбраны достаточно малыми и спектры соседних периодов не перекрываются:

S_a(?_x,?_y) = 0 при |?_x|?/?x, |?_y|?/?x,

то, как и в одномерном случае, сигнал s_a(x,y) может быть восстановлен по дискретному сигналу с использованием двумерного аналога ряда Котельникова-Шеннона:

s_a(x,y) = ?_n ?_m s(n,m). (4.5)

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

Произвольный растр дискретизации. Понятие прямоугольной дискретизации обобщается на произвольный растр дискретизации с линейно независимыми векторами v₁ = (v₁₁,v₂₁)^T и v₂ = (v₁₂,v₂₂)^T, где T - индекс транспонирования (рис. 4.1). Координаты двумерного периодического множества отсчетов на плоскости (x,y):

x = v₁₁n + v₁₂m,

y = v₂₁n + v₂₂m.

С использованием векторных обозначений:

=

где = (x,y)^T, =(n,m)^T, =(v₁|v₂)- матрица дискретизации. Определитель матрицы не равен нулю, если вектора v₁ и v₂ линейно независимы. При дискретизации непрерывного сигнала s_a(x,y) матрицей формируется дискретный сигнал:

s() s_a( ).

Двумерное интегральное преобразование Фурье непрерывного сигнала по непрерывному вектору = (?₁,?₂)^T:

S_a() =s_a() exp(-j^T) d, (4.6)

s_a() =S_a() exp(j^T) d, (4.7)

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

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

S() = ?_n s() exp(-j^T), (4.8)

s() =S() exp(j^T) d. (4.9)

где: = (?_х,?_у)^T .

s() = s_a() =S_a() exp(j^T) d.

После подстановки в это выражение значения = ^T, получаем:

s() = S_a(/^T) exp(j^Т) d.

Или, с учетом периодичности по квадратным областям плоскости:

s() = S_a((-2?)/^T) exp(j^Т) d, (4.10)

где - вектор целочисленных значений периодов дискретизированной функции по осям ?_х и ?_у. Сравнивая последнее выражение с выражением (4.9), получаем:



S() =S_a((-2??/^T),

S(^T) = S_a(-), (4.11)

где - матрица периодичности:

^Т = 2?, (4.12)

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

Как и в одномерном случае, интервалы дискретизации ?x и ?y определяют главный период двумерного спектра соответственно по осям ?_x и ?_y и частоты Найквиста: ?_xN = ?/?x и ?_yN = ?/?y. Спектр дискретного сигнала также является периодическим продолжением спектра аналогового сигнала. Для исключения искажений спектра (наложения спектров боковых периодов на главный период) предельные частоты сигнала должны быть меньше частот Найквиста.



На рис. приведен пример центральной части спектра дискретного сигнала при ?x=1 и ?y=1.

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

, det = ?х?у, (4.13)

. (4.14)

Интерполяция дискретных сигналов. Для сигнала с ограниченным спектром изменением матрицы дискретизации можно подобрать матрицу периодичности таким образом, чтобы в правой части выражения (4.11) не было перекрытия спектров. Тогда для значений по точкам ^T области С главного периода спектра выражение (4.11) упрощается:

S(^T) = S_a() / |det |. (4.15)

S_a() = |det | S(^T) = |det | S(), С. (4.16)

Из выражения (4.16) следует, что при корректной дискретизации непрерывной двумерной функции ее спектр с точностью до нормировочного множителя |det | может быть восстановлен по спектру дискретной функции. Соответственно, выполнив обратное преобразование Фурье левой и правой части равенства (4.16), получим уравнение восстановления непрерывной функции по ее дискретному варианту (многомерный аналог интерполяционного ряда Котельникова-Шеннона):

s_a() = s()exp(j^T (-)) d.

s_a() = s() f(-), (4.17)

где f(..) - интерполяционная функция:

f(-) = exp(j^T (-)) d. (4.18)

Все приведенные векторные уравнения могут быть обобщены на Р-мерные функции с заменой константы 4?² там, где она встречается, на (2?)^P.

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

На практике для двумерных сигналов используют, как правило, только два варианта растров дискретизации - прямоугольный и гексагональный. Прямоугольному варианту соответствуют диагональные матрицы дискретизации и периодичности (4.13-14). Для гексагональной дискретизации, пример которой приведен на рис. 4.1, в частном случае при ?t = ?х каждый отсчет располагается на равном расстоянии от шести ближайших отсчетов, при этом матрицы дискретизации:

, .

Допустим, имеем сигнал с частотным спектром, ограниченным круговой областью частот ?_r:

S_a(?_х,?_у) = 0 при ?_х²+?_у² > ?_r².



Круговая область частот вписывается без перекрытий в квадрат со стороной 2?_r или в шестиугольник со стороной 2?_r/. Матрицы дискретизации:

_пр =, det = ?²/?_r²,

_гекс =, det = 2?²/(?_r²).

Поскольку плотность отсчетов пропорциональна 1/|det|, то отсюда следует, что для представления одного и того же сигнала гексагональный растр дискретизации требует на 4% меньше отсчетов по сравнению с прямоугольным. Эффективность "гексагональной" матрицы возрастает при увеличении размерности сигнала. Так, при 4-мерном сигнале для "гексагональной" матрицы требуется в 2 раза меньше отсчетов, чем для "прямоугольной".

5. Многомерный спектральный анализ

Периодические последовательности. Двумерная последовательность s(n,m) прямоугольно периодична, если

s_п(n,m) = s(n,n+M) = s(n+N,m)

для всех (n,m) при целочисленных значениях N и M. Минимальные значения N и M, при которых выполняется данное равенство, называют горизонтальным и вертикальным периодами функции s_п(..), которыми ограничивается прямоугольная область R_N,M главного периода, содержащая NM независимых отсчетов: 0nN-1, 0mM-1.

Последовательность s_п(n,m) с периодами N и M можно представить в виде конечной суммы (ряда Фурье) комплексных синусоид с кратными частотами:

s_п(n,m) =(1/NM)S_п(k,l) exp(j2?nk/N+j2?ml/M). (5.1)

S_п(k,l) =s_п(n,m) exp(-j2?nk/N-j2?ml/M). (5.2)

Конечные последовательности. Если s(n,m) представляет собой последовательность конечной протяженности, имеющей опорную область R_N,M, то периодическую последовательность s_п(n,m) с главным периодом R_N,M можно сформировать периодическим продолжением s(n,m):

s_п(n,m) = s(n-aN,m-bM),

s(n, m) = s_п(n, m), (n, m) < R_N,M.

в остальных случаях.

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

Аналогично можно записать и для частотной области:

S_п(k, l) = ?_a ?_b S(k-aN, l-bM).

S(k, l) = S_п(k, l), 0kN, 0lM

Отсюда значения s(n,m) и S(k,l) можно вычислить с использованием выражений (5.1-2) путем последовательности операций:

s(n,m) ?? s_п(n,m) ? S_п(k,l) ? S(k,l).



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

Таким образом:

1. Дискретизация сигнала в пространственной области вызывает периодизацию частотного спектра сигнала.

2. Дискретизация частотного спектра сигнала вызывает периодизацию его пространственного представления.

3. Прямое и обратное ДПФ сигнала ограниченной протяженности автоматически означает периодизацию как его спектра, так и его пространственного представления.

4. Сигналы, ограниченные в пространстве, можно точно отобразить отсчетами их фурье-преобразования.

5. Частотное представление сигнала с ограниченным спектром обратным фурье-преобразованием может быть точно переведено в пространственную область.

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

Многомерные последовательности. Определение ДПФ для Р-мерной последовательности с опорной областью R_P = {: 0n_iN_i-1, i=1,2,3, ... ,P} производится введением диагональной матрицы значений N_i:

при этом P-мерное ДПФ записывается в виде:

S() =s() exp(-j^T2?/). (5.3)

s() =S() exp(j^T2?/). (5.4)

Кратко рассмотрим особенности многомерных ДПФ (на примере двумерных последовательностей).

ДПФ суммы двух последовательностей с опорной областью на R_N,M равно сумме их ДПФ:

аs(n,m)+bz(n,m) ? aS(k,l)+bZ(k,l),

но при этом все ДПФ должны быть одного размера и этот размер должен быть достаточным, чтобы включить всю опорную область суммарной последовательности аs(n,m)+bz(n,m). Практически это означает, что х(..) и z(..) должны иметь одну и ту же опорную область. Опорная область каждой последовательности при необходимости дополняется нулями.

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

Допустим, что s(n,m) имеет опорную область R_P1,P2, a h(n,m) - R_Q1,Q2. Результат линейной свертки:

s(n,m) = ?_k ?_l h(k,l) s(n-k,m-l).

Опорная область последовательности s(n,m):

0nP1+Q1-1, 0mP2+Q2-1.

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

NP1+Q1-1, MP2+Q2-1.

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

Вейвлетное преобразование сигналов является обобщением спектрального анализа. Термин "вейвлет" (wavelet) в переводе с английского означает "маленькая (короткая) волна". Вейвлеты - это обобщенное название семейств математических функций, которые локальны во времени и по частоте, и в которых все функции получаются из одной базовой (порождающей) посредством ее сдвигов и растяжений по оси времени. Вейвлет-преобразования (WT) подразделяют на дискретное (DWT) и непрерывное (CWT).

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

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

истоки Вейвлет - преобразования

Историческая справка. История спектрального анализа восходит к И. Бернулли, Эйлеру и Фурье, который впервые построил теорию разложения функций в тригонометрические ряды. Однако это разложение долгое время применялось как математический прием и не связывалось с какими-либо физическими понятиями. Однако, начиная с 20-х годов прошлого века, по мере развития радиотехники и акустики, спектральные разложения приобрели физический смысл и практическое применение. Основным средством анализа реальных физических процессов стал гармонический анализ, а математической основой анализа - преобразование Фурье. Преобразование Фурье разлагает произвольный процесс на элементарные гармонические колебания с различными частотами с помощью одной базисной функции exp(j?t) или двух действительных функций sin(?t) и cos(?t). Гармонические колебания имеют широкое распространение в природе, и поэтому смысл преобразования Фурье интуитивно понятен независимо от математической аналитики.

Вейвлет-анализ является разновидностью спектрального анализа, в котором роль простых колебаний играют вейвлеты. Базисная вейвлетная функция - это некоторое "короткое" колебание, но не только. Понятие частоты спектрального анализа здесь заменено масштабом, а чтобы перекрыть "короткими волнами" всю временную ось, введен сдвиг функций во времени. Базис вейвлетов - это функции типа ?((t-b)/a), где b - сдвиг, а - масштаб. Функция ?(t) должна иметь нулевую площадь, а Фурье-образ таких функций равен нулю при ?=0 и имеет вид полосового фильтра. При различных значениях масштабного параметра 'a' это будет набор полосовых фильтров. Семейства вейвлетов используются для представления сигналов и функций в виде суперпозиций вейвлетов на разных масштабных уровнях декомпозиции (разложения) сигналов.

Первое упоминание о подобных функциях появилось в работах Хаара (Haar) еще в начале прошлого века. Вейвлет Хаара - это короткое прямоугольное колебание на интервале [0,1], показанное на рис. 1.1. Сам термин "вейвлет" ввели J. Morlet и A. Grossman в 1984 г. Они занимались исследованиями сейсмических сигналов с помощью базиса, который и назвали вейвлетом. Весомый вклад в теорию вейвлетов внесли Ингрид Добеши, разработавшая ортогональные вейвлеты (1988), Натали Делпрат, создавшая время-частотную интерпретацию CWT (1991), и многие другие. В настоящее время пакеты расширений по вейвлетам присутствуют в основных системах компьютерной математики (Matlab, Mathematica, Mathcad, и др.), а вейвлет-преобразования и вейвлетный анализ используются во многих областях науки и техники для самых различных задач. Вейвлет-анализ называют "математическим микроскопом" для точного изучения внутреннего состава и структур неоднородных сигналов и функций. Вейвлеты позволяют существенно расширить инструментальную базу информационных технологий обработки данных.

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

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

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

||s(t)||² = s(t), s(t) =s(t) s^*(t) dt,

где s^*(t) - функция, комплексно сопряженная с s(t). Если норма конечна (интеграл сходится), то функция принадлежит пространству функций L²[R], R=[0,T], интегрируемых с квадратом (пространство Гильберта), и имеет конечную энергию. В этом пространстве на основе совокупности ортогональных функций с нулевым скалярным произведением

v(t,?), v(t,?) =v(t,?) v^*(t,?) dt = 0

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

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



S_n = (1/T)s(t) exp(-jn??t) dt, n (-?, ?),

где ??=2?/T - частотный аргумент векторов. Тем самым, сигнал s(t) может быть представлен в виде ряда Фурье и однозначно определяется совокупностью спектральных коэффициентов S_n этого ряда:

s(t) =S_n exp(jn??t).

Таким образом, ряд Фурье - это разложение сигнала s(t) по базису пространства L²(0,T) ортонормированных гармонических функций exp(jn??t) с изменением частоты, кратным частоте первой гармоники ?₁=??. Ортонормированный базис пространства L²(0,T) построен из одной функции exp(j??t) = cos(??t)+j·sin(??t) с помощью кратного масштабного преобразования независимой переменной. Как правило, Ряд Фурье ограничивается определенным количеством членов N и равномерно сходится к s(t) при N>?.

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

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

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

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

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

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

S(?,b_k) = s(t) w*(t-b_k) exp(-j?t) dt.



Функция w*(t-b) представляет собой функцию окна сдвига преобразования по координате t, где параметром b задаются фиксированные значения сдвига. При сдвиге окон с равномерным шагом значения b_k принимаются равными k?b. В качестве окна преобразования может использоваться как простейшее прямоугольное окно, так и специальные весовые окна (Бартлетта, Гаусса, и пр.), обеспечивающие малые искажения спектра при вырезке оконных отрезков сигналов (нейтрализация явления Гиббса).

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

S(t,?) = s(t-?) w(?) exp(-j??) d?.

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

На рис. приведен пример вычисления и представления (модуль правой части главного диапазона спектра) результатов спектрограммы при дискретном задании зашумленного входного сигнала sq(n). Сигнал представляет собой сумму трех последовательных радиоимпульсов с разными частотами без пауз, с отношением сигнал/шум, близким к 1. Оконная функция w_i задана в одностороннем варианте с эффективной шириной окна b 34 и полным размером М = 50. Установленный для результатов шаг по частоте ?? = 0.1 несколько выше фактической разрешающей способности 2?/M = 0.126. Для обеспечения работы оконной функции по всему интервалу сигнала задавались начальные и конечные условия вычислений (продление на M точек обоих концов сигнала нулевыми значениями).

Как видно из приведенных примеров, оконное преобразование позволяет выделить информативные особенности сигнала и по времени, и по частоте. Разрешающая способность локализации по координатам и по частоте определяется принципом неопределенности Гейзенберга. В силу этого принципа невозможно получить произвольно точное частотно-временное представление сигнала. На рис. 1.4 приведен пример частотно-временного оконного преобразования сигнала, состоящего из 4-х непересекающихся интервалов, в каждом из которых сумма двух гармоник разной частоты. В качестве окна применена гауссова функция разной ширины. Узкое окно обеспечивает лучшее временное разрешение и четкую фиксацию границ интервалов, но широкие пики частот в пределах интервалов. Широкое окно напротив - четко отмечает частоты интервалов, но с перекрытием границ временных интервалов. При решении практических задач приходится выбирать окно для анализа всего сигнала, тогда как разные его участки могут требовать применения разных окон.

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



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