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

Следующий рисунок показывает дерево 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$$ где Фс является частотой дискретизации.

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

Продемонстрируйте, что эти два компонента локализуются временем 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: 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

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