Документация

От вейвлетов к вейвлетным пакетам

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

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

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

Вейвлет-пакеты в действии: введение

Следующие простые примеры иллюстрируют некоторые различия между вейвлет-анализом и вейвлет-анализом пакетов.

fs = 1000; % sampling rate
t = 0:1/fs:2; % 2 secs at 1kHz sample rate
y = sin(256*pi*t); % sine of period 128
level = 6;
wpt = wpdec(y,level,'sym8');
[Spec,Time,Freq] = wpspectrum(wpt,fs,'plot');

figure;
windowsize = 128;
window = hanning(windowsize);
nfft = windowsize;
noverlap = windowsize-1;
[S,F,T] = spectrogram(y,window,noverlap,nfft,fs);
imagesc(T,F,log10(abs(S)))
set(gca,'YDir','Normal')
xlabel('Time (secs)')
ylabel('Freq (Hz)')
title('Short-time Fourier Transform spectrum')

fs = 1000;
t = 0:1/fs:2;
y = sin(128*pi*t) + sin(256*pi*t); % sine of periods 64 and 128.
level = 6;
wpt = wpdec(y,level,'sym8');
[Spec,Time,Freq] = wpspectrum(wpt,fs,'plot');

figure;
windowsize = 128;
window = hanning(windowsize);
nfft = windowsize;
noverlap = windowsize-1;
[S,F,T] = spectrogram(y,window,noverlap,nfft,fs);
imagesc(T,F,log10(abs(S)))
set(gca,'YDir','Normal')
xlabel('Time (secs)')
ylabel('Freq (Hz)')
title('Short-time Fourier Transform spectrum')

fs = 500;
t = 0:1/fs:4;
y = sin(32*pi*t).*(t<2) + sin(128*pi*t).*(t>=2);
level = 6;
wpt = wpdec(y,level,'sym8');
[Spec,Time,Freq] = wpspectrum(wpt,fs,'plot');

figure;
windowsize = 128;
window = hanning(windowsize);
nfft = windowsize;
noverlap = windowsize-1;
[S,F,T] = spectrogram(y,window,noverlap,nfft,fs);
imagesc(T,F,log10(abs(S)))
set(gca,'YDir','Normal')
xlabel('Time (secs)')
ylabel('Freq (Hz)')
title('Short-time Fourier Transform spectrum')

fs = 1000;
t = 0:1/fs:2;
y = sin(256*pi*t.^2);
level = 6;
wpt = wpdec(y,level,'sym8');
[Spec,Time,Freq] = wpspectrum(wpt,fs,'plot');

figure;
windowsize = 128;
window = hanning(windowsize);
nfft = windowsize;
noverlap = windowsize-1;
[S,F,T] = spectrogram(y,window,noverlap,nfft,fs);
imagesc(T,F,log10(abs(S)))
set(gca,'YDir','Normal')
xlabel('Time (secs)')
ylabel('Freq (Hz)')
title('Short-time Fourier Transform spectrum')

y = wnoise('quadchirp',10);
len = length(y);
t = linspace(0,5,len);
fs = 1/t(2);
level = 6;
wpt = wpdec(y,level,'sym8');
[Spec,Time,Freq] = wpspectrum(wpt,fs,'plot');

windowsize = 128;
window = hanning(windowsize);
nfft = windowsize;
noverlap = windowsize-1;
imagesc(T,F,log10(abs(S)))
set(gca,'YDir','Normal')
xlabel('Time (secs)')
ylabel('Freq (Hz)')
title('Short-time Fourier Transform spectrum')

Построение вейвлет-пакетов

Схема вычисления для генерации вейвлет-пакетов проста при использовании ортогонального вейвлета. Начнем с двух фильтров длины 2N, где h(n) и g(n), соответствуют вейвлету.

Теперь по индукции определим следующую последовательность функций:

около

где W0(x) = φ(x) является функцией масштабирования и W1(x) = ψ(x) - вейвлет-функция.

Например, для вейвлета Хаара, который у нас есть

и

Уравнения становятся

и

W0(x) = φ(x) является функцией масштабирования Хаара и W1(x) = ψ(x) является вейвлетом Хаара, поддерживаемым в [0, 1]. Затем мы можем получить W2n путем добавления двух 1/2 масштабированных версий _Wn с различными опорами [0,1/2] и [1/2,1] и получить W2n ₊ 1 путем вычитания тех же версий Wn.

Для n = от 0 до 7 мы имеем W-функции, показанные на рисунке Пакеты вейвлета Haar.

Это можно получить с помощью следующей команды:

[wfun,xgrid] = wpfun('db1',7,5);

который возвращается в wfun приблизительные значения _Wn для n = от 0 до 7, вычисленные на 1/25 сетке опоры xgrid.

Начиная с более регулярных оригинальных вейвлетов и используя аналогичную конструкцию, получаем сглаженные версии этой системы W-функций, все с поддержкой в интервале [0, 2N-1]. На рисунке db2 Vavelet Packets представлена система W-функций для исходного db2 вейвлет.

Атомы вейвлет-пакетов

Начиная с функций $$({}_{\mathrm{Wn}}(x),$$n∈N) и следуя той же строке, ведущей к ортогональным вейвлетам, рассмотрим трёхиндексное семейство анализирующих функций (формы волн):

где n∊N и (j, k) ^∊Z2.

Как и в вейвлет-фреймворке, k можно интерпретировать как параметр локализации времени, а j как параметр масштаба. Так какова интерпретация n?

Основная идея вейвлет-пакетов состоит в том, что для фиксированных значений j и k Wj, _n, k анализирует флуктуации сигнала приблизительно вокруг положения 2j· k, на шкале 2j и на различных частотах для различных допустимых значений последнего параметра n.

На самом деле, тщательно исследуя вейвлет-пакеты, отображаемые в пакетах вейвлета Haar и пакетах вейвлета db2, естественно упорядоченный _Wn для n = 0, 1,..., 7, не совпадает точно с порядком, определенным числом колебаний. Точнее, подсчет количества нулевых переходов (пересечений вверх и вниз) для db1 вейвлет-пакеты, у нас есть следующее.

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

Как видно на предыдущих чертежах, _{Wr (n})(x) «колеблется» приблизительно n раз.

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

При построении графиков коэффициентов различные опции, связанные с выбором порядка «Частота» или «Естественный», доступны с помощью приложения Wavelet Analyzer.

Эти параметры также доступны в режиме командной строки при использовании wpviewcf функция.

Организация вейвлет-пакетов

Набор функций Wj_, n =(Wj,n,k(x),k∊Z) является (j, n) вейвлет-пакетом. Для положительных значений целых чисел j и n вейвлет-пакеты организованы в деревьях. Дерево на рисунке Vavelet Packets Organized in a Tree; Масштаб j Определяет глубину и частоту n Определяет положение в дереве, чтобы получить максимальное разложение уровня, равное 3. Для каждой шкалы j возможные значения параметра n равны 0, 1,..., ^2j-1.

Обозначение _Wj, n, где j обозначает параметр масштаба и n параметр частоты, согласуется с обычной маркировкой дерева положения глубины.

У нас $${}_{\mathrm{}\mathrm{есть\; W0,0}}=((x-k)$$, k∈Z) и $${}_{\mathrm{}\mathrm{}}\frac{}{\mathrm{W1,1}}=((\mathrm{x2}$$− k), k∈Z).

Оказывается, что библиотека баз вейвлет-пакетов содержит базис вейвлета, а также несколько других баз. Давайте посмотрим на некоторые из этих баз. Точнее, пусть V0 обозначает пространство (охватываемое W0,0 семейства), в котором лежит анализируемый сигнал; затем (Wd,1; d ≥ 1) - ортогональный базис V0.

Для каждого строго положительного целого числа D, (WD,0, (Wd,1; 1 ≤ d ≤ D)) - ортогональный базис V0.

Известно также, что семейство функций {(_{Wj + 1,} 2n), _{(Wj + 1}, 2n + 1)} является ортогональным базисом пространства, _{охватываемого} Wj, n, которое разделено на два _{подпространства}: Wj + 1, 2n охватывает первое _{подпространство, и} Wj + 1, 2n + 1 второе.

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

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

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

Выбор оптимальной декомпозиции

На основе организации библиотеки вейвлет-пакетов естественно подсчитывать разложения, выданные из данного ортогонального вейвлета.

Сигнал длиной N = ^2L может быть расширен α различными способами, где α - число двоичных поддеревьев полного двоичного дерева глубины L. В результате $${\mathrm{}}^{\mathrm{\alpha \ge 2N/2}}$$ (см. [Mal98] стр. 323).

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

Функции проверки свойства типа аддитивности хорошо подходят для эффективного поиска структур двоичного дерева и фундаментального разделения. Классические критерии, основанные на энтропии, соответствуют этим условиям и описывают связанные с информацией свойства для точного представления данного сигнала. Энтропия - распространённое понятие во многих областях, главным образом в обработке сигналов. Приведем четыре различных критерия энтропии (см. [CoiW92]); многие другие доступны и могут быть легко интегрированы (тип help wentropy). В следующих выражениях s - сигнал и (_si) - коэффициенты s в ортонормированном базисе.

Энтропия E должна быть функцией аддитивной стоимости, так что E (0) = 0 и

Эти функции энтропии доступны с помощью wentropy файл.

Некоторые интересные поддеревы

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

Полное двоичное дерево глубины D, соответствующее дереву разложения вейвлет-пакетов, разработанному на уровне D, обозначается WPT.

У нас есть следующие интересные поддеревы.

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

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

с естественным исходным условием на ссылочном дереве, Eopt(t) = E(t) для каждого терминального узла t.

load noisdopp

dwtmode('per');
T = wpdec(noisdopp,5,'sym4');
plot(T)

wpc = wpcoef(T,16);
% wpc is length 64

rwpc = wprcoef(T,16);
% rwpc is length 1024
plot(noisdopp,'k'); hold on;
plot(rwpc,'b','linewidth',2);
axis tight;

Topt = besttree(T);
% plot the best tree
plot(Topt)

rsig = wprcoef(Topt,7);
% rsig is length 1024
plot(noisdopp,'k'); hold on;
plot(rsig,'b','linewidth',2);
axis tight;

Node = depo2ind(2,[3 0]);

Вейвлет-пакеты 2-D структура декомпозиции

Точно так же, как в случае вейвлет-разложения, предыдущая структура 1-D может быть расширена для анализа изображения. Незначительные прямые изменения приводят к определениям, связанным с четвертичным деревом. На следующем рисунке показан пример глубины 2.

Вейвлет-пакеты для сжатия и деноизирования

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