Практическое введение в мультирезолюционный анализ

Этот пример показывает, как выполнить и интерпретировать основное мультиразрешение сигнала (MRA). Пример использует как моделируемые, так и реальные данные, чтобы ответить на такие вопросы, как: Что означает мультирезолюционный анализ? Какие сведения о моем сигнале я могу получить, выполнив мультирезолюционный анализ? Каковы некоторые преимущества и недостатки различных методов MRA? Многие из анализов, представленных здесь, могут быть реплицированы с помощью приложения Signal Multiresolution Analyzer.

Что такое мультирезолюционный анализ?

Сигналы часто состоят из нескольких физически значимых компонентов. Довольно часто вы хотите изучить один или несколько из этих компонентов в изоляции на той же временной шкале, что и исходные данные. Мультирезолюционный анализ относится к разбиению сигнала на компоненты, которые формируют исходный сигнал именно при сложении обратно вместе. Чтобы быть полезным для анализа данных, важно, как сигнал разлагается. Компоненты идеально разлагают изменчивость данных на физически значимые и интерпретируемые части. Термин мультиразрешения анализ часто ассоциируется с вейвлетами или вейвлет-пакетами, но существуют невейвлет методов, которые также дают полезные MRA.

Как мотивирующий пример понимания, которое вы можете получить от MRA, рассмотрите следующий синтетический сигнал. Дискретизация сигнала производится на частоте 1000 Гц в течение одной секунды.

Fs = 1e3;
t = 0:1/Fs:1-1/Fs;
comp1 = cos(2*pi*200*t).*(t>0.7);
comp2 = cos(2*pi*60*t).*(t>=0.1 & t<0.3);
trend = sin(2*pi*1/2*t);
rng default
wgnNoise = 0.4*randn(size(t));
x = comp1+comp2+trend+wgnNoise;
plot(t,x)
xlabel('Seconds')
ylabel('Amplitude')
title('Synthetic Signal')

Figure contains an axes. The axes with title Synthetic Signal contains an object of type line.

Сигнал явно составлен из трех основных компонентов: локализованное во времени колебание с частотой 60 циклов/секунду, локализованное во времени колебание с частотой 200 циклов/секунду и термин тренда. Термин тренда здесь также синусоидален, но имеет частоту 1/2 цикла в секунду, поэтому он завершает только 1/2 цикла за одну секунду интервала. Колебание 60 оборотов в секунду или 60 Гц происходит между 0,1 и 0,3 секундами, в то время как колебание 200 Гц происходит между 0,7 и 1 секундой.

Не все это сразу видно из графика необработанных данных, потому что эти компоненты смешаны.

Теперь постройте график сигнала с частотной точки зрения.

xdft = fft(x);
N = numel(x);
xdft = xdft(1:numel(xdft)/2+1);
freq = 0:Fs/N:Fs/2;
plot(freq,20*log10(abs(xdft)))
xlabel('Cycles/second')
ylabel('dB')
grid on

Figure contains an axes. The axes contains an object of type line.

Из частотного анализа нам гораздо проще различить частоты колебательных компонентов, но мы потеряли их локализованный во времени характер. Также трудно визуализировать тренд в этом представлении.

Чтобы получить некоторую одновременную информацию о времени и частоте, мы можем использовать метод частотно-временного анализа, такой как непрерывное вейвлет (cwt).

cwt(x,Fs)

Figure contains an axes. The axes with title Magnitude Scalogram contains 3 objects of type image, line, area.

Теперь вы видите временные границы компонентов 60 Гц и 200 Гц. Однако у нас до сих пор нет никакой полезной визуализации тренда.

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

Мультирезолюционный анализ достигает этого. На самом деле, полезным способом думать о мультирезолюционном анализе является то, что он предоставляет способ избежать необходимости частотно-временного анализа, позволяя вам работать непосредственно во временном интервале.

Разделение компонентов сигнала во времени

