Exponenta Banner
exponenta event banner

vmd

Разложение в вариационном режиме

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

Синтаксис

мвф = vmd (x)
[мвф, остаток] = vmd (x)
[МВФ, остаток, информация] = vmd (x)
[___] = vmd (x, имя, значение)
vmd (___)

Описание

пример

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

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

пример

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

пример

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

пример

vmd(___) готовит оригинальный сигнал, IMFs и остаточный сигнал как побочные сюжетные линии в том же числе.

Примеры

свернуть все

Открыть сценарий в реальном времени

Создайте сигнал с частотой 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. Axes 1 contains an object of type line. This object represents data. Axes 2 contains an object of type line. This object represents data. Axes 3 contains an object of type line. This object represents data. Axes 4 contains an object of type line. This object represents data. Axes 5 contains an object of type line. This object represents data.

Открыть сценарий в реальном времени

Генерация многокомпонентного сигнала, состоящего из трех синусоид частот 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));

Вычислите значения В шумного сигнала и визуализируйте их на 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. The axes contains 5 objects of type line.

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

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

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

Открыть сценарий в реальном времени

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

x (t) = 6t2 + cos (4.dt + 10.dt2) + {cos (60.dt), t≤0.5,cos (100.dt-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. Axes 1 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.

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

[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. Axes 1 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 2 objects of type line. These objects represent Reconstructed signal, Original signal.

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

norm(x-sig',Inf)
ans = 0
Открыть сценарий в реальном времени

Сигналы, отмеченные в этом примере, поступают из базы данных аритмии MIT-BIH [3] (панель инструментов обработки сигналов). Сигнал в базе данных отбирали с частотой 360 Гц.

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

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

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

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

Вычислите девять функций внутреннего режима оконного сигнала. Визуализируйте 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. Axes 1 with title IMF1 contains an object of type line. Axes 2 with title IMF2 contains an object of type line. Axes 3 with title IMF3 contains an object of type line. Axes 4 with title IMF4 contains an object of type line. Axes 5 with title IMF5 contains an object of type line. Axes 6 with title IMF6 contains an object of type line. Axes 7 with title IMF7 contains an object of type line. Axes 8 with title IMF8 contains an object of type line. Axes 9 with title IMF9 contains an object of type line.

Первый режим содержит наибольшее количество шумов, а второй режим колеблется на частоте пульса. Построение чистого ЭКГ-сигнала путем суммирования всех, кроме первого и последнего режимов 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. The axes contains 2 objects of type line. These objects represent Original ECG, Clean ECG.

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

свернуть все

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

Примечание

Если расписание имеет отсутствующие или повторяющиеся моменты времени, его можно исправить с помощью подсказок в «Чистом расписании» с «Отсутствующим», «Повторяющимся» или «Неуниформным временем».

Аргументы пары «имя-значение»

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

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

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

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

Примечание

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

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

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

«Количество извлекаемых файлов, указанных как пара, разделенная запятыми, состоящая из»NumIMF"и положительное скалярное целое число.

Начальные центральные частоты МВФ, указанные как разделенная запятыми пара, состоящая из 'CentralFrequencies«и вектор длины» NumIMFs. Значения вектора должны находиться в пределах [0, 0,5] циклов/выборок, что указывает, что истинные частоты находятся в пределах [0, fs/2], где fs - частота выборок.

Начальный 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 возвращается как:

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

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

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

residual возвращается как:

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

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

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

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

  • CentralFrequencies - центральные частоты ХХ.

  • NumIterations - общее число итераций.

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

  • RelativeImprovement - Среднее относительное улучшение к сходимости IMFs между заключительными двумя повторениями.

  • LagrangeMultiplier - умножитель Лагранжа в частотной области на последней итерации.

Подробнее

свернуть все

Функции внутреннего режима

vmd функция разлагает сигнал x (t) на небольшое число K узкополосных функций внутреннего режима (ВС):

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

Стативы имеют следующие характеристики:

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

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

    где β k (t) - фаза режима, а Ak (t) - его огибающая.

  2. Режимы имеют положительные и медленно изменяющиеся огибающие.

  3. Каждый режим имеет мгновенную частоту start' 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. Оцените полосу пропускания путем вычисления квадрата 2-нормы градиента демодулированного аналитического сигнала.

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

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

Алгоритмы

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

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

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

Для (n + 1) -й итерации алгоритм выполняет следующие шаги:

  1. Выполнить итерацию по режимам К сигнала (указанным с помощью '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-го режима, вычисленное в (n + 1) -й итерации.

    2. k-я центральная частота 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] Бойд, Стивен, Нил Парих, Эрик Чу, Борха Пелеато и Джонатан Экштейн. «Распределенная оптимизация и статистическое обучение с помощью метода чередующегося направления множителей». Основы и тенденции ® в машинном обучении. Том 3, номер 1, 2011, стр. 1-122.

[2] Драгомирецкий, Константин и Доминик Зоссо. «Декомпозиция в вариационном режиме». Транзакции IEEE ® при обработке сигналов. Том 62, номер 3, 2014, стр. 531-534.

[3] Муди, Джордж Б. и Роджер Г. Марк. «Влияние базы данных аритмии MIT-BIH». IEEE Engineering in Medicine and Biology Magazine. Том 20, № 3, май-июнь 2001 года, стр. 45-50.

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

.

См. также

|

Представлен в R2020a

Документация инструментария вейвлет

Поддержка