Вейвлеты в СКМ Mathematica 8

Вейвлеты в Mathematica 8

Для предшествующих версий СКМ Mathematica 4/5/6/7 разработчиком — фирмой Wolfram Research (США) — был создан пакет расширения Wavelet Explorer, второй по значимости после Wavelet Toolbox для MATLAB. Однако в составе самой СКМ он не поставлялся, а в России был большой редкостью. Во время стажировки автора этой статьи в США пакет Wavelet Explorer был любезно предоставлен автору и описан в книгах [1, 4, 5].

В появившейся около года назад новейшей версии СКМ Mathematica 8 основные, тщательно отобранные и переработанные средства работы с вейвлетами были впервые введены в откомпилированное ядро системы и существенно расширили возможности обработки сигналов. Они стали поставляться в составе этой СКМ и доступны всем ее пользователям. Эти средства достаточно полно отражены в обширной англоязычной справке по этой системе (рис. 1).

Рис. 1. Окно СКМ Mathematica 8 с данными о вейвлетах

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

Немного теории

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

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

В основе вейвлет-преобразований лежит использование двух непрерывных и интегрируемых по всей оси t (или x) функций:

• материнская вейвлет-функция Psi ψ(t)
  
• отцовская масштабирующая или скейлинг-функция Phi φ(t) с
  

Phi-функции, задающие грубое приближение (аппроксимацию) сигнала, присущи только ортогональным вейвлетам. А детализирующие Psi-функции создаются на основе той или иной базисной функции ψ₀(t), которая и определяет тип вейвлета. Часто Psi-функции создаются рекурсивными равенствами. Обе эти функции характерны только для ортогональных вейвлетов. Есть небольшое число вейвлетов (как правило, непрерывные), у которых Phi-функция отсутствует. Разработка вейвлетов ведется опытными математиками и отражена в их книгах [2, 3].

Psi-функция должна иметь свойства смещения во времени и масштабирования:

где a — параметр масштабирования; b — параметр сдвига.

С помощью вейвлетов возможно приближение произвольных зависимостей и сигналов s(t) с помощью ряда:

где k задает номер вейвлет-коэффициента в ходе вейвлет-приближения.

Прямое непрерывное вейвлет-преобра-зование сигнала s(t) задается вычислением вейвлет-коэффициентов по формуле:

С учетом обычно ограниченной области определения сигналов и a, b ∈ R, a ≠ 0 (4) можно представить в виде:

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

где

В настоящее время созданы десятки семейств, содержащих множество вейвлетов.

Дискретные вейвлеты

Часто используются дискретные вейвлеты, у которых в ходе преобразований, в процессе итерационных приближений заданного сигнала, ширина вейвлета на каждом шаге итераций меняется вдвое при переходе на очередной уровень L преобразований, так что a = 2^j и b = k2^j. Такое преобразование называют диадическим. Существуют алгоритмы быстрого дискретного вейвлет-преобразования как во временной, так и в частотной области.

В Mathematica 8 вейвлеты имеют имена, обычно заданные по имени ученого или ученых, впервые их описавших. Имена дискретных вейвлетов и характерный вид их Psi- (желтый цвет) и Phi-функций (красный цвет) показаны на рис. 2. Есть и частные типы вейвлетов, например с комплексными функциями, стационарные и пакетные [1–3].

Рис. 2. Семейство дискретных вейвлетов в СКМ Mathematica 8

Задающие вейвлеты функции СКМ Mathematica 8 состоят из имени вейвлета и приставки Wavelet. Например, простейшие дискретные вейвлеты Хаара (Haar) задаются функцией HaarWavelet[n], где n — порядок вейвлета. Поскольку Mathematica изначально была системой символьной (аналитической) математики, она позволяет легко находить аналитические представления того или иного вейвлета, если оно существует. Так, используя функцию для вейвлета Хаара, найдем для него аналитические выражения (рис. 3).

Рис. 3. Получение аналитических выражений для вейвлета Хaaра

Дискретное прямое вейвлет-преобразование сводится к вычислению:

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

где Z = (…–1, 0, 1…) — множество целых чисел.

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

Функции Phi и Psi, их графики и фильтры

Функции WaveletPhi и WaveletPsi позволяют в составе графической функции Plot построить графики функций Phi и Psi. Для этого используются фрагменты программ:

Их графики для вейвлета Хаара показаны на рис. 4. На рис. 2 эти графики даны для всех типов дискретных вейвлетов, определяемых по умолчанию. Большинство вейвлетов могут иметь разный порядок и образуют семейства.

Рис. 4. Графики функций Phi и Psi для вейвлета Хаара

В частотной области для характеристик вейвлетов используется фильтровый подход. Между коэффициентами вейвлетов и их фильтрами имеется однозначная связь [1–3]. Фильтры характеризуются своими коэффициентами, которые выводит функция

в формате

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

Далее мы ограничимся рассмотрением вейвлетов во временной области, опуская описанные в [1–3] детали математического характера.

