авторефераты диссертаций БЕСПЛАТНАЯ  БИБЛИОТЕКА

АВТОРЕФЕРАТЫ КАНДИДАТСКИХ, ДОКТОРСКИХ ДИССЕРТАЦИЙ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:   || 2 |
-- [ Страница 1 ] --

Министерство образования Российской Федерации

НОВОСИБИРСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ

УНИВЕРСИТЕТ

А.Н. ЯКОВЛЕВ

ВВЕДЕНИЕ

В ВЕЙВЛЕТ-ПРЕОБРАЗОВАНИЯ

Утверждено Редакционно-издательским советом

в качестве учебного пособия

НОВОСИБИРСК

2003

УДК 621.372(075.8)+004.92(075.8)

Я 474

Рецензенты:

Кафедра радиотехнических систем Сибирского государственного ун-та телекоммуникаций и информатики (СибГУТИ) Д-р техн. наук, проф. В.Н. Васюков (НГТУ), канд. техн. наук, доц. Г.Х. Гарсков (СибГУТИ) Работа выполнена на кафедре теоретических основ радиотехники (ТОР) Яковлев А.Н.

Я 474 Введение в вейвлет-преобразования: Учеб. пособие. – Новосибирск: Изд-во НГТУ, 2003. – 104 с.

ISBN 5-7782-0405- В сжатой и доступной форме рассмотрена новейшая технология об работки информации – непрерывное (глава 1), дискретное и быстрое вейвлет-преобразования (глава 2) сигналов, а также двумерное вейвлет преобразование и обработка изображений (глава 3). Приведены приме ры применения ВП для анализа, фильтрации и сжатия сигналов и изо бражений. В приложениях дан «инструментарий» по ВП в системе ком пьютерной математики Wavelet Toolbox для MATLAB 6.

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

УДК 621.372(075.8)+004.92(075.8) ISBN 5-7782-0405-1 © Новосибирский государственный технический университет, © Яковлев А. Н., 50-летию НГТУ (НЭТИ) посвящаю Автор Введение В конце прошлого века возникло и успешно развивается но вое и важное направление в теории и технике обработки сигна лов, изображений и временных рядов, получившее название вейвлет-преобразование (ВП), которое хорошо приспособлено для изучения структуры неоднородных процессов.

Термин вейвлет (wavelet) ввели в своей статье Гроссманн (Grossmann) и Морле (Morlet) в середине 80-х годов XX века в связи с анализом свойств сейсмических и акустических сигна лов. Их работа послужила началом интенсивного исследования вейвлетов в последующее десятилетие рядом ученых таких, как Добеши (Dobechies), Мейер (Meyer), Малл (Mallat), Фарж (Farge), Чуи (Chui) и др.

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

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

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

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

Подтверждением значимости ВП является и тот факт, что алгоритмы ВП представлены в широко распространенных сис темах компьютерной математики (СКМ), таких как Mathcad, MATLAB, Mathematica. Международные стандарты JPEC-200, MPEG-4 и графические программные средства Corel DRAW 9/10 широко используют ВП для обработки изображений и, в частности, для сжатия изображений для каналов с ограниченной пропускной способностью, например, для Интернет. Кроме то го, фирмой Analog Devices разработаны и выпускаются одно кристальные дешевые микропроцессоры ADV6xx (ADV601, ADV601LC, ADV611, ADV612), основанные на ВП и предна значенные для сжатия и восстановления видеоинформации в реальном масштабе времени.

В России первые немногочисленные работы по применению ВП были опубликованы примерно с десятилетней задержкой. В основном они носили обзорный характер по материалам зару бежных публикаций [15, 19, 20, 26].

В последние несколько лет интерес к ВП у нас резко возрос.

Появились учебные пособия [1–3], монография [4], переведены на русский фундаментальные теоретические книги И. Добеши [5] и Т. Чуи [6], вейвлетам посвящены разделы и главы в учеб никах [10–14]. Несмотря на приличную стоимость (от 100 до 300 рублей), книги мгновенно исчезли с прилавков магазинов.

Кроме малого тиража и высокой стоимости книги по вейвле там имеют существенные недостатки – слабую практическую реализацию ВП, отсутствие практических примеров и описания современного компьютерного «инструментария» для практи ческого использования ВП.

В справочной книге [7], написанной В.П. Дьяконовым и И.В. Абраменковой и изданной в 2002 г. (тираж 5000 экз., 38 п.л., цена 160 руб.), авторами сделана попытка преодолеть последние из указанных выше недостатков. Справочник содер жит сведения по современным средствам обработки и фильтра ции сигналов и изображений, входящим в пакет MATLAB 6.0/6.1. В последующей книге [8] (тираж 3000, 28 п.л., цена 270 руб.) впервые, наряду с теоретическими сведениями о вейв летах, детально описаны известные пакеты по вейвлетам – Wavelet Toolbox, Wavelet Extension Pack и Wavelet Explorer, ис пользуемые новейшими массовыми системами компьютерной математики (СКМ) MATLAB 6.0/6.1, Mathcad 2001 и Mathe matica 4. Однако пользователь должен иметь сведения и навыки по этим СКМ хотя бы в рамках книг, опубликованных В.П. Дья коновым в 2001–2002 гг.

В Санкт-Петербурге благодаря энтузиазму профессорско преподавательского состава СПбГУ, СПбГТУ и других вузов регулярно действует городской семинар «Wavelet Seminar – вейвлеты (всплески) и их приложение». На Интернет-сайте это го семинара можно получить сведения о русскоязычных публи кациях по вейвлетам, о зарубежных вейвлет-ресурсах, о новых книгах и просто полезные советы.

Защищен ряд диссертаций по теории и практике ВП и, в ча стности [22, 25].

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

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

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

Автор выражает глубокую благодарность рецензентам д-ру техн. наук проф. В.Н. Васюкову и канд. техн. наук. доц.

Г.Х. Гарскову за полезные советы и критические замечания, а также специалисту по программированию Ю.В. Морозову за консультации при запуске программ в системе MATLAB.

Глава Непрерывное вейвлет-преобразование 1.1. Обобщенный ряд Фурье Известно [13, 14], что произвольный сигнал S (t ), для ко торого выполняется условие t [ S (t )] dt, (1.1) t может быть представлен ортогональной системой функций {n (t )}, т.е.

Cnn (t ), S (t ) = C0 0 (t ) +.... + Cn n (t ) +... = (1.2) n= коэффициенты которой определяются из соотношения t S (t )n (t )dt, Cn = (1.3) n t t || n ||2 = n (t )dt где (1.4) t – квадрат нормы, или энергия базисной функции n (t ). При этом предполагается, что никакая из базисных функций не равна тождественно нулю и на интервале ортогональности [t1, t2 ] вы полняется условие t || n ||2, k = n, k (t )n (t )dt = 0, k n. (1.5) t Базисная функция n (t ), для которой квадрат нормы равен единице ( || n ||2 = 1 ), называется нормированной (нормальной), а вся система функций {n (t )} – ортонормированной или орто нормальной. В этом случае говорят, что задан ортонормирован ный базис.

Ряд (1.2), в котором коэффициенты Cn определяются по формуле (1.3), называется обобщенным рядом Фурье.

Произведение вида Cn n (t ), входящее в ряд (1.2), представ ляет собой спектральную составляющую сигнала S (t ), а сово купность коэффициентов {C0,.., Cn,..} называется с п е к т р о м сигнала. Графическое изображение {C0,.., Cn,..} в виде верти кальных отрезков, называемое спектральной диаграммой, дает наглядное представление о спектре сигна Cn ла (рис. 1.1).

Суть спектрального анализа сигнала S (t ) состоит в определении коэффициен- C C1 C2 Ci тов Cn (экспериментально или аналитиче ски) в соответствии с (1.3). 0 i 12 n На основе ряда (1.2) возможен синтез (аппроксимация) сигналов при фиксиро- Рис. 1. ванном числе N ряда N І Cn n (t ).

S (t ) = C0 0 (t ) +... + C N N (t ) = (1.2’) n = При этом обобщенный ряд Фурье обладает следующим важным свойством: при заданной системе базисных функций {n (t )} и числе слагаемых N он обеспечивает наилучший синтез (ап проксимацию), давая минимум среднеквадратической ошибки, под которой понимается величина t2 t N І = S (t ) S (t ) dt = S (t ) Cn n (t ) dt. (1.6) t1 n = t Ортогональная система называется полной, если увеличени ем N можно сделать сколь угодно малой. Ряд (1.2) называет ся в этом случае сходящимся в среднеквадратическом.

Относительная ошибка синтеза определяется по формуле = /Э, (1.7) где Э – энергия сигнала (на сопротивлении 1 Ом), численно рав ная квадрату нормы сигнала S (t ), т.е.

t = [ S (t )]2 dt.

Э= S (1.8) t Формула (1.8) с учетом ряда (1.2) может быть записана t2 [ S (t )] dt = Cn 2 Э= n, (1.9) n = t а при использовании ортонормированной системы функций {n (t )} Cn.

Э= n= Очевидно, что средняя за период T = t2 t1 мощность сигнала t Э [ S (t )] 2 Cn n Pсp = = dt =. (1.10) ТТ Т n = t Выражение вида (1.9) или (1.10) называется равенством Парсеваля.

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

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

изменяются лишь амплитуда и начальная фаза.

Во-вторых, широко используется хорошо разработанный в тео рии цепей символический метод. Представление сигналов в ба зисе гармонических функций традиционно рассматривается в радиотехнических курсах, например в [13, 14].

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

Для представления непрерывных сигналов применяются поли номы и функции Лагерра, Лежандра, Чебышева, Эрмита и др.

Для представления сигналов с точками разрыва используются кусочно-постоянные функции Уолша, Хаара, Радемахера [14].

Для дискретизации непрерывных сигналов во времени исполь зуется ортогональный ряд Котельникова [13, 14].

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

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

Вейвлет-преобразование (ВП) одномерного сигнала – это его представление в виде обобщенного ряда или интеграла Фурье по системе базисных функций 1 t b ab (t ) =, (1.11) a a сконструированных из материнского (исходного) вейвлета (t ), обладающего определенными свойствами за счет операций сдвига во времени ( b ) и изменения временного масштаба ( a ) (рис. 1.2). Множитель 1/ a обеспечивает независимость нор мы этих функций от масштабирующего числа a.

Итак, для заданных значений параметров a и b функция ab (t ) и есть вейвлет, порождаемый материнским вейвлетом (t ).

На рис. 1.2 в качестве примера приведены вейвлет «мекси канская шляпа» (а) и модуль его спектральной плотности (б).

Малые значения а соответствуют мелкому масштабу ab (t ) или высоким частотам ( ~ 1/ a ), большие параметры a – а б Рис. 1. крупному масштабу ab (t ), т.е. растяжению материнского вейвлета (t ) и сжатию его спектра.

Таким образом, в частотной области спектры вейвлетов по хожи на всплески (волночки) с пиком на частоте 0 и полосой, т.е. имеют вид полосового фильтра;

при этом 0 и уменьшаются с ростом параметра a.

Следовательно, вейвлеты локализованы как во временной, так и частотной областях.

а а а2 э э а b b Рис. 1. В соответствии с принципом неопределенности произведе ние эффективной длительности ( э ) и эффективной ширины спектра ( э ) функции ab (t ) (площадь прямоугольников на рис. 1.3) остается неизменным. Кроме того, из-за масштабирова ния и временного сдвига ( b / a = = const ) сохраняется относи тельная «плотность» расположения базисных функций по оси t.

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

1.3. Главные признаки вейвлета В качестве базисных функций, образующих ортого нальный базис, можно использовать широкий набор вейвлетов.

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

Ограниченность. Квадрат нормы функции должен быть ко нечным:

2 = (t ) dt. (1.12) Локализация. ВП в отличие от преобразования Фурье ис пользует локализованную исходную функцию и во времени, и по частоте. Для этого достаточно, чтобы выполнялись условия:

(t ) C (1 + t )1 и S () C (1 + ) 1, при 0. (1.13) Например, дельта-функция (t ) и гармоническая функция не удовлетворяют необходимому условию одновременной локали зации во временной и частотной областях.

Нулевое среднее. График исходной функции должен осцил лировать (быть знакопеременным) вокруг нуля на оси времени (см. рис. 1.2) и иметь нулевую площадь (t )dt = 0. (1.14) Из этого условия становится понятным выбор названия «вейв лет» – маленькая волна.

Равенство нулю площади функции (t ), т.е. нулевого мо мента, приводит к тому, что фурье-преобразование S () этой функции равно нулю при = 0 и имеет вид полосового фильт ра. При различных значениях a это будет набор полосовых фильтров.

Часто для приложений бывает необходимо, чтобы не только нулевой, но и все первые n моментов были равны нулю n t (t ) dt = 0. (1.15) Вейвлеты n -го порядка позволяют анализировать более тон кую (высокочастотную) структуру сигнала, подавляя медленно изменяющиеся его составляющие.

Автомодельность. Характерным признаком ВП является его самоподобие. Все вейвлеты конкретного семейства ab (t ) имеют то же число осцилляций, что и материнский вейвлет (t ), поскольку получены из него посредством масштабных преобразований ( a ) и сдвига ( b ).

1.4. Примеры материнских вейвлетов Основные вейвлетообразующие функции, или материн ские вейвлеты, приведены в табл. 1.1.

Наиболее распространенные вещественные базисы конст руируются на основе производных функции Гаусса ( g 0 (t ) = exp( t 2 / 2)). Это обусловлено тем обстоятельством, что функция Гаусса имеет наилучшие показатели локализации как во временной, так и в частотной областях.

На рис. 1.4 показаны вейвлеты первых четырех порядков и модули их спектральной плотности. При n = 1 получаем вейвлет первого порядка, называемый WAVE-вейвлетом с равным нулю нулевым моментом. При n = 2 получаем MHAT-вейвлет, назы ваемый «мексиканская шляпа» (mexican hat – похож на сомбре ро). У него нулевой и первый моменты равны нулю. Он имеет лучшее разрешение, чем WAVE-вейвлет.

Свойства гауссовых вейвлетов подробно описаны в [13]. В работе [30] показано, что совместное использование вейвлетов g1 g 4 для ВП существенно повышает точность вейвлет анализа.

Т а б л и ц а 1. Аналитическая запись Спектральная плотность Вейвлеты () (t ) Вещественные непрерывные базисы Гауссовы:

– первого порядка, или t exp(t 2 / 2) (i) 2 exp(2 / 2) WAVE-вейвлет, – второго порядка, или MHAT-вейвлет «мек (1 t 2 )exp(t 2 / 2) (i)2 2 exp( 2 / 2) сиканcкая шляпа» – mexican hat), dn exp(t 2 / 2) (1) n (i) n 2 exp(2 / 2) (1) n – n-го порядка, dt n DOG – difference of 2 2 2 (e e 2 ) e t 0,5e t / 2 / gaussians (2)1/ 2, t 2, (t ) 1 (sin 2t sin t ) LP-Littlewood & Paley 0, в противном случае Вещественные дискретные 1, 0 t 1/ 2, sin 2 / iei / 1, 1/ 2 t 1, HAAR-вейвлет / 0, t 0, t 0.

1, t 1/ 3, FHAT-вейвлет, или 4 sin 3 / «французская шляпа» 1/ 2, 1/ 3 t 1, (French hat – похож на 3 / 0, t 1.

цилиндр) Комплексные () 2 e(0 ) / ei0t e t Морле (Morlet) / Пауля (Paul) (чем in больше n, тем больше () 2 () n e Г (n + 1) (1 n) n + нулевых моментов име ет вейвлет) Наиболее простой пример дискретного вейвлета – это HAAR-вейвлет. Недостатком его являются несимметричность формы и негладкость – резкие границы в t-области, вследствие 5. 1 S () g (t) 4. 1 g S () g (t) 0 3. 2 g S () g (t) 2. 3 g S () g (t) 1. 4 g 2 0. 3 432101234 0 0.5 1 1.5 2 2.5 3 3.5 t t Рис. 1. чего возникает бесконечное чередование «лепестков» в частот ной области, хотя и убывающих как 1/.

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

Среди комплексных вейвлетов наиболее часто использует ся базис, основанный на хорошо локализованном и во вре менной и в частотной областях вейвлете Морле. Характерный параметр 0 позволяет изменять избирательность базиса. Ве щественная и мнимая части (t ) – это амплитудно-модулиро ванные колебания.

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

В настоящее время выбор вейвлетов довольно обширен.

Только в пакете Wavelet Toolbox 2.0/2.1 (MATLAB 6) пред ставлено полтора десятка материнских вейвлетов;

при этом для ряда из них дано ещё множество вариантов. Для получения справки по какому-либо типу вейвлета при работе в командном режиме MATLAB достаточно исполнить команду waveinfo (‘type’), указав тип вейвлета. Для просмотра же вейвлетов доста точно исполнить команду wavemenu и в появившемся окне со списком разделов ВП нажать кнопку Wavelet Display. Нажатие этой кнопки выводит окно просмотра вейвлетов, в котором име ется возможность просмотра: общей информации о вейвлетах, выбранного вейвлета (с именем ‘Name’) и информации о нем.

На рис. П.1 дано окно просмотра Wavelet Display с данными о вейвлете Добеши db4.

Сведения по сравнению вейвлетов различного типа приведе ны в [8, разд. 2.9].

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

1.5. Непрерывное вейвлет-преобразование Непрерывное (интегральное) вейвлет-преобразование (НВП или СWT – continuous wavelet transform).

Сконструируем базис ab (t ) с помощью непрерывных мас штабных преобразований ( a ) и переносов ( b ) материнского вейвлета (t ) с произвольными значениями базисных парамет ров a и b в формуле (1.11).

Тогда по определению прямое (анализ) и обратное (синтез) HВП (т.е. ПНВП и ОНВП) сигнала S (t ) запишутся:

