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

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

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

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

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

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

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

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

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

Насколько хорошо это полосовое приближение, зависит от того, насколько частотно локализованы фильтры. Для фильтров Фейера-Коровкина 'fk18' (18 коэффициентов), аппроксимация довольно хорошая. Для такого фильтра, как Haar ('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 записей с 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 перечислены в приложении. Для запуска классификатора необходимо иметь Toolbox™ статистики и машинного обучения и 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 уровня-3 и спектров мощности возвращает один кластер в качестве оптимального решения.

Если повторить вышесказанное с уровнем, равным 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