Широкое применение получило обширное семейство ортогональных вейвлетов Добеши [2]. Такие вейвлеты высокого порядка имеют достаточно плавные функции Phi и Psi. Получить графики функций Phi и Psi вейвлетов Добеши порядка 8 можно с помощью следующих фрагментов программ:

Сменив цифру 8, можно получить графики для других вейвлетов этого семейства. Для изучения поведения этих функций иногда полезно их представление в виде графика, представленного на рис. 5, а иногда — в виде трехмерного графика (рис. 6).

Рис. 5. Графики функций Phi и Psi вейвлета Добеши 8

Рис. 6. Представление вейвлет-функции вейвлета Добеши в виде поверхности

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

Дискретное вейвлет-преобразование (DWT)

Дискретное вейвлет-преобразование по умолчанию выполняется с применением вейвлета Хаара. Ниже дан пример прямого дискретного преобразования для вектора из пяти элементов — чисел:

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

Восстановление происходит с точностью до машинного нуля. Команда dwd[“TreeView”] выводит дерево дискретных вейвлет-преобразований. Для приведенного выше примера оно показано на рис. 7. На графике обозначены узлы и ветви вейвлет-преобразования заданного уровня L.

Рис. 7. Дерево прямого дискретного преобразования вектора из пяти элементов

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

Рис. 8. Пример прямого и обратного дискретного вейвлет-преобразования вектора с комбинированными (символьными и численными) элементами

Как уже отмечалось, коэффициенты вейвлет-преобразований являются функциями времени t (или x). Ниже показан пример построения графика вейвлет-коэффициентов прямого преобразования функции cos(x²):

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

Рис. 9. Графики вейвлет-коэффициентов прямого преобразования функции cos (x²)

Дискретные вейвлет-преобразования звуковых сигналов и изображений

Дискретные прямые и обратные вейвлет-преобразования в СКМ Mathematica 8 возможны для звуковых сигналов. На рис. 10 сверху показан случай синтеза достаточно сложного звукового сигнала и пример его воспроизведения с помощью виртуального проигрывателя, входящего в состав СКМ Mathematica 8. А ниже показан пример прямого и обратного вейвлет-преобразования с заданным типом вейвлета — Добеши 4. Восстановленный сигнал воспроизводится вторым проигрывателем. Осциллограммы сигналов проигрывателей и создаваемые ими звуки абсолютно идентичны — так же, как и звуки, создаваемые проигрывателями. С вейвлетами связаны надежды на появление новых стандартов высококачественного воспроизведения звуков.

Рис. 10. Примеры прямого и обратного вейвлет-преобразований синтезированного звукового сигнала при использовании вейвлета Добеши

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

Рис. 11. Вейвлет-преобразования монохромного изображения

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

Рис. 12. Вейвлет-преобразования цветного изображения

Компрессия и очистка от шума при дискретном вейвлет-преобразовании

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

Рис. 13. Пример очистки от шума ступенчатого сигнала (вейвлет Хаара)

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

Рис. 14. Пример очистки от шума сложного сигнала (вейвлет Добеши 4)

Рис. 15. Пример очистки от шума изображения (часть диска планеты)

Вейвлет-сглаживание и удаление шума с изменением порога

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

Параметр tspec задается в форме tfun, {tfun, pars}, например:

На рис. 16 показан пример задания зашумленного синусоидального сигнала с линейно нарастающей амплитудой и построен его график. А рис. 17 иллюстрирует очистку этого сигнала функцией WaveletThreshold с разными типами порогов, определяющих условия удаления вейвлет-коэффициентов. Эти пороги вводятся по их имени-идентификатору и характеризуются видом их передаточной характеристики. На рис. 18 показан фрагмент программы, строящей графики передаточной характеристики различных типов порогов. Они дают хорошее представление о характере влияния на вейвлет-коэффициенты, которые, напомним, являются функциями времени.

Рис. 16. Пример задания зашумленного сигнала xsin(x)

Рис. 17. Пример очистки от шума сигнала, показанного на рис. 16, с разными типами порогов

Рис. 18. Примеры передаточных характеристик для разных типов порогов

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

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

Рис. 19. Пример задания сигнала с высокочастотной помехой и очистки от нее с помощью дискретного стационарного вейвлет-преобразования

Вейвлет-очистка изображений от шума и их компрессия

Описанные выше методы очистки от шума одномерных сигналов легко распространяются на двумерные сигналы, например монохромные изображения. Особенно эффективны методы сглаживания и компрессии изображений с адаптивным заданием порогов (рис. 20). Это возможно и в отношении трехмерных изображений, например образованных ступеньками (рис. 21). На рис. 22 представлены примеры для работы с простым трехмерным изображением, которые достаточно хорошо характеризуют реальную степень очистки и влияние остаточного шума на изображение. Пример очистки от шума изображения космонавта на поверхности Луны с помощью вейвлет-преобразования с биортогональным сплайновым вейвлетом показан на рис. 23. Функции обработки сигналов изображений в Mathematica 8 позволяют обрабатывать и цветные изображения (рис. 24).

Рис. 20. Пример задания зашумленного синусоидального сигнала и очистки его от шума с помощью адаптивного порога