Реальные сигналы являются смесью различных компонентов. Часто вас интересует только подмножество этих компонентов. Мультирезолюционный анализ позволяет вам сузить анализ путем разделения сигнала на компоненты при разных разрешениях.

Извлечение компонентов сигнала в различных разрешениях сводится к разложению изменений в данных по различным временным шкалам или эквивалентно в разных полосах (различные скорости колебаний). Соответственно, вы можете визуализировать изменчивость сигнала в разных шкалах или полосах одновременно.

Анализируйте и стройте график синтетического сигнала с помощью вейвлет. Сигнал анализируется в восьми разрешениях или уровнях.

mra = modwtmra(modwt(x,8));
helperMRAPlot(x,mra,t,'wavelet','Wavelet MRA',[2 3 4 9])

Figure contains 10 axes. Axes 1 with title Wavelet MRA contains an object of type line. Axes 2 contains an object of type line. Axes 3 contains an object of type line. Axes 4 contains an object of type line. Axes 5 contains an object of type line. Axes 6 contains an object of type line. Axes 7 contains an object of type line. Axes 8 contains an object of type line. Axes 9 contains an object of type line. Axes 10 contains an object of type line.

Не объясняя, что означают обозначения на графике, давайте используем наши знания о сигнале и попробуем понять, что показывает нам этот вейвлет MRA. Если вы начнете с самого верхнего графика и продолжите вниз, пока не достигнете графика исходных данных, то увидите, что компоненты становятся постепенно более гладкими. Если Вы предпочитаете думать о данных с точки зрения частоты, частоты, содержавшиеся в компонентах, становятся ниже. Напомним, что исходный сигнал имел три основных компонента, высокую частоту колебания на 200 Гц, более низкие частотные колебания на 60 Гц и тренд срок, все испорченные аддитивным шумом.

Если вы посмотрите на D2 График можно увидеть, что локализованный по времени высокая частота компонента изолирована там. Вы можете увидеть и исследовать эту важную функцию сигнала, по существу, изолированно. Следующие два графика содержат более низкие частотные колебания. Это важный аспект мультиразрешения, а именно важные компоненты сигнала могут не оказаться изолированными в одном компоненте MRA, но они редко расположены более чем в двух. Наконец, мы видим S8 график содержит термин trend. Для удобства, цвет осей в этих компонентах был изменен, чтобы подсветить их в MRA. Если вы предпочитаете визуализировать этот график или последующие таковые без подсветки, оставьте последний числовой вход в helperMRAPlot.

MRA вейвлет использует фиксированные функции, называемые вейвлетами, чтобы разделить компоненты сигнала. k-й вейвлета MRA компонента, обозначаемый Dk на предыдущем графике может рассматриваться как фильтрация сигнала в полосы частот вида [12k+1Δt,12kΔt] где Δt - период дискретизации или интервал дискретизации. Конечный гладкий компонент, обозначенная на графике SL, где L - уровень MRA, захватывает полосу [0,12L+1Δt]. Точность этого приближения зависит от вейвлета, используемого в MRA. Подробное описание вейвлет и вейвлет смотрите в разделе [5].

Однако существуют другие методы MRA.

Эмпирическое разложение моды (EMD) является адаптивным к данным методом мультиразрешения. EMD рекурсивно извлекает различные разрешения из данных без использования фиксированных функций или фильтров. EMD рассматривает сигнал как состоящий из быстрого колебания, наложенного на более медленный. После извлечения быстрого колебания процесс обрабатывает оставшийся более медленный компонент как новый сигнал и снова рассматривает ее как быстрое колебание, наложенное на более медленную. Процесс продолжается до достижения некоторого критерия остановки. В то время как EMD не использует фиксированные функции, подобные вейвлетам, для извлечения информации, EMD-подход концептуально очень похож на вейвлет-метод разделения сигнала на детали и приближения и затем снова разделения аппроксимации на детали и аппроксимацию. Компоненты MRA в EMD называются функциями внутреннего режима (МВФ). Подробную информацию о ЭМД см. в разделе [4].

