Пакеты вейвлета: разложение деталей

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

Это делает пакеты вейвлета привлекательной альтернативой DWT во многих приложениях. Двумя примерами, представленными здесь, является частотно-временной анализ и классификация сигнала. У вас должны быть Statistics and Machine Learning Toolbox™ и Signal Processing Toolbox™, чтобы выполнить классификацию.

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

Особое внимание в качестве примера на 1D случае, но многие концепции расширяет непосредственно к пакетному преобразованию вейвлета изображений.

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

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

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

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

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

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

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

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

[nFs2j+1,(n+1)Fs2j+1)n=0,1,2j-1 где Fs частота дискретизации.

То, насколько хороший это полосовое приближение, зависит от того, насколько локализованный частотой фильтры. Поскольку 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')

Figure contains an axes object. The axes object with title Time-Frequency Analysis -- Undecimated Wavelet Packet Transform contains an object of type contour.

Обратите внимание на то, что пакет вейвлета преобразовывает, может разделить компоненты на 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')

Figure contains an axes object. The axes object with title Time-Frequency Analysis -- Undecimated Wavelet Transform contains an object of type contour.

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

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

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

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

Figure contains an axes object. The axes object with title Kobe Earthquake Data contains an object of type line.

Получите и подкошенные и неподкошенные пакетные преобразования вейвлета данных вниз к уровню 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

Figure contains 2 axes objects. Axes object 1 with title Earthquake Recording contains an object of type line. Axes object 2 with title Explosion Recording contains an object of type line.

Сформируйте пакетный характеристический вектор вейвлета путем разложения каждых временных рядов вниз, чтобы выровнять три использования '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.

Смотрите также

| | |

Для просмотра документации необходимо авторизоваться на сайте