ДВУХЭТАПНЫЙ КВАЗИОПТИМАЛЬНЫЙ АЛГОРИТМ ВЕЙВЛЕТ-ФИЛЬТРАЦИИ СИГНАЛОВ И ИЗОБРАЖЕНИЙ Текст научной статьи по специальности «Математика»

УДК 519.2

Двухэтапный квазиоптимальный алгоритм вейвлет-фильтрации сигналов и изображений

Ю.Е. Воскобойников ФГБОУВПО НГАСУ (Сибстрин), Новосибирск, Россия

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

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

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

ВВЕДЕНИЕ

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

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

• обработка «зашумленных» коэффициентов разложения (устранение шумовой составляющей);

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

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

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

используемых алгоритмов носят пороговый характер: коэффициент разложения меньший по абсолютной величине некоторой пороговой величины зануляется, в противном случае коэффициент сохраняется или подвергается некоторому (в общем случае нелинейному) преобразованию. Распространение на практике получили «жесткая» и «мягкая» пороговые функции, зависящие только от одного параметра - величина порога. При этом выбор пороговой величины существенно влияет на ошибку фильтрации и эта величина, по сути, является управляющим параметром алгоритмов вейвлет - фильтрации [1-3].

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

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

Поэтому в данной работе рассматривается задача построения алгоритма вейвлет-фильтрации, реализуемого на практике, но имеющего СКО незначительно превышающего СКО винеровского алгоритма фильтрации (не более чем на 10 %).

1. ТЕОРЕТИЧЕСКИЕ ОСНОВЫ МЕТОДА

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

Удобной формой вейвлет - представления является многомасштабное (multiresolution) представление вида [4-5]:

ДО = X а/0+з к(г)+ + X Xа],к¥],к(о

/=/о +1 к

которое можно интерпретировать как восстановление исходного сигнала /(г) по его коэффициентам разложения на 3 -ом уровне. Функции {ф/ к (г)} называют масштабирующими

(или отцовскими), а функции {у/ к (г)} - вейвлет

- функциями (или материнскими). Коэффициенты

разложения а.к называют аппроксимирующими, ё/ к - детализирующими и они определяются выражениями:

а/,к = 1/(?)фк (№ , / = | /(г)/ (г)ёг, (6)

к к

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

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

Напомним, что системы функций {ф/к (?)}, {У/,к (?)} составляют ортогональный

базис пространств у/, -~< / <~, и определяются соотношениями [2-3]:

Ф/к (г) = 2-//2 ф(2-1 г - к), уик (г) = 2-//2 у(2-1 г - к), где функции {ф(г - к)}, {у(г - к)} -ортогональный базис «нулевых» пространств у0,ш0. При этом пространство ортогонально

пространству V/ и У/_1 = V/ + Ш/ . Чем меньше

номер /, тем более «мелкие» структуры исходной функции / (г) могут быть

представлены в базисах {ф,к (г)}, {У/ к (г)} и

тем ближе реконструированный сигнал /(г) к исходному. Вейвлеты, базисные функции, которых удовлетворяют условиям ортонормиро-ванности, называются ортогональными и к ним относятся вейвлеты Добеши, Симлета, Кайфлета [4,5]. Ниже рассматриваются только ортогональные вейвлеты.

На практике для вычисления коэффициентов а/к, ё/ к вместо интегрирования (2) используется пирамидальный алгоритм Малла [4-5], требующий порядка N операций, что делает его очень эффективным инструментом в вейвлет - технологиях. В этом алгоритме порядок вычисления коэффициентов следующий: первоначально (точка старта алгоритма) задаются аппроксимирующие коэффициенты на некотором «нулевом» уровне

/0: а /о к, к = 1,..., N . Величина N

удовлетворяет условию N = 2т, т > 1 . Затем

по этим коэффициентам вычисляются аппроксимирующие и детализирующие коэффициенты вычисляются а/ +1к , ё +1к, по

коэффициентам а/ +1к находятся а +2,к,

ё/0 +2,к ; по а/о +2,к вычисляются а+з,к , ^ +зД . Таким образом, по аппроксимирующим коэффициентам а/ + /-1,к определяются

коэффициенты а](+ + /^ , ^ , / = 1,...,3 ,

N

к = 1,2,...^~~. В алгоритме Малла для

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

Для задания коэффициентов а/^к, к = 1,...,N , используется несколько

способов. Наиболее частым и простым является задание а/0,к = /^), к = 1,...,N , где /^) -

значение функции в равноотстоящих узлах сетки. Такое задание означает, что исходная функция / (г) аппроксимируется функцией

N

ц(г) = X / (гк )Ф/0,к (г). В дальнейшем будем

к=1

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

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

среднее и дисперсию а . Это означает, что коэффициенты ^/^к,к = 1,...,N , задаются с

такой же погрешностью и коэффициенты разложения <Л/0 +/к , й/0+/к , / = 1,...,3 , также

вычисляются со случайными ошибками. Для упрощения записи уровней разложения в дальнейшем полагаем /0 = 0. С учетом этого допущения на 3 - ом уровне разложения имеем набор коэффициентов:

Мё2,к К...^/,к К-.Д^ ,к },{аз ,к }, N

к = 1,2,...,— , / = 1,...,3 . (7) 2 /

По этим коэффициентам вычисляются оценки {ё/,к },{<53 ,к} для точных коэффициентов {ё/к },{а/к } и строится «сглаженная» функция

/(г) вида (1). Определим

среднеквадратическую ошибку (СКО) фильтрации выражением:

Д(/) = M

J (ДО - f (t))2 dt

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

Д( /)=M

+ M

V (а - а )2

J ,k J ,k'

VV (dj,k - djk )2 ,j=1 k

(5)

Величины aJ к, с1] к , доставляющие минимум

функционалу Д(7"), назовем оптимальными оценками коэффициентов вейвлет-разложения. Оптимальные оценки будем искать в виде:

а-- =<k • äJk, d k -wj,k • dj,k ■ (6)

J,к ^ ,к Jк ],к

Из необходимых и достаточных условий

минимума Д( /) получаем оптимальные множители [2,3]:

w;

1

oPtj ,k

1 +

wo

1

°Pti,k

J ,k

1 + 5

и оптимальные оценки 1

oPj ,k

1 + 5

J ,k

а ; d„„t

J k oPtjk

1 + 5

j,k

-г-d

d J ,k j,k

(7)

(8)

где

отношения

с<а _

SJ ,k -

s

aJk

od S j = d2

можно трактовать как отношение «шум/сигнала» соответствующего

коэффициента разложения. Можно показать,

что СКО Д(/) при оптимальных множителях (7) определяется выражением:

Д( Ipt) - V

S'

J ,k

k

1 + S

0/,k

J ,k

+VV

j-1 k

S

d j,k

1 + S

■ d

jk

2},k

(9)

из которого видно, что при достаточно малых а2 величина Д(/ор1:) имеет порядок а2, т.е.

Л 2

Д()~ о . В дальнейшем алгоритм вейвлет-

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

Для удобства дальнейшего изложения

введем:

• вектор «точных» коэффициентов, вычисленных по не зашумлённым значениям функции

I • • • • Т

0 = +1,к :А}0 +2,к :"':¿¡0 +J,к :а}0 +J,к\ ;

• вектор «зашумленных» коэффициентов, вычисленных по зашумлённым значениям функции

0 = \<*]0+1,к : ¿¡0 +2,к '': ¿¡0 +J,к : а]0 +J,к | . При этом I -ую проекцию этих векторов будем

