emd

Эмпирическое разложение моды

Описание

пример

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

пример

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

пример

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

пример

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

Примеры

свернуть все

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

load('sinusoidalSignalExampleData.mat','X','fs')
t = (0:length(X)-1)/fs;
plot(t,X)
xlabel('Time(s)')

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

Чтобы создать Гильбертов график спектра, вам нужны внутренние функции режима (IMFs) сигнала. Выполните эмпирическое разложение моды, чтобы вычислить IMFs и остаточные значения сигнала. Поскольку сигнал не является гладким, задайте 'pchip'как метод интерполяции.

[imf,residual,info] = emd(X,'Interpolation','pchip');
Current IMF  |  #Sift Iter  |  Relative Tol  |  Stop Criterion Hit  
      1      |        2     |     0.026352   |  SiftMaxRelativeTolerance
      2      |        2     |    0.0039573   |  SiftMaxRelativeTolerance
      3      |        1     |     0.024838   |  SiftMaxRelativeTolerance
      4      |        2     |      0.05929   |  SiftMaxRelativeTolerance
      5      |        2     |      0.11317   |  SiftMaxRelativeTolerance
      6      |        2     |      0.12599   |  SiftMaxRelativeTolerance
      7      |        2     |      0.13802   |  SiftMaxRelativeTolerance
      8      |        3     |      0.15937   |  SiftMaxRelativeTolerance
      9      |        2     |      0.15923   |  SiftMaxRelativeTolerance
The decomposition stopped because the number of extrema of the residual signal is less than 'MaxNumExtrema'.

Таблица, сгенерированная в командном окне, показывает, что количество отсеивает итерации, относительный допуск, и отсеять критерий остановки каждого сгенерировал МВФ. Эта информация также содержится в info. Можно скрыть таблицу путем добавления 'Display',0 пара значение-имя.

Создайте Гильбертов график спектра с помощью imf компоненты, полученные с помощью эмпирического разложения моды.

hht(imf,fs)

Частота по сравнению с графиком временной зависимости является разреженным графиком с вертикальной цветной полосой, указывающей на мгновенную энергию в каждой точке в МВФ. График представляет мгновенный спектр частоты каждого компонента, анализируемого от исходного смешанного сигнала. Три IMFs появляются в графике с отличным изменением в частоте в 1 секунду.

Эта тригонометрическая идентичность представляет два других взгляда на тот же физический сигнал:

52потому что 2πf1t+14(потому что 2π(f1+f2)t+потому что 2π(f1-f2)t)=(2+потому что2 πf2 t)потому что 2πf1 t.

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

t = 0:1e-3:10;
omega1 = 2*pi*100;
omega2 = 2*pi*20;
s = 0.25*cos((omega1-omega2)*t) + 2.5*cos(omega1*t) + 0.25*cos((omega1+omega2)*t);
z = (2+cos(omega2/2*t).^2).*cos(omega1*t);

norm(s-z,Inf) 
ans = 3.2729e-13

Постройте синусоиды и выберите 1 второй интервал, запускающийся в 2 секунды.

plot(t,[s' z'])
xlim([2 3])
xlabel('Time (s)')
ylabel('Signal')

Получите спектрограмму сигнала. Спектрограмма показывает три отличных синусоидальных компонента. Анализ Фурье рассматривает сигналы как суперпозицию синусоид.

pspectrum(s,1000,'spectrogram','TimeResolution',4)

Используйте emd вычислить внутренние функции режима (IMFs) и дополнительной диагностической информации сигнала. Функция выходными параметрами по умолчанию таблица, которая указывает на количество отсеивания итераций, относительного допуска и критерия остановки отсеивания каждого МВФ. Эмпирическое разложение моды рассматривает сигнал как z.

[imf,~,info] = emd(s);
Current IMF  |  #Sift Iter  |  Relative Tol  |  Stop Criterion Hit  
      1      |        1     |   1.8025e-06   |  SiftMaxRelativeTolerance
The decomposition stopped because the current energy ratio is greater than 'MaxEnergyRatio'.

Количество нулевых пересечений и локальных экстремальных значений отличается самое большее один. Это удовлетворяет необходимому условию для сигнала быть МВФ.

info.NumZerocrossing - info.NumExtrema
ans = 1