t b Ws (a, b) = ( S (t ), ab (t ) ) = S (t ) dt, (1.16) a a 1 dadb Ws (a, b) ab (t ) S (t ) =, (1.17) a C где C – нормирующий коэффициент С = () d, (, ) – скалярное произведение соответствующих сомножителей, () – фурье-преобразование вейвлета (t ). Для ортонорми рованных вейвлетов C = 1.

Из (1.16) следует, что вейвлет-спектр Ws (a, b) (wavelet spec trum, или time-scale-spectrum – масштабно-временной спектр) в отличие от фурье-спектра (single spectrum) является функцией двух аргументов: первый аргумент а (временной масштаб) ана логичен периоду осцилляций, т.е. обратен частоте, а второй b – аналогичен смещению сигнала по оси времени.

Следует отметить, что Ws (b, a0 ) характеризует временную зависимость (при a = a0 ), тогда как зависимости Ws (a, b0 ) мож но поставить в соответствие частотную зависимость (при b = b0 ).

Если исследуемый сигнал S (t ) представляет собой одиноч ный импульс длительностью u, сосредоточенный в окрестно сти t = t0, то его вейвлет-спектр будет иметь наибольшее значе ние в окрестности точки с координатами a = u, b = t0.

Способы представления (визуализации) Ws (a, b) могут быть различными. Спектр Ws (a, b) является поверхностью в трех мерном пространстве (см. рис. 1.5). Однако часто вместо изо бражения поверхности представляют её проекцию на плоскость ab с изоуровнями (рис. 1.6), позволяющими проследить изме нение интенсивности амплитуд ВП на разных масштабах (а) и во времени ( b ). Кроме того, изображают картины линий ло кальных экстремумов этих поверхностей, так называемый ске летон (sceleton), который выявляет структуру анализируемого сигнала.

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

Линейность. Она следует из скалярного произведения (1.16):

W [S1 (t ) + S2 (t )] = W1 (a, b) + W2 (a, b).

Сдвиг. Смещение сигнала во временной области на b0 ведет к сдвигу вейвлет-образа также на b0 :

W [ S (t b0 )] = W [a, b b0 ].

Масштабирование. Растяжение (сжатие) сигнала приводит также к растяжению (сжатию) его в области W (a, b) :

1 a b W [ S (t a0 )] = W,.

a0 a0 a Дифференцирование:

W [dtm S ] = (1) m S (t )dtm [ ab (t )]dt, где dtm = d m [...]/ dt m, m 1. Из этого свойства следует, что про игнорировать, например, крупномасштабные составляющие и проанализировать особенности высокого порядка или мелко масштабные вариации сигнала S (t ) можно дифференцировани ем нужное число раз либо вейвлета, либо самого сигнала. Если учесть, что часто сигнал задан цифровым рядом, а анализирую щий вейвлет–формулой, то это свойство весьма полезное.

Масштабно-временная локализация. Она обусловлена тем, что элементы базиса ВП хорошо локализованы и обладают под вижным частотно-временным окном.

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

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

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

Рассмотрим эту возможность на некоторых простейших примерах.

1.7.1. Определение вейвлет-спектра на основе «мексиканской шляпы»

в системе Mathcad Программирование ВП базируется на соотношениях (1.11), (1.16), (1.17).

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

Следует отметить, что на практике невозможно проводить вычисления с непрерывными значениями a и b ;

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

Пример 1.1. Гармоническое колебание s (t ) = A sin(t ), где A = 1 В, = 2 / T = 2 / 50, = 0.

Вейвлетобразующая функция:

d exp(t 2 / 2).

N := 256, MHAT (t ) := dt t b ( a, b, t ) = Вейвлеты:.

MHAT a a Вейвлет-спектр:

N (a, b, t )s(t )dt, a := 1..30, b := 0..50, W (a, b) := N ab := W (a, b) N График двухпараметрического спектра N ab = W ( a, b) выведен на рис. 1. в виде поверхности в трехмерном пространстве, а на рис. 1.6 – в виде привычных для ВП изоуровней на плоскости (a,b). Следует отметить, что сечение W (a, b) для временного масштаба a = a0 характеризует исходное колебание s(t ) ;

при этом амплитуда его максимальна при a0 : 1/. Зависимости W (a, b0 ) можно поставить в соответствие текущий спектр колебания при b = b0.

Рис. 1. Рис. 1. Пример 1.2. Сумма двух гармонических колебаний Сигнал имеет вид: s (t ) := A1 sin(1t ) + A2 sin(2t ), где A1 = A2 = 1 В, 1 = 2 /T1, 2 = 2 /T2, T1 = 50 с, T2 = 10 с.

1. 0. s ( t) 0. 1. 2 0 25 50 75 100 125 150 175 200 225 0 t d 2 t 2 / MHAT (t ) :=, N:=256, e dt N t b (a, b, t ) f (t )dt, (a, b, t ) := MHAT W (a, b) :=, a N N ab := W (a, b).

a := 1…30, b : = 0…50, График двухпараметрического спектра W (a, b) выведен на рис. 1.7 в виде поверхности в трехмерном пространстве, а на рис. 1.8 в виде изоуровней на плоскости (a,b).

Рис. 1. Рис. 1. На рис. 1.9, а приведены «сечения» вейвлет-спектра для двух значений параметра а. При относительно небольшом параметре временного масштаба a, т.е. при a1 = 2 ( a1 : 1/ 2 ), сечение спектра несет информацию только о высокочастотной составляющей сигнала, отфильтровывая (подавляя) его низ кочастотный компонент. С ростом a происходит растяжение базисной функ ции [(t b) / a ] и, следовательно, сужение ее спектра, т.е. уменьшение поло сы частотного «окна». В результате при a2 = 15 ( a2 : 1/ 1 ) сечение спектра представляет собой лишь низкочастотный компонент сигнала. При дальней шем увеличении a полоса окна еще уменьшается и уровень этого низкочас тотного компонента убывает до постоянной составляющей (при a 25 ), рав ной нулю для анализируемого сигнала.

а б Рис. 1. На рис 1.9, б даны сечения вейвлет-спектра W (a, b), характеризующие текущий спектр сигнала при b1 = 13 и b2 = 17.

Итак, очевидно, что спектр рассматриваемого сигнала в отличие от гармо нического (см. пример 1.1) содержит еще высокочастотный компонент в об ласти малых значений временного масштаба a ( a : 1..3 ), который соответст вует второй составляющей сигнала, т.е. A2 sin(2t ).

Пример 1.3. Прямоугольный импульс U := 5 t0 := 20 := s ( t) := U if t0 t t0 + 0 otherwise t d2 exp MHAT ( t) := s ( t) dt N := 0 20 40 60 80 100 t N t b (a, b, t ) f (t )dt, (a, b, t ) := MHAT W (a, b) :=, a N Nba := W (a, b).

a := 1..50, b := 0..100, Рис. 1. Вейвлет-спектры приведены на рис. 1.10.

Как видно из рис. 1.10, вейвлет-спектр хорошо передает тонкие особенно сти сигнала – его скачки на отсчетах b = 20 и b = 80 (при a : 1..5 ).

Следует отметить, что выполнение непрерывного ВП в Mathcad требует больших вычислительных затрат. Это обуслов лено довольно медленным методом интегрирования. На ПК с процессором Intel Celeron (667 МГц и оперативной памятью в 128 Мбайт) время вычислений составляет до 5 минут. Гораздо более эффективен ВА в системе MATLAB.

1.7.2. Вейвлет-анализ в системе MATLAB Пакет Wavelet Toolbox системы MATLAB располагает средствами для построения вейвлет-спектров сигналов с улуч шенной визуализацией. Спектрограммы представляют значения вейвлет-коэффициентов в плоскости масштаб–время;

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

В разделах «Continuous Wavelet 1-D» и «Complex Continuous Wavelet 1-D» пакета Wavelet Toolbox даны примеры на выявле ние и анализ тонкой структуры сигналов.

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

В приводимых ниже примерах прямое ВП реализуется функцией cwt, которая может быть представлена в нескольких формах (см.П.2.2), например:

COEF = cwt(S, SCALES, ‘wname’ PLOTMODE, XLIM), где строка S – сигнал, строка SKALES – задание диапазона и ша га изменения параметра a, строка ‘wname’ – имя (тип) вейвлета, Рис. 1. строка PLOTMODE – настройка цвета: ‘lvl’ – окраска шаг за ша гом, ‘glb’ – окраска с учетом всех коэффициентов, ‘absvil’ или ‘lvlabs’ – окраска шаг за шагом с использованием абсолютных значений коэффициентов, строка XLIM – строка переменных настройки.

Пример 1.4. Гармоническое колебание function garm t = 0:0.00001:0.0004;

A1 = 1;

F1 = 10000;

a1 = 45;

s = A1*soc(2*pi*F1*t+a1);

figure (1);

plot(t,s);

axis([0 0.0004 -3 3]);

grid on;

subplot(211), plot(t,s);

title('Сигнал S(t)') subplot(212), c = cwt(s, 1:2:32, 'mexh', 'abs1vl', [0 10]);

title('Вейвлет-спектр ');

