В этом примере показано, как пакеты вейвлета отличаются от дискретного вейвлета преобразовывает (DWT). Пример показывает, как пакет вейвлета преобразовывает результаты в фильтрацию поддиапазона равной ширины сигналов в противоположность более грубой фильтрации полосы октавы, найденной в DWT.
Это делает пакеты вейвлета привлекательной альтернативой DWT во многих приложениях. Двумя примерами, представленными здесь, является частотно-временной анализ и классификация сигнала. У вас должны быть Statistics and Machine Learning Toolbox™ и Signal Processing Toolbox™, чтобы выполнить классификацию.
Если вы используете ортогональный вейвлет с пакетом вейвлета, преобразовывают, вы дополнительно заканчиваете с разделением энергии сигнала среди поддиапазонов равной ширины.
Особое внимание в качестве примера на 1D случае, но многие концепции расширяет непосредственно к пакетному преобразованию вейвлета изображений.
Следующий рисунок показывает дерево DWT для 1D входного сигнала.
Рисунок 1: Дерево DWT вниз к уровню 3 для 1D входного сигнала.
масштабирующийся (lowpass) аналитический фильтр и представляет вейвлет (highpass) аналитический фильтр. Метки в нижней части показывают раздел оси частоты [0,1/2] в поддиапазоны.
Рисунок показывает, что последующие уровни DWT управляют только на выходных параметрах lowpass (масштабирование) фильтром. На каждом уровне DWT делит сигнал на полосы октавы. В критически произведенном DWT выходные параметры полосовых фильтров прорежены два на каждом уровне. В неподкошенном дискретном вейвлете преобразовывают, выходные параметры не прорежены.
Сравните дерево DWT с полным пакетным деревом вейвлета.
Рисунок 2: Полное пакетное дерево вейвлета вниз к уровню 3.
В вейвлете пакет преобразовывает, операции фильтрации также применяются к вейвлету, или детали, коэффициентам. Результат состоит в том, что пакеты вейвлета обеспечивают фильтрацию поддиапазона входного сигнала на прогрессивно более прекрасные интервалы равной ширины. На каждом уровне, , ось частоты [0,1/2] разделена на поддиапазоны. Поддиапазоны в герц на уровне приблизительно
где частота дискретизации.
То, насколько хороший это полосовое приближение, зависит от того, насколько локализованный частотой фильтры. Поскольку Fejér-Korovkin фильтрует как 'fk18'
(18 коэффициентов), приближение довольно хорошо. Для фильтра как Хаар ('haar'
), приближение не точно.
В критически произведенном вейвлете пакет преобразовывает выходные параметры полосовых фильтров, прорежены два. В неподкошенном вейвлете пакет преобразовывает, выходные параметры не прорежены.
Поскольку пакеты вейвлета делят ось частоты на более прекрасные интервалы, чем 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 (в Гц)
[0,31.25)
[31.25,62.5)
[62.5,125)
[125,250)
[250,500)
Продемонстрируйте, что эти два компонента локализуются временем 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') axis tight
Получите и подкошенные и неподкошенные пакетные преобразования вейвлета данных вниз к уровню 3. Чтобы гарантировать сопоставимые результаты для подкошенного пакета вейвлета преобразовывают, устанавливают граничный дополнительный режим на 'periodic'
.
wptreeDecimated = dwpt(kobe,'Level',3,'Boundary','periodic'); wptUndecimated = modwpt(kobe,3);
Вычислите полную энергию и на подкошенном и на неподкошенном пакетном уровне вейвлета три коэффициента и сравните с энергией в исходном сигнале.
decimatedEnergy = sum(cell2mat(cellfun(@(x) sum(abs(x).^2),wptreeDecimated,'UniformOutput',false)))
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 Стоффер [4], который добавляет Шумвея и Стоффера [3]. Автор любезно позволил его использование в этом примере.
Постройте временные ряды от каждой группы для сравнения.
load EarthQuakeExplosionData subplot(2,1,1) plot(EarthQuakeExplosionData(:,3)) xlabel('Time') title('Earthquake Recording') grid on axis tight subplot(2,1,2) plot(EarthQuakeExplosionData(:,9)) xlabel('Time') title('Explosion Recording') grid on axis tight
Сформируйте пакетный характеристический вектор вейвлета путем разложения каждых временных рядов вниз, чтобы выровнять три использования '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 матрица. Исходный код для helperEarthQuakeExplosionClassifier
перечислен в приложении. У вас должны быть 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: 2
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
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.
helperEarthQuakeExplosionClassifier
function [WavPacketCluster,WtCluster,PspecCluster] = helperEarthQuakeExplosionClassifier(Data,Level) % This function is provided solely in support of the featured example % waveletpacketsdemo.mlx. % This function may be changed or removed in a future release. % Data - The data matrix with individual time series as column vectors. % Level - The level of the wavelet packet and wavelet transforms NP = 2^Level; NW = Level+1; WpMatrix = zeros(16,NP); WtMatrix = zeros(16,NW); for jj = 1:8 [~,~,~,~,RE] = modwpt(Data(:,jj),Level,'fk6'); wt = modwt(Data(:,jj),Level,'fk6'); WtMatrix(jj,:) = sum(abs(wt).^2,2)./norm(Data(:,jj),2)^2; WpMatrix(jj,:) = RE; end for kk = 9:16 [~,~,~,~,RE] = modwpt(Data(:,kk),Level,'fk6'); wt = modwt(Data(:,kk),Level,'fk6'); WtMatrix(kk,:) = sum(abs(wt).^2,2)./norm(Data(:,kk),2)^2; WpMatrix(kk,:) = RE; end % For repeatability rng('default'); WavPacketCluster = evalclusters(WpMatrix,'kmeans','gap','KList',1:6,... 'Distance','cityblock'); WtCluster = evalclusters(WtMatrix,'kmeans','gap','KList',1:6,... 'Distance','cityblock'); Pxx = periodogram(Data,hamming(numel(Data(:,1))),[],1,'power'); PspecCluster = evalclusters(Pxx','kmeans','gap','KList',1:6,... 'Distance','cityblock'); end
[1] Wickerhauser, Младен Виктор. Адаптированный анализ вейвлета от теории до программного обеспечения. Веллесли, MA: А.К. Питерс, 1994.
[2] Персиваль, Дональд Б. и Эндрю Т. Уолден. Методы вейвлета для анализа временных рядов. Кембриджский ряд в статистической и вероятностной математике. Кембридж ; Нью-Йорк: Издательство Кембриджского университета, 2000.
[3] Shumway, Роберт Х. и Дэвид С. Стоффер. Анализ Временных рядов и Его Приложения: С Примерами R. 3-и тексты редактора Спрингера в Статистике. Нью-Йорк: Спрингер, 2011.
[4] Стоффер, D. H. "asta: Прикладной Статистический Анализ Временных рядов". R версия пакета 1.3. https://CRAN.R-project.org/package=astsa
, 2014.