Постройте МВФ и выберите 0,5 вторых интервала, запускающиеся в 2 секунды. МВФ является сигналом AM потому что emd просматривает сигнал как модулируемую амплитуду.

plot(t,imf)
xlim([2 2.5])
xlabel('Time (s)')
ylabel('IMF')

Симулируйте сигнал вибрации от поврежденного подшипника. Выполните эмпирическое разложение моды, чтобы визуализировать IMFs сигнала и искать дефекты.

Терпение диаметра подачи 12 см имеет восемь прокручивающихся элементов. Каждый элемент прокрутки имеет диаметр 2 см. Внешняя гонка остается стационарной, когда внутренняя гонка управляется в 25 циклах в секунду. Акселерометр производит колебания подшипника на уровне 10 кГц.

fs = 10000;
f0 = 25;
n = 8;
d = 0.02;
p = 0.12;

Сигнал вибрации от здорового подшипника включает несколько порядков ведущей частоты.

t = 0:1/fs:10-1/fs;
yHealthy = [1 0.5 0.2 0.1 0.05]*sin(2*pi*f0*[1 2 3 4 5]'.*t)/5;

Резонанс взволнован в вибрации подшипника на полпути посредством процесса измерения.

yHealthy = (1+1./(1+linspace(-10,10,length(yHealthy)).^4)).*yHealthy;

Резонанс вводит дефект во внешней гонке подшипника, который приводит к прогрессивному износу. Дефект вызывает ряд ударов, которые повторяются на частоте передачи шара внешняя гонка (BPFO) подшипника:

BPFO=12nf0[1-dpпотому чтоθ],

где f0 пропускная способность, n количество прокручивающихся элементов, d диаметр прокручивающихся элементов, p диаметр подачи подшипника, и θ угол контакта подшипника. Примите угол контакта 15 ° и вычислите BPFO.

ca = 15;
bpfo = n*f0/2*(1-d/p*cosd(ca));

Используйте pulstran функционируйте, чтобы смоделировать удары, когда периодическое обучается синусоид с 5 миллисекундами. Каждая синусоида на 3 кГц является оконной окном с плоской вершиной. Используйте закон о степени, чтобы ввести прогрессивный износ в сигнале вибрации подшипника.

fImpact = 3000;
tImpact = 0:1/fs:5e-3-1/fs;
wImpact = flattopwin(length(tImpact))'/10;
xImpact = sin(2*pi*fImpact*tImpact).*wImpact;

tx = 0:1/bpfo:t(end);
tx = [tx; 1.3.^tx-2];

nWear = 49000;
nSamples = 100000;
yImpact = pulstran(t,tx',xImpact,fs)/5;
yImpact = [zeros(1,nWear) yImpact(1,(nWear+1):nSamples)];

Сгенерируйте сигнал вибрации BPFO путем добавления ударов на здоровый сигнал. Постройте сигнал и выберите 0,3 вторых интервала, запускающиеся в 5,0 секунд.

yBPFO = yImpact + yHealthy;

xLimLeft = 5.0;
xLimRight = 5.3;
yMin = -0.6;
yMax = 0.6;

plot(t,yBPFO)

hold on
[limLeft,limRight] = meshgrid([xLimLeft xLimRight],[yMin yMax]);
plot(limLeft,limRight,'--')
hold off

Увеличьте масштаб выбранного интервала, чтобы визуализировать эффект ударов.

xlim([xLimLeft xLimRight])

Добавьте белый Гауссов шум в сигналы. Задайте шумовое отклонение 1/1502.

rn = 150;
yGood = yHealthy + randn(size(yHealthy))/rn;
yBad = yBPFO + randn(size(yHealthy))/rn;

plot(t,yGood,t,yBad)
xlim([xLimLeft xLimRight])
legend('Healthy','Damaged')

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

imfGood = emd(yGood,'MaxNumIMF',5);
Current IMF  |  #Sift Iter  |  Relative Tol  |  Stop Criterion Hit  
      1      |        3     |     0.017132   |  SiftMaxRelativeTolerance
      2      |        3     |      0.12694   |  SiftMaxRelativeTolerance
      3      |        6     |      0.14582   |  SiftMaxRelativeTolerance
      4      |        1     |     0.011082   |  SiftMaxRelativeTolerance
      5      |        2     |      0.03463   |  SiftMaxRelativeTolerance
The decomposition stopped because 'MaxNumIMF' was reached.

Используйте emd без выходных аргументов, чтобы визуализировать первые три режима и невязку. Установите 'Display' к 0 скрыть таблицу.

emd(yGood,'MaxNumIMF',5,'Display',0)

Вычислите и визуализируйте IMFs дефектного сигнала подшипника. Первый эмпирический режим показывает высокочастотные удары. Этот высокочастотный режим увеличения энергии как износ прогрессирует. Третий режим показывает резонанс в сигнале вибрации.

imfBad = emd(yBad,'MaxNumIMF',5);
Current IMF  |  #Sift Iter  |  Relative Tol  |  Stop Criterion Hit  
      1      |        2     |     0.041274   |  SiftMaxRelativeTolerance
      2      |        3     |      0.16695   |  SiftMaxRelativeTolerance
      3      |        3     |      0.18428   |  SiftMaxRelativeTolerance
      4      |        1     |     0.037177   |  SiftMaxRelativeTolerance
      5      |        2     |     0.095861   |  SiftMaxRelativeTolerance
The decomposition stopped because 'MaxNumIMF' was reached.
emd(yBad,'MaxNumIMF',5,'Display',0)

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

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

load('sinusoidalSignalExampleData.mat','X','fs')
t = (0:length(X)-1)/fs;
plot(t,X)
xlabel('Time(s)')

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

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

emd(X,'Interpolation','pchip')
Current IMF  |  #Sift Iter  |  Relative Tol  |  Stop Criterion Hit  
      1      |        2     |     0.026352   |  SiftMaxRelativeTolerance
      2      |        2     |    0.0039573   |  SiftMaxRelativeTolerance
      3      |        1     |     0.024838   |  SiftMaxRelativeTolerance
      4      |        2     |      0.05929   |  SiftMaxRelativeTolerance
      5      |        2     |      0.11317   |  SiftMaxRelativeTolerance
      6      |        2     |      0.12599   |  SiftMaxRelativeTolerance
      7      |        2     |      0.13802   |  SiftMaxRelativeTolerance
      8      |        3     |      0.15937   |  SiftMaxRelativeTolerance
      9      |        2     |      0.15923   |  SiftMaxRelativeTolerance
The decomposition stopped because the number of extrema of the residual signal is less than 'MaxNumExtrema'.

Сгенерирован интерактивный график с исходным сигналом, первыми 3 IMFs и невязкой. Таблица, сгенерированная в командном окне, показывает, что количество отсеивает итерации, относительный допуск, и отсеять критерий остановки каждого сгенерировал МВФ. Можно скрыть таблицу путем определения Display как 0.

Щелкните правой кнопкой по пробелу в графике открыть окно селектора МВФ. Используйте селектор МВФ, чтобы выборочно просмотреть сгенерированный IMFs, исходный сигнал и невязку.

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

Выбранные IMFs теперь отображены на графике.

Используйте график визуализировать отдельные компоненты, анализируемые из исходного сигнала наряду с невязкой. Обратите внимание на то, что невязка вычисляется для общего количества IMFs и не изменяется на основе IMFs, выбранного в окне селектора МВФ.

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

свернуть все

Однородно произведенный сигнал временной области, заданный или как векторное или как одно расписание столбца данных.

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

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

Пример: 'MaxNumIMF',5

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

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

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

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

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

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

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

Сигнал к остаточному энергетическому отношению, заданному как разделенная запятой пара, состоящая из 'MaxEnergyRatio'и скаляр. MaxEnergyRatio отношение энергии сигнала в начале отсеивания и средней энергии конверта. MaxEnergyRatio один из критериев остановки разложения, то есть, остановки разложения, когда текущее энергетическое отношение больше, чем MaxEnergyRatio.

Метод интерполяции для конструкции конверта, заданной как разделенная запятой пара, состоящая из 'Interpolation'и любой 'spline' или 'pchip'.

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

  • 'spline', если X сглаженный сигнал

  • 'pchip', если X несглаженный сигнал

сплайн'метод интерполяции использует кубический сплайн, в то время как 'pchip' использует кусочный кубический метод интерполяционного полинома Эрмита.

Переключите отображение информации в командном окне, заданном как разделенная запятой пара, состоящая из 'Display'и или 1 или 0. Таблица, сгенерированная в командном окне, показывает, что количество отсеивает итерации, относительный допуск, и отсеять критерий остановки каждого сгенерировал МВФ. Задайте Display как 1, чтобы показать таблицу или 0, чтобы скрыть таблицу.

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

свернуть все

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

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

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

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

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

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

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

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

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

  • NumIMF - Количество IMFs извлечено из сигнала

  • NumExtrema - Количество экстремального значения в каждом МВФ

  • NumZeroCrossing - Количество нулевых пересечений для каждого МВФ

  • NumSifting - Количество отсеиваний выполняется для каждого МВФ

  • MeanEnvelopeEnergy - Энергия среднего значения верхнего и более низкого конверта получена в каждом расчете МВФ

  • RelativeTolerance - Относительный допуск достигнут для каждого МВФ

Алгоритмы

Empirical Mode Decomposition

emd разлагает X сигнала (t) на количество k внутренних функций режима (IMF) и остаточный rk (t) с помощью процесса отсеивания. Краткий обзор процесса отсеивания, перечисленного в [1] и [2], следующие:

  1. Найдите, что локальные максимумы и минимумы для X(t) сигнала создают верхний конверт s + (t) и более низкий конверт s (t).

  2. Вычислите средний конверт для i th итерация, mk,i (t),

    mk,i(t)=12[s+(t)+s(t)]

  3. С ck (t) = X (t) для первой итерации, вычтите средний конверт из остаточного сигнала,

    ck(t)=ck(t)mk,i(t)

    Если ck (t) не совпадает с критериями МВФ, шаги 4 и 5 пропущены. Процедура выполнена с помощью итераций снова на шаге 1 с новым значением ck (t).

  4. Если ck(t) совпадает с критериями МВФ, новая невязка вычисляется. Чтобы обновить остаточный сигнал, вычтите k th МВФ от предыдущего остаточного сигнала,

    rk(t)=rk1(t)ck(t)

  5. Затем начните с шага 1, с помощью невязки, полученной как новый rk сигнала (t), и сохраните ck (t) как внутренняя функция режима.

Для функций режима внутреннего параметра N исходный сигнал представлен как,

X(t)=i=1Nci(t)+rN(t)

Для получения дополнительной информации о процессе отсеивания, см. [1] и [2].

SiftRelativeTolerance

SiftRelativeTolerance критерий остановки типа Коши, предложенный в [4]. Отсеивание остановок, когда текущий относительный допуск меньше SiftRelativeTolerance. Текущий относительный допуск задан как

RelativeTolerancec(t)предыдущийc(t)текущий2c(t)текущий2.

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

MaxEnergyRatio

Энергетическое отношение является отношением энергии сигнала в начале отсеивания и средней энергии конверта [3]. Разложение останавливается, когда текущее энергетическое отношение больше, чем MaxEnergyRatio. Для k IMFs, EnergyRatio задан как

EnergyRatio10журнал10(X(t)2rk(t)2).

Ссылки

[1] Хуан, Норден Э., Чжен Шен, Стивен Р. Лонг, Манли К. Ву, Син Х. Ши, Куэнэн Чжен, Най-Чюань Янь, Ши Чао Туньг и Генри Х. Лю. "Эмпирическое разложение моды и Гильбертов спектр для нелинейного и неустановившегося анализа временных рядов". Продолжения Королевского общества Лондона. Серии А: Математические, Физические и Технические науки. Издание 454, 1998, стр 903–995.

[2] Вытекание струей, G., Патрик Фландрен и Паулу Гонзальвес. "На эмпирическом разложении моды и его алгоритмах". Семинар IEEE-EURASIP по Нелинейной Обработке сигналов и Обработке изображений 2003. NSIP-03. Градо, Италия. 8–11.

[3] Rato, R.T., Мануэль Ортигуейра и Арнальдо Батиста. "На HHT, его проблемах и некоторых решениях". Механические Системы и Издание 22, 2008 Обработки сигналов, стр 1374–1394.

[4] Ван, Банда, Сиань-Yao Чен, Литий Клыка Цяо, Жаохуа Ву и Норден Хуан. "На Внутренней Функции Режима". Усовершенствования в Адаптивном Анализе данных. Издание 2, Номер 3, 2010, стр 277–293.

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

Генерация кода C/C++
Генерация кода C и C++ с помощью MATLAB® Coder™.

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

Введенный в R2018a