xlabel('Временной сдвиг, b');

ylabel('Временной масштаб, a');

end И хотя спектрограмма гармонического колебания, представленная на рис. 1.12, особой выразительностью не отличается, на ней отчетливо видны переходы сигнала через нуль (темный тон) и экстремальные точки (светлый тон).

Рис. 1. Пример 1.5. Сумма двух гармонических колебаний Сигнал S (t ) представляет собой сумму двух гармонических колебаний с кратными частотами:

function binar t = 0:0.000001:0.0004;

A1 = 1;

A2 = 1;

F1 = 10000;

F2 = 2*F1;

a = 90;

b = 90;

a1 = a*0.0174533;

a2 = b*0.0174533;

s = A1*sin(2*pi*F1*t-a1) + A2*sin(2*pi*F2*t-a2);

figure (1);

plot(t,s);

axis([0 0.0004 -3 3]);

grid on;

subplot(211), plot(t,s);

title('Сигнал S(t)') Рис. 1. subplot(212), c = cwt(s,1:2:50,'mexh','abs1v',[0 1]);

title('Вейвлет-спектр сигнала S(t)');

xlabel('Временной сдвиг, b');

ylabel('Временной масштаб, a');

end Временная диаграмма S (t ) и спектральная диаграмма сигнала W (a, b) приведены на рис. 1.13. В нижней части спектрограммы хорошо просматрива ется структура второй гармоники, а с ростом a – первой;

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

Пример 1.6. Бигармонический сигнал с шумом Сигнал x(t ) представляет собой сумму бигармонического сигнала S (t ) и белого гауссова шума n(t ) с математическим ожиданием m = 0 и среднеквад ратическим отклонением g.

function bigarm_rauch t = 0:0.000001:0.001;

A1 = 1;

A2 = 1;

F1 = 10000;

F2 = 2*F1;

a = 90;

b = 90;

a1 = a*0.0174533;

a2 = b*0.0174533;

s1(1:200) = 0;

t2 = 0.0002:0.000001:0.0007;

s2 = A1*sin(2*pi*F1*t2-a1) + A2*sin(2*pi*F2*t2-a2);

s3(1:300) = 0;

s = [s1 s2 s3];

randn('state',0);

g = 0.5;

n =g *randn(size(t));

x = s+n;

figure (1);

subplot(211), plot(t,x,'k');

t itle('Сигнал x(t)');

grid on;

gtext('F=10кГц, А1=А2=1В, g=0.5 B');

Рис. 1. subplot(212), c = cwt(x,1:124,'mexh','absglb',[0 50]);

title('Вейвлет-спектр W(a,b)');

xlabel('Временной сдвиг, b');

ylabel('Временной масштаб, a');

end На рис. 1.14 приведены диаграмма сигнала и его спектрограмма. В ниж ней части спектрограммы хорошо видна сложная структура вейвлет-спектра шума ( a : 1 5 ), с ростом параметра a просматривается структура второй гармоники, а затем – первой;

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

Пример 1.7. Прямоугольный импульс с шумом function pr_rauch_wav t = 0:0.000001:0.000300;

A1 = 2;

F1 = 0;

s1(1:75) = 0;

t2 = 0.000075:0.000001:0.000175;

s2 = A1*cos(2*pi*F1*t2);

s3(1:125) = 0;

s = [s1 s2 s3];

randn('state',0);

g = 0.5;

n = g*randn(size(t));

x = s + n;

figure (1);

subplot(211), plot(t,x,'k');

title('Сигнал x(t)');

grid on;

subplot(212), c = cwt(x,1:27,'mexh','absglb',[0 10]);

title('Вейвлет-спектр');

xlabel('Временной сдвиг, b');

ylabel('Временной масштаб, a');

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

Рис. 1. Пример 1.8. Звуковой сигнал Загрузим звуковой сигнал из файла mtlb с выборкой в 200 отсчетов:

function ss load mtlb;

v=mtlb(1:200);

lv = length(v);

subplot(211), plot(v);

title('Звуковой сигнал');

set(gca, 'Xlim',[0 200]);

[c,l] = wavedec(v,5,'sym2');

cfd = zeros(5,lv);

subplot(212) ccfs = cwt(v,1:128,'sym4','plot');

title('Вейвлет-спектр') colormap(pink(32));

xlabel('Временной сдвиг, b');

ylabel('Временной масштаб,a');

end Рис. 1. Вейвлет-спектрограмма (рис. 1.16) демонстрирует мельчайшие детали частотного образа сигнала: в нижней части отчетливо видны высокочастотные компоненты, а в верхней – низкочастотные (где изменения яркости менее час тые, чем в нижней части).

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

1.8. Сопоставление с преобразованием Фурье Классическое преобразование Фурье (ПФ) является тра диционным математическим аппаратом для анализа стационар ных процессов. При этом сигналы разлагаются в базисе косину сов и синусов или комплексных экспонент. Эти базисные функции простираются вдоль всей оси времени.

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

Часть указанных трудностей преодолевается при использо вании оконного ПФ:

S (t ) w(t b)e jt dt, S (, b) = в котором применяется предварительная операция умножения сигнала S (t ) на «окно» w(t b) ;

при этом окном является ло кальная во времени функция (например, прямоугольная, т.е.

w(t ) = 1 при 0 t и w(t ) = 0 при t 0 и t ), перемещае мая вдоль оси времени t (рис. 1.17) для вычисления ПФ в раз ных позициях b. В результате получается текущий спектр, т.е.

частотно-временное описание сигнала.

s+t W(t–b) окно 0 t b+ b Рис. 1. Если окно, показанное на рис. 1.17, перемещать скачками (через ) вдоль всего времени существования сигнала S (t ), то за некоторое число таких перемещений возможен «просмотр»

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

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

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

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

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

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

ВП просто иное и позволяет посмотреть на исследуемый про цесс с иной точки зрения.

Глава Разложение сигналов в ряд по вейвлетам 2.1. Диадное вейвлет-преобразование При непрерывном изменении параметров a и b для расче та вейвлет-спектра необходимы большие вычислительные за траты. Множество функций ab (t ) избыточно. Необходима дискретизация этих параметров при сохранении возможности восстановления сигнала из его преобразования. Дискретизация, как правило, осуществляется через степени двойки [1–3]:

1 t b 1 m a = 2m, b = k 2m, mk (t ) = = m (2 t k ), (2.1) a a где m и k – целые числа. В этом случае плоскость a, b превра щается в соответствующую сетку m, k. Параметр m называется параметром масштаба.

Рассмотренная дискретизация наиболее распространена.

Сетка дискретизации называется диадной и соответственное преобразование – диадным (dyadic) ВП.

Рис. 2.1 на примере вейвлета Хаара иллюстрирует дискрети зацию a, b. При фиксированном параметре m вейвлеты имеют одинаковые масштабы и лишь смещаются во времени. При уве личении параметра m на 1 масштаб увеличивается вдвое и вейвлеты вдвое растягиваются. Для различных значений m ши рина mk (t ) различна и выбор b = k 2m гарантирует, что растя нутые вейвлеты на уровне m «покрывают» ось времени так же, как это делают исходные вейвлеты на уровне m = 0.

30(t) 0 2 4 6 t 21(t) 20(t) 0 2 4 6 t 10(t) 11(t) 13(t) 0 2 4 6 8 t 07(t) (t) 01(t) 1.

0 2 4 6 8 t 1.

Рис. 2. Прямое и обратное диадное ВП непрерывных сигналов за пишутся в виде:

cmk = ( S (t ), mk (t ) ) = S (t ) mk (t )dt, (2.2) S (t ) = cmk mk (t ). (2.3) m,k Проводя аналогию с преобразованием Фурье, коэффициенты cmk разложения (2.3) можно определить через непрерывное ВП Ws (a, b) ( ) cmk = W 2m, k 2m. (2.4) Обращаясь к (2.2) и (2.4), видим, что вейвлет-спектр cmk можно представить как «лес» из вертикальных отрезков, разме щенных над m, k – плоскостью (сеткой);

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

Из (2.3) следует, что сигнал S (t ) может быть представлен суммой «вейвлетных волн» с коэффициентами cmk. Формаль но обобщенный ряд Фурье (2.3) отличается от традиционного тем, что суммирование проводится не по одному, а по двум индексам. Однако это несущественно, так как обе системы ин дексации принадлежат одному классу бесконечных счетных множеств.

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

a, b и базисные Примечаниe. В некоторых публикациях параметры функции задаются в виде a = 1/ 2 j, b = k / 2 j, jk (t ) = 2 j (2 j t k ), (2.1’) т.е. с ростом j параметр a уменьшается, что соответствует сжатию функции jk (t ). Согласно формулам (2.1) с ростом m увеличивается и коэффициент a, т.е. функция mk (t ) растягивается.

Фреймы. Это особый вид вейвлетов, занимающих промежуточное поло жение между непрерывным и диадным ВП. Вейвлет-фреймы используют кратное двум масштабирование ( a = 2m ), но непрерывный сдвиг. Следова тельно, они сохраняют избыточность, которая присуща непрерывному ВП, но в гораздо меньшей мере по сравнению с ним. Они не входят в пакеты расши рения систем компьютерной математики (СКМ). Но если необходимо, то соот ветствующие им инструментальные средства легко получить незначительной модификацией средств непрерывного ВП.

2.2. Дискретное преобразование Статьи, касающиеся практического использования ВП, содержат в основной своей массе результаты компьютерных расчетов, в которых использовано дискретное вейвлет-преобра зование (ДВП или DWT). При этом не только параметры a и b, но и сигналы также дискретизируются во времени.

На основании теоремы Котельникова (теоремы отсчетов) непрерывный сигнал S (t ), спектр которого не содержит частот выше f m, полностью определяется дискретной последователь ностью своих мгновенных значений {Si }, i = 0,1,..., N 1, от считываемых через интервалы времени t :

t = 1/ 2 f m, f д = 1/ t = 2 f m, (2.5) где t и f д – интервал (шаг) и частота дискретизации.

Таким образом, дискретизированный с шагом t сигнал можно определить выражением:

N S (it )(t it ), Sд (t ) = {Si } = (2.6) i = где (t ) – дельта-функция.

Если число отсчетов составляет N = 2n0, то максимальное значение m в формулах (2.1) будет равно n0 1. Наибольшее значение k для текущего m определяется: k = 2n0 m 1. В ча стности, для m = 0 (т.е. a = 1) число сдвигов k базисного вейв лета составит 2n0 1 = N 1 ;

с каждым последующим значени ем m (1, 2, …) вейвлет mk (t ) расширяется в два раза, а число сдвигов k уменьшается в два раза. Для максимального значе ния m = mmax, равного n0 1, k = 0, т.е. один вейвлет mmax 0 (t ) «накрывает» весь интервал сигнала (рис. 2.1;

N = 8 ).

Вейвлет-коэффициенты cmk (или c jk ) можно вычислить с помощью итерационной процедуры, известной под названием быстрого вейвлет-преобразования БВП [1–3, 19, 25, 29]. Алго ритм БВП приведен в п. 2.4. При этом, если необходимо, можно сжать полученные данные, отбросив некоторую несуществен ную часть закодированной таким образом информации. Осуще ствляется это квантованием, в процессе которого приписывают ся разные весовые множители различным вейвлет-коэффициен там. Аккуратно проведенная процедура позволяет не только удалить некоторые статистические флюктуации и повысить роль динамических характеристик сигнала, но и существенно сократить компьютерную память и требования к передаче ин формации и, следовательно, снизить расходы.

2.3. Примеры дискретного вейвлет-преобразования (ДВП) 2.3.1. ДВП в Mathcad Системы компьютерной математики Mathcad первыми использовали прямое и обратное дискретное ВП. В ядро систем (начиная с версии Mathcad 8) встроен единственный вейвлет – Добеши db4 (или DB4). При этом реализация ВП происходит с большой скоростью (т.е. эффективностью) и можно осуществ лять практическое исследование различных сигналов и времен ных рядов на выявление как их свойств, так и свойств ВП.

Ядро систем Mathcad содержит две следующие функции ВП:

wave(x) – вектор прямого ВП;

iwave(w) – вектор обратного ВП.

Вектор данных x и вектор вейвлет-спектра w должны иметь ровно N = 2n0 элементов ( n0 – целое число). Результатом функ ции wave(x) является вектор, скомпонованный из коэффициен тов двухпараметрического вейвлет-спектра cmk.

Пример 2.1. Прямоугольный импульс с шумом Исследуемый сигнал x(t ) представляет собой аддитивную смесь x(t ) = S (t ) + n(t ) прямоугольного видеоимпульса S (t ) и белого нормального шума n(t ) (рис. 2.2):

U = 5 (В), t0 = 40 (мкс), = 60 (мкс) s(t):= U ift 0 t t 0 +.

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

xi 0 20 40 60 80 100 120 140 160 180 200 220 240 0 i Рис. 2. Представление сигнала и шума в дискретном виде:

n0 = 8, N = 2n0, N = 256, i := 0..N 1, si := s (i ) ni := 2ln(rnd (1)) sin(2 rnd (1)), xi := si + ni := 0, Вейвлет-анализ, т.е. прямое ДВП:

i := 0..N 1 y := x w := wave( y ) z := n0 1 z := 7 m := 1, 2..z level + level coeffs (level ) := submatrix( w, 2 1,0,0) ci, z m := coeffs (m),2 i flor N m Семейства коэффициентов вычисленного вейвлет-спектра показаны на рис. 2.3, а весь спектр – на рис. 2.4.

m Примечание. У коэффициентов (c )i нижний индекс i означает номер текущего отсчета времени и принимает N значений от 0 до N–1, а верхний m имеет тот же смысл, что и у вейвлет-коэффициентов cmk, определяемых по формуле (2.2). Напомним, что параметры m и k (которым соответствуют индексы вейвлет-коэффициентов) характеризуют дискретные изменения вре менного масштаба ( a = 2m ) вейвлета и его сдвига (b = k 2m ) во времени. Для текущего масштаба m параметр k имеет 2n0 m значений от 0 до 2n0 m 1. В частности, для m = 0 ( a = 1 ) вейвлет 0 k ( x) смещается N раз (включая нуль), т.е. индекс k в cmk и индекс i в (c )i совпадают. При m = 1 вейвлет 1k ( x) расширяется по сравнению с вейвлетом 0 k ( x) в два раза и общее число сдвигов будет в два раза меньше;

при этом значение k будет изменять ся через два отсчета i. Для наибольшего временного масштаба, когда m = n0 1 (в данном случае 7), k = 0 и один вейвлет 7,0 ( x) «накроет» весь временной интервал;

при этом значение (c )i будет постоянным и равным c7,0 при всех значениях i от 0 до N – 1.

2. ( c 0 ) i ( c 1 ) i 1.

0 50 100 150 200 0 i (c )i ( c 3 ) i ( c 4 ) i 6.

0 50 100 150 200 0 i (c )i ( c 6 ) i ( c 7 ) i 0 50 100 150 200 0 i Рис. 2. c Рис. 2. m := 7 xi x1i 2 0 20 40 60 80 100 120 140 160 180 200 220 240 0 i m := 7 x1i 2 0 20 40 60 80 100 120 140 160 180 200 220 240 0 i m := 7 x1i 2 0 20 40 60 80 100 120 140 160 180 200 220 240 0 i m := 7 x1i 2 0 20 40 60 80 100 120 140 160 180 200 220 240 0 i Рис. 2. Вейвлет-синтез, т.е. обратное ДВП. Синтезируемый сигнал:

x1i := iwave( w).

Осуществим синтезирование сигнала с подавлением коэффициентов cmk при быстрых (высокочастотных) слагаемых обобщенного ряда (2.3):

j := 2 z m... N 1 w j := 0.

Результаты представлены на рис. 2.5. Очевидно, что при m = 0 синтез проис ходит без подавления составляющих и исследуемый xi и синтезируемый x1i сигналы полностью совпадают.

С увеличением параметра m расширяется полоса подавления составляю щих в вейвлет-спектре, что эквивалентно пропусканию сигнала через фильтр низких частот с уменьшающейся полосой пропускания фильтра и, следова тельно, росту подавления шума и относительно высокочастотных компонентов сигнала;

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

Пример 2.2. Вейвлет-фильтрация бигармонического сигнала с шумом На вход низкочастотного приемного тракта фазового параметрического гидролокатора [38] поступает сигнал x(t ) = n(t ) + sэ (t ), где n(t ) – помеха в виде белого гауссова шума с математическим ожиданием mn = 0 и среднеквадратическим отклонением ;

sэ (t ) – бигармонический (двухкомпонентный) сигнал sэ (t ) = A1 sin[2F1t 1(t )] + A2 sin[2F2t 2 (t )], t з t t з + и, где A1, A2 – амплитуды компонентов эхосигнала на частотах F1 = F и F2 = 2 F, 1 (t ), 2 (t ) – фазовые сдвиги, зависящие от акустической жестко сти и структуры подводных объектов, tз и и – время задержки и длитель ность эхоимпульса. Измеритель фазового сдвига приемника дает на своем выходе: во время действия сигнала sэ (t ) напряжение, пропорциональное фа зовому сдвигу (t ) = 21 (t ) 2 (t ), и равномерный шум вне интервала tз t tз + и. Присутствие шума n(t ) привносит случайную ошибку в изме рения (t ).

Задача ВП – это осуществление фильтрации sэ (t ) из шума n(t ). Отфильт рованный сигнал подается на измеритель фазового сдвига (ИФС на рис. 2.6).

(t ) y (t ) x(t ) ИФС ВП Рис. 2. Моделирование сигнала x(t ), его ВП, алгоритма измерения фазового сдвига (t ) и оценки параметров входного и выходного сигналов осуществле но в пакете Mathcad (2001).

Дискретное ВП осуществлено на основе встроенной в пакет Mathcad ба зисной функции Добеши db 4 : wave(x) – вектор прямого ВП, iwave(w) – вектор обратного ВП. Сигнал x(t ) был подвергнут (как и в примере 2.1) прямому ВП, и по найденному вектору коэффициентов cmk осуществлено обратное ВП с подавлением коэффициентов при быстрых (высокочастотных) слагаемых ряда (2.3): yi := iwave( w), j := 2 z m..N 1, w j := 0.

На рис. 2.7 представлены входной сигнал x(t ) и результат его вейвлет фильтрации y (t ). Здесь A1 = A2 = 1 В, F = 10 кГц, 1 (t ) = 2 (t ) = = 90°, = 0,5 В, t з = 0, 2 мс, = 0,6 мс, n0 = 10 (z = 9), N = 1024.

В ходе исследований изменялся параметр m и измерялись: среднеквадра % тическое значение шума y после ВП, фазовый сдвиг, отклонение % % измеренного фазового сдвига от истинного = 2 = и = 1 среднеквадратическое отклонение фазового сдвига. Результаты сведены в табл. 2.1.

При m = 0 восстановление (синтез) сигнала происходит без подавления шума и синтезируемый y1 и исследуемый xi полностью совпадают. С увели чением m (на рис. 2.7 m = 3 ) расширяется полоса подавления спектра, сужа ется полоса пропускания низкочастотной части спектра, что уменьшает уро вень шума ( y ).

Рис. 2. Т а б л и ц а 2. m 0 1 2 3 4 5 y, В 0.493 0.327 0.198 0.125 0.105 0.078 0., град –5.73 –5.35 –5.01 –4.24 –4.84 –6.91 120.

, град 14.87 14.81 18.47 14.08 14.54 18.8 Установлено, что существует оптимальное mopt, при котором отклонение, или среднеквадратическая ошибка, в измерении фазового сдвига эхо сигнала минимально. При m mopt наряду с дальнейшим подавлением шума происходит искажение формы бигармонического импульсного сигнала, что приводит к росту и. Отметим, что при отсутствии шума ( = 0 ) из мерение фазового сдвига происходит без ошибок (т.е. = 0 и = 0 ).

2.3.2. ДВП в Matlab Изучение ДВП, как и непрерывного ВП, лучше осущест влять с помощью графического интерфейса GUI (прил.1). Для вызова меню следует исполнить команду wavemenu. В появив шемся окне со списком разделов ВП активизировать позицию Wavelet 1-D, а далее установить File Demo Analys и выбрать один из 32 примеров применения вейвлет-технологии.

Пакет расширения систем MATLAB 6.0/6.1 Wavelet Toolbox 2/2.1 содержит несколько функций нахождения вейвлет коэффициентов (прил. 2, П.2.2), например, coef = detcoef(C,L,M).

Эта функция возвращает коэффициенты на уровне М из струк туры wavelet разложения [C, L];

при этом уровень М должен быть целым числом, таким, что 1 M MMAX, где MMAX = length( L) 2.Функция [С, L] = wavedec(S, M, 'wname') возвращает векторы wavelet разложения сигнала Х на уровне М, используя выбранный вейвлет (с именем ‘wname’).

Пример 2.3. Бигармонический сигнал с шумом Модель такого сигнала приведена в примере 2.2. Найдем его непрерывный и дискретный вейвлет-спектры:

function binar_rauch_wav_ t = 0:0.000001:0.001;

A1 = 1;

A2 = 1;

F1 = 10000;

F2 = 2*F1;

a1 = 0;

a2 = 0;

s1(1:200) = 0;

t2 = 0.0002:0.000001:0.0008;

s2 = A1*sin(2*pi*F1*t2 + a1) + A2*sin(2*pi*F2*t2+a2);

s3(1:200) = 0;

s = [s1 s2 s3];

randn('state',0);

g = 0.5;

n = g*randn(size(t));

x = s + n;

figure (1);

subplot(311), plot(t,x,'k');

title('Сигнал x(t)');

grid on;

gtext('F=10кГц, А1=А2=1В, g=0.5 B');

subplot(312), c = cwt(x,1:64,'mexh','absglb',[0 400]);

title('Вейвлет-спектр');

xlabel('Временной сдвиг, b');

ylabel('Временной масштаб,a');

set(gca,'Xlim',[0 1000]);

[c,l] = wavedec(s,6,'db4');

for m = 1: d = detcoef(c,l,m);

d = d(ones(1,2^m),:);

cfd(m,:) = wkeep(d(:)',1000);

end cfd = cfd(:);

I = find(abs(cfd)sqrt(eps));

cfd(I) = zeros(size(I));

cfd = reshape(cfd,6,1000);

subplot(313), colormap(pink(16));

img = image(flipud(wcodemat(cfd,64,'row')));

set(get(img,'parent'), 'YtickLabel',[]);

title('Дискретное преобразование');

ylabel('Уровень, m');

xlabel('Временной сдвиг, b');

end Рис. 2. На рис. 2.8 приведены диаграмма сигнала и его спектрограммы. Особен ности спектрограммы непрерывного ВП этого сигнала обсуждены в примере 1.6. Очевидно, что детали сигнала просматриваются и на спектрограмме дис кретного ВП, но с худшим разрешением.

Пример 2.4. Звуковой сигнал Загрузим звуковой сигнал из файла mtlb с выборкой в 200 отсчетов (см.

пример 1.7) и построим его график и две спектрограммы – непрерывного и дискретного ВП:

function ss_cd load mtlb;

v = mtlb(1:200)', lv = length(v);

subplot(311), plot(v);

title('Звуковой сигнал');

set(gca, 'Xlim',[0 200]);

[c,l] = wavedec(v,6,'sym2');

cfd = zeros(6,lv);

subplot(312), ccfs = cwt(v,4:127,'sym2','plot');

title('Непрерывное преобразование'), colormap(pink(32));

ylabel('Временной масштаб,a');

xlabel('Временной сдвиг, b');

for m = 1: d = detcoef(c,l,m);

d = d(ones(1,2^m),:);

cfd(m,:) = wkeep(d(:)',lv);

end cfd=cfd(:);

I = find(abs(cfd)sqrt(eps));

cfd(I) = zeros(size(I));

cfd = reshape(cfd,6,lv);

subplot(313), colormap(pink(32));

img = image(flipud(wcodemat(cfd,64,'row')));

set(get(img,'parent'), 'YtickLabel',[]);

title('Дискретное преобразование');

ylabel('Уровень');

xlabel('Временной сдвиг, b');

end Полученные в результате выполнения программы графики приведены на рис. 2.9.

Рис. 2. Очевидно, что все мельчайшие детали сложной временной зависимости S (t ) отчетливо просматриваются на спектрограмме как непрерывного, так и дискретного ВП. Однако последнее выполняется значительно быстрее, хотя по детальности представления уступает непрерывному ВП.

Примечание. Обычное дискретное ВП (DWT) осуществляется исходя из предположения нестационарности сигнала. Если сигнал стационарный, то в этом частном случае можно использовать стационарное ВП [8]. Это преобра зование, применяемое для очистки сигнала от шума, имеет ряд форм записи (см. прил. 2.6).

2.4. Быстрое вейвлет-преобразование При исследовании сигналов полезно их представление в виде совокупности последовательных приближений грубой (ап проксимирующей) Am (t ) и уточненной (детализирующей) Dm (t ) составляющих m S (t ) = Am ( t ) + D j ( t ), (2.7) j = с последующим их уточнением итерационным методом. Каж дый шаг уточнения соответствует определенному масштабу a m (т.е. уровню m ) анализа (декомпозиции) и синтеза (реконст рукции) сигнала. Такое представление каждой составляющей сигнала вейвлетами можно рассматривать как во временной, так и в частотной областях. В этом суть кратномасштабного анализа (КМА). В прил. 3 описан КМА для непрерывных сигналов.

Как уже отмечалось в разд. 2.3, в практике ВП в большинст ве случаев мы имеем дело с дискретными сигналами. Однако формулы для ВП дискретных сигналов не могут быть получены простой дискретизацией формул диадного ВП для непрерывно го сигнала. Найдем их из предпосылок КМА.

Пусть имеется непрерывный сигнал S (t ) V0. Дискретный сигнал Sд интерпретируем как последовательность коэффици ентов ak, полученную в ходе КМА сигнала S (t ) при масштаби рующих функциях 0 k (t ) :

S (t ) = A0 (t ) = a0 k 0 k (t ), (2.8) k где a0 k = ak = ( S (t ), 0 k (t )) – коэффициенты аппроксимации на уровне m = 0.

По концепции КМА сигнал S (t ) декомпозируется на две со ставляющие (принадлежащие подпространствам V1 и W1 ):

S (t ) = A1 (t ) + D1 (t ) = a1k 1k (t ) + d1k 1k (t ). (2.9) k k Следовательно, получены две новые последовательности a1k и d1k. Отметим, что последовательности a1k и d1k имеют поло винную длину по сравнению с a0k. Далее процесс декомпози ции может быть продолжен по A1 (t ) (подпространства V2 и W2 ). Сигнал S (t ) на уровне декомпозиции m будет представ лен совокупностью коэффициентов amk и d mk.

Однако вычисления amk и d mk по-прежнему зависят от не прерывных базисных функций (t ) и (t ). Как показано в прил. 3, эти функции однозначно определяются коэффициента ми hl :

(t ) = 2 hl (2t l ), (2.10) l (t ) = 2 (1)l h1l (2t l ) = 2 gl (2t l ), (2.11) l l hl = ((t ), (2t l ), (2.12) gl = (1)l h2 n 1l, (2.13) где l = 0,1,..., lo = 2n 1, n – порядок вейвлета. Вейвлеты n-го порядка существуют только на интервале длиной 2n 1 и име ют 2n отличающихся от нуля коэффициентов hl.

Из (2.10) и (2.11) можно получить следующие соотношения:

amk = ( S (t ), mk (t )) = hl 2 k ((t ), m 1,l (t )) = hl 2 k al,m1, l l (2.14) d mk = ( S (t ), mk (t )) = gl 2k ((t ), m 1,l (t )) = gl 2k al,m1.

l l (2.15) Итерационная процедура быстрого вейвлет-анализа получи ла название анализа от «тонкого» к «грубому» масштабу.

На практике наименьший возможный масштаб (наибольший уровень разрешения n0 ) определяется числом N дискретных значений сигнала ( N = 2n0 ). На самом «тонком» значении мас штаба ( m = 0, a = 2m = 1 ) за аппроксимирующие коэффициенты a0k принимаются сами отсчеты Si сигнала S (t ), т.е. a0k = Si, k = i, i = 0, 1,..., N 1. При переходе от текущего масштаба m к следующему m + 1 число вейвлет-коэффициентов уменьшается в два раза и они определяются по рекуррентным соотношениям:

am +1,k = hl 2k aml, d m +1,k = gl 2 k aml. (2.16) l l Процесс останавливается после конечного числа уровней m = MMAX, которое зависит от протяженности сигнала ( N ) и порядка ( l ) фильтра hl.

При восстановлении (реконструкции) сигнала по его вейв лет-коэффициентам процесс идет от крупных масштабов к мел ким и на каждом шаге описывается выражением am 1,k = (hk 2l aml + g k 2l d ml ), (2.17) l которое получается из соотношений (2.10) и (2.11).

Число операций умножения при прямом быстром ВП (БВП) будет 2LN, где L = 2n [29]. Столько же операций необходимо и для реконструкции сигнала. Таким образом, для анализа синтеза сигнала в базисе вейвлетов необходимо выполнить 4LN операций, что не превышает ( и даже меньше) числа опе раций для быстрого преобразования Фурье ( N log 2 N ).

Таким образом, в практических приложениях с применением БВП используются только коэффициенты hl, сами же вейвлеты не вычисляются и в расчетах не используются.

Пакет Wavelet Toolbox позволяет осуществлять КМА с ис пользованием БВП. При этом порядок следования коэффициен тов – «дерево» коэффициентов приведено на рис. 2.10: декомпо зиция сигнала – сверху-вниз и реконструкция – снизу-вверх.

S m=0 S = {Si} m=1 S = A1 + D A1 D D A2 m=2 S = A2 + D2+D m=M S = AM + DM+DM–1,…,+D Рис. 2. На рис. 2.11 приведен пример КМА, взятый из раздела де монстрационных примеров «wavedemo». Слева под сигналом представлены аппроксимирующие коэффициенты am, а справа – детализирующие d m ( m от 1 до 5). Очевидно, что коэффици енты аппроксимации являются грубыми копиями сигнала, а де тализирующие коэффициенты выделяют локальные особенно сти и свойства сигнала. Справа сверху приведен также вейвлет спектр сигнала – cfs.

Функции для нахождения этих коэффициентов имеют ряд форм и, в частности:

A = appcoef (C, L,' wname ', M ), (2.18) D = det coef (C, L, M ) – возвращают аппроксимирующие и детализирующие коэффи циенты из структуры вейвлет-разложения [C, L], ' wname ' – строка, содержащая имя вейвлета, уровень M должен быть 1 M MMAX, целым числом, таким, что где MMAX = length( L) 2.

Рис. 2. Структура вейвлет-разложения [C, L] = wavedec( X, N, ' wname ') (2.19) возвращает векторы вейвлет-разложения сигнала X на уровне N ;

выходная структура разложения содержит векторы C и L.

Реконструкцию (восстановление) сигнала S с многоуровне вой структурой разложения [C, L] осуществляет функция wa verec:

S = waverec(C, L,' wname '). (2.20) 2.5. О вейвлетах для БВП Большинство используемых вейвлетов не имеют, к сожа лению, аналитического выражения. Однако из предыдущего рассмотрения следует, что для практических расчетов исполь зуются не сами вейвлеты, а их коэффициенты hl.

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

Следует отметить, что процесс определения коэффициентов hl, т.е. конструирования вейвлетов, достаточно сложен для пользователя. Да в этом и нет особой необходимости, так как уже создано большое число вейвлетов, в том числе входящих в пакет расширения Wavelet Toolbox, например, вейвлеты Добе ши (dbN), Симплета (sumN), Койфлета (coifN), Хаара (haar) и др.;

их подробное описание приведено в [7, 8].

Особо следует отметить вейвлеты Добеши. Это один из са мых известных и используемых во многих практических прило жениях типов. Вейвлеты порядка N (dbN) отличны от нуля лишь на интервале длиной 2N–1 и имеют 2N отличных от нуля коэф фициентов фильтров hl и gl. Исключая случай N = 1 (а это есть базис Хаара, который не регулярен), функции (t ) и (t ) не прерывны и дифференцируемы.

Решение уравнения (2.11) для этих вейвлетов дает [5, 8]:

для n = 2 (четырехточечный фильтр Добеши):

h0 = (1 + 3) /(4 2) = 0.482963, h1 = (3 + 3) /(4 2) = 0.836516, h2 = (3 3) /(4 2) = 0.224144, h3 = (1 3) /(4 2) = – 0.129409, g0 = h3, g1 = h2, g 2 = h1, g3 = h0 ;

для n = 3 (шеститочечный фильтр):

h0 = 0.332670, h1 = 0.806891, h2 = 0.459877, h3 = – 0.135011, h4 = – 0.085441, h5 = 0.035227.

для n = 4 (восьмиточечный фильтр):

h0 = 0.230377, h1 = 0.714847, h2 = 0.630881, h3 = – 0.027984, h4 = – 0.187035, h5 = 0.030841, h6 = 0.032883, h7 = – 0.010597.

На рис. 2.12 приведены отцовский (t ) (сплошной линией) и материнский (t ) вейвлеты второго, третьего и четвертого по рядков, которые задаются приведенными выше коэффициента ми фильтров.

2.0 0.2 0. n= n=2 n= 1. 1.0 0.1 0. 0. 0 -0. -1.0 -0. -0. -1. 0 1 2 3 4 6 2 0 Рис. 2. Очевидно, что вейвлеты высокого порядка ( n = 3 и n = 4) более гладкие по сравнению с db 2 ;

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

В вейвлет-преобразованиях, осуществляемых системой Mathcad, используется вейвлет db 4.

2.6. Частотный подход к ВП До сих пор рассмотрение ВП базировалось на временном подходе. Однако также плодотворна трактовка ВП в частотной области на базе частотной фильтрации. В этом случае КМА сиг нала рассматривается как поэтапная процедура фильтрации.

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

Обратимся к схеме рис. 2.13. Сигнал S подается на низко частотный (нижняя часть схемы) и высокочастотный фильтры декомпозиции Lo _ D и Hi _ D. В них вычисляется свертка (цифровая фильтрация) по формуле:

2 n S (l ) q ( k l ), y (k ) = (2.21) l = где 2n – число отсчетов импульсной характеристики q () фильтра.

В соответствии с (2.21) и (2.7) на выходе фильтров будут ВЧ и НЧ компоненты сигнала:

D1 = y H (k ), A1 = yL ( k ).

Секция анализа Секция синтеза (декомпозиции) (реконструкции) % % cD D y H (k ) S (n) S (k ) Hi _ R Hi _ D SH ( z) S ( z) % A c A yL ( k ) 2 Lo _ R Lo _ D SL ( z) Рис. 2. Из сопоставления (2.21) и (2.16) следует, что для вычисления коэффициентов amk и d mk (на первом этапе m = 1 ) аргументы весовых коэффициентов фильтров h(l ) = hl и g (l ) = gl должны быть взяты с обратным знаком (порядком следования), т.е. hl и g l. Такие фильтры называются транспонированными.

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

cD1 = d1k, cA1 = a1k.

Правая часть схемы рис. 2.13 осуществляет вейвлет реконструкцию сигнала. Эта процедура использует операции интерполяции и фильтрации фильтрами реконструкции Lo _ R и Hi _ R. Операция интерполяции 2, обратная децимации 2, осуществляется путем увеличения в два раза числа составляю щих добавлением нулевых компонентов вперемежку с имею % % щимися. При сложении сигналов ( A1 и D1 ), полученных на вы % ходе фильтров Lo _ R и Hi _ R, будем иметь сигнал S (k ), близкий к исходному S (k ), т.е. произойдет его реконструкция на начальном уровне.

Для последующей итерации ( m = 2) используются значения a1k с предыдущей и т.д. до m = MMAX.

Схема многошаговой итерационной процедуры анализа син теза показана на рис. 2.14, где представлены диаграммы (а) и структура (б) многошагового алгоритма декомпозиции и рекон струкции сигнала, называемого алгоритмом Малла (Mallat).

Здесь для наглядности сигнал представлен 512 отсчетами ( N = 2n0, n0 = 9 ).

cDm 2 Hi _ R Hi _ D 2 cDm (+) % S S m Lo _ D 2 cAm cAm 2 Lo _ R а Вейвлет- Секция синтеза Секция анализа коэффициенты (реконструкции) (декомпозиции) 2 H H 2 % S (k 512 CD CD1 D L 2 H H 2 2 L A1 cA 2 L L б Рис. 2. Таким образом, БВП во временной и частотной областях – это две стороны единой многошаговой структуры, позволяющей быстро осуществить как декомпозицию, так и реконструкцию сигнала.

Следовательно, при частотном подходе к дискретному ВП можно использовать прежние функции, например, (2.18)–(2.20), но вместо имени вейвлета (‘wname’) в качестве входного аргу мента должны быть заданы соответствующие НЧ и ВЧ фильтры разложения и восстановления: Lo _ D, Hi _ D, Lo _ R, Hi _ R.

Основные функции дискретного ВП в пакете Wavelet Tool box приведены в прил. 2. Для просмотра коэффициентов фильтров (совместно с вейвлетом) достаточно исполнить ко манду «wavemenu» и в появившемся окне с описанием разделов ВП нажать кнопку «Wavelet Display». Выводится следующее ок но, в котором, выбрав имя, «wname», можно просмотреть коэф фициенты фильтров декомпозиции ( Lo _ D – low-pass, Hi _ D – high-pass) и реконструкции ( Lo _ R, Hi _ R ). На рис. П.1 дан пример для вейвлета db 4. Количественные данные о вейвлет фильтрах можно получить в командном режиме с помощью простых команд, например load, wname, length (длина вектора коэффициентов), sum (сумма коэффициентов), norm (норма вектора коэффициентов) и др. Например, загружая командой load фильтр Добеши db 4, получим:

load db db db4 = 0.1629 0.5055 0.4461 – 0.0198 – 0.1323 0.0218 0.0233 – 0. length (db4) ans = sum(db4) ans = 1. norm(db4) ans = 0. hl даны с учетом норми Примечание. Коэффициенты вейвлет-фильтра ровки – множителя 1/ 2.

Следующий пример командой load загружает вейвлет db 4, строит его график и графики вейвлет-коэффициентов hl и коэф фициентов фильтров декомпозиции и реконструкции (рис. 2.16):

function db load db4;

w = db4;

iter = 10;

wav = 'db4';

[phi,psi,xval] = wavefun(wav,iter);

subplot(321);

plot(xval,psi);

title('Wavelet');

Рис. 2. subplot(322);

stem(w);

title('Original scaling filter');

[Lo_D, Hi_D, Lo_R, Hi_R] = orthfilt(w);

subplot(323);

stem(Lo_D);

title('Decomposition low-pass filter');

subplot(324);

stem(Hi_D);

title('Decomposition high-pass filter');

subplot(325);

stem(Lo_R);

title('Reconstruction low-pass filter');

subplot(326);

stem(Hi_R);

title('Reconstruction high-pass filter');

end 2.7. Пакетные вейвлеты и вейвлет-алгоритмы При рассмотрении БВП по алгоритму Малла на каждом шаге происходит октавополосное «расщепление» (splitting) сиг нала на ВЧ и НЧ составляющие и «отсечение» ВЧ составляю щей. Причина такого подхода заключена в неявном предполо жении, что НЧ область содержит больше информации об исходном сигнале, чем ВЧ область. В результате получается «однобокое» дерево (рис. 2.10 ). Такое предположение оправда но для многих реальных сигналов, однако для некоторых оно не выполняется.

Р. Койфман и М. Викерхаузер усовершенствовали алгоритм Малла, предложив применить процесс расщепления как для НЧ, так и ВЧ составляющих сигнала. В результате получается «пол ное» (бинарное или сбалансированное) дерево, представленное на рис. 2.16, а.

H эквивалентно L S S (0,0) (0,0) D A A1 D (1,1) (1,0) (1,0) AA2 DA AD2 DD2 AA DA (2,1) (2,1) (2,2) (2,0) (2,3) DDA ADA AAA DDD DAA а б Рис. 2. Ветвям дерева будет соответствовать набор подпространств сигнала с базисами, построенными, как и для однобокого дере ва, согласно КМА. Функции и фильтры, порождающие эти бази сы, называются соответственно вейвлет-пакетами и пакетными фильтрами.

На рис. 2.17 в качестве примера приведены начальные па кетные вейвлеты для функции Добеши db 4.

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

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

Разработан ряд методов для выбора оптимального или квазиоп тимального дерева.

Рис. 2. В качестве функции стоимости используется энтропия. Су ществует много разных определений энтропии, например, опре деление по Шеннону (Shannon) E = Si2 log( Si2 ), (2.22) i или через логарифм энергии E = log( Si2 ).

i Однако суть у них одна – большая энтропия свидетельствует о «размазанности» сигнала по базисным функциям;

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

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

Рис. 2. Рис. 2.18 иллюстрирует наилучшее дерево (слева) по крите рию энтропии, выведенное на экран при анализе звукового сиг нала (см. пример 2.4):

Function ss_tree load mtlb;

x = mtlb(1:200);

wpt = wpdec(x,3,'db1');

wpt = wpsplt(wpt, [3 0]);

plot(wpt) bst = besttree(wpt);

plot(bst);

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

Возможен упрощенный вариант [7, 8], состоящий в подборе оптимальной высоты (уровня) полного дерева, при которой эн тропия минимальна. Известны и другие алгоритмы ВП с ис пользованием вейвлет-пакетов [1, 3, 7, 8].

Пакетные вейвлет-алгоритмы встроены в Wavelet Toolbox.

Функции одно– и двумерного пакетного вейвлет-анализа и син теза, вычисления энтропии, определения наилучшего дерева по критерию энтропии и другим приведены в прил. 2 (П.2.5).

2.8. Удаление шумов и компрессия сигналов Традиционно для решения этих задач применяется из вестный из практики фильтрации метод подавления высокочас тотных составляющих спектра. Этот метод был использован в примерах 2.2 и 2.3.

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

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

Подробно познакомиться с этими проблемами и приобрести определенные навыки читатель может с помощью демонстраци онных примеров, работая с подразделами De-noise и Compres sion интерфейса GUI.

Для практических применений Wavelet Toolbox имеет ряд соответствующих функций (см. прил. П.2.7). Использование не которых из них проиллюстрировано в приводимых ниже при мерах.

Пример 2.5. Установка мягкого или жесткого порогов Осуществляется с помощью функции y = wtresh( X, SORN, T ). Она воз вращает жесткий ythard ( SORN = ' h ') или мягкий ytsoft ( SORN = ' s ') порог ( thersold ) T для входного вектора или матрицы X.

y = linspace(-1,1,100);

thr = 0.4;

ythard = wthresh(y,'h',thr);

ytsoft = wthresh(y,'s',thr);

subplot(131), plot(y);

title('No thersold') subplot(132), plot(ythard);

title('thard thersold') subplot(133), plot(ytsoft);

title('tsoft thersold') На рис. 2.19 приведен вид зависимости y = f ( x), полученной с помощью функции y = wtresh( X, SORN, T ), как при отсутствии порога, так и при жест ком и мягком порогах.

Рис. 2. Пример 2.6. Очистка от шума тестового сигнала с использованием функции wdencmp В памяти компьютера имеется шесть тестовых зашумленных сигналов. За грузим один из них из файла noismima и очистим его с помощью функции удаления шума и сжатия wdencmp (прил. П.2.7). Функция [xd, cxd, lxd, perf0, perfl2] = wdencmp ('lvd',c,l,wname,lev,thr,'h') возвращает очищенный и сжатый сигнал XC, полученный из исходного сигнала X (как одномерного, так и двумерного) с использованием глобального порога THR. Выходные аргумен ты [CXD, LXD] – структура вейвлет-разложения вектора XC. PERFO и L2 – нормы восстановления и сжатия в процентах.

PERFL 2 – PERFL 2 = 100 * ( norm(C C ) / norm(C ))2, где norm – норма вектора;

для од 2 номерного сигнала PERFL 2 = 100 XC / X. N – уровень вейвлет-разложе ния. SORH ( ' s ' или ' h ' ) – установка гибкого или жесткого порога.

Процедура удаления шума и сжатия включает в себя три шага:

1) Разложение для уровня lev. Выполняется функцией wavedec(x, lev, wname);

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



Pages:   || 2 |
 














 
2013 www.netess.ru - «Бесплатная библиотека авторефератов кандидатских и докторских диссертаций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.