vmd

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

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

Синтаксис

imf = vmd(x)
[imf,residual] = vmd(x)
[imf,residual,info] = vmd(x)
[___] = vmd(x,Name,Value)
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(___) строит графики исходного сигнала, МВФ и остаточного сигнала как подграфиков на том же рисунке.

Примеры

свернуть все

Открыть 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. 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.

Открыть 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));

Вычислите МВФ сигнала с шумом и визуализируйте их на 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.

Открыть 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. 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

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

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
Открыть 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. The axes contains an object of type line.

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

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

[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.

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

свернуть все

Равномерно дискретизированный сигнал временной области, заданный как вектор или timetable. Если 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, f s/2], где f s - частота дискретизации.

Первоначальные МВФ, заданные как разделенная разделенными запятой парами, состоящая из '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 возвращается как:

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

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

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

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

  • CentralFrequencies - Центральные частоты МВФ.

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

  • AbsoluteImprovement - Среднее квадратное абсолютное улучшение к сходимости МВФ между двумя последними итерациями.

  • RelativeImprovement - Среднее относительное улучшение к сходимости МВФ между двумя последними итерациями.

  • LagrangeMultiplier - Частотный множитель Лагранжа при последней итерации.

Подробнее о

свернуть все

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

vmd функция разлагает сигнал <reservedrangesplaceholder3> (<reservedrangesplaceholder2>) на небольшое число K узкополосных intrinsic mode functions (IMFs):

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

МВФ обладает следующими характеристиками:

  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) путем наложения квадратичного штрафа и включения множителя Лагранжа. The 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', МВФ инициализируются в нуле. Инициализация '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) -й итерации алгоритм выполняет следующие шаги:

  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-го режима, вычисленное в (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)), her- частоты обновления множителя Лагранжа, заданная с помощью 'LMUpdateRate'.

Ссылки

[1] Бойд, Стивен, Нил Парих, Эрик Чу, Боря Пелеато и Джонатан Экштейн. Распределенная оптимизация и статистическое обучение через метод переменного направления умножителей. Основы и тренды® в машинном обучении. Том 3, Номер 1, 2011, стр. 1-122.

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

[3] Moody, George B., and Roger G. Mark. «The влияния of the MIT-BIH Arrhythmia Database». IEEE Engineering in Medicine and Biology Magazine. Том 20, № 3, май-июнь 2001 года, стр. 45-50.

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

.

См. также

|

Введенный в R2020a

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

Поддержка

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