Метод вейвлета пакета является обобщением разложения вейвлета, которое предлагает более богатый анализ сигнала.
Атомы вейвлета пакета являются формами волны, индексируемыми тремя естественно интерпретируемыми параметрами: положением, масштабом (как в вейвлет разложении) и частотой.
Для заданной функции ортогонального вейвлета мы генерируем библиотеку основ, называемую вейвлетом packet основ. Каждая из этих основ предлагает конкретный способ кодирования сигналов, сохранения глобальной энергии и восстановления точных функций. Вейвлет могут использоваться для многочисленных расширений заданного сигнала. Затем мы выбираем наиболее подходящее разложение данного сигнала относительно основанного на энтропии критерия.
Существуют простые и эффективные алгоритмы как для вейвлета разложения пакетов, так и для оптимального выбора разложения. Затем мы можем создать адаптивные алгоритмы фильтрации с прямыми приложениями в оптимальном кодировании сигнала и сжатии данных.
В процедуре ортогонального разложения вейвлетов общий шаг разделяет коэффициенты приближения на две части. После разделения мы получаем вектор приближения коэффициентов и вектор детальных коэффициентов, оба в более грубой шкале. Информация, потерянная между двумя последовательными приближениями, получена в коэффициентах детализации. Затем следующий шаг состоит из разделения нового вектора коэффициента приближения; последующие детали никогда не повторяются.
В соответствующей ситуации с вейвлетом пакетом каждый вектор коэффициента детализации также разлагается на две части с использованием того же подхода, что и при приближении векторном разделении. Это предлагает богатейший анализ: полное двоичное дерево получается как показано на следующем рисунке.
Вейвлет разложения вейвлет-пакетов на уровне 3
Идея этого разложения состоит в том, чтобы начать с масштабно-ориентированного разложения, а затем анализировать полученные сигналы на частотных поддиапазонах.
Следующие простые примеры иллюстрируют определенные различия между вейвлет и вейвлет-пакетным анализом.
Спектральный анализ широкополосных стационарных сигналов с помощью преобразования Фурье хорошо установлен. Для нестационарных сигналов существуют локальные методы Фурье, такие как краткосрочное преобразование Фурье (STFT). Краткое описание смотрите в Кратковременном Преобразовании Фурье.
Поскольку вейвлеты локализованы по времени и частоте, возможно использовать основанные на вейвлетах аналоги STFT для частотно-временного анализа нестационарных сигналов. Для примера возможно создать скалограмму (wscalogram
) на основе непрерывного вейвлет (CWT). Однако потенциальным недостатком использования CWT является то, что он является в вычислительном отношении дорогим.
Дискретное вейвлет (DWT) позволяет частотно-временное разложение входного сигнала, но степень разрешения частоты в DWT обычно рассматривается как слишком грубая для практического частотно-временного анализа.
В качестве компромисса между методами, основанными на DWT- и CWT, вейвлет обеспечивают вычислительно эффективную альтернативу с достаточным разрешением частоты. Вы можете использовать wpspectrum
для выполнения частотно-временного анализа вашего сигнала с помощью вейвлет.
Следующие примеры иллюстрируют использование вейвлета пакетов для выполнения локального спектрального анализа. Следующие примеры также используют spectrogram
(Signal Processing Toolbox) из программного обеспечения Signal Processing Toolbox™ в качестве бенчмарка для сравнения со спектром вейвлет-пакета. Если у вас нет программного обеспечения Signal Processing Toolbox, можно просто запустить вейвлет примеры спектра пакетов.
Вейвлет пакетный спектр синусоиды.
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');
Если у вас есть программное обеспечение Signal Processing Toolbox, можно вычислить краткосрочное преобразование Фурье.
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')
Сумма двух синусоид с частотами 64 и 128 герц.
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');
Если у вас есть программное обеспечение Signal Processing Toolbox, можно вычислить краткосрочное преобразование Фурье.
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')
Сигнал с резким изменением частоты от 16 до 64 герц на две секунды.
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');
Если у вас есть программное обеспечение Signal Processing Toolbox, можно вычислить краткосрочное преобразование Фурье.
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');
Если у вас есть программное обеспечение Signal Processing Toolbox, можно вычислить краткосрочное преобразование Фурье.
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');
Если у вас есть программное обеспечение Signal Processing Toolbox, можно вычислить краткосрочное преобразование Фурье.
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')
Расчет для генерации вейвлет-пакетов прост при использовании ортогонального вейвлета. Начнем с двух фильтров длины 2 N, где h (n)
и g (n)
, соответствуют вейвлет.
Теперь по индукции давайте зададим следующую последовательность функций:
(W n (x), n = 0, 1, 2 ,...)
около
где W 0 (x)
=, (x)
- функция масштабирования и W 1 (x)
=, (x)
является вейвлет.
Для примера для вейвлета Haar мы имеем
и
Уравнения становятся
и
<reservedrangesplaceholder1> 0 <reservedrangesplaceholder0> =,
(x)
- функция масштабирования Haar и W 1 (x)
=, (x)
- вейвлет Haar, оба поддерживаемые в [0, 1]. Затем мы можем получить W 2 n путем добавления двух 1/2-масштабированных версий Wn с различными поддержкой [0,1/2] и [1/2,1] и получить W 2 n + 1 путем вычитания тех же версий Wn.
Для n = 0 до 7 у нас есть W-функции, показанные на рисунке Haar Wavelet Packets.
Пакеты Вейвлет
Это можно получить с помощью следующей команды:
[wfun,xgrid] = wpfun('db1',7,5);
который возвращается в wfun
приблизительные значения Wn для n = 0 до 7, вычисленные на 1/25 сетка опорной xgrid
.
Начиная с более регулярных исходных вейвлетов и используя подобную конструкцию, мы получаем сглаженные версии этой системы W-функций, все с поддержкой в интервале [0, 2 N -1]. Рисунок db2 Wavelet Packets представляет систему W-функций для исходного db2
вейвлет.
db2 Вейвлет
Начиная с функций и следуя той же линии, ведущей к ортогональным вейвлетам, рассмотрим трехиндексированное семейство функций анализа (формы волны):
где <reservedrangesplaceholder4> <reservedrangesplaceholder3> и (j, k) <reservedrangesplaceholder0>2.
Как и в вейвлет среды, k можно интерпретировать как параметр локализации во времени и j как параметр шкалы. Так какова же интерпретация n?
Основная идея вейвлет заключается в том, что для фиксированных значений j и k Wj,n,k анализирует колебания сигнала примерно вокруг положения 2j· k, в шкале 2j и на различных частотах для различных допустимых значений последнего n параметра.
Фактически, внимательно исследуя вейвлет, отображенные в Haar Wavelet Packets и db2 Wavelet Packets, естественно упорядоченное Wn для n = 0, 1,..., 7, не соответствует точно порядку, заданному количеством колебаний. Точнее, подсчет количества нулевых переходов (вверх-переходы и вниз-переходы) для db1
вейвлет, у нас есть следующее.
Так, чтобы восстановить свойство, что основная частота увеличивается монотонно с порядком, удобно задать частотный порядок, полученный из естественной рекурсивно.
Как видно из предыдущих рисунков, Wr (n) (x)
«колеблется» примерно n раз.
Чтобы анализировать сигнал (щебет примера 2 для образца), лучше построить график коэффициентов пакета вейвлетов, следующих за порядком частот от низких частот в нижней части до высоких частот в верхней части, чем естественно упорядоченные коэффициенты.
При построении графика коэффициентов с помощью приложения Wavelet Analyzer доступны различные опции, связанные с выбором порядка «Частота» или «Естественный».
Эти опции также доступны в командной строке при использовании wpviewcf
функция.
Набор функций W j, n = (Wj,n,k(x),k∊Z)
- пакет (j, n) вейвлет. Для положительных значений целых и n целых чисел j вейвлет организованы в деревьях. Дерево на рисунке Wavelet Packets Organized in a Tree; Шкала j Задаёт Глубину и Частота n Задает Положение в дереве Создается, чтобы получить максимальное разложение уровня, равное 3. Для каждого масштабного j возможные значения параметра n равны 0, 1,..., 2j–1.
Вейвлет, организованные в дереве; Шкала j Определение глубины и частоты n определяет положение в дереве
Обозначение W j,n, где j обозначает параметр шкалы и n частотный параметр, согласуется с обычной меткой дерева положения глубины.
У нас есть , и .
Оказывается, библиотека вейвлета пакетных основ содержит вейвлет базиса а также несколько другие основы. Давайте рассмотрим некоторые из этих основ. Точнее, позвольте V 0 обозначить пространство (охватываемое семейством W 0,0), в котором лежит сигнал, подлежащий анализу; затем (W d, 1; d ≥ 1) является ортогональным базисом V 0.
Для каждого строго положительного целого числа D (W D, 0, (W d, 1; 1 ≤ d ≤ D)) является ортогональным базисом V 0.
Мы также знаем, что семейство функций {(Wj + 1,2 n), (Wj + 1,2 n + 1)} является ортогональным базисом пространства, охватываемого Wj,n, которое разделено на два подпространства: Wj + 1,2 n охватывает первое подпространство, а Wj + 1,2 n + 1 второе.
Это последнее свойство дает точную интерпретацию разделения в дереве организации вейвлет, потому что все разработанные узлы имеют форму, показанную на рисунке Wavelet Packet Tree: Split и Merge.
Из этого следует, что листья каждого связанного двоичного поддерева полного дерева соответствуют ортогональному базисный исходного пространства.
Для сигнала с конечной энергией, принадлежащего V 0, любой базис пакета вейвлет обеспечит точную реконструкцию и предложит конкретный способ кодирования сигнала, используя распределение информации в поддиапазонах шкалы частот.
Основываясь на организации библиотеки вейвлет-пакетов, естественно считать декомпозиции, выданные из данного ортогонального вейвлета.
Сигнал длины N = 2L может быть расширен α различными способами, где α - количество двоичных подтревов полного двоичного дерева глубин L. В результате, (см. [Mal98] стр. 323).
Поскольку это число может быть очень большим, и поскольку явное перечисление в целом неуправляемо, интересно найти оптимальное разложение относительно удобного критерия, вычисляемого эффективным алгоритмом. Мы ищем минимум критерия.
Функции, проверяющие свойство типа аддитивности, хорошо подходят для эффективного поиска двоичных древовидных структур и фундаментального разделения. Классические основанные на энтропии критерии соответствуют этим условиям и описывают связанные с информацией свойства для точного представления данного сигнала. Энтропия является общей концепцией во многих областях, в основном в обработке сигналов. Давайте перечислим четыре различных критерия энтропии (см. [CoiW92]); многие другие доступны и могут быть легко интегрированы (тип help
wentropy
). В следующих выражениях s является сигналом и (si) являются коэффициентами s в ортонормированном базисе.
Энтропийная E должна быть аддитивной функцией стоимости, такой что E (0) = 0 и
Энтропия (ненормализованная) Шеннона
так
с условием 0log (0) = 0.
Концентрация в лp норма с 1 ℜ ≤ p
так
Логарифм «энергетической» энтропии
так
с журналом соглашения (0) = 0.
Пороговая энтропия
если и 0 в другом месте, так E 4 (s) = # {i, что } - количество моментов времени, когда сигнал больше порога
Эти функции энтропии доступны с помощью wentropy
файл.
Сгенерируйте сигнал энергии, равный 1.
s = ones(1,16)*0.25;
Вычислите энтропию Шеннона s.
e1 = wentropy(s,'shannon') e1 = 2.7726
Вычислите l1.5 энтропия s, эквивалентная norm(s,1.5)
1.5.
e2 = wentropy(s,'norm',1.5) e2 = 2
Вычислите энтропию «логарифмической энергии» s.
e3 = wentropy(s,'log energy') e3 = -44.3614
Вычислите пороговую энтропию s с помощью порогового значения 0,24.
e4 = wentropy(s,'threshold', 0.24) e4 = 16
Этот простой пример иллюстрирует использование энтропии, чтобы определить, представляет ли интерес новое расщепление для получения минимального энтропийного разложения.
Начнем с постоянного исходного сигнала. Двух частей информации достаточно для определения и восстановления сигнала (то есть длины и постоянного значения).
w00 = ones(1,16)*0.25;
Вычислите энтропию исходного сигнала.
e00 = wentropy(w00,'shannon') e00 = 2.7726
Затем разделите w00
использование воловьего вейвлета.
[w10,w11] = dwt(w00,'db1');
Вычислите энтропию приближения на уровне 1.
e10 = wentropy(w10,'shannon') e10 = 2.0794
Детали уровня 1, w11
, равен нулю; энтропия e11
равен нулю. Из-за свойства аддитивности энтропия разложения определяется e10+e11=2.0794
. Это необходимо сравнить с начальной энтропией e00=2.7726
. У нас есть e10 + e11 < e00
, так что расщепление интересно.
Теперь разделите w10
(не w11
потому что разделение вектора null не имеет интереса, поскольку энтропия равна нулю).
[w20,w21] = dwt(w10,'db1');
У нас есть w20=0.5*ones(1,4)
и w21
равен нулю. Энтропия уровня приближения 2
e20 = wentropy(w20,'shannon') e20 = 1.3863
Снова у нас есть e20 + 0 < e10
таким образом, расщепление приводит к уменьшению энтропии.
Тогда
[w30,w31] = dwt(w20,'db1'); e30 = wentropy(w30,'shannon') e30 = 0.6931 [w40,w41] = dwt(w30,'db1') w40 = 1.0000 w41 = 0 e40 = wentropy(w40,'shannon') e40 = 0
В последней операции разделения мы обнаруживаем, что для восстановления исходного сигнала необходима только одна часть информации. Вейвлет на уровне 4 является лучшим базисом согласно энтропии Шеннона (с нуль оптимальной энтропии с e40+e41+e31+e21+e11 = 0
).
Выполните вейвлет пакетов разложение сигнального s, заданного в примере 1.
t = wpdec(s,4,'haar','shannon');
Дерево вейвлета пакетов в Значения Entropy показывает узлы, маркированные исходными номерами энтропии.
Значения энтропии
Вычислите лучшее дерево.
bt = besttree(t);
Лучшее дерево показано на следующем рисунке. В этом случае лучшее дерево соответствует вейвлет. Узлы маркированы оптимальной энтропией.
Оптимальные значения энтропии
Использование вейвлет требует древовидных действий и маркировки. Реализация пользовательского интерфейса построена вокруг этого фактора. Для получения дополнительной информации о технических деталях см. страницы с описанием.
Полное двоичное дерево глубин D соответствующий вейвлет дереву разложения пакетов, разработанному на уровне D, обозначено WPT.
У нас есть следующие интересные поддеревья.
Дерево разложения | Поддерево, такое, что набор листьев является базисом |
---|---|
Вейвлет разложения вейвлет-пакетов | Полное двоичное дерево: WPT глубины D |
Оптимальное дерево разложения вейвлет-пакетов | Двоичный поддерево WPT |
Вейвлет дерева наилучшего уровня | Полный двоичный поддерево WPT |
Дерево разложения вейвлетов | Левый односторонний двоичный поддерево WPT глубины D |
Вейвлет лучшее базисное дерево | Левый односторонний двоичный поддерево WPT |
Мы выводим следующие определения оптимальных разложений, относительно E критерия энтропии.
Разложения | Оптимальное разложение | |
---|---|---|
Вейвлет вейвлет-пакетов | Поиск среди 2D деревья | Поиск среди D деревьев |
Вейвлет | Поиск среди D деревьев | Поиск среди D деревьев |
Для любого нетерминального узла мы используем следующий основной шаг, чтобы найти оптимальное поддерево относительно заданного критерия энтропии E (где Eopt обозначает оптимальное значение энтропии).
Энтропийные Условия | Действие с древовидной и энтропийной маркировкой |
---|---|
| Если (node≠root), объедините и установите Eopt (node) = E (node) |
Разделение и установка |
с естественным начальным условием на ссылку дереве, Eopt (t)
= <reservedrangesplaceholder1> <reservedrangesplaceholder0> для каждого терминального узла t.
Вы можете использовать функцию wprcoef
чтобы восстановить приближение вашего сигнала от любого узла в дереве вейвлет. Это верно независимо от того, работаете ли вы с полным вейвлетом деревом пакетов или поддеревом, определяемым критерием оптимальности. Использовать wpcoef
если вы хотите извлечь вейвлет пакетные коэффициенты из узла, не восстанавливая приближения к сигналу.
Загрузите шумный сигнал Доплера.
load noisdopp
Вычислите разложение вейвлета пакета до уровня 5 с помощью sym4
вейвлет. Используйте режим периодизации.
dwtmode('per'); T = wpdec(noisdopp,5,'sym4'); plot(T)
Постройте график дерева двоичных вейвлетов пакетов и щелкните по (4,1) дуплету (узел 16).
Извлеките вейвлет пакетные коэффициенты из узла 16.
wpc = wpcoef(T,16); % wpc is length 64
Получите приближение к сигналу от узла 16.
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)
Восстановите приближение к сигналу от (3,0) дублета (узел 7).
rsig = wprcoef(Topt,7); % rsig is length 1024 plot(noisdopp,'k'); hold on; plot(rsig,'b','linewidth',2); axis tight;
Если вы знаете, какой дуплет в двоичном дереве вейвлет вы хотите извлечь, можно определить узел, относящийся к этому дуплету с depo2ind
.
Для примера, чтобы определить узел, относящийся к дублету (3,0), введите:
Node = depo2ind(2,[3 0]);
Точно так же, как в случае разложения вейвлет, предыдущая 1-D среда может быть расширена до анализа изображений. Незначительные прямые изменения приводят к четвертичным древовидным определениям. Пример показан на следующем рисунке для глубины 2.
Четвертичное дерево глубины 2
В вейвлет пакетной среде идеи сжатия и шумоподавления идентичны идеям, разработанным в вейвлет среде. Единственной новой возможностью является более полный анализ, который обеспечивает повышенную гибкость. Одно разложение с использованием вейвлет генерирует большое количество основ. Затем можно искать лучшее представление относительно цели проекта, используя besttree
с функцией энтропии.