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

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

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

ТЕМА

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

Содержание

Введение

1. Непрерывное вейвлет - преобразование. CWT типовыми средствами Mathcad. Логарифмическая шкала масштабов. Связь шкалы масштабов с частотой. CDWT типовыми средствами Mathcad. Размер вейвлета. Процесс преобразования

2. Диадное вейвлет - преобразование. Уровень декомпозиции. Структура записи спектра. Визуализация спектра. Форма вейвлета Добеши db4.

3. Вейвлетная очистка сигналов от шумов. Подготовка преобразования. Анализ шумов по вейвлетному спектру

Введение

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

В настоящей теме рассматриваются непрерывное вейвлет-преобразование, выполняемое типовыми средствами в системе Mathcad вейвлетами произвольной формы, задаваемой пользователем, а также простое диадное вейвлет-преобразование вейвлетом Добеши db4.

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

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

CWT типовыми средствами Mathcad. Качественный анализ сигналов, особенно модельных сигналов на основе гладких финитных функций, удобно производить симметричными аналитическими вейвлетами. Так, например, базовый вейвлет MXAT конструируется на основе второй производной функции Гаусса:

Отсюда:

, (1.1)

где значения 'a' и 'b' - параметры масштаба вейвлета и сдвига, а значение коэффициента К определяется приведением нормы вейвлетной функции ||y(t,a,b)|| к 1 при нулевом среднем значении:

y(t,a,b) dt 0, y(t,a,b)2 dt 1.

Рис. 1.1.

Для MXAT-вейвлета значение К = 1.0314. Форма базового вейвлета при a=1 и b=0 приведена на рис. 1.1. Аналогично могут задаваться вейвлеты нечетного типа и вейвлеты любой другой формы.

Для демонстрации CWT зададим модельный сигнал в виде функции Гаусса на интервале 0-Т, Т=100 (рис. 1.2): s(t) = 4 exp(-(t-50)2/10).

Рис. 1.2.

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

