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

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

Фигура 1: Дерево DWT до уровня 3 для 1-D входного сигнала.

$$\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$$ поддиапазоны. Поддиапазоны в hertz на уровне j приблизительно

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

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

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

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

Поскольку вейвлет делят ось частоты на более мелкие интервалы, чем 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 уровня 4 являются (в Гц)

Продемонстрировать, что эти два компонента локализованы по времени 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 записей с 2048 выборками. Первые 8 записей (столбцов) - сейсмические записи землетрясений, последние 8 записей (столбцов) - сейсмические записи взрывов. Загрузите данные и постройте график регистрации землетрясений и взрывов для сравнения. Данные взяты из пакета R astsa Stoffer [4], который дополняет Shumway и Stoffer [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. Классификатор использует статистический критерий Gap, чтобы определить оптимальное количество кластеров для векторов функций между 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 - выход кластеризации для неразрешенных векторов вейвлета packet функции. 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 классифицируется неправильно. Классификация, основанная на undecimated 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