vmd

Вариационное разложение режима

свернуть все на странице

Синтаксис

imf = vmd(x)
[imf,residual] = vmd(x)
[imf,residual,info] = vmd(x)
[___] = vmd(x,Name,Value)
vmd(___)

Описание

пример

imf = vmd(x) возвращает вариационное разложение режима xИспользование vmd анализировать и упростить сложные сигналы в конечное число внутренних функций режима (IMFs), требуемый выполнять Гильбертов спектральный анализ.

[imf,residual] = vmd(x) также возвращает остаточный residual сигнала соответствие вариационному разложению режима x.

пример

[imf,residual,info] = vmd(x) возвращает дополнительную информацию info на IMFs и остаточном сигнале в диагностических целях.

пример

[___] = vmd(x,Name,Value) выполняет вариационное разложение режима с дополнительными опциями, заданными одним или несколькими Name,Value парные аргументы.

пример

vmd(___) строит исходный сигнал, IMFs и остаточный сигнал как подграфики на том же рисунке.

Примеры

свернуть все

Скрипт Open Live Script

Создайте сигнал, произведенный на уровне 4 кГц, который напоминает набор номера всех ключей цифрового телефона. Сохраните сигнал как расписание MATLAB®.

fs = 4e3;
t = 0:1/fs:0.5-1/fs;

ver = [697 770 852 941];
hor = [1209 1336 1477];

tones = [];

for k = 1:length(ver)
    for l = 1:length(hor)
        tone = sum(sin(2*pi*[ver(k);hor(l)].*t))';
        tones = [tones;tone;zeros(size(tone))];
    end
end

% To hear, type soundsc(tones,fs)

S = timetable(tones,'SampleRate',fs);

Постройте вариационное разложение режима расписания.

vmd(S)

Figure contains 5 axes objects. Axes object 1 contains an object of type line. This object represents data. Axes object 2 contains an object of type line. This object represents data. Axes object 3 contains an object of type line. This object represents data. Axes object 4 contains an object of type line. This object represents data. Axes object 5 contains an object of type line. This object represents data.

Скрипт Open Live Script

Сгенерируйте многокомпонентный сигнал, состоящий из трех синусоид частот 2 Гц, 10 Гц и 30 Гц. Синусоиды производятся на уровне 1 кГц в течение 2 секунд. Встройте сигнал в белый Гауссов шум отклонения 0,01 ².

fs = 1e3;
t = 1:1/fs:2-1/fs;
x = cos(2*pi*2*t) + 2*cos(2*pi*10*t) + 4*cos(2*pi*30*t) + 0.01*randn(1,length(t));

Вычислите IMFs сигнала с шумом и визуализируйте их в 3-D графике.

imf = vmd(x);
[p,q] = ndgrid(t,1:size(imf,2));
plot3(p,q,imf)
grid on
xlabel('Time Values')
ylabel('Mode Number')
zlabel('Mode Amplitude')

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

Используйте вычисленный IMFs, чтобы построить Гильбертов спектр многокомпонентного сигнала. Ограничьте частотный диапазон [0, 40] Гц.

hht(imf,fs,'FrequencyLimits',[0,40])

Figure contains an axes object. The axes object with title Hilbert Spectrum contains 5 objects of type patch.

Скрипт Open Live Script

Сгенерируйте кусочный составной сигнал, состоящий из квадратичного тренда, щебета и косинуса с резким переходом между двумя постоянными частотами в t = 0.5.

x(t)=6t2 +cos(4πt+10πt2)+{cos(60πt),t≤0.5,cos(100πt-10π),t>0.5.

Сигнал производится на уровне 1 кГц в течение 1 секунды. Постройте каждый отдельный компонент и составной сигнал. 

fs = 1e3;
t = 0:1/fs:1-1/fs;

x = 6*t.^2 + cos(4*pi*t+10*pi*t.^2) + ...
    [cos(60*pi*(t(t<=0.5))) cos(100*pi*(t(t>0.5)-10*pi))];

tiledlayout('flow')
nexttile
plot(t,[zeros(1,length(t)/2) cos(100*pi*(t(length(t)/2+1:end))-10*pi)])
xlabel('Time (s)')
ylabel('Cosine')

nexttile
plot(t,[cos(60*pi*(t(1:length(t)/2))) zeros(1,length(t)/2)])
xlabel('Time (s)')
ylabel('Cosine')

nexttile
plot(t,cos(4*pi*t+10*pi*t.^2))
xlabel('Time (s)')
ylabel('Chirp')


nexttile
plot(t,6*t.^2)
xlabel('Time (s)')
ylabel('Quadratic trend')

nexttile(5,[1 2])
plot(t,x)
xlabel('Time (s)')
ylabel('Signal')

Figure contains 5 axes objects. Axes object 1 contains an object of type line. Axes object 2 contains an object of type line. Axes object 3 contains an object of type line. Axes object 4 contains an object of type line. Axes object 5 contains an object of type line.

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

[imf,res] = vmd(x,'NumIMFs',4);

tiledlayout('flow')

for i = 1:4
    nexttile
    plot(t,imf(:,i))
    txt = ['IMF',num2str(i)];
    ylabel(txt)
    xlabel('Time (s)')
    grid on
end

Восстановите сигнал путем добавления функций режима и невязки. Постройте и сравните исходные и восстановленные сигналы.

sig = sum(imf,2) + res;

nexttile(5,[1 2])
plot(t,sig,'LineWidth',1.5)

hold on

plot(t,x,':','LineWidth',2)
xlabel('Time (s)')
ylabel('Signal')
hold off
legend('Reconstructed signal','Original signal', ...
       'Location','northwest')

Figure contains 5 axes objects. Axes object 1 contains an object of type line. Axes object 2 contains an object of type line. Axes object 3 contains an object of type line. Axes object 4 contains an object of type line. Axes object 5 contains 2 objects of type line. These objects represent Reconstructed signal, Original signal.

Вычислите норму различия между оригиналом и восстановленными сигналами.

norm(x-sig',Inf)
ans = 0
Скрипт Open Live Script

Сигналы, помеченные в этом примере, от Базы данных Аритмии MIT-BIH [3] (Signal Processing Toolbox). Сигнал в базе данных был произведен на уровне 360 Гц.

Загрузите сигналы базы данных MIT, соответствующие, чтобы записать 200 и построить сигнал.

load mit200
Fs = 360;
plot(tm,ecgsig)
ylabel('Time (s)')
xlabel('Signal')

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

Сигнал ECG содержит скачки, управляемые ритмом heartbeat и колеблющегося низкочастотного шаблона. Отличные спицы ECG создают важные гармоники высшего порядка.

Вычислите девять внутренних функций режима оконного сигнала. Визуализируйте IMFs. 

[imf,residual] = vmd(ecgsig,'NumIMF',9);
t = tiledlayout(3,3,'TileSpacing','compact','Padding','compact');
for n = 1:9
    ax(n) = nexttile(t);
    plot(tm,imf(:,n)')
    xlim([tm(1) tm(end)])
    txt = ['IMF',num2str(n)];
    title(txt)
    xlabel('Time (s)')
end
title(t,'Variational Mode Decomposition')

Figure contains 9 axes objects. Axes object 1 with title IMF1 contains an object of type line. Axes object 2 with title IMF2 contains an object of type line. Axes object 3 with title IMF3 contains an object of type line. Axes object 4 with title IMF4 contains an object of type line. Axes object 5 with title IMF5 contains an object of type line. Axes object 6 with title IMF6 contains an object of type line. Axes object 7 with title IMF7 contains an object of type line. Axes object 8 with title IMF8 contains an object of type line. Axes object 9 with title IMF9 contains an object of type line.

Первый режим содержит самое шумовое, и второй режим колеблется на частоте heartbeat. Создайте чистый сигнал ECG путем подведения итогов всех кроме первых и последних режимов VMD, таким образом отбрасывания низкочастотного базового колебания и большей части высокочастотного шума.

cleanECG = sum(imf(:,2:8),2);
figure
plot(tm,ecgsig,tm,cleanECG)
legend('Original ECG','Clean ECG')
ylabel('Time (s)')
xlabel('Signal')

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent Original ECG, Clean ECG.

Входные параметры

свернуть все

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

Примечание

Если расписание имеет пропавших без вести или дублирующиеся моменты времени, можно зафиксировать его с помощью советов в Чистом Расписании с Пропавшими без вести, Копией, или Неоднородные Времена.

Аргументы name-value

Задайте дополнительные разделенные запятой пары Name,Value аргументы. Name имя аргумента и Value соответствующее значение. Name должен появиться в кавычках. Вы можете задать несколько аргументов в виде пар имен и значений в любом порядке, например: Name1, Value1, ..., NameN, ValueN.

Пример: 'NumIMF',10

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

Относительная погрешность сходимости режима в виде разделенной запятой пары, состоящей из 'RelativeTolerance'и положительный действительный скаляр. RelativeTolerance один из критерия остановки для оптимизации, то есть, остановки оптимизации, когда среднее относительное улучшение к сходимости IMFs, в двух последовательных итерациях, меньше RelativeTolerance.

Примечание

Процесс оптимизации останавливается когда 'AbsoluteTolerance' и 'RelativeTolerance' совместно удовлетворены.

Максимальное количество итераций оптимизации в виде разделенной запятой пары, состоящей из 'MaxIterations'и положительное скалярное целое число. MaxIterations один из критерия остановки для оптимизации, то есть, остановки оптимизации, когда количество итераций больше MaxIterations.

MaxIterations может быть задан с помощью только положительные целые числа.

Количество IMFs, извлеченного в виде разделенной запятой пары, состоящей из 'NumIMF'и положительное скалярное целое число.

Начальные центральные частоты МВФ в виде разделенной запятой пары, состоящей из 'CentralFrequencies'и вектор из длины NumIMFs. Векторные значения должны быть в [0, 0.5] циклы/выборка, который указывает, что истинные частоты в [0, f s/2], где f s является частотой дискретизации.

Начальный IMFs в виде разделенной запятой пары, состоящей из 'InitialIMFs'и действительная матрица. Строки соответствуют выборкам времени, и столбцы соответствуют режимам.

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

Начальный множитель Лагранжа в виде разделенной запятой пары, состоящей из 'InitialLM'и комплексный вектор. Область значений начального множителя Лагранжа в частотном диапазоне является [0, 0.5] циклами/выборкой. Множитель осуществляет ограничение реконструкции. Длина множителя зависит от входного размера.

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

Метод, чтобы инициализировать центральные частоты в виде разделенной запятой пары, состоящей из 'InitializeMethod'и любой 'peaks', 'random', или 'grid'.

Задайте InitializeMethod как:

  • 'peaks' инициализировать центральные частоты как пиковые местоположения сигнала в частотном диапазоне (значение по умолчанию).

  • 'random' инициализировать центральные частоты как случайные числа, распределенные однородно в интервале [0,0.5] циклы/выборка.

  • 'grid' инициализировать центральные частоты как однородно произведенную сетку в интервале [0,0.5] циклы/выборка.

Переключите отображение прогресса в командном окне в виде разделенной запятой пары, состоящей из 'Display'и любой 'true' (или 1) или 'false' (или 0). Если вы задаете 'true', функция отображает среднее абсолютное и относительное улучшение режимов и центральных частот каждые 20 итераций, и покажите итоговую информацию об остановке.

Задайте Display как 1, чтобы показать таблицу или 0, чтобы скрыть таблицу.

Выходные аргументы

свернуть все

Внутренние функции режима, возвращенные как матрица или расписание. Каждый imf амплитуда и частота модулируемый сигнал с положительными и медленно различными конвертами. Каждый режим имеет мгновенную частоту, которая не уменьшается, медленно варьируется и сконцентрирована вокруг центрального значения. Используйте imf применять преобразование Гильберта-Хуанга, чтобы выполнить спектральный анализ на сигнале.

imf возвращен как:

  • Матрица A, где каждым столбцом является imf, когда x вектор

  • Расписание, когда x одно расписание столбца данных

Остаточный сигнал, возвращенный как вектор-столбец или одно расписание столбца данных. residual представляет фрагмент исходного x сигнала не анализируемый vmd.

residual возвращен как:

  • Вектор-столбец, когда x вектор.

  • Одно расписание столбца данных, когда x одно расписание столбца данных.

Дополнительная информация для диагностики, возвращенной как структура с этими полями:

  • ExitFlag – Флаг Termination. Значение 0 указывает на алгоритм, остановленный, когда он достиг максимального количества итераций. Значение 1 указывает на алгоритм, остановленный, когда он соответствовал абсолютным и относительным погрешностям.

  • CentralFrequencies – Центральные частоты IMFs.

  • NumIterations – Общее количество итераций.

  • AbsoluteImprovement – Среднее значение придало абсолютному улучшению квадратную форму к сходимости IMFs между итоговыми двумя итерациями.

  • RelativeImprovement – Среднее относительное улучшение к сходимости IMFs между итоговыми двумя итерациями.

  • LagrangeMultiplier – Множитель Лагранжа частотного диапазона в последней итерации.

Больше о

свернуть все

Внутренние функции режима

vmd функция разлагает x сигнала (t) на небольшое число K узкополосного intrinsic mode functions (IMFs):

x(t)=∑k=1Kuk(t).

IMFs имеют эти характеристики:

  1. Каждый режим uk является амплитудой, и частота модулировала сигнал формы

    uk(t)=Ak(t)cos(ϕk(t)),

    где ϕk (t) является фазой режима, и Ak (t) является своим конвертом.

  2. Режимы имеют положительные и медленно различные конверты.

  3. Каждый режим имеет мгновенную частоту ϕ'k (t), который не уменьшается, медленно варьируется и сконцентрирован вокруг центрального значения fk.

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

Оптимизация

Чтобы вычислить uk и fk, процедура находит оптимум увеличенной функции Лагранжа

L(uk(t),fk,λ(t))=α∑k=1K‖ddt[(δ(t)+jπt)∗uk(t)]e−j2πfkt‖22+‖x(t)−∑k=1Kuk(t) ‖22+〈λ(t),x(t)−∑k=1Kuk(t)〉(i)(ii)(iii),

где скалярное произведение 〈p(t),q(t)〉=∫−∞∞p∗(t) q(t) dt и 2-норма ‖p(t)‖22=〈p(t),p(t)〉. Термин регуляризации (i) включает эти шаги:

  1. Используйте преобразование Гильберта, чтобы вычислить аналитический сигнал, сопоставленный с каждым режимом, где * обозначает свертку. Это приводит к каждому режиму, имеющему чисто положительный спектр.

  2. Демодулируйте аналитический сигнал к основной полосе путем умножения его с комплексной экпонентой.

  3. Оцените полосу пропускания путем вычисления квадратичной нормы градиента демодулируемого аналитического сигнала.

Условия (ii) и (iii) осуществляют ограничение x(t)=∑k=1Kuk(t) путем наложения квадратичного штрафа и слияния множителя Лагранжа. PenaltyFactor α измеряет относительную важность (i) по сравнению с (ii) и (iii).

Алгоритм решает задачу оптимизации с помощью переменного метода направления множителей, описанных в [1] (Signal Processing Toolbox).

Алгоритмы

vmd функция вычисляет IMFs в частотном диапазоне, восстанавливая X (f) = ДПФ {x (t)} в терминах Uk (f) = ДПФ {uk (t)}. Чтобы удалить краевые эффекты, алгоритм расширяет сигнал путем зеркального отражения половины его длины с обеих сторон.

Множитель Лагранжа, введенный в Оптимизации (Signal Processing Toolbox), имеет преобразование Фурье Ʌ (f). Длина вектора множителей Лагранжа является длиной расширенного сигнала.

Если в противном случае не задано в 'InitialIMFs', IMFs инициализируются в нуле. Инициализируйте 'CentralFrequencies' использование одного из методов задано в 'InitializeMethod'. vmd итеративно обновляет режимы, пока одно из этих условий не соблюдают:

  • ∑k‖ukn+1(t)−ukn(t)‖22/‖ukn(t)‖22<εr и ∑k‖ukn+1(t)−ukn(t)‖22<εa совместно удовлетворены, где εr и εa заданы с помощью 'RelativeTolerance' и 'AbsoluteTolerance', соответственно.

  • Алгоритм превышает максимальное количество итераций, заданных в 'MaxIterations'.

Для (n + 1)-th итерация, алгоритм выполняет эти шаги:

  1. Выполните итерации по режимам K сигнала (заданное использование 'NumIMF') вычислить:

    1. Формы волны частотного диапазона для каждого использования режима

      Ukn+1(f)=X(f)−∑i<kUkn+1(f)−∑i>kUkn(f)+Λn2(f)1+2α{2π(f−fkn)}2,

      где Ukn+1(f) преобразование Фурье k th режим, вычисленный в (n + 1)-th итерация.

    2. k th центральная частота fkn+1 использование

      fkn+1=∫0∞|Ukn+1(f)|2 f df∫0∞|Ukn+1(f)|2 df≈∑f|Ukn+1(f)|2 ∑|Ukn+1(f)|2 .

  2. Обновите использование множителя Лагранжа Λn+1(f)=Λn(f)+τ(X(f)−∑kUkn+1(f)), где τ является частотой обновления множителя Лагранжа, заданное использование 'LMUpdateRate'.

Ссылки

[1] Бойд, Стивен, Нил Пэрих, Эрик Чу, Боря Пелеато и Йонатан Экштайн. "Распределенная Оптимизация и Статистическое Изучение с помощью Переменного Метода Направления Множителей". Основы и Тренды® в Машинном обучении. Vol 3, Номер 1, 2011, стр 1–122.

[2] Dragomiretskiy, Константин и Доминик Зоссо. "Вариационное разложение режима". IEEE® Транзакции на Обработке сигналов. Издание 62, Номер 3, 2014, стр 531–534.

[3] Капризный, Джордж Б. и Роджер Г. Марк. "Удар Базы данных Аритмии MIT-BIH". Разработка IEEE в Журнале Медицины и Биологии. Издание 20, № 3, мочь-июнь 2001, стр 45–50.

Расширенные возможности

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

|

Введенный в R2020a

Документация Wavelet Toolbox

Поддержка

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