S(a,b) = s(t), y(t,a,b = s(t) y(t,a,b) dt. (1.2)

Если сигнал s(t) является финитным и выходит практически на нулевые значения по концам интервала задания, то пределы интегрирования в (1.2) можно задавать непосредственно по интервалу задания сигнала.

Для графического просмотра результатов преобразования следует задать шаг дискретизации функции S(a,b) по параметрам 'a', 'b' и выполнить дискретизацию. Например, для заданного интервала 0-Т функции s(t):

Db = 1, i = 0..T, Da = 0.5, j = 1..200, WSi,j = S(jDa, iDb).

Графическое представление вейвлетного спектра сигнала в двух вариантах представлено на рис. 1.3.

Рис. 1.3.

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

Спектр сигнала в логарифмической шкале масштабов приведен на рис. 2.1.4.

Рис. 1.4.

На рис. 1.5 приведен пример преобразования более сложного ЛЧМ сигнала.

Рис. 1.5.

Рис. 1.6.

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

непрерывное диадное вейвлет преобразование

p(f) = y(t,1,0) cos(2pft) dt.

После определения значения f0 (рис. 1.6) частотная шкала вычисляется по масштабной шкале трансформацией f0/a f.

Рис. 1.7.

CDWT типовыми средствами Mathcad. При вычислении вейвлетных спектров в дискретной форме существенное значение для наглядности и выразительности спектров имеет шаг дискретизации параметров 'a' и 'b'. Шаг дискретизации Db, как правило, принимается равным Dt анализируемых сигналов, т.е. равен 1, как и условное значение Dt, при этом временной масштаб вейвлетного спектра соответствует временному масштабу сигнала и удобен для локализации особенностей в сигнале.

Для обеспечения возможности изменения базового размера вейвлета в его формулу множителем к параметру 'a' вводится постоянный масштабный коэффициент 'd', как это показано на рис. 1.7 для вейвлета MXAT. Средняя частота спектра базового вейвлета при а=1 не должна превышать частоты Найквиста, но должна быть больше максимальной частоты в спектре сигнала.

Размер вейвлета в единицах Dt (целочисленные значения) задается параметром L. Так как значение L = 1 (размер вейвлета 2Dt) соответствует частоте Найквиста, а частота «короткой волны» в пределах этого размера будет больше частоты Найквиста, то минимальное значение L равно 2, что обеспечивает максимальную разрешающую способность по времени. Под это значение L и рекомендуется подобрать такое максимальное значение коэффициента `d', при котором обеспечивается (с определенной точностью) выполнение условий:

y(t,1,0) dt 0, y(t,1,0)2 dt 1. (2.1.3)

Сохранение найденного значения const = L•d обеспечивает выполнение условий (2.1.3) для любой пары значений L и d. Минимальное значение L=2 действительно для CWT, но не рекомендуется для CDWT, т.к. дискретизация вейвлета порождает периодизацию его частотного спектра и перекрытие главного диапазона с боковыми, что наглядно видно на рис. 1.8.

Рис. 1.8.

Для дискретного вейвлета МХАТ при L=2 и 3 форма спектра вейвлета в области высоких частот существенно искажена, что приведет к искажению строк вейвлет-спектра. Только начиная с L=4, при средней частоте вейвлета порядка p/2 перекрытие с боковыми диапазонами становится не столь значимым.

Процесс преобразования начинается с задания интервалов и шага дискретизации по масштабу и сдвигу вейвлетов. Для наглядного сравнения с CWT повторим преобразование сигнала, приведенного на рис. 1.5, методом CDWT с заданием параметров дискретизации:

T = 200, B = T, b = 0..B, Db =1, A = 60, a = 1..A, Da = 1.

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

(2.1.4)

Рис. 1.9.

Рис. 1.10.

Примеры результатов преобразования приведены на рис. 1.9 при различных значениях параметра L для наглядного представления искажений, которые вносятся в дискретный спектр при неправильном выборе размеров базового вейвлета. При L=4 сравнением с рис. 1.5 можно убедиться в практической идентичности спектров CWT и CDWT.

На рис. 1.10 приведен пример преобразования ЛЧМ сигнала (рис. 1.5), на который наложен статистический шум, мощность которого в 3 раза больше мощности самого сигнала. Однако вид спектра хорошо фиксирует как общую картину спектра ЛЧМ сигнала, так и распределение во времени локальных неоднородностей.

2. Диадное вейвлет - преобразование

Диадное вейвлет-преобразование (DWT) одномерных сигналов в системе Mathcad выполняется вейвлетом Добеши четвертого порядка db4 функциями прямого преобразования S:= wave(s) и обратного s:= iwave(S). Значения параметров а и b задаются в виде степенных функций:

a = ао-m, b = k·ао-m, ao = 2, m, k I,

где I - пространство целых чисел, m - параметр масштаба, k - параметр сдвига. Базис пространства L2(R) в дискретном представлении:

ymk(t) = |ао|-m/2y(ао-mt-k), m,k I, y(t) L2(R). (2.1)

Вейвлет-коэффициенты прямого преобразования:

Cmk =s(t)--ymk(t) dt. (2.2)

Вейвлет Добеши относится к числу ортогональных вейвлетов и обратное преобразование для непрерывных сигналов при нормированном ортогональном вейвлетном базисе пространства:

s(t) = Cmk--ymk(t). (2.3)

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

Учитывая диадность преобразования, размер Т (в единицах Dt, t:= 0..T-1) дискретного массива st должен задаваться кратным 2M, где M - число уровней декомпозиции сигнала (m:= 1..М), при этом отсчет уровней детализирующих коэффициентов спектра начинается с m=1 (минимальный масштабный коэффициент вейвлета а=2), до m=M (a=2M). На последнем M-уровне записывается 1 детализирующий коэффициент и дополнительно конечный аппроксимирующий коэффициент. Общее количество отсчетов спектра равно количеству отсчетов сигнала. Если количество отсчетов сигнала на удовлетворяет условию 2M, то выполняется либо передискретизация сигнала с уменьшением шага Dt, либо дополнение интервала Т сигнала до ближайшего большего значения 2N. Обычно используется второй способ, т.к. с учетом односторонней формы вейвлета Добеши для исключения искажений спектра требуется задание начальных условий. Методы дополнения сигналов не отличаются от задания начальных условий при выполнении свертки сигналов.

Структура записи спектра. Все коэффициенты спектра записываются в массив спектра в следующем порядке. Размер каждой строки (уровня декомпозиции m = 1..M) занимает 2M/2m отсчетов. Каждая m - строка коэффициентов уровня начинается в массиве спектра с номера 2M-m.

Рис. 2.1.

На рис. 2.1. приведен график вейвлетного разложения импульса Кронекера в точке t=3 массива размером 128 отсчетов.

На рис. 2.2 приведен график модельного сигнала st с двумя перекрывающимися во времени синусоидами. Для просмотра коэффициентов вейвлет-преобразования S сигнала раздельно по уровням декомпозиции m можно использовать функцию формирования субматрицы c(m) = submatrix(S, r1, r2, p1, p2), которая записывает в строки m из вектора S отсчеты с номера r1 по r2 (из столбцов с p1 по p2, в данном случае столбцов нет). Пример применения функции можно видеть на рисунке.

Рис. 2.2.

Рис. 2.3.

Визуализация спектра. При визуализации картины коэффициентов субматрица переводится в двумерный массив с приведением к единой числовой оси входного сигнала (растягивание коэффициентов по оси сдвигов без изменения их значений), как это показано на рис. 2.3. Это позволяет выводить в графической форме только наиболее информативные уровни разложения, а также представлять вейвлетный спектр в 2D и в 3D форме (рис. 2.4).

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

Рис. 2.4.

Качество визуализации спектра может быть улучшено при переводе субматрицы на единую временную ось с применением методов интерполяции. В качестве примера приводится листинг подпрограммы интерполяции коэффициентов спектра с использованием кубического сплайна. Результаты применения интерполяции можно наглядно видеть на рис. 2.5 при сопоставлении с рисунками 2.3 и 2.4.

Рис. 2.5.

Рис. 2.6.

Форму вейвлета db4 и средние частоты уровней декомпозиции (масштабов вейвлета) можно определить обратным преобразованием импульса Кронекера, задаваемого на какой-либо уровень декомпозиции спектра, при обнулении всех других уровней спектра. Пример операции приведен на рис. 2.2.

Для более точного вычисления формы вейвлета число уровней декомпозиции следует устанавливать достаточно большим (N = 8-12), а импульс Кронекера в массиве спектра задавать на 2-3 уровня меньше М и на 2-3 точку интервала отсчетов этого уровня. Восстановленный при обратном преобразовании массив отсчетов вейвлета можно с использованием БПФ перевести в спектральную область и определить среднюю частоту вейвлета (с учетом вывода вейвлета на K/2 точек временной оси). При известной средней частоте вейвлета на m-уровне декомпозиции wm, частотная шкала разложения при переходе на 1 уровень ниже возрастает в 2 раза, а на 1 уровень выше - в 2 раза уменьшается. Соответственно, частотная шкала 1-го уровня (максимальная частота декомпозиции) по измерениям средней частоты вейвлета на m-уровне определяется выражением: w1 = wm 2m-1, и для вейвлета db4 составляет 2.356 радиан. Соответственно, на любом другом i-уровне: wi = w1/2i.

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

3. Вейвлетная очистка сигналов от шумов

Рис. 3.1.

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

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

Рис.3.2.

Для оценки количества возможных уровней декомпозиции, в которых будет преобладать шум, можно выполнить БПФ массива. На рис. 3.2 представлен спектр массива в интервале 0-p (256 отсчетов), т.е. на первой половине главного диапазона БПФ. При диадном делении спектра на каждом уровне декомпозиции, первый уровень детализирующих коэффициентов будет сформирован из высокочастотной части спектра сигнала от p/2 до p (в односторонней физической шкале частот). Вторая часть спектра от 0 до p/2 конвертируется в аппроксимирующие коэффициенты. На втором уровне декомпозиции аппроксимирующие коэффициенты диапазона 0-p/2 также будут разделены пополам с преобразованием диапазона p/4-p/2 в детализирующие, а диапазона 0-p/4 в аппроксимирующие коэффициенты 2-го уровня декомпозиции, и т.д. Это позволяет непосредственно по частотному спектру сигнала установить ориентировочную границу шумов и, соответственно уровни декомпозиции, в которых мощность шумов соизмерима и выше мощности сигнала.

Рис. 3.3.

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

Рис. 3.4.

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

Рис. 3.5.

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

После формирования новой строки qu данного уровня декомпозиции вейвлетного спектра можно заменить этой строкой соответствующий уровень разложения в полном массиве коэффициентов и визуально оценить результаты операции. Пример приведен на рис. 3.6.

Рис. 3.6.

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

Рис. 3.7.

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

Рис. 3.8.

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

Литература

2. Дремин И.Л. и др. Вейвлеты и их использование. / Успехи физических наук, 2001, т.171, № 5, стр. 465-501.

3. Дьяконов В., Абраменкова И. MATLAB. Обработка сигналов и изображений. Специальный справочник. - СПб.: Питер, 2002, 608 с.

5. Петухов А.П. Введение в теорию базисов всплесков. - СПб.: Изд. СПбГТУ, 1999, 132 с.

12. Смоленцев Н.К. Основы теории вейвлетов. Вейвлеты в Matlab. М.: LVR Пресс, 2005. - 304 с.

13. Переберин А.В. О систематизации вейвлет-преобразований. - Вычислительные методы и программирование, 2001, т.2.

14. Воробьев В.И., Грибунин В.Г. Теория и практика вейвлет-преобразования. - СПб, ВУС, 1999. 204 с.

39. Дьяконов В.П. Вейвлеты. От теории к практике. - М.: СОЛОН-Р, 2002. - 448 с.

Размещено на Allbest.ru