CC BY
30ISSN 0868-5886
НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2009, том 19, № 1, c. 67-79
МАТЕМАТИЧЕСКИЕ МОДЕЛИ
УДК 621.391; 519.21
© В. М. Малыхин, Г. Ф. Малыхина, А. В. Меркушева
МОДИФИКАЦИЯ ВЕЙВЛЕТ-ПРЕОБРАЗОВАНИЯ: КОНЦЕПЦИЯ И ПРИЛОЖЕНИЕ
Описана модификация вейвлет-преобразования (ВП). Построение основано на аналогии с обобщенной формой преобразования Фурье, использовании аппарата преобразования Радона и гибком способе масштабирования полученной модификации ВП. Полученная модификация вейвлет-преобразования (и особенно ее вариант с дополнительным масштабированием) имеет определенное преимущество перед традиционной формой ВП. Это относится к задаче увеличения информационнозначимой компоненты сигнала за счет снижения уровня шума, включая задачу анализа двумерного сигнала — изображения, которое искажено сопутствующим шумом (наличием так называемого блера).
Кл. сл.: вейвлет-преобразование, преобразование Радона, модификация, операторы, локальная и расширенная формы, снижение шума.
ВВЕДЕНИЕ
Вейвлет-преобразования (ВП) находят применение главным образом в задачах анализа нестационарных сигналов информационно-измерительных систем (ИИС)1-1. Класс ВП является расширением и некоторым обобщением класса время-частотных преобразований (ВЧП)2). Различные типы преобразований двух классов достаточно часто используются в таких приложениях, как анализ речевых сигналов [5], скремблирование (как средства скрытой передачи информации) [6], идентификация типа интервала нефть/газ в системах контроля расхода товарной нефти [7], согласованное многоканальное разделение сигнала [8], получение оценок динамического спектра нестационарных сигналов [9]3).
Параллельно с практическим использованием идет развитие новых возможностей ВЧП и ВП.
- Использование концепции фреймов породило новую форму представления ВЧП и ВП [10].
1) Вейвлет-преобразование (появившееся в 80-90-х годах) имеет фундаментальную теоретическую основу [1-3]. В [4] приведена краткая (достаточная для практического применения) сводка основных соотношений по ВП и расширенный перечень библиографии.
2) При представлении динамического спектра нестационарного сигнала на время-частотной плоскости ВП (в отличие от ВЧП) вместо частоты т использует параметр масштаба a — a = c/т, — который (с точностью до единиц измерения) эквивалентен величине обратной частоты.
3) Ссылки [5-9] отражают разработки авторов статьи, но
они не только позволяют познакомиться с типами при-
ложений, связанных с применением ВП и ВЧП, но и
содержат обширную библиографию по теории и приложениям этих преобразований.
- Усовершенствование структуры ядра повысило разрешение ВЧП [11].
- Для преобразования Бертрана установлена двойственность свойств ВЧП и ВП [12].
- Разработан метод модификации ВЧП с повышенным разрешением в области текущей частоты сигнала [13].
- Получен метод снижения уровня шума на базе фильтрации сигнала в области его ВП [14, 15] (метод усовершенствован в [5]).
- Предложен метод построения новых видов ВЧП на основе унитарных преобразований [16].
Наиболее продуктивной оказалась модификация преобразований, связанная с получением их обобщенной структуры, которая при определенном задании параметров переходит в традиционный вид преобразования. На этой основе получено обобщенное преобразование Фурье (ОПФ), построенное с использованием вращаемой системы координат. Аналогичный подход целесообразно применить для получения модифицированного ВП5).
4) См., например, [17] со ссылками на другие работы, связанные с модификацией преобразования Фурье (ПФ).
5) Можно отметить использованный в [18] термин "ап-паратурно-ориентированные вейвлеты" (АОВ) наряду с применявшимся в [19] выражением "модифицированные вейвлеты". Однако суть конструкции заключается во включении передаточной функции прибора в состав вейвлет-функции, используемой для получения представления сигнала. Конструкция эмулирует основные соотношения ВП и достаточно громоздкая, однако передаточная функция прибора (или каналов ИИС) обычно неизвестна и именно поэтому обычно ВП применяется к сигналу на выходе сенсоров или отдельного прибора.
МОДИФИЦИРОВАННАЯ ФОРМА ВЕЙВЛЕТ-ПРЕОБРАЗОВАНИЯ (ОБЩАЯ КОНЦЕПЦИЯ)
Прототипом модифицированного ВП (МВП) является обычное двухмерное ВП, причем МВП также служит двухмерным отображением
Я2 ^ Я2 и позволяет анализировать вектор-сиг-
() ^
. Для построения МВП использу-
^))
ется гладкая одномерная и достаточно быстро убывающая функция ), для частотной формы которой, определяемой преобразованием Фурье (ПФ) — $(£), выполняется условие6):
Г И^Г г
интеграл I-1-■ йд ограничен, т.е. <г (1)
нал s(t) =
Ш
При этом можно считать, что нормировка функции y/(t) такова, что ПФ ¥(Ш) удовлетворяет ус-
?|¥(Ш)Г аШ ,
ловию J J-— • ад =1 .
- Ш2
МВП определяется выражением (2), где а > 0 и b — параметры масштаба и смещения ВП, а ве [0,2п]— угол вращения вектор-сигнала s с компонентами s1, s2:
Yabe(s) = a~V2 ¥ (s1cose+s2sine-b)/a ). (2)
Функция y/abe (s) постоянна на прямой линии
s1cose+s2sine= const, а при значениях в= 0 и в = п / 2 обращается в вейвлет-функцию одной из компонент вектор-сигнала s:
¥a,b,0(s) = a ~V2 •¥{(. S1 - b)/ a );
-1/2
¥a,b,n/2 (s) = a •¥{(S2 - b)/ a )•
Если задана функция f (s), то могут быть определены коэффициенты Rf (a, b,e) разложения ее по МПФ:
Rf (a,b,e) = J J ¥a,b,e(s) •f (s)dsi d^2,
или в более компактной форме
(3)
R
■ (a,Ь,0) = \ ь ^/) .
6) Условие (1), в частности, выполняется, если среднее значение равно нулю, т.е. | \(£) ■ йд = 0.
При этом набор этих коэффициентов оказывается достаточным для восстановления самой функ-
7)
ции :
/(8) = | ] |Я/ (а,Ь, в) ■ \аДв(8) • Т~(а,Ь, в). (4)
п „„ п а 4п
Определение коэффициентов Я/ (а,Ь,в) разложения функции и представления ее с помощью обобщенного вейвлет-преобразования производится с использованием преобразования Радона. Преобразование Радона (ПР) осуществляет представление функции двух переменных набором (множеством) двухпараметрических коэффициентов { Я/(в,?) }(в, ()е[0 2п]хК — набором интегралов
этой функции вдоль линий, проходящих через начало координат с различным углом наклона8). В нашем случае это функция / (8) двухмерного вектор-сигнала 8 с компонентами 51 и s2. Таким образом, преобразование Радона Я определяется соотношением
Я / (в, Г) =
= | / (51, s2) ■ 3^1СО§в + s2sinв-(5)
где 3 — распределение Дирака, т. е. функция со свойствами
П при х = 0, Г „
8( х) = 1п П 1 8( х)йх = 1.
[0 прих ф 0; -Г
Коэффициенты МВП функции / (8) выражаются в виде вейвлет-преобразования коэффициентов Я / (в, I) преобразования Радона:
Яу (а, Ь,в) = 1 Я /(в, 0 ■ а~у\((г - Ь)/а (6)
7) Как и у вейвлет-преобразования (или преобразования Фурье), соотношению (4) соответствует условие Парсе-валя [20, 21]. Для МВП условие Парсеваля имеет вид:
11/(8)|2а^2 = У 1 0 | (а,Ь,в)|2П и так же, как
у ВП, это соотношение свидетельствует об устойчивости (4) и о возможности представить произвольную функцию в виде разложения по ОВП. Дискретный аналог соотношения (4) и соответствующего ему условия Парсеваля проанализированы Донохо (БоиоИо) [22].
8) Математическая форма преобразования Радона разработана в 30-х годах, но только с 80-90-х годов оно нашло приложение в области эмиссионной (гамма- и так называемой позитронной) томографии [23].
МОДИФИЦИРОВАННАЯ ФОРМА ВЕЙВЛЕТ-ПРЕОБРАЗОВАНИЯ С ПЕРЕМЕННЫМ МАСШТАБИРОВАНИЕМ
Введенная выше форма МВП предусматривает локальное использование этого преобразования с фиксированной длиной, представление о которой (или ее грубую оценку) дает 2 5(5 — фиксированное). Теперь проанализируем, как ввести меняющийся масштаб у МВП и получить новую форму набора элементов (базиса) этого преобразования.
Элементы формализации
Если ввести в ^рассмотрение квадрат Q с размером стороны 1/2 ("с масштабом 5м), для которого целые числа к1 и к2 фиксируют его положение:
Q = [к /25,(к +1) /25]х[к2/Т,(к +1) /25], (7)
то удобно обозначить набор (множество) всех таких квадратов масштаба 5 как Qs, т. е.
Q(
(при любых к! и к2) ^
Q5.
Используем некоторую весовую функцию удовлетворяющую соотношению9)
w,
X w2(5l -к1,52 - к) = 1. (8)
к1, к2
Перенесем ее на стороны квадрата Q е Qs и масштабируем в соответствии с их размером (т. е. сожмем функцию w в 25 раз), вследствие чего она будет удовлетворять условию, эквивалентному (8), на сторонах такого квадрата и являться функцией распределения (единицы) энергии. Для этой масштабируемой весовой функции используем обозначение WQ, где Qе Qs 10).
В качестве следующего элемента формализации примем оператор TQ О е Qs) сдвига и масштабирования. Действие оператора TQ на функцию £(8) определяется в виде
TQg 00 = (Тэ § )(51, 52) = =(2> -к,2Ч -к).
(9)
Результат выполнения оператора TQ аналогичен
действию процедуры сдвига и растяжения (при 5 < 0) в вейвлет-преобразовании.
Локальное модифицированное вейвлет-преобразование
При использовании обозначений, отраженных в (7), (9), и компактного изображения интеграла двух функций в виде скалярного произведения (*,*), как это показано в (3), можно выразить представление функции / • WQ через ОВП в виде
/(Э) • WQ =
= ПК/(8')• w^ • ъ,в(*'))х
„ , Аа ли dв х Те^а, ъ,в00—М —• ^ а 4п
(10)
Суммирование равенства (10) с весом Wд по квадратам масштаба 5 приводит к выражению (11), поскольку весовая функция wQ, согласно (8), представляет разложение энергии:
/(*) = X / • wQ =
QеQ5
= ХШК/(*')• we • ^ъА*')>
Я. I о — о
„ ^а dв
• Тд¥а,ъdЪ —
и ^ • • а 4п
х
(11)
Из соотношения (11) видно, что любая функция может быть представлена в виде суперпозиции элементов, имеющих форму wQ • Т0 у/а ъ 6(*), т. е. в виде элементов, локализованных вблизи квадратов Q. Преобразования wQ • Т^аЪ 6(в) являются МВП
¥аеъеве (в) с параметрами а,ъ,в, величина которых видоизменена действием оператора TQ:
ад = 2-5 а, вд = 6,
ЪQ = ъ + к12-5 008 в + к2 2-5 8т 6.
(12)
9) В приложениях по обработке сигналов ее иногда называют весовой функцией распределения (единицы) энергии.
10) В соотношении (8) и далее по умолчанию предполагается, что Q е Qs, т. е. Q принадлежит множеству квадратов типа (7) — квадратов масштаба 5 (размера 1/ 25 х1/25), левая верхняя вершина которых имеет координаты (к / 25; к21/ 25) при различных целых значениях к и к2 .
Таким образом, wQ • ТдУаЪв(*) — это МВП,
сформированное весовой функцией ("окном") WQ
и имеющее область определения вблизи квадрата Q. Такую форму МВП будем называть локальным МВП (Л_МВП).
Рассмотренная форма локального МВП имеет определенную "длину" (область определения), которая характеризуется значением масштабного
коэффициента 2 я (при фиксированной величине 5). Если, оставляя локализацию вблизи квадрата Q, предусмотреть изменение масштаба 5 и, значит, изменение эффективной длины ЛВП, то можно получить набор (множество) ЛВП с различными масштабами:
{\м: 5 > 50,Qе Q,, а > 0, Ь е Я, ве[0,2п)}, (13) где \Ь в определено соотношением в =
= WQ ■ ,Ь,в..
Набор Л_МВП является своеобразной пирамидальной структурой, компонентами которой служат локальные (локализованные) МВП, характеризуемые различным расположением и длиной. С точки зрения представления простых достаточно гладких функций набор локализованных МВП со многими различными значениями масштаба является избыточным, поэтому разложение типа (6) с использованием Л_МВП нецелесообразно. Однако, как показано ниже, для двухмерных сигналов — изображений использование Л_МВП способствует более точному воспроизведению прямолинейных элементов повышенной интенсивности.
Вычисление МВП
Вычислительные процедуры, связанные с получением МВП, целесообразно рассмотреть сначала для непрерывной формы МВП (проанализированной выше), а затем для дискретной формы МВП, включая особенности ее получения.
Получение непрерывной формы МВП. Для преобразования Радона функции /(г,в) двух переменных справедливо соотношение [23]
/(Д cos sin =
1
>/2Л
] R f (t в)
e-at dt,
(14)
где f — ПФ, которое выражает двухмерное ПФ функции f через ПФ преобразования Радона R f (t,в) этой функции. Следовательно, R f (t,в) может быть выражено обратным (одномерным) ПФ для /(Д cos в, Д sin в) при фиксированном значении 9, т. е. на линии, проходящей через центр под углом 9 к оси х:
R f (t, в) = Í /(Д cos в, Д sin в) • eatdt. (15) V2n jL
Согласно (6), коэффициенты Rf(a,b,e) разложения f по МВП получаются вейвлет-преобразова-нием R f (t,в).
Таким образом, вычислительная процедура получения непрерывного МВП функции f включает:
1) получение преобразования Радона Я / (г ,в) с использованием двухмерного ПФ функции /;
2) вычисление коэффициентов Я/(а,Ь,в) представления функции / в виде разложения ее по МПФ путем выполнения ВП Я /(г, в).
При обработке цифровых данных приближенное значение преобразования Радона основывается на использовании алгоритма быстрого ПФ (БПФ). Для функции /, заданной значениями в дискретных точках квадрата размера пхп, т. е. /(¡1,/2), 0 < /1, ?2 < п — 1, последовательно выполняются следующие шаги.
- Двухмерное БПФ функции / которое дает
массив ( /(к^к2)), — п /2 <к1,к2 < п /2—1.
- Переход от прямоугольной системы координат к полярной. Величины БПФ(/), полученные на квадратной сетке, подставляются в качестве значений / на полярной сетке, т. е. на сетке, точки которой лежат на линиях, проходящих через начало координат.
- Обратное (одномерное) БПФ на каждой линии, т. е. для каждой величины углового параметр а в.
Цифровые данные в полярных (квазиполярных) координатах. Для реализации превращения прямоугольной системы в полярную оказывается достаточным использовать квазиполярную сетку. В этой "прямоугольно-полярной" сетке квазирадиальная переменная имеет набор уровней, которые определяются квадратами, а не окружностями (рис. 1). Формирование "прямоугольно-полярной" сетки и соответствующих ей координат может иллюстрироваться следующим образом. Выбирается в частотной плоскости изображения (размера пхп ) несколько радиальных линий, которые получаются связью начала координат с вертикалями, проходящими через точку (к1, к2) и лежащими на границе двухмерного массива (кь к2}, т. е. такими, что к\ или к2 е{ — п / 2, п / 2}. Квазиполярная сетка д1т, где I — индекс данной радиальной линии, а т отмечает положение точки (сетки) на этой линии, образуется точками пересечения между набором радиальных линий и линий, параллельных осям прямоугольной системы координат. Происходит это так, что точки отсчетов вдоль радиальной линии £, угол которой с вертикальной осью не более, чем п/ 4 (т. е. <п/ 4), получаются пересечением £ с набором горизонтальных линий {х2 = к2, к2 =— п/2, — п/2 +1,..., п/2}. Аналогично, пересечение с вертикальными линиями {х1 = к1, к1 =— п/2, — п/2 +1,..., п/2} определяет точки отсчетов в том месте, где угол между £ и горизонтальной осью не более, чем п / 4 .
Рис. 1. Квазиполярная сетка в частотной области для двухкомпо-нентного сигнала — изображения размером пхп (п = 8). Набор проходящих через центр линий (направленных по радиус-вектору) соединяет пары симметричных точек, лежащих на границе квадрата. Квазиполярная сетка — это набор точек (отмеченных кружочками) на пересечении "радиальных линий" с линиями, которые параллельны прямоугольным осям
-4 -3
2
0
2
З
4
Размер квазиполярной сетки равен 2n2, поскольку получается 2n радиальных линий L и n величин отсчетов на каждой из этих линий. В результате структуры данных, связанных с такой сеткой, будут иметь прямоугольный формат. Нетрудно заметить, что сделанный выбор квазиполярной сетки соответствует неравномерным интервалам между значениями угловой переменной 9.
Интерполирование координат к квазиполярной сетке. Для получения координат в квазиполярной сетке следует осуществить интерполяцию по отсчетам в прямоугольной системе координат. Интерполирование на уровне ПФ затруднительно, т. к. для двухмерного сигнала, в том числе для изображения, ПФ имеет сильно осциллирующий характер. Альтернативой такого интерполирования является "грубое" интерполирование, которое
основано на замене отсчетов f(^lm) ПФ сигнала в каждой точке полярной сетки на значение, которое является ближайшим к точке в точке прямоугольной сетке. Старком, Кандесом и Донохо (Stark, Candes, Donoho [24]) показано, что обратное "грубое" интерполирование (восстановление дискретных отсчетов ПФ на прямоугольной сетке по их значениям на полярной сетке) устойчиво и
обеспечивает точность не более половины размера прямоугольной сетки.
Несмотря на упрощенность "грубого" интерполирования оно дает достаточно хорошие результаты при анализе изображения с использованием МВП. Некоторые искажения получаются только в области высокочастотных компонент изображения (на его границах и в области резких перепадов интенсивности). Ниже показано, что такие искажения устраняются при использовании расширенной модификации ВП.
Одномерное вейвлет-преобразование. Для завершения МВП требуется выполнить ВП вдоль радиальной переменной в пространстве преобразования Радона. Выбор типа ВП, подходящего для этой цели, не произволен. Использование ВП, базис которого имеет компактное ПФ, ведет к заметным искажениям изображения, особенно при нелинейной форме обработки сигнала. Кроме того, это ВП обладает плохой временной локализацией, а флюктуации вейвлет-коэффициентов грубого масштаба порождают флюктуации коэффициентов мелкого масштаба, которые нежелательны в анализируемой задаче.
Однако в задаче рассматривается частотный подход, где дискретное ПФ восстанавливается из обратного преобразования Радона. Это обстоя-
тельство приводит к использованию вейвлет-базиса, ограниченного в частотной области, т. е. с компактной областью определения ПФ, а не временной. Поэтому в соответствии с анализом Старка и Муртафа ^агок, Murtagh) [25] для нашей задачи использовано ВП на основе масштабирующей функции р с ПФ, равным нулю вне интервала [-Ус ус ], т. е. (р = 0 при У£[-Ус ус ]. ПФ масштабирующей функции (р определено в виде (нормализованного) сплайна В3 [26]
(((V) = 3 Вз(4к), (16)
а ПФ вейвлет-функции у/ определяется как разность двух масштабирующих функций последовательных уровней разрешения:
у(2у) = р(у)-р(2у). (17)
Поскольку у/у) имеет компактную область определения, то можно построить набор 2п элементов (п + п/2 + ... +1) пирамидальной структуры [25]11). Такое вейвлет-преобразование имеет довольно привлекательные характеристики.
- Вейвлет-коэффициенты вычисляются непосредственно для частотной формы, т. е. в области ФП. Применительно к МВП это позволяет избежать вычисления (одномерного) обратного ПФ вдоль каждой радиальной линии квазиполярной сетки.
- Отсчеты в каждом частотном диапазоне получаются с частотой выше частоты Найквиста— Котельникова (ЧНК), поэтому практически исчезает наложение частот (так называемое явление элайзинга) — явление, которое встречается при ВП, когда частота дискретизации не превышает ЧНК [27].
- Процедура восстановления сигнала проста: чтобы восстановить входной сигнал в данной точке достаточно просто совместно сложить вейвлет-коэффициенты. В рассматриваемой нами концепции для того, чтобы восстановить коэффициенты ФП требуется совместно сложить коэффициенты МВП.
Анализируемая форма ВП вносит фактор значительной избыточности. С точки зрения привычных требований ортогональности ВП и использования дискретизации на уровне минимально приемлемой частоты — ЧНК — эта избыточность может трактоваться как некоторый недостаток. Однако в рассматриваемом приложении целью является не сжатие данных или эффективное коди-
11) В [25] авторами обоснован выбор ВП с аналогичными частотными свойствами в задаче для улучшения изображения с использованием интерферометрического метода.
рование (для которых минимальная частота дискретизации желательна), а удаление шума, для которого избыточность может обеспечить существенное преимущество [15].
Промежуточные итоги и представление МВП потоковым графом. Поскольку процедура получения МВП состоит из последовательности этапов — элементов преобразования, каждое из которых обратимо, точно также обратимо и МВП. При этом обратное (относительно МВП) преобразование обладает свойством точного восстановления сигнала, подлежащего МВП, и процедура восстановления устойчива относительно небольших искажений коэффициентов МВП.
На рис. 2 показан потоковый граф, отражающий этапы реализации МВП. МВП изображения размера п х п является изображением 2п х 2п , так что вводится избыточность с коэффициентом 4. С вычислительной точки зрения алгоритм получения МВП предъявляет достаточно умеренные требования. При увеличении области обрабатываемого изображения (или области определения другого двухмерного сигнала) до размера п х п количество арифметических действий, необходимых для получения МВП, возрастает пропорционально величине п 2^(п).
Цифровая и дискретная форма МВП. Описанные выше процедуры, связанные с численным выполнением МВП, показывают некоторое расхождение с общей теоретической концепцией этого преобразования. Этому целесообразно дать расширенную трактовку.
В пирамидальной структуре МВП количество ориентаций радиальных линий постоянно (не зависит от масштаба), в то время как, согласно теоретической концепции (ТК), количество ориентаций обратно пропорционально масштабу [28]. Так что по ТК на наборе ориентаций требуется уменьшение частоты отсчетов вдвое каждый раз при переходе к большему, более грубому, масштабу. При используемой процедуре выполнения МВП число отсчетов угловой переменной 9 увеличивается вдвое (сравнительно с ЧНК) каждый раз при переходе к более грубому масштабу, и это обеспечивает соответствующую избыточность в численном представлении МВП.
По этому поводу уместно отметить аналогию с алгоритмом ВП для удаления шума, сопутствующего сигналу. По принятой методологии для большинства задач в пирамидальной структуре ортонормированных вейвлетов число отсчетов для определенного масштаба и локализации (вейвлета) является фиксированным на уровне ЧНК и не зависящим от масштаба. Но для ВП, применяемого для удаления шума, этот принцип не используется: при переходе к большему значению параметра
Изображение Двухмерное преобразование Фурье П ' А // ФТУ/у'
УС/Х/1 / ху-л-д \| \\\\ "Т V уХ
Обратное ПФ (одномерное)
МВП
Рис. 2. Потоковый граф МВП. Каждая из радиальных линий в области ПФ обрабатывается отдельно. Одномерное обратное ПФ вычисляется вдоль каждой радиальной линии и после этого следует (одномерное) неортогональное ВП. Практически коэффициенты одномерного ВП могут вычисляться непосредственно в области ПФ
-►
Частота
масштаба получается возрастающее количество отсчетов для определенного масштаба и локализации. При этом (для задачи удаления шума) показано преимущество такой схемы использования ВП.
Поскольку основным приложением МВП является обработка двухмерного сигнала (и особенно изображений) с целью удаления шума и сохранения всех особенностей (двухмерного) спектрального рельефа сигнала или особенностей контраста изображения, включая контуры ("хребты" контрастности), то есть определенное преимущество в том, что при возрастании параметра масштаба получается большее число радиальных ориентаций, чем это необходимо.
Локальные МВП (гладко распределенные элементы МВП). С помощью рассмотренного выше цифрового варианта МВП исходное изображение размера п х п представляется гладко накладывающимися блоками с длиной стороны Ь пикселей. При этом наложение двух соседних вертикальных блоков образует прямоугольник размера Ь х Ь /2 12). Поэтому для изображения размера п х п имеется 2п / Ь таких блоков в каждом из двух направлений.
12) Такое наложение устраняет возможные искажения,
которые могли бы появиться при непосредственной стыковке соседних блоков.
Такое разделение изображения на блоки вносит избыточность, т. к. каждый пиксель изображения принадлежит четырем соседним блокам.
Возможны два варианта выполнения анализа и синтеза изображения.
- На стадии анализа величинам (значений интенсивности в пикселях) придаются определенные веса. Эти веса таковы, что совместное сложение взвешенных отсчетов во всех блоках (на стадии синтеза) точно восстанавливает значения в пикселях исходного изображения.
- Величины в блоках соответствуют величинам в пикселях исходного изображения (стадия анализа), но затем для восстановления изображения (на стадии синтеза) они взвешиваются.
К анализу и синтезу могут быть применены и промежуточные подходы, например с использованием скользящего окна, упоминавшегося в п. "Элементы формализации". В первом подходе данные сглаживаются окном на стадии синтеза. К преимуществам этого подхода относится устранение возможных искажений, обычно появляющихся в области границы изображения; к недостаткам относится небольшое снижение чувствительности. Так, например, если вертикальная линия с интенсивностью Ь пересекает блок размера Ь и дисперсия шума равна 1, то когда угловой параметр 9 преобразования Радона соответствует положению линии Ь, будет получено измерение с интенсив-
ностью сигнала ъ • Ь , СКО шума равно л/Ь , а отношение сигнал/шум (ОСШ) будет Если взвешивание сглаживающим окном производится на стадии анализа, то ОСШ будет
Ь Xъ=1 .
Потеря чувствительности может быть существенной особенно в приложениях по фильтрации, поэтому второй подход лучше подходит к задаче восстановления изображения с подавлением шума.
Для вычисления величины пикселей /(у) из значений, находящихся в смежных четырех блоках половинного размера 1/2 (а именно в блоках
^^ jl), B2(i2,jl), Вз(^Л) и В4(^2,32), где ^/1 > > ъ /2 и i2 = \ -1, /2 = /1 -1) могут быть использованы соотношения:
/1 = w(i2/1) • В1 (il, /1) + w(1 - ц/1) • В1 (^, /1);
/2 = w(i2/1) • Въ(^, /2) + w (1 - 1) • В4О2, /2); (18)
/ 0', /) = w(i2/1) • /1 + w(1 - /2/1) • /2,
где w( х) = 0082 (п • х /2).
Можно было бы выбрать и другую невозрас-тающую функцию, удовлетворяющую условию w(0) = 1, w(1) = 0, w'(0) и подчиняющуюся свойству симметрии w( х) + w(1 - х) = 1.13)
Как отмечалось выше, используемый способ пространственного блочного разделения изображения вводит избыточность с коэффициентом, равным четырем.
РАСШИРЕННАЯ МОДИФИКАЦИЯ ВП И ЕЕ ДИСКРЕТНАЯ ФОРМА
Расширенная модификация ВП (РМВП) ориентирована на получение наиболее точного восстановления изображения после реализации процедуры снижения уровня шума (искажений, присущих первичной регистрации изображения). При этом, по методологии, описанной в [5], снижение шума
13) Для того чтобы было более полное соответствие с теоретической формой концепции МВП, следовало бы произвести нормализацию с коэффициентом, зависящим от размера блоков. Однако при использовании МВП для задачи снижения уровня шума и применения фильтрации в пространстве МВП (аналогично фильтрации в области ВП [5]) путем установления уровней дискриминации коэффициентов МВП нормализация не является отдельной задачей. Введение нормализации коэффициентов МВП автоматически ведет к точно такой же нормализации соответствующих порогов дискриминации коэффициентов, используемой при фильтрации изображения в области МВП. Дополнительно этот вопрос рассмотрен в следующем разделе.
осуществляется методом фильтрации в области преобразования, т. е. дискриминацией коэффициентов (с использованием жесткого или гибкого порогов) в разложении сигнала по РМВП14).
Как отмечено выше, рассмотренная пирамидальная структура МВП с различными масштабами является избыточной формой представления сигнала. В этом случае трудно получить методом пороговой фильтрации разложение с малым количеством коэффициентов, отражающих информационно значимую часть сигнала с исключенным шумом. РМВП отличается от МВП устранением избыточности в разложении сигнала за счет уменьшения количества масштабов. Для этого вводится ортогональность между элементами МВП с различными масштабами, которая достигается субполосной фильтрацией. Уровни с различными масштабами в пирамидальной структуре МВП используются для представления субполосных входов системы согласованного многоканального разделения (системы СМР [29]). В то же время это СМР налагает определенное соотношение между шириной и длиной элементов системы МВП так, что они становятся неоднородными и удовлетворяют соотношению ширина = (длина)2.
Дискретное РМВП функции /(х1, х2 ) использует двоичную последовательность значений масштаба и систему СМР (Р0/,Д1 /,А2/,...) со свойством, что полосовой фильтр А5 концентрируется в области интервала частот [225, 225+2]. Это свойство выражается соотношением
А5 = У */ , у5(4) =у(2-254). (19)
В (обычном) ВП используется разложение по элементам базиса, у которых частотное представление принадлежит субполосам с "двоично орга-низиванными диапазонами" (ДОД) [25 , 25+1]. В дискретной форме РМВП используются субполосы с нестандартными ДОД: [225,225+2]. Эту особенность следует учитывать при рассмотрении этапов формирования РМВП.
При обозначениях, которые были использованы выше (в разделе "Вычисление МВП"), алгоритмическая схема построения РМВП для функции / включает следующие этапы.
- Разложение / на субполосы:
/ Ы-(Р0/,А/, Д2/,...).
- Гладкое разделение: каждая субполоса гладко разделяется окнами на "квадраты" подходящего масштаба (со стороной 2-5)
14) В [5] развит метод фильтрации в области ВП, и он использован для эффективного снижения уровня шума для нестационарного сигнала (речи).
A J ^ ( WQ A J W •
- Ренормализация: каждый полученный квадрат ренормализуется к единичному масштабу
gQ = (Tq )-1( Wq A J) QeQs .
- Представление функции на каждом квадрате разложением в формате дискретного МВП.
При введенной последовательности процедур, определяющих концепцию РМВП, до выполнения дискретного МВП две субполосы [ 22s, 22s+1 ] и [22s+1, 22s+2] объединяются.
Дискретная форма РМВП и ее реализация. При рассмотренном ранее выполнении МВП для функции f, заданной в дискретной форме на квадрате размера n х n, осуществлялась последовательная (на каждом этапе) замена непрерывной концепции на дискретную концепцию. Аналогичную логику перехода к дискретной форме целесообразно использовать также и для РМВП. В основном этот переход воспроизводит этапы, связанные с получением дискретного МВП. Однако предварительное опробование РМВП показывает, что слияние двух субполос [ 22s, 22s+1 ] и [22s+1, 22s+2] (которое предусмотрено в предложенной выше последовательности этапов формирования расширенного МВП) дает несколько худшее восстановление изображения, чем то которое получается при раздельном использовании субполос, т. е. при делении на квадраты и применении МВП в каждой субполосе отдельно.
При формировании РМВП используется субполосная фильтрация по методу СМР [8], хорошо адаптирующемуся к схеме алгоритма получения дискретного вида преобразования РМВП. С помощью алгоритма изображение I размера n х n (или другого двухмерного сигнала) представляется совокупностью выражений вида (20), каждое из которых справедливо для отдельного пикселя:
j
I(X, y) = Cj (x, y) + £ Wj (x, y), (20)
j=i
где Wj представляет парциальную часть изображения I при масштабе 21 , а cJ — остаточная часть, не вошедшая в первые J масштабов разложения изображения [25].
Таким образом, выходом алгоритма являются массивы размерности n х n, каждый из которых соответствует одной из субполос (j = 1, ..., J) или остаточной части cJ (x, y) изображения, не вошедшей в первые J масштабов.
Структура алгоритма получения дискретной формы РМВП схематично может быть представлена следующим образом.
- Применяется процедура субполосной фильтрации (по методу СМР) с J масштабами.
- Устанавливается Я = B„t„.
1 мин
- Для значений индекса j = 1,...,J выполняются операции:
1) субполоса Wj разделяется на блоки размера Bj, и к каждому блоку применяется МВП;
2) для четных значений j полагается Bj+1 = Bj,
а для нечетных значений j полагается Bj+1 = 2 Bj .15)
Длина стороны локализующих окон удваивается при каждой следующей субполосе, и таким образом поддерживается основное свойство РМВП: элементы с величиной длины 2-3 2 служат для анализа и синтеза j-й субполосы [ 2j, 2j+1 ]. При этом в разложении (20) начальные значения индекса j соответствуют более тонкому масштабу, а "остаток" разложения cJ не обрабатывается. Схема организации описанного алгоритма показана на рис. 3.
Выполнение РМВП обладает некоторой избыточностью. Если при разложении функции f (i, j) используется J масштабов, то коэффициент избыточности составляет 16 J+1. Вместе с тем, достоинством метода РМВП является устойчивость преобразования и возможность точного восстановления изображения. Это обусловлено обратимостью всех процедур, участвующих в алгоритме формирования РМВП.
Возможности использования РМВП для удаления шума изображения
Применительно к изображению шумом являются искажения, связанные с несовершенством регистрации изображения (неоднородностью и частичной непрозрачностью среды, микродвижением объекта, погрешностью оптической фокусировки). Эффективное снижение уровня шума возможно только при его низком уровне. Вследствие многофакторного формирования шума без преобладания вклада отдельных его компонент в качестве модели зашумленного изображения можно использовать соотношение
xhJ = f(i, j) + zhJ, (21)
где f — изображение, которое требуется восстановить; z,j — шум с независимыми нормально-распределенными отсчетами (zi j е N(0,1) ); а — параметр, определяющий уровень шума.
15) Индексирование удобно производить так, что } = 1
соответствует наиболее "тонкому" масштабу — высоким частотам (т. е. большому значению параметра масштаба 5).
Частота
Рис. 3. Потоковый граф, схематично представляющий алгоритм формирования РМВП
В отличие от ПФ и ВП, построенные выше дискретные формы преобразований (МВП и РМВП) не сохраняют норму, первоначально введенную для функции / . Поэтому дисперсия коэффициентов разложения зашумленной функции / зависит от индекса X у МВП (и соответственно у РМВП). Так, если F означает матрицу дискретной формы РМВП, то после преобразования компонента шума в изображении Fz будет иметь нормальное (гауссово) распределение со средним значением, равным FFT, т. е. F ге N (0,FFT)1б).
Удаление шума в зарегистрированном изображении производится с помощью процедуры дискриминации коэффициентов РМВП в соответствующем представлении первоначально полученного изображения 17). Если ул — коэффициенты РМВП (у = F х , а оператором F для краткости символически обозначено дискретное РМВП зарегистрированного изображения х, определяемого соотношением (21)), то процедура дискриминации коэффициентов "с жестким" порогом определяется выражением:
16) Поскольку произведение двух матриц ПФ при высокой их размерности (в случае изображения) представляет определенные трудности, то можно просто использовать приближенные величины дисперсий о\, получаемее моделированием (по методу Монте-Карло), причем диагональные элементы матрицы FFT оцениваются с помощью вычисления РМВП нескольких изображений, состоящих из стандартного белого шума.
17) Процедура дискриминации коэффициентов разложения функции, которое основано на вейвлет-преобра-зовании (ВП), введена Донохо и Джонстоном [14], [31]. В таком первоначальном виде она вошла в Toolbox системы MATLAB. Позднее эта процедура получила существенное развитие в [5], где даны уточненный вариант процедуры Донохо и расширенное "семейство" вероятностных распределений для вейвлет-коэффици-ентов — коэффициентов разложения функции при ВП.
Ул
|y,|
у л пРи — ^ к ал; а
0 при
|Ул
(22)
<к а.
а
где ал — оценка дисперсии коэффициентов РМВП изображения, которое модифицируется шумом.
Процедура дискриминации (коэффициентов РМВП изображения) "с гибким порогом" совпадает с описанной в [5] процедурой дискриминации коэффициентов при дискретном ВП.
По предварительным оценкам в рассматриваемой методике снижения шума в зарегистрированном изображении целесообразно использовать для к в (22) некоторую зависимость от масштаба. Удовлетворительные результаты дает использование к = 4 для первого масштаба (j = 1) и к = 3 — для других масштабов (j > 1).
Для оценки степени снижения шума изображения может быть рекомендован метод, который включает несколько этапов.
- Получение цифровой формы 2-3 изображений, представленных на сайтах [31] или в системе MATLAB [32]. Размер изображений выбирается в соответствии с необходимым количеством масштабов, используемых в применяемом преобразовании — МВП, РМВП или ВП. Размер изображения 512х512 позволяет при анализе использовать 4 масштаба (размера блоков) — 8, 16, 32 и 64. (Изображения Lena и Barbara используются наиболее часто и служат своеобразным стандартом при сопоставлении методов сжатия и фильтрации).
- Введение зашумления изображения: к изображению добавляется гауссов шум с величиной среднеквадратичного отклонения 10-20.
- Введение (в выбранные для анализа изображения) гауссова шума с двумя или тремя значениями интенсивности. Среднеквадратичное значение шума выбирается в пределах 10-20.
- Выполнение РМВП изображения для четырех масштабов (для размера блоков 8, 16, 32, 64) и осуществление фильтрации в области РМВП: применение описанной выше схемы дискриминации коэффициентов в представлении изображения разложением по (дискретному) РМВП.
- Для получения сравнительной оценки снижения шума изображения и качества изображения, полученного после фильтрации с "жестким" (или "гибким") порогом дискриминации коэффициентов, кроме РМВП следует выполнить МВП и ВП.
- В качестве объективной меры снижения шума изображения может служить отношение сигнал/шум (по мощности) на выходе процедуры фильтрации (дискриминации) коэффициентов разложения изображения по соответствующему пре-
образованию. Основой вспомогательного критерия для сопоставления качества восстановления изображения после фильтрации шума является способность зрительной идентификации остаточных артефактов, которые не имеют достаточно хорошего количественного отражения в величине отношения сигнал/шум. Так, ВП (без прореживания по частоте) после фильтрации дает остаточное отношение сигнал/шум, близкое к значению этого отношения после фильтрации изображения в области РМВП, однако появление зрительно заметных искажений контуров при использовании ВП и потеря некоторых деталей свидетельствует о преимуществе РМВП.
Таким образом, предварительное изучение возможностей РМВП для снижения шума изображения позволяет считать, что это преобразование имеет определенное преимущество как перед используемым для этой цели ВП, так и МВП. При восстановлении изображения после фильтрации шума ВП (как отмечалось и ранее [25], [33]) имеются искажения в области границ изображения: гладкие элементы кривых (контуры) часто превращаются в прерывистые плохо определенные фрагменты, а некоторые хорошо восстанавливаемые участки пограничных контуров приобретают колебательную структуру. Дополнительно появляются достаточно многочисленные вкрапления ("пятнышки") малого масштаба, а отсечение их повышением порога дискриминации коэффициентов приводит к потере еще большего количества фрагментов изображения.
Преимущество РМВП перед ВП и МВП наблюдается при восстановлении изображения после фильтрации шума различного уровня (со среднеквадратичным значением шума в пределах 10-20).
ЗАКЛЮЧЕНИЕ
Рассмотрено модифицированное вейвлет-пре-образование (МВП), построение которого основано на некоторой аналогии с обобщенной формой преобразования Фурье (ПФ), использовании аппарата преобразования Радона и гибкого способа масштабирования полученной МВП. Модификация вейвлет-преобразования (и особенно ее вариант с дополнительным масштабированием — расширенная модификация вейвлет-преобразования (РМВП)), имеет определенное преимущество перед традиционной формой вейвлет-преобразо-вания (ВП). Это относится к задаче увеличения информационно значимой компоненты сигнала за счет снижения шума применительно к анализу двумерного сигнала — изображения, которое искажено сопутствующим шумом (наличием так называемого блера).
- Введена концепция модифицированной фор-
мы ВП. Прототип МВП — ВП, у которого аргумент вейвлет-функции суммируется из проекций компонент (51, 52) двухмерного сигнала на полярную ось с углом в, т. е. ¥аьв(%) =
= а-12-^(^0080 +s2sinв-Ь)/a ). При значении углов в = 0 и в = п /2 у/а ь 0 (8) обращается в
вейвлет-функцию одной из компонент вектор-сигнала 8.
- Двухмерная функция (или изображение) представима разложением по МВП, которое обратимо и позволяет восстановить без искажений начальную функцию. Коэффициенты МВП выражаются в виде ВП коэффициентов преобразования Радона.
- Вычислительные процедуры, связанные с получением МВП, проанализированы для непрерывной и дискретной формы. Для дискретной формы рассмотрен переход к квазиполярной системе координат (ее отсчеты принадлежат точкам пересечения полярной оси с концентрическими квадратами).
- Рассмотрена методология построения МВП с изменяющимся параметром масштаба — расширенной модификации вейвлет-преобразования (РМВП).
- Основным направлением приложений МВП и РМВП является решение задачи удаления шума, присущего первично зарегистрированному изображению. Концепция этих преобразований предусматривает использование набора (пирамидальной структуры) масштабов аналогично ВП, но разложение функции (изображения) с помощью МВП и РМВП дает существенно меньшее число коэффициентов. Поэтому при фильтрации шума изображения методом дискриминации коэффициентов разложения МВП и РМВП имеют значительное преимущество перед ВП.
- Преимущество РМВП перед ВП и МВП наблюдается при восстановлении изображения после фильтрации шума различного уровня (со средне -квадратичным значением шума в пределах 10-20).
- Идея новых преобразований, эффективно повышающих качество изображения, и метод перехода к полярной системе координат с неравномерной градацией углового параметра разрабатывались Донохо, Старком и Кандесом, Мухтафой, Койфманом и Джонстном.
СПИСОК ЛИТЕРАТУРЫ
1. Воробьев В.И., Грибунин В.Г. Теория и практика вейвлет. СПб.: Изд. ВУС, 1999. 203 с.
2. Добеши И. Десять лекций по вейвлетам. Москва—Ижевск: РХД, 2001. 460 с.
3. Фрейзер М. Введение в вейвлеты в свете линейной алгебры. М.: Бином, 2008. 487 с.
4. Меркушева А. В. Классы преобразований нестационарного сигнала в информационно-измерительных системах. III. Время-масштабные (вейвлет-) преобразования для спектрально-временного анализа // Научное приборостроение. 2002. Т. 12, № 3. С. 68-82.
5. Меркушева А.В. Фильтрация нестационарного сигнала (речи) в вейвлет-области с адаптацией к виду и динамике шума // Научное приборостроение. 2003. Т. 13, № 2. С. 73-87.
6. Меркушева А.В. Скремблирование на основе преобразований вейвлет и Адамара // Проблемы информационной безопасности. 2000. № 1. С.86-92.
7. Кратиров Д.В., Меркушева А.В. Алгоритм, основанный на вейвлет-преобразовании и нейронной сети, для бесконтактного измерения параметров газожидкостного потока // Датчики и системы: Сборник докладов международной конференции, том 3. СПб.: Изд. СПб ГПУ, 2002. С. 51-55.
8. Малыхина Г.Ф., Меркушева А.В. Согласованное многоканальное разделение сигнала с вейвлет-преобразованием и модифицированной полифазной структурой для спектрального вычитания фона // Научное приборостроение. 2008. Т. 18, № 1. С. 110-118.
9. Малыхина Г.Ф., Меркушева А.В. Расширенная характеристика динамического спектра нестационарного сигнала // Научное приборостроение. 2008. Т. 18, № 2. С. 66-72.
10. Pey S.Ch., Yeh M.-H. An Introduction to Discrete Finite Frames // IEEE Signal Processing Magazine. 1997.N 11. P.84-92.
11. Choi H., Williams W.J. Time-Frequency Transformation with Reduced Interference on the Basis of Exponential Kernel // Proceedings of IEEE. 1996. V. 84, N 9. P. 1264-1271.
12. Consalves P., Baraniuk R.G. Bertran's Bilinear Time-Frequency Transformations // Signal Processing Letters. 1996. V. 3, N 3. P. 82-84.
13. Stankovic L.J. Time-Frequency Distribution Concentrated along the Instantaneous Frequency // Signal Processing Letters. 1996. V. 3, N 3. P. 8992.
14. Donoho D.L. De-noising by Soft-Thresholding // IEEE Transactions on Information Theory. 1995. V. 41, N 3. P. 613-627.
15. Coifman R.R., Donoho D.L. Translation Invariant De-noising // Wavelets and Statistics / A. Anto-niodis, G. Oppenheimer ^ds.). N.Y.: Springer Verlag, 1995. P. 125-150.
16. Baraniuk R.J., Jones D.L. Unitary Equivalence: a New Twist on Signal Processing // IEEE Transactions on Signal Processing. 1995. V. 43, N 10. P.1251-1261.
17. Меркушева А.В. Обобщенная модификация преобразования Фурье для обработки неста-
ционарных сигналов в ИИС // Информационные технологии. 2006. № 5. С. 62-72.
18. Новиков Л.В. Аппаратно-ориентированные вейвлеты и их применение для обработки данных. Автореф. дис. ... д-ра физ.-мат. наук. СПб.: ИАнП РАН, 2007. 31 с.
19. Новиков Л. В. Модифицированные вейвлеты в обработке данных аналитических приборов // Научное приборостроение. 2006. Т. 16, № 1. С. 3-14.
20. Gaunaurd G.C., Stifors H.C. Signal Analysis by Means of Time-Frequency Transformation of Wigner Type // Proceedings of IEEE. 1996. V. 84, N 9. P. 1231-1247.
21. Меркушева А.В. Классы преобразований нестационарного сигнала в информационно-измерительных системах. I. Элементы теории // Научное приборостроение. 2002. Т. 12, № 2. С.50-58.
22. Donoho D.L. Orthogonal Ridgelets and Linear Singularities // Journal on Mathematical Analysis (SIAM). 2000. V. 31, N 5. P. 1062-1089
23. Федоров Г.А., Терещенко С.А. Вычислительная эмиссионная томография. М.: Энергоатом-издат, 1990. 182 с.
24. Starck J.-L., Candes E.J., Donoho D.L. The Curvelet Transform for Image De-noising // IEEE Transactions on Image Processing. 2002. V. 11, N 6. P. 669-680.
25. Starck J.-L., Murtagh F. A New Entropy Measure Based on Wavelet Transform and Noise Modeling // IEEE Transactions on Circuits and Systems II. 1998. V. 45, N 8. P. 1118-1124.
26. Чуи К.С. Введение в вейвлеты. М.: Мир, 2001. 412 с.
27. Simoncelli E.P., Freeman W.T., Adelson E.H. Shiftable Multiscale Transforms (What's wrong with orthogonal wavelets) // IEEE Transactions on Information Theory. 1992. V. 38, N 3. P. 587607.
28. Candes E.J. Harmonic Analysis of Nuclear Networks // Applied Computations and Harmonic Analysis. 1999. V. 6. P. 197-218.
29. Меркушева А.В., Малыхина Г.Ф. Согласованное многоканальное разделение сигнала: фильтрация и мультиплексирование // Научное приборостроение. 2008. Т. 18, № 1. С. 98-109.
30. Donoho D.L., Johnstone I.M. Ideal Spatial Adaptation by Wavelet Shrinkage // Biometrika. 1994. V. 81. P. 425-455.
31. (http://vismod.www.media.mit.rdu/vismod/imagin ery/Vision_Texture/vistex.html); (http://compression.graphicon.ru).
32. Дьяконов В., Круглов В. Математическое расширение MATLAB. СПб.: Питер, 2001. 480 с.
33. Crouze M., Novak R., Baraniuk R. Wavelet-Based Statistical Signal Prosessing Using Hidden Markov Models // IEEE Transactions on Signal Processing. 1998. V. 46. P. 886-902.
Санкт-Петербург
Материал поступил в редакцию 23.10.2008.
MODIFIED WAVELET-TRANSFORMATION. CONCEPTION AND APPLICATION
V. M. Malychin, G. F. Malychina, A. V. Merkusheva
Saint-Petersburg
The conception of modification for wavelet transformation (WT) is given. The construction is based on analogy with generalized form of Fourier transformation, on using the Radon transformation method and on flexible form of scaling for designed modified wavelet transformation (MWT). The MWT (and especially its mode with additional scaling) have superior performance over tradition WT. This relate to the problem of increasing informational-significant component of signal by means of lowering the noise level, including the problem of two dimensional signal analysis, — reconstructing the image that is distorted by concomitant noise (by existence of so termed blur).
Keywords: wavelet transformation, Radon transformation, modification, operators, local and extended forms, noise reduction.