обозначать как: в(0, в(1) ■ (8) Тогда выражения

переписать в виде:

1

(8), (9) можно

dopt(i) -

1 + S (i)

в(о -

q2(i)

q2(i)+s2

D(fopt) = V I^« = V

S

e2(i) + s2

в(о;

q2(i) >

где

S (i) -

- s

q2(i)

(10) (11)

(12)

- отношение «шум/сигнал», которое имеет следующий предел:

2

lim S(i) -¥ , если s или в(0 ® 0. (13)

При этом оптимальный фильтрующий

множитель

Wopt(i) -

1

q2(i)

1 + S (i)

q2(i)+s2

имеет предел:

lim W;pt(i) - 0.

S (i

В матричном виде оптимальный можно записать в виде:

1r

(14)

(15) алгоритм

(16)

/—1

/ор{ = Ж" 'ПорЖГ ,

Где Ж, Ж_1 - матричные обозначения прямого и обратного дискретного вейвлет-преобразования,

= ¿т\Wopt (1), Wo.pt (2),..., Фор? (#)} -

диагональная матрица оптимальных

фильтрующих множителей.

К сожалению, построенные

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

Эмпирический винеровский

алгоритм вейвлет-фильтрации. В этом алгоритме вместо неизвестного «точного» коэффициента 0(1) берется «зашумленный»

коэффициент разложения 0(0 и тогда

получаем следующие характеристики алгоритма:

Л/ .... 1

S(i) -

- s

■~г и w(i) = ■ (17) в (i) 1 + S(i)

В работах [2,3] показано, что коэффициент 0(0 можно представить в виде:

0(0 = 0(0 + £(0,

я

+

k

1

где

- случайная величина с нулевым

средним и дисперсией s2 . Тогда из (17) видно,

2 уже является

что

s (0 =

(0(i) + e(i) )2

,2

случайной величиной и при а ® ¥ среднее значение Ё(') не будет стремиться к ¥, а, следовательно, и д') также не будет стремиться к 0.

Для иллюстрации этого важного факта на Рис. 1 приведены кривые: 1 - значения оптимального фильтрующего множителя (14); 2 - значения фильтрующего множителя д(/) (17) при разных значениях отношения

«шум/сигнал» Ё(').

Рис. 1. Значения фильтрующих множителей

Видно, что среднее значение д(/) (вычисленное по выборке объемом 400) при увеличении Ё(') стремится к постоянной величине (примерно 0.35), в то время как оптимальный фильтрующий множитель стремиться к нулю . Следовательно, эмпирический винеровский алгоритм вейвлет-фильтрации имеет существенный недостаток: при отношении «шум/сигнал» больше 1 фильтрующие множители д') не осуществляют

эффективную фильтрацию шума, и результат такой фильтрации имеет СКО фильтрации значительно большую по сравнению с оптимальным (винеровским) алгоритмом вейвлет-фильтрации.

Предположим, что вместо точного (но неизвестного) коэффициента 0(1) в (12) будет

использоваться некоторая оценка 0(1) этого коэффициента. Возникает вопрос: насколько увеличиться СКО фильтрации при использовании эмпирического отношения «шум/сигнал»:

„2

S(i) =

s

q2(i)

(18)

Можно показать, что СКО фильтрации (4) (обозначим как А) при использовании (18) можно представить в виде:

(19)

А = Aopt + Aemp ,

где

А = ^ S(i)

Хopt = L

г 1+s (i)

q2(i)

(20)

СКО фильтрации оптимального фильтра;

1

1

(21)

Ае = Yq2(i) • (1+s(i)) • .

emp w v yjj 1 + S(i) 1 + s(i) ^

- дополнительная (эмпирическая) СКО, обусловленная заменой точного отношения S(i) на эмпирическое S(i). Очевидно, что, если S(i) = S (i), то Aemp = 0.

Вывод: чем точнее будет оценка 0(i), тем меньше будет величина Aemp и меньше будет суммарная СКО фильтрации (19).

Двухэтапный алгоритм вейвлет-фильтрации. Как способ получения более точной оценки 0(i) (по сравнению с 0(i)) для выражения (18) в работе [6] был предложен алгоритм (названный эмпирическим

винеровским фильтром), суть которого можно выразить следующими этапами:

Этап 1. С использованием дискретного ортогонального вейвлет-преобразования

(матричное обозначение W) вычисляются

коэффициенты разложения 0, т.е. 0 = Wx f. Затем, с использованием порогового алгоритма обработки Ti(0) вычисляются «сглаженные» коэффициенты разложения 0 = Т^(0). После этого берется обратное вейвлет-преобразование W"1 и получается вектор f отфильтрованных значений, вычисление которого в матричном виде представимо как:

Л = W-1 • Ti(0) -WJ = W{X0X. (22)

Напомним основные понятия пороговых алгоритмов фильтрации (подробнее см. [1, 3]). В этих алгоритмах каждый коэффициент разложения 0(i) подвергается индивидуальной обработке с использованием пороговых функций (в зарубежной литературе такая обработка получили название thresholding). Наиболее часто используются две пороговые функции:

• «жесткая» пороговая функция вида: 0(i), если 0(i) > l,

Thi (0(i)) = ■

0, если 0(i) < l;

«мягкая» пороговая функция вида:

2

П (0(1)) = ^п (0(1)) • [|0(/)| - 2]+ , (24)

где [*1+ =

обозначение г, если г > 0; 0, если г < 0.

[

определяется как

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

• из-за уменьшения амплитуды коэффициента разложения на величину 2 в

функции ТЦ (в(1)) возможно сглаживание

(размытие) контрастных элементов обрабатываемого сигнала;

• наличие в функции Т^ (в (г)) разрыва в

окрестности точки 2 в определенной степени обуславливает появление осцилляций (эффекта Гиббса) в особых точках фильтруемого сигнала.

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

переглаженным.

В работе [6] для получения оценок 0(1) обрабатываются только детализирующие коэффициенты {¿¡,к} (аппроксимирующие {aJ к} остаются без изменения из-за их

небольшой относительной погрешности (см.

[3])) с использованием функция (23) с так

называемым универсальным порогом (для каждого ¡ - го уровня разложения):

1=а 21п

(25)

где

N ,■

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

Этап 2. К вектору / применяется дискретное вейвлет преобразование Ж2, которое отличается от Ж . Например, если в качестве Ж использовались вейвлеты Добеши

БаиЬ(8), то

в качестве

Ж

можно

использовать ОаыЬ(\6) или БаиЬ(4). Таким образом, вычисляется вектор коэффициентов разложения

0ц = Ж 71. (26)

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

Ф(0 =

^О)

1

1

42о)+о2

1+

о2

42о)

1 + ^(О

(27)

и находятся разложения

ад =

сглаженный коэффициенты 1

42(0 + о2

ад=-

- ^(О, (28)

1 + ¿2(0

где 62(/) - коэффициенты разложения

62 = Ж /. В завершении этапа 2 вычисляется

вектор /2 = Ж2-1 #2, содержащий результат фильтрации описанным алгоритмом.

Если обозначить через Nт - число

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

что дискретное преобразование Ж21 = Ж2Ж-1 трансформирует («растягивает») Nт

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

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

II" II2

ИЛ- л\ г.. Ц/ -/||

: [2.2 - 2.8]

7ор? Л

/ор? /

: [1.4 +1.9]. (29)

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

Квазиоптимальный двухэтапный

алгоритм вейвлет-фильтрации. По-

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

низкая точность оценки " = т1 (в),

вычисленные с помощью пороговой функции

(23) с пороговой величиной (25). Показано [3], что эта пороговая величина не является

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

Критерий оптимальности был предложен в работе [7] и широко используется для выбора параметра регуляризации линейных регуляризирующих алгоритмов решения некорректных задач [8,9]. Алгоритмы, построенные на основе этого критерия, позволяют оценить оптимальный (с хорошей точностью) параметр регуляризации, не используя при этом количественную информацию об искомом решении задачи. Для выбора параметра сглаживания в линейных алгоритмах Фурье - фильтрации этот критерий использовался в работе [3, с. 37-45].

Заметим, что все выше названные алгоритмы являются линейными (при заданной величине параметра). Пороговые алгоритмы вейвлет -фильтрации являются нелинейными. Однако идея критерия оптимальности позволила и в этом случае эффективно оценить оптимальную пороговую величину.

Следуя работам [2, 3] величину порога на каждом j-ом уровне разложения будем задавать выражением

1j (b) = b2 • ln(Nj), (30)

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

фильтрации при использовании пороговых величин (30). Для этого введем критерий pW (a), определяемый выражением

m eo• f

Pw (b) = I^y-, (31)

i=1 s

где eь = fi - f д. - компоненты вектора невязки,

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

при параметре b. Можно показать, что при использовании ортогональных вейвлетов, критерий pW (b) вычисляются через коэффициенты вейвлет-разложения:

• для функции soft

1 j0 +J

Pw (b) =—■ I Idj,k• (djk-soft(djk,ij(b)));

s j=J0 +1 k

• для функции hard

1 j0 +J

Pw (b)--2 • I Idj,k • (dj,k - hard(djk ;1j (b)))

s j=J0 +1 k

В качестве оценки bW оптимального параметра b значение bW , для которого выполняется неравенство

Jm,y/2 £ PW (aw ) £ Jm,1-g2 , (32)

где Цп,г/2, «т,1-^2 квантили Ст -распределения с т степенями свободы уровней у/2 , 1 -у/2 соответственно, у - вероятность ошибки первого рода при проверке статистической гипотезы об оптимальности параметра сглаживания, т = N - число значений фильтруемого сигнала. Тогда в качестве пороговой величины, вычисленной на основе критерия оптимальности, принимается значение

^^ГП^). (зз)

Для вычисления Ьш используется любая итерационная процедура решения уравнения

рш (а) = т, (34)

которая завершается, как только очередное

приближение Ь(и) удовлетворяет условию (32).

2

Если число степеней свободы т > 30, то %т -распределение хорошо аппроксимируется нормальным распределением с математическим ожиданием т и дисперсией 2т . Тогда, приняв вероятность ошибки первого рода у= 0.05, получаем формулы для вычисления квантилей, входящих в неравенство (32)

J

m,0.025

J

m,0.975

= m - 1.9^V5ot ; l.9^yf2m .

= m +

(35)

Существенной чертой приведенного

алгоритма вычисления Ьш является

2

использование дисперсии шума а . На практике, как правило, эта величина неизвестна, и в этом случае можно использовать оценку „ те<И<2п(\ ё!,к |)

а =-!—, (36)

0.6745

где оператор тв<1ап(\ к \) вычисляет медиану абсолютных величин детализирующих коэф-фиентов уровня разложения /0 + 1 (объем выборки равен N/2). Эта оценка широко используется в робастных алгоритмах регрессионного анализа. Применительно к алгоритмам вейвлет - фильтрации эта оценка исследовалась в работе [2], где была показана приемлемая точность оценки (36). Так для

2 -1 заданной дисперсии а = 0.91-10 значения

оценки (36), вычисленные по 30 реализациям

длиной (N/2=1024) находились в интервале

[0.88 10-1,0.97-10-1].

Заметим, что проведенные в работах [2, 3] исследования (с применением и других алгоритмов выбора пороговой величины), показали, что в качестве наилучшей пороговой величины в пороговых алгоритмах обработки коэффициентов разложения зашумленных сигналов и изображений, целесообразно

#

использовать пороговую величину Лу ,

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

Проиллюстрируем это только одним рисунком. На Рис. 2 как функции параметра Д показаны:

• кривая 1 - значения относительной

f - f

I = 0.12). Выбор такой формы сигнала

ошибки фильтрации Д(Д) =

fWp f

• кривая 2 - значения критерия рш (Д);

• кривые 3,4 - границы интервала

[^2048,0.025, ^2048,0.975 ] , ЮТОрьге в масштабе рисунка практически совпадают.

Видно, что значения параметра Д, при которых рш (Р) попадают в интервал

[^2048,0.025, ^048,0.975], т. е. значения Д (на рис. обозначены РЖ ) соответствуют минимуму относительной ошибки фильтрации. Поэтому на первом этапе описанного двухэтапного алгоритма будем использовать пороговые

величины

в эффективности которых

убедимся на результатах эксперимента.

следующего

Рис. 2. Выбор параметра Д на основе критерия оптимальности

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

которого (вектор /) приведены на рис. 3 (относительный уровень шума

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

Рис. 3. Зашумленные значения сигнала

По зашумленным значениям функции (вектор f) были вычислены следующие векторы:

• fiu - вейвлет-фильтрация

(Wi = Daub(8)) с пороговой функцией (23) и пороговыми величинами (25);

• f1W - вейвлет-фильтрация (W1 = Daub(8)) с пороговой функцией (23) и пороговыми величинами (33);

• fopt - оптимальная (винеровская) вейвлет-фильтрация Daub(8) с фильтрующими множителями (7);

• f2U - результат двухэтапной вейвлет-фильтрация (в качестве оценки f1 брался вектор /ш), W2 = Daub(4));

• f2W - результат двухэтапной вейвлет-фильтрация (в качестве оценки f1 брался вектор f1W ), W2 = Daub(4)).

На Рис. 4 показаны «точные» значения функции (сплошная кривая) и значения проекций вектора f2W, из сравнения которых видно хорошее качество фильтрации предложенным двухэтапным алгоритмом.

2

2

4.25

2.5

0.75

500

lxlO3

5x1 ff*

Рис. 4. Отфильтрованные значения функции

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

1|2

f1U f

2

fopt f

= 2.65:

= 1.38:

f1W f

= 1.41. (37)

= 1.08. (38)

fopt f

Сравнивая отношения (37), видим, что использование пороговой величины (33) на первом этапе существенно уменьшает ошибку фильтрации и это обуславливает «вторичное» уменьшение ошибки фильтрации на втором этапе - превышение ошибки фильтрации над минимальным значением составляет всего 8 %. Это позволяет назвать предложенный двухэтапный алгоритм с пороговыми величинами (33) квазиоптимальным.

ЗАКЛЮЧЕНИЕ

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

ЛИТЕРАТУРА

[1] Fodor, I. K., Kamath C. Denoising through Wavelet Shrinkage: An Empirical Study. SPIE Journal on Electronic Imaging. 2003. Vol.12, N 1, pp. 151-160.

[2] Воскобойников Ю.Е., Гочаков А.В. Оценивания оптимальных пороговых величин в алгоритмах вейвлет-фильтрации изображений. Автометрия. 2011. т. 47. №2. С. 3-12.

[3] Воскобойников Ю.Е., Гочаков А.В., Колкер А.Б. Фильтрации сигналов и изображений: Фурье и вейвлет алгоритмы (с примерами в Mathcad) Новосибирск: НГАСУ (Сибстрин), 2010. 188 с.

[4] Mallat S. Multiresolution approximation and wavelet orthonormal bases of LA2(R). Trans. AMS. 1989. v.315. N1. P. 69-87.

[5] Mallat S. А theory of multiresolution signal decomposition: the wavelet representation. IEEE Trans. Pattern Anal. Machine Intell. 1989. v.11. N 9. P. 674-693.

[6] Ghael S.P., Sayeed A.M., Baraniuk R.G. Improved wavelet denoising via empirical Wiener filtering. Proceedings of SPIE, 1997. Vol.3169. P.389-399.

[7] Воскобойников Ю. Е. Оценивание оптимального параметра регуляризирующего алгоритма восстановления изображений. Автометрия. 1995. № 3. С. 64-72.

[8] Воскобойников Ю. Е. Устойчивые методы и алгоритмы параметрической идентификации: Монография Новосибирск: НГАСУ (Сибстрин), 2006.-186 с.

[9] Воскобойников Ю.Е. Устойчивые алгоритмы решения обратных измерительных задач: Монография. Новосибирск: НГАСУ (Сибстрин), НГАСУ. 2007. - 184 с.

Воскобойников Юрий Евгеньевич, доктор физ.-мат. наук, профессор, Заслуженный работник Высшей школы РФ, Соросовский профессор,

действительный член МАИ, РАЕ, МАН ВШ, заведующий кафедрой прикладной

математики Новосибирского государственного архитектурно-строительного университета (Сибстрин), профессор кафедры автоматики НГТУ. Автор более 286 публикаций, 5 монографий, посвященных решению некорректных задач интерпретации данных и обработке сигналов и изображений и большого числа учебных пособий. E-mail: voscob@mail.ru

Two-Step Quasi-Optimal Algorithm of Wavelet Filtering of Signals and Images

Yu.E. VOSKOBOYNIKOV

Abstract: This paper considers the construction of quasi-optimal algorithm, similar in structure to the Wiener optimal algorithm. In the proposed algorithm the evaluation of unknown noise-to-signal relationship is performed for each coefficient in the second step on the base of signal obtained in the first step using a thresholding wavelet filtration algorithms. It is shown that with the choice of threshold values by the optimality criterion with the use of such filtering algorithm the error is significantly reduced as compared with the conventional wavelet-filtering algorithms.

Keywords: wavelet filtering algorithms, optimal algorithm, quasi-optimal algorithms, two-stage

2

2

algorithms of wavelet filtering.

REFERENCES

[1] Fodor, I. K., Kamath C. Denoising through Wavelet Shrinkage: An Empirical Study//SPIE Journal on Electronic Imaging. 2003. Vol.12, N 1, pp. 151-160.

[2] Voskobojnikov Ju.E., Gochakov A.V. Ocenivanija optimal'nyh porogovyh velichin v algoritmah vejvlet-fil'tracii izobrazhenij //Avtometrija. 2011. t. 47. №2. S. 3-12.

[3] Voskobojnikov Ju.E., Gochakov A.V., Kolker A.B. Fil'tracii signalov i izobrazhenij: Fur'e i vejvlet algoritmy (s primerami v Mathcad) Novosibirsk: NGASU (Sibstrin), 2010. 188 s.

[4] Mallat S. Multiresolution approximation and wavelet orthonormal bases of LA2(R) // Trans. AMS. 1989. v.315. N1. P. 69-87.

[5] Mallat S. A theory of multiresolution signal decomposition: the wavelet representation // IEEE Trans. Pattern Anal. Machine Intell. 1989. v.11. N 9. P. 674-693.

[6] Ghael S.P., Sayeed A.M., Baraniuk R.G. Improved wavelet denoising via empirical Wiener filtering//Proceedings of SPIE, 1997. Vol.3169. P.389-399.

[7] Voskobojnikov Ju. E. Ocenivanie optimal'nogo parametra reguljarizirujushhego algoritma vosstanovlenija izobrazhenij. Avtometrija. 1995. № 3. S. 64-72.

[8] Voskobojnikov Ju. E. Ustojchivye metody i algoritmy parametricheskoj identifikacii: Monografija / Ju. E. Voskobojnikov. - Novosibirsk: NGASU (Sibstrin), 2006.-186 s.

[9] Voskobojnikov Ju.E. Ustojchivye algoritmy reshenija obratnyh izmeritel'nyh zadach : Monografija. Novosibirsk: NGASU (Sibstrin), NGASU. 2007. -184 s.