Рис. 21. Пример создания простого трехмерного изображения (сверху оригинал, снизу варианты зашумленного изображения с разным уровнем шума)

Рис. 22. Примеры очистки простого трехмерного изображения от шума

Рис. 23. Пример очистки от шума реального изображения

Рис. 24. Пример обработки цветного изображения цветка

Вейвлет-компрессия и очистка сигналов от шума уже нашли массовое применение, например в стандарте сжатия фотографических изображений JPEG 2000. При большой степени сжатия (более чем в 10–15 раз) они обеспечивают лучшее качество изображений после декомпрессии и меньшее проявление артефактов. Им не присущ распад изображения на квадратики, характерный для косинусных Фурье-преобразований. Это особенно важно при сжатии телевизионных изображений в системах охранной сигнализации, при передаче сообщений со спутников и межпланетных космических станций. Для таких применений уже выпускаются специальные интегральные микросхемы, например видеокодек ADV601 фирмы Advanced Device.

Непрерывные вейвлеты и вейвлет-преобразования

Непрерывными являются вейвлеты, у которых параметры сдвига и масштабирования (Psi-функции) могут принимать произвольные значения. Семейство таких вейвлетов в СКМ Mathematica 8 наглядно показано на рис. 25. На нем указаны названия вейвлетов и их вид, определяемый видом Psi-функции (Phi-функция у этих вейвлетов отсутствует).

Рис. 25. Семейство непрерывных вейвлетов

Итак, в Mathematica 8 заданы пять типов непрерывных вейвлетов:

• MexicanHatWavelet[σ] — вейвлет «мексиканская шляпа»;
• GaborWavelet[w] — вейвлет Габора;
• DGaussianWavelet[n] — производная функции Гаусса;
• MorletWavelet[] — вейвлет Морлета;
• PaulWavelet[n] — вейвлет Паула.

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

Рис. 26. Аналитическое представление вейвлета «мексиканская шляпа» первого и второго порядков и график этого вейвлета первого порядка

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

Рис. 27. Пример прямого и обратного преобразования для нестационарного сигнала с применением непрерывного вейвлета «мексиканская шляпа»

Скейлограммы и их построение

Результаты прямого вейвлет-преобразования часто отображают с помощью скейлограмм, создаваемых функцией Scalogram. Она впервые введена в СКМ Mathematica 8. Скейлограмма строится в плоскости масштаб-время (или номер отсчета), причем уровень задается функциональной окраской с применением операции интерполяции. При построении скейлограммы задается ее цветовая схема. Пример без явного задания цветовой схемы показан на рис. 27.

К более редким вейвлетам относится вейвлет Паула порядка n (по умолчанию n = 6). Аналитическое представление и вывод графика приведены на рис. 28.

Рис. 28. Аналитическое представление и вывод графика вейвлета Паула

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

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

Рис. 29. Пример прямого и обратного вейвлет-преобразования цепочки из четырех синусоид с разной частотой (вейвлет Паула)

Рис. 30. Скейлограммы цепочки из четырех синусоид с разной частотой при двух цветовых схемах

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

Рис. 31. Задание двух пачек синусоид с колоколообразной формой их огибаюшей и построение скейлограммы такого сигнала

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

Рассмотрим применение различных непрерывных вейвлетов. На рис. 32 приведен пример создания сложного нестационарного сигнала и построение его скейлограмм с применением комплексного вейвлета Габора и разной математической обработкой его коэффициентов. Левая скейлограмма получена с использованием действительных вейвлет-коэффициентов, а правая — с вычисленным абсолютным значением каждой пары вейвлет-коэффициентов.

Рис. 32. Создание сложного нестационарного сигнала и построение его скейлограмм с применением вейвлета Габора

Создание сложного нестационарного сигнала и построение его скейлограммы с применением вейвлета Морле показано на рис. 33. Одна из составляющих сигнала — синусоида с растущей частотой. На нее наложена пачка синусоид. Скейлограмма (без указания цветовой схемы) прекрасно справляется с идентификацией этого сигнала. На ней видны обе составляющие сигнала и зависимости их частоты от времени.

Рис. 33. Создание сложного нестационарного сигнала и построение его скейлограммы с применением вейвлета Морле

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

Рис. 34. Построение скейлограммы синусоиды с коротким пиком на ней и разрывом по вертикали

Скейлограммы при дискретном вейвлет-преобразовании

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

Рис. 35. Построение скейлограммы прямоугольных импульсов с частотной модуляцией (дискретное вейвлет-преобразование)

Трехмерное представление результатов вейвлет-преобразований

В некоторых случаях полезно трехмерное представление вейвлетов. На рис. 36 приведен пример такого представления для вейвлет-преобразования с применением вейвлета «мексиканская шляпа». Подобное представление с применением вейвлета Паула показано на рис. 37.

Рис. 36. Трехмерное представление непрерывного вейвлет-преобразования с вейвлетом «мексиканская шляпа»

Рис. 37. Трехмерное представление непрерывного вейвлет-преобразования с вейвлетом Паула

Заключение

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