Постройте график EMD-анализа того же сигнала.

[imf_emd,resid_emd] = emd(x);
helperMRAPlot(x,imf_emd,t,'emd','Empirical Mode Decomposition',[1 2 3 6])

Figure contains 7 axes. Axes 1 with title Empirical Mode Decomposition contains an object of type line. Axes 2 contains an object of type line. Axes 3 contains an object of type line. Axes 4 contains an object of type line. Axes 5 contains an object of type line. Axes 6 contains an object of type line. Axes 7 contains an object of type line.

Хотя количество компонентов MRA отличается, EMD и вейвлет MRA создают аналогичное изображение сигнала. Это не случайно. Смотрите [2] для описания сходства между вейвлет и EMD.

При разложении EMD высокочастотные колебания локализуются до первой функции внутреннего режима (МВФ 1). Более низкие колебания частоты локализованы в основном на МВФ 2, но вы можете увидеть некоторый эффект также в МВФ 3. Компонент тренда в МВФ 6 очень похож на компонент тренда, извлеченный вейвлетом метода.

Еще одним методом для адаптивного мультирезолюционного анализа является вариационное разложение режима (VMD). Как и EMD, VMD пытается извлечь функции внутреннего режима или режимы колебаний из сигнала, не используя фиксированные функции для анализа. Но EMD и VMD определяют режимы очень по-разному. EMD работает рекурсивно на сигнале временной области, чтобы извлечь постепенно более низкую частоту МВФ. VMD запускается путем идентификации peaks сигнала в частотный диапазон и извлекает все режимы одновременно. См. [1] для лечения VMD.

[imf_vmd,resid_vmd] = vmd(x);
helperMRAPlot(x,imf_vmd,t,'vmd','Variational Mode Decomposition',[2 4 5])

Figure contains 6 axes. Axes 1 with title Variational Mode Decomposition contains an object of type line. Axes 2 contains an object of type line. Axes 3 contains an object of type line. Axes 4 contains an object of type line. Axes 5 contains an object of type line. Axes 6 contains an object of type line.

Ключевым моментом является то, что, подобно разложениям вейвлет и EMD, VMD разделяет три интересующих компонента на совершенно отдельные режимы или на небольшое количество смежных режимов. Все три техники позволяют вам визуализировать компоненты сигнала в то же время, шкала и исходный сигнал.

Существует адаптивный к данным метод, который фактически создает вейвлет на основе частотного содержимого данных. Этот метод является эмпирическим вейвлет (EWT) [3]. Одним из основных критических замечаний EMD является то, что его определение является чисто алгоритмическим. В результате он не легко поддается математическому анализу. EWT, с другой стороны, фактически создает вейвлеты Мейера на основе частотного содержимого анализируемого сигнала. Соответственно, результаты EWT поддаются математическому анализу, потому что фильтры, используемые в анализе, доступны пользователю. Повторите анализ синтетического сигнала с помощью EWT.

[mra_ewt,cfs,wfb,info] = ewt(x,'MaxNumPeaks',5);
helperMRAPlot(x,mra_ewt,t,'EWT','Empirical Wavelet Transform',[1 2 3 5])

Figure contains 6 axes. Axes 1 with title Empirical Wavelet Transform contains an object of type line. Axes 2 contains an object of type line. Axes 3 contains an object of type line. Axes 4 contains an object of type line. Axes 5 contains an object of type line. Axes 6 contains an object of type line.

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

Nf = length(x)/2+1;
cf = mean(info.FilterBank.Passbands,2);
f = 0:Fs/length(x):Fs-Fs/length(x);
clf
ax = newplot;
plot(ax,f(1:Nf),wfb(1:Nf,:))
ax.XTick = sort(cf.*Fs);
ax.XTickLabel = ax.XTick;
xlabel('Hz')
grid on
title('Empirical Meyer Wavelets')

Figure contains an axes. The axes with title Empirical Meyer Wavelets contains 5 objects of type line.

Другим преимуществом метода EWT является то, что коэффициенты анализа сохраняют энергию исходного сигнала. Это свойство разделяется неадаптивными методами вейвлета, но не найдено в невейвлете адаптивных методах.

EWTEnergy = sum(vecnorm(cfs,2).^2)
EWTEnergy = 875.5768
sigEnergy = norm(x,2)^2
sigEnergy = 875.5768

Для реального примера полезного разделения компонентов рассмотрим сейсмограф (вертикальное ускорение, nm/s2) землетрясения в Кобе, зафиксированного в Университете Тасмании, Хобарт, Австралия, 16 января 1995 года, начиная с 20:56:51 (GMT) и продолжаясь в течение 51 минут с 1 секундными интервалами.

load KobeTimeTable
T = KobeTimeTable.t;
kobe = KobeTimeTable.kobe;
figure
plot(T,kobe)
title('Kobe Earthquake Seismograph')
ylabel('Vertical Acceleration (nm/s^2)')
xlabel('Time')
axis tight
grid on

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

Получите и постройте график вейвлет данных.

mraKobe = modwtmra(modwt(kobe,8));
figure
helperMRAPlot(kobe,mraKobe,T,'Wavelet','Wavelet MRA Kobe Earthquake',[4 5])

Figure contains 10 axes. Axes 1 with title Wavelet MRA Kobe Earthquake contains an object of type line. Axes 2 contains an object of type line. Axes 3 contains an object of type line. Axes 4 contains an object of type line. Axes 5 contains an object of type line. Axes 6 contains an object of type line. Axes 7 contains an object of type line. Axes 8 contains an object of type line. Axes 9 contains an object of type line. Axes 10 contains an object of type line.

График показывает разделение первичных и отложенных компонентов вторичной волны в компонентах MRA D4 и D5. Компоненты в сейсмической волне перемещаются с различными скоростями, первичные (сжимающие) волны перемещаются быстрее, чем вторичные (сдвигающие) волны. Методы MRA могут позволить вам исследовать эти компоненты изолированно на исходной временной шкале.

Восстановление сигналов от MRA

Точка разделения сигналов на компоненты часто состоит в том, чтобы удалить определенные компоненты или уменьшить их эффект на сигнал. Решающее значение для методов MRA имеет способность восстанавливать исходный сигнал.

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

sigrec_wavelet = sum(mra);
sigrec_emd = sum(imf_emd,2)+resid_emd;
sigrec_vmd = sum(imf_vmd,2)+resid_vmd;
figure
subplot(3,1,1)
plot(t,sigrec_wavelet); title('Wavelet reconstruction');
set(gca,'XTickLabel',[]);
ylabel('Amplitude');
subplot(3,1,2);
plot(t,sigrec_emd); title('EMD reconstruction');
set(gca,'XTickLabel',[]);
ylabel('Amplitude');
subplot(3,1,3)
plot(t,sigrec_vmd); title('VMD reconstruction');
ylabel('Amplitude');
xlabel('Time');

Figure contains 3 axes. Axes 1 with title Wavelet reconstruction contains an object of type line. Axes 2 with title EMD reconstruction contains an object of type line. Axes 3 with title VMD reconstruction contains an object of type line.

Максимальная ошибка реконструкции на основе выборки по выборкам для каждого из методов имеет порядок 10-12 или меньше, что указывает на то, что они являются совершенными методами реконструкции. Поскольку сумма компонентов MRA восстанавливает исходный сигнал, это означает, что включение или исключение подмножества компонентов может привести к полезному приближению.

Вернитесь к нашему исходному вейвлет синтетического сигнала и предположим, что вы не интересовались терминами тренда. Поскольку термин тренда локализован в последнем компоненте MRA, просто исключить этот компонент из реконструкции.

sigWOtrend = sum(mra(1:end-1,:));
figure
plot(t,sigWOtrend)
xlabel('Time (s)')
ylabel('Amplitude')
title('Trend Term Removed')

Figure contains an axes. The axes with title Trend Term Removed contains an object of type line.

Чтобы удалить другие компоненты, можно создать логический вектор с false значения в компонентах, которые вы не хотите включать. Здесь мы удаляем тренд и самую высокую частотную составляющую вместе с первым компонентом MRA (который выглядит в значительной степени шумом). Постройте график фактического второго компонента сигнала (60 Гц) вместе с реконструкцией для сравнения.

include = true(size(mra,1),1);
include([1 2 9]) = false;
ts = sum(mra(include,:));
plot(t,comp2,'b')
hold on
plot(t,ts,'r')
title('Trend Term and Highest Frequency Component Removed')
xlabel('Time (s)')
ylabel('Amplitude')
legend('Component 2','Partial Reconstruction')
xlim([0.0 0.4])

Figure contains an axes. The axes with title Trend Term and Highest Frequency Component Removed contains 2 objects of type line. These objects represent Component 2, Partial Reconstruction.

В предыдущем примере мы рассматривали термин тренда как неприятный компонент, который будет удален. Существует ряд приложений, в которых тренд может быть основным компонентом интереса. Давайте визуализируем условия тренда, извлеченные нашими тремя примерами MRA.

figure
plot(t,trend,'LineWidth',2)
hold on
plot(t,[mra(end,:)' imf_vmd(:,end) imf_emd(:,end) mra_ewt(:,end)])
grid on
legend('Trend','Wavelet','VMD','EMD','EWT')
ylabel('Amplitude')
xlabel('Time (s)')
title('Trend in three MRAs')

Figure contains an axes. The axes with title Trend in three MRAs contains 5 objects of type line. These objects represent Trend, Wavelet, VMD, EMD, EWT.

Обратите внимание, что тренд более плавный и наиболее точно захвачен вейвлетом метода. EMD находит плавный термин тренда, но он смещен относительно истинной амплитуды тренда, в то время как метод VMD кажется по своей сути более смещенным к нахождению колебаний, чем методы вейвлет и EMD. Последствия этого далее обсуждаются в разделе MRA Technologies - Преимущества и недостатки.

Обнаружение переходных изменений с помощью MRA

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

Чтобы проиллюстрировать это, рассмотрим ежеквартальные взвешенные по цепочке данные о реальном валовом внутреннем продукте (ВВП) США для 1947Q1- 2011Q4. Ежеквартальные выборки соответствуют частоте дискретизации 4 пробы в год. Вертикальная черная линия знаменует начало «Великой модерации», означающей период пониженной макроэкономической волатильности в США, начиная с середины 1980-х годов. Обратите внимание, что это трудно различить, рассматривая необработанные данные.

load GDPData
figure
plot(year,realgdp)
hold on
plot([year(146) year(146)],[-0.06 0.14],'k')
title('GDP Data')
xlabel('Year')
hold off

Figure contains an axes. The axes with title GDP Data contains 2 objects of type line.

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

mra = modwtmra(modwt(realgdp,'db2'));
figure
plot(year,mra(1,:))
hold on
plot([year(146) year(146)],[-0.015 0.015],'k')
title('Wavelet MRA - Quarter-to-Quarter Changes')
xlabel('Year')
hold off

Figure contains an axes. The axes with title Wavelet MRA - Quarter-to-Quarter Changes contains 2 objects of type line.

Уменьшение изменчивости или экономической волатильности гораздо более очевидно в компоненте MRA с самым высоким разрешением, чем в необработанных данных. Методы для обнаружения изменений отклонения, таких как MATLAB findchangepts (Signal Processing Toolbox) часто работает лучше с компонентами MRA, чем с необработанными данными.

Методов MRA - Преимущества и недостатки

В этом примере мы обсудили вейвлет и адаптивные к данным методы для мультирезолюционного анализа. Каковы некоторые преимущества и недостатки различных методов? Другими словами, для каких приложений я могу выбрать один над другим? Начнем с вейвлетов. Вейвлеты методов в этом примере используют фиксированные фильтры для получения MRA. Это означает, что вейвлет MRA имеет четко определенное математическое объяснение, и мы можем предсказать поведение MRA. Мы также можем связать события в MRA с конкретными временными шкалами в данных, как это было сделано в примере ВВП. Недостатком является то, что преобразование вейвлета делит сигнал на октавные полосы (уменьшение центральной частоты на 1/2 в каждом компоненте), так что при высоких центральных частотах полосы пропускания намного больше, чем при более низких частотах. Это означает, что два близко расположенных высоких частот колебания могут легко оказаться в одном и том же компоненте MRA с вейвлетом метода.

Повторите первый синтетический пример, но перемещайте два колебательных компонента в пределах одной октавы друг от друга.

Fs = 1e3;
t = 0:1/Fs:1-1/Fs;
comp1 = cos(2*pi*150*t).*(t>=0.1 & t<0.3);
comp2 = cos(2*pi*200*t).*(t>0.7);
trend = sin(2*pi*1/2*t);
rng default;
wgnNoise = 0.4*randn(size(t));
x = comp1+comp2+trend+wgnNoise;
plot(t,x)
xlabel('Seconds')
ylabel('Amplitude')

Figure contains an axes. The axes contains an object of type line.

Повторите и постройте график вейвлет.

mra = modwtmra(modwt(x,8));
helperMRAPlot(x,mra,t,'wavelet','Wavelet MRA',[2 9])

Figure contains 10 axes. Axes 1 with title Wavelet MRA contains an object of type line. Axes 2 contains an object of type line. Axes 3 contains an object of type line. Axes 4 contains an object of type line. Axes 5 contains an object of type line. Axes 6 contains an object of type line. Axes 7 contains an object of type line. Axes 8 contains an object of type line. Axes 9 contains an object of type line. Axes 10 contains an object of type line.

Теперь мы видим, что D2 содержит как компоненты 150 Гц, так и компоненты 200 Гц. Если вы повторяете этот анализ с помощью EMD, вы видите тот же результат.

Теперь давайте используем VMD.

[imf_vmd,~,info_vmd] = vmd(x);
helperMRAPlot(x,imf_vmd,t,'vmd','VMD',[1 2 3 5]);

Figure contains 6 axes. Axes 1 with title VMD contains an object of type line. Axes 2 contains an object of type line. Axes 3 contains an object of type line. Axes 4 contains an object of type line. Axes 5 contains an object of type line. Axes 6 contains an object of type line.

VMD может разделить эти два компонента. Колебания высокой частоты локализованы на МВФ 1, в то время как второй компонент распределен по двум соседним МВФ. Если вы посмотрите на предполагаемые центральные частоты режимов VMD, метод локализовал первые два компонента около 200 и 150 Гц. Третий МВФ имеет центральную частоту, близкую к 150 Гц, из-за чего мы видим второй компонент в двух компонентах MRA.

info_vmd.CentralFrequencies*Fs
ans = 5×1

  202.7204
  153.3278
  148.8022
   84.2802
    0.2667

VMD способен сделать это, потому что он начинается с определения центральных частот кандидатов для МВФ, рассматривая анализ данных в частотном диапазоне.

Хотя вейвлет MRA не может разделить две высокочастотные компоненты, существует дополнительный вейвлет-пакет MRA, который может.

wpt = modwptdetails(x,3);
helperMRAPlot(x,flip(wpt),t,'wavelet','Wavelet Packet MRA',[5 6 8]);

Figure contains 9 axes. Axes 1 with title Wavelet Packet MRA contains an object of type line. Axes 2 contains an object of type line. Axes 3 contains an object of type line. Axes 4 contains an object of type line. Axes 5 contains an object of type line. Axes 6 contains an object of type line. Axes 7 contains an object of type line. Axes 8 contains an object of type line. Axes 9 contains an object of type line.

Теперь вы видите, что два колебания разделены в D5 и D6. Из этого примера мы можем извлечь общее правило. Если начальное разложение вейвлета или EMD показывает компоненты с, по-видимому, различными скоростями колебаний в том же компоненте, рассмотрите VMD или MRA вейвлет-пакета. Если вы подозреваете, что ваши данные имеют высокую частоту компоненты, близкие друг к другу по частоте, VMD или подход с вейвлет-пакетом обычно работают лучше, чем вейвлет или EMD подход.

Напомним, проблема извлечения плавного тренда. Повторите как вейвлет MRA, так и EMD.

mra = modwtmra(modwt(x,8));
helperMRAPlot(x,mra,t,'wavelet','Wavelet MRA',[2 9])

Figure contains 10 axes. Axes 1 with title Wavelet MRA contains an object of type line. Axes 2 contains an object of type line. Axes 3 contains an object of type line. Axes 4 contains an object of type line. Axes 5 contains an object of type line. Axes 6 contains an object of type line. Axes 7 contains an object of type line. Axes 8 contains an object of type line. Axes 9 contains an object of type line. Axes 10 contains an object of type line.

imf_emd = emd(x);
helperMRAPlot(x,imf_emd,t,'EMD','Empirical Mode Decomposition',[1 2 6])

Figure contains 7 axes. Axes 1 with title Empirical Mode Decomposition contains an object of type line. Axes 2 contains an object of type line. Axes 3 contains an object of type line. Axes 4 contains an object of type line. Axes 5 contains an object of type line. Axes 6 contains an object of type line. Axes 7 contains an object of type line.

Тренды, извлеченные вейвлетом и методами EMD, ближе к истинному тренду, чем извлеченные с помощью VMD и вейвлета метода пакета. VMD по своей сути смещен в сторону нахождения узкополосных колебательных компонентов. Это - сила VMD в обнаружении тесно расположенных колебаний, но недостаток при извлечении плавных трендов в данных. То же самое относится и к MRA вейвлет, когда разложение выходит за пределы нескольких уровней. Это приводит ко второй общей рекомендации. Если вы заинтересованы в характеристике плавного тренда в ваших данных для идентификации или удаления, попробуйте вейвлет или метод EMD.

Как насчет обнаружения переходных изменений, как мы видели в данных по ВВП? Повторим анализ ВВП с помощью VMD.

[imf_vmd,~,vmd_info] = vmd(realgdp);
figure
subplot(2,1,1)
plot(year,realgdp)
title('Real GDP');
hold on
plot([year(146) year(146)],[-0.1 0.15],'k')
hold off
subplot(2,1,2)
plot(year,imf_vmd(:,1));
title('First VMD IMF');
xlabel('Year')
hold on
plot([year(146) year(146)],[-0.02 0.02],'k')
hold off

Figure contains 2 axes. Axes 1 with title Real GDP contains 2 objects of type line. Axes 2 with title First VMD IMF contains 2 objects of type line.

Хотя компонент VMD с самой высокой частотой также, по-видимому, показывает некоторое снижение изменчивости, начиная с середины 1980-х годов, это не так легко очевидно, как в вейвлет. Поскольку метод VMD смещен в сторону нахождения колебаний, в первом МВФ VMD существует существенный звон, который заслоняет изменения волатильности.

Повторите этот анализ с помощью EMD.

imf_emd = emd(realgdp);
figure
subplot(2,1,1)
plot(year,realgdp)
title('Real GDP');
hold on
plot([year(146) year(146)],[-0.1 0.15],'k')
hold off
subplot(2,1,2)
plot(year,imf_emd(:,1));
title('First EMD IMF')
xlabel('Year')
hold on
plot([year(146) year(146)],[-0.06 0.05],'k')
hold off

Figure contains 2 axes. Axes 1 with title Real GDP contains 2 objects of type line. Axes 2 with title First EMD IMF contains 2 objects of type line.

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

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

Заключения

Этот пример показал, как методы мультирезолюционного разложения, такие как вейвлет, вейвлет-пакет, эмпирическое разложение моды, эмпирическое вейвлет и вариационное разложение режима позволяют вам изучать компоненты сигнала в относительной изоляции в той же временной шкале, что и исходные данные. Каждый метод зарекомендовала себя мощной в ряде применений. Пример дал несколько правил большого пальца, чтобы начать, но они не должны рассматриваться как абсолютные. В следующей таблице приведены свойства методов MRA, представленных здесь, вместе с некоторыми общими правилами большого пальца. Две порции обозначают конкретную прочность, одна плюс указывает, что метод применим, но не является особой прочностью. Для двоичного свойства, такого как сохранение энергии в анализе, флажок указывает, что метод имеет это свойство, и «x» указывает, что свойство отсутствует.

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

  • Если данные, по-видимому, содержат тесно расположенные высокочастотные колебательные компоненты, попробуйте VMD или вейвлет. VMD определяет важные центральные частоты непосредственно от данных, в то время как вейвлет используют анализ фиксированной частоты, который может быть менее гибким, чем VMD. Вейвлет и эмпирические методы являются постоянными Q-фильтрами сигнала, что означает, что полосы пропускания пропорциональны центральной частоте. На высоких частотах это затрудняет разделение тесно расположенных компонентов.

  • Если вы заинтересованы в переходных событиях в ваших данных, таких как импульсивные события или переходные сокращения и увеличения изменчивости, попробуйте вейвлет или эмпирический вейвлет MRA. В любом MRA эти события обычно локализуются до компонентов MRA с самой высокой шкалой (самая высокая центральная частота).

  • Если вы заинтересованы в представлении термина плавного тренда в данных, рассмотрите EMD или вейвлет MRA.

С помощью Signal Multiresolution Analyzer можно выполнить мультиразрешение анализ по сигналу, получить метрики на различных компонентах MRA, экспериментировать с частичными реконструкциями и сгенерировать скрипты MATLAB, чтобы воспроизвести анализ в командной строке.

Ссылки

[1] Драгомирецкий, Константин, Доминик Зоссо. «Разложение вариационного режима». Транзакции IEEE по обработке сигналов 62, № 3 (февраль 2014): 531-44. https://doi.org/10.1109/TSP.2013.2288675.

[2] Flandrin, P., G. Rilling и P. Goncalves. Эмпирическое разложение моды как банк фильтров. IEEE Signal Processing Letters 11, No. 2 (February 2004): 112-14. https://doi.org/10.1109/LSP.2003.821662.

[3] Giles, J. «Empirical wavelet Transform», IEEE Transactions on Signal Processing, vol. 61, no. 16 (May 2013): 3999 - 4010.

[4] Huang, Norden E., Zheng Shen, Steven R. Long, Manli C. Wu, Hsing H. Shih, Quanan Zheng, Nai-Chyuan Yen, Chi Chao Tung, и Henry H. Liu. Эмпирическое разложение моды и спектр Гильберта для нелинейного и нестационарного анализа временных рядов. Материалы Лондонского королевского общества. Серия A: Математические, физические и инженерные науки 454, № 1971 (8 марта 1998): 903-95. https://doi.org/10.1098/rspa.1998.0193.

[5] Персиваль, Дональд Б. и Эндрю Т. Уолден. Вейвлет для анализа временных рядов. Кембриджская серия в статистической и вероятностной математике. Кембридж; Нью-Йорк: Cambridge University Press, 2000.

См. также

Функции

Приложения

Похожие темы