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

Дискретный вейвлет и дискретный пакет вейвлета преобразовывают

Следующий рисунок показывает дерево DWT для 1D входного сигнала.

Рисунок 1: Дерево DWT вниз к уровню 3 для 1D входного сигнала.

$$\underset{}{\overset{\sim }{G}}(f)$$ масштабирующийся (lowpass) аналитический фильтр и $$\underset{}{\overset{\sim }{H}}(f)$$ представляет вейвлет (highpass) аналитический фильтр. Метки в нижней части показывают раздел оси частоты [0,1/2] в поддиапазоны.

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

Сравните дерево DWT с полным пакетным деревом вейвлета.

Рисунок 2: Полное пакетное дерево вейвлета вниз к уровню 3.

В вейвлете пакет преобразовывает, операции фильтрации также применяются к вейвлету, или детали, коэффициентам. Результат состоит в том, что пакеты вейвлета обеспечивают фильтрацию поддиапазона входного сигнала на прогрессивно более прекрасные интервалы равной ширины. На каждом уровне, $$j$$, ось частоты [0,1/2] разделена на $$2^j$$ поддиапазоны. Поддиапазоны в герц на уровне j приблизительно

$$[\frac{nFs}{2^{j+1}},\frac{(n+1)Fs}{2^{j+1}})n=0,1,\dots 2^j-1$$ где Фс является частотой дискретизации.

То, насколько хороший это полосовое приближение, зависит от того, насколько локализованный частотой фильтры. Для фильтров Фейера-Коровкина как 'fk18' (18 коэффициентов), приближение довольно хорошо. Для фильтра как Хаар ('Хаар') приближение не точно.

В критически произведенном вейвлете пакет преобразовывает выходные параметры полосовых фильтров, прорежены два. В неподкошенном вейвлете пакет преобразовывает, выходные параметры не прорежены.

Частотно-временной анализ с пакетами вейвлета

Поскольку пакеты вейвлета делят ось частоты на более прекрасные интервалы, чем DWT, пакеты вейвлета выше в частотно-временном анализе. Как пример, рассмотрите две неустойчивых синусоиды с частотами 150 и 200 Гц в аддитивном шуме. Данные производятся на уровне 1 кГц. Чтобы предотвратить потерю разрешения времени, свойственного от критически произведенного пакета вейвлета, преобразовывают, используют неподкошенный, преобразовывают.

dt = 0.001;
t = 0:dt:1-dt;
x = ...
cos(2*pi*150*t).*(t>=0.2 & t<0.4)+sin(2*pi*200*t).*(t>0.6 & t<0.9);
y = x+0.05*randn(size(t));
[wpt,~,F] = modwpt(x,'TimeAlign',true);
contour(t,F.*(1/dt),abs(wpt).^2)
grid on
xlabel('Time (secs)')
ylabel('Hz')
title('Time-Frequency Analysis -- Undecimated Wavelet Packet Transform')

Обратите внимание на то, что пакет вейвлета преобразовывает, может разделить компоненты на 150 и 200 Гц. Это не верно для DWT, потому что падение на 150 и 200 Гц той же полосы октавы. Полосы октавы для уровня четыре DWT (в Гц)

Продемонстрируйте, что эти два компонента локализуются временем DWT, но не разделяются в частоте.

mra  = modwtmra(modwt(y,'fk18',4),'fk18');
J = 5:-1:1;
freq = 1/2*(1000./2.^J);
contour(t,freq,flipud(abs(mra).^2))
grid on
ylim([0 500])
xlabel('Time (secs)')
ylabel('Hz')
title('Time-Frequency Analysis -- Undecimated Wavelet Transform')

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

Энергетическое сохранение в пакете вейвлета преобразовывает

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

load kobe
plot(kobe)
grid on
xlabel('Seconds')
title('Kobe Earthquake Data')

Получите и подкошенные и неподкошенные пакетные преобразования вейвлета данных вниз к уровню 3. Чтобы гарантировать сопоставимые результаты для подкошенного пакета вейвлета преобразовывают, пример устанавливает dwtmode к periodization и возвращает его в вашу исходную установку в конце примера.

st = dwtmode('status','nodisplay');
dwtmode('per','nodisplay');
wptreeDecimated = wpdec(kobe,3,'fk18');
wptUndecimated = modwpt(kobe,3);

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

wpcfs = zeros(8,381);
for kk = 7:14
    wpcfs(kk-6,:) = wpcoef(wptreeDecimated,kk);
end

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

decimatedEnergy = sum(sum(abs(wpcfs).^2,2))

decimatedEnergy = 2.0189e+11

undecimatedEnergy = sum(sum(abs(wptUndecimated).^2,2))

undecimatedEnergy = 2.0189e+11

signalEnergy = norm(kobe,2)^2

signalEnergy = 2.0189e+11

Доли DWT это важное свойство с пакетом вейвлета преобразовывают.

wt = modwt(kobe,'fk18',3);
undecimatedWTEnergy = sum(sum(abs(wt).^2,2))

undecimatedWTEnergy = 2.0189e+11

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

Классификация пакетов вейвлета - землетрясение или взрыв?

Сейсмические записи поднимают действие со многих источников. Важно смочь классифицировать это действие на основе своего источника. Например, землетрясения и взрывы могут представить подобные подписи временного интервала, но очень важно дифференцироваться между этими двумя событиями. Набор данных EarthQuakeExplosionData содержит 16 записей с 2 048 выборками. Первые 8 записей (столбцы) являются сейсмическими записями землетрясений, последние 8 записей (столбцы) являются сейсмическими записями взрывов. Загрузите данные и постройте землетрясение и запись взрыва для сравнения. Данные взяты из пакета R astsa Стоффер (2014), который добавляет Шумвея и Стоффера (2011). Автор любезно позволил его использование в этом примере.

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

load EarthQuakeExplosionData
subplot(2,1,1)
plot(EarthQuakeExplosionData(:,3))
xlabel('Time')
title('Earthquake Recording')
grid on
subplot(2,1,2)
plot(EarthQuakeExplosionData(:,9))
xlabel('Time')
title('Explosion Recording')
grid on

Сформируйте пакетный характеристический вектор вейвлета путем разложения каждых временных рядов вниз, чтобы выровняться, три использования 'fk6' вейвлета с неподкошенным пакетом вейвлета преобразовывают. Это приводит к 8 поддиапазонам с аппроксимированной шириной 1/16 циклов/выборки. Используйте относительную энергию в каждом поддиапазоне, чтобы создать характеристический вектор. Как пример, получите пакетный энергетический характеристический вектор родственника вейвлета для первой записи землетрясения.

[wpt,~,F,E,RE] = modwpt(EarthQuakeExplosionData(:,1),3,'fk6');

RE содержит относительную энергию в каждом из этих 8 поддиапазонов. Если вы суммируете все элементы в RE, это равно 1. Обратите внимание на то, что это - значительное сокращение данных. Временные ряды длины 2048 уменьшаются до вектора с 8 элементами, каждым элементом, представляющим относительную энергию в пакетных узлах вейвлета на уровне 3. Функция помощника helperEarthQuakeExplosionClassifier получает относительные энергии для каждой из этих 16 записей на уровне 3 с помощью 'fk6' вейвлета. Это приводит к 16 8 матрица и использует эти характеристические векторы в качестве входных параметров к kmeans классификатору. Классификатор использует критерий статистической величины Разрыва, чтобы определить оптимальное количество кластеров для характеристических векторов между 1 и 6 и классифицирует каждый вектор. Кроме того, классификатор выполняет ту же самую классификацию на неподкошенном вейвлете, преобразовывают коэффициенты на уровне 3, полученном с 'fk6' вейвлетом и спектрами мощности для каждых из временных рядов. Неподкошенный вейвлет преобразовывает результаты в 16 4 матрица (3 поддиапазона вейвлета и 1 масштабирующийся поддиапазон). Спектры мощности приводят к 16 1025 матрица. У вас должны быть Statistics and Machine Learning Toolbox™ и Signal Processing Toolbox™, чтобы запустить классификатор.

Level = 3;
[WavPacketCluster,WtCluster,PspecCluster] = ...
helperEarthQuakeExplosionClassifier(EarthQuakeExplosionData,Level)

WavPacketCluster = 
  GapEvaluation with properties:

    NumObservations: 16
         InspectedK: [1 2 3 4 5 6]
    CriterionValues: [-0.0039 0.0779 0.1314 0.0862 0.0296 -0.0096]
           OptimalK: 2

WtCluster = 
  GapEvaluation with properties:

    NumObservations: 16
         InspectedK: [1 2 3 4 5 6]
    CriterionValues: [-0.0095 0.0474 0.0706 0.0260 -0.0280 -0.0346]
           OptimalK: 1

PspecCluster = 
  GapEvaluation with properties:

    NumObservations: 16
         InspectedK: [1 2 3 4 5 6]
    CriterionValues: [0.3804 0.3574 0.3466 0.2879 0.3256 0.3326]
           OptimalK: 1

dwtmode(st,'nodisplay')

WavPacketCluster кластеризация выход для неподкошенных пакетных характеристических векторов вейвлета. WtCluster кластеризация выход с помощью неподкошенных характеристических векторов DWT и PspecCluster выход для кластеризации на основе спектров мощности. Классификация пакетов вейвлета идентифицировала два кластера (OptimalK: 2) как оптимальный номер. Исследуйте результаты пакетной кластеризации вейвлета.

res = [ WavPacketCluster.OptimalY(1:8)' ; WavPacketCluster.OptimalY(9:end)']

res = 2×8

     1     2     2     2     2     2     2     2
     1     1     1     1     1     1     1     2

Вы видите, что были совершены только две ошибки. Из восьми записей землетрясения, семь классифицируются вместе (группа 2), и каждый по ошибке классифицируется как принадлежащий группе 1. То же самое верно для 8 записей взрыва - 7, классифицируются вместе, и 1 неправильно классифицируется. Классификация на основе уровня три неподкошенных DWT и спектры мощности возвращает один кластер как оптимальное решение.

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

Заключения

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

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

Ссылки

Wickerhauser, Т-х "Адаптировал анализ вейвлета от теории до программного обеспечения". Нажатие IEEE, 1994.

Персиваль, Д.Б. и Уолден, A.T. "Методы вейвлета для анализа временных рядов". Издательство Кембриджского университета, 2000.

Shumway, Р.Х. и Стоффер, D.H. "Анализ временных рядов и его приложения с примерами R". Спрингер, 2011.

Стоффер, D.H. "astsa: Прикладной Статистический Анализ Временных рядов". R версия пакета 1.3. https://CRAN.R-project.org/package=astsa, 2014.

Приложение

Следующая функция помощника используется в этом примере: