hht

Преобразование Гильберта-Хуанга

Описание

пример

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

пример

hs = hht(imf,fs) возвращает спектр Гильберта hs сигнала, дискретизированного со скоростью fs.

пример

[hs,f,t] = hht(___) возвращает вектор частоты f и временной вектор t в дополнение к hs. Эти выходные аргументы могут использоваться с любым из предыдущих входных синтаксисов.

пример

[hs,f,t,imfinsf,imfinse] = hht(___) также возвращает мгновенные частоты imfinsf и мгновенные энергии imfinse функций внутреннего режима для диагностики сигналов.

пример

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

пример

hht(___) без выходных аргументов строит графики спектра Гильберта в текущую фигуру окне. Можно использовать этот синтаксис с любым из входных параметров в предыдущих синтаксисах.

hht(___,freqlocation) строит графики спектра Гильберта с помощью необязательных freqlocation аргумент для определения местоположения оси частоты. Частота представлена на оси y по умолчанию.

Примеры

свернуть все

Сгенерируйте Гауссов-модулированный квадратичный щебет. Задайте частоту дискретизации 2 кГц и длительность сигнала 2 секунды.

fs = 2000;
t = 0:1/fs:2-1/fs;
q = chirp(t-2,4,1/2,6,'quadratic',100,'convex').*exp(-4*(t-1).^2);
plot(t,q)

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

Использование emd визуализация функций внутреннего режима (МВФ) и невязки.

emd(q)

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.

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

imf = emd(q,'Display',1);
Current IMF  |  #Sift Iter  |  Relative Tol  |  Stop Criterion Hit  
      1      |        2     |    0.0063952   |  SiftMaxRelativeTolerance
      2      |        2     |       0.1007   |  SiftMaxRelativeTolerance
      3      |        2     |      0.01189   |  SiftMaxRelativeTolerance
      4      |        2     |    0.0075124   |  SiftMaxRelativeTolerance
Decomposition stopped because the number of extrema in the residual signal is less than the 'MaxNumExtrema' value.

Используйте вычисленные МВФ, чтобы построить график спектра квадратичного щебета. Ограничьте частотную область значений от 0 Гц до 20 Гц.

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

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

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

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

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

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

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

[imf,residual,info] = emd(X,'Interpolation','pchip');

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

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

hht(imf,fs)

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

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

Загрузите файл, содержащий аудио данных от тихоокеанского синего кита, дискретизированный с частотой дискретизации 4 кГц. Файл получен из библиотеки вокализаций животных, поддерживаемой Программой исследований биоакустики Корнеллского университета. Шкала времени в данных сжимается в 10 раз, чтобы поднять тангаж и сделать вызовы более слышимыми. Преобразуйте сигнал в расписание MATLAB ® и постройте график. Из шума в сигнале выделяются четыре функции. Первый известен как трель, а другие три известны как стоны.

[w,fs] = audioread('bluewhale.wav');
whale = timetable(w,'SampleRate',fs);
stackedplot(whale);

Figure contains an object of type stackedplot.

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

emd(whale,'MaxNumIMF',3)

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.

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

imf = emd(whale,'MaxNumIMF',3,'Display',1);
Current IMF  |  #Sift Iter  |  Relative Tol  |  Stop Criterion Hit  
      1      |        1     |      0.13523   |  SiftMaxRelativeTolerance
      2      |        2     |     0.030198   |  SiftMaxRelativeTolerance
      3      |        2     |      0.01908   |  SiftMaxRelativeTolerance
Decomposition stopped because maximum number of intrinsic mode functions was extracted.

Используйте вычисленные МВФ, чтобы построить график Гильберта спектра сигнала. Ограничьте частотную область значений от 0 Гц до 1400 Гц.

hht(imf,'FrequencyLimits',[0 1400])

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

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

[hs,f,t] = hht(imf,'FrequencyLimits',[0 1400]);

mesh(seconds(t),f,hs,'EdgeColor','none','FaceColor','interp')
xlabel('Time (s)')
ylabel('Frequency (Hz)')
zlabel('Instantaneous Energy')

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

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

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

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

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

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

[imf,residual,info] = emd(X,'Interpolation','pchip');

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

Вычислите параметры спектра Гильберта: hs спектра Гильберта, вектор частоты f, временной вектор t, мгновенная частота imfinsf, и мгновенная энергия imfinse.

[hs,f,t,imfinsf,imfinse] = hht(imf,fs);

Используйте вычисленные параметры спектра Гильберта для частотно-временного анализа и диагностики сигналов.

Сгенерируйте многокомпонентный сигнал, состоящий из трех синусоидов частот 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.

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

Подшипник тангажа диаметром 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-dpcosθ],

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

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

Используйте pulstran функция для моделирования влияний как периодического train 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

Figure contains an axes. The axes contains 3 objects of type line.

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

xlim([xLimLeft xLimRight])

Figure contains an axes. The axes contains 3 objects of type line.

Добавьте белый Гауссов шум к сигналам. Задайте отклонение шума 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')

Figure contains an axes. The axes contains 2 objects of type line. These objects represent Healthy, Damaged.

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

imfGood = emd(yGood,'MaxNumIMF',5,'Display',1);
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
Decomposition stopped because maximum number of intrinsic mode functions was extracted.

Использование emd без выходных аргументов для визуализации первых трех МВФ и невязки.

emd(yGood,'MaxNumIMF',5)

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.

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

imfBad = emd(yBad,'MaxNumIMF',5,'Display',1);
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
Decomposition stopped because maximum number of intrinsic mode functions was extracted.
emd(yBad,'MaxNumIMF',5)

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.

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

figure
hht(imfBad(:,1),fs)

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

Спектр Гильберта третьего режима показывает резонанс в сигнале вибрации. Ограничьте частотную область значений от 0 Гц до 100 Гц.

hht(imfBad(:,3),fs,'FrequencyLimits',[0 100])

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

Для сравнения постройте график спектров Гильберта первого и третьего режимов исправного сигнала подшипника.

subplot(2,1,1)
hht(imfGood(:,1),fs)
subplot(2,1,2)
hht(imfGood(:,3),fs,'FrequencyLimits',[0 100])

Figure contains 2 axes. Axes 1 with title Hilbert Spectrum contains an object of type patch. Axes 2 with title Hilbert Spectrum contains an object of type patch.

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

свернуть все

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

hht обрабатывает каждый столбец в imf как внутреннюю функцию mode. Для получения дополнительной информации об вычислениях imf, см. emd.

Частота дискретизации, заданная как положительная скалярная величина. Если fs не подается, нормированная частота используется для вычисления спектра Гильберта. Если imf задается как расписание, из него выводится частота дискретизации.

Расположение оси частоты на графике, заданное как 'yaxis' или 'xaxis'. Чтобы отобразить данные о частоте на y -оси или x -оси графика, задайте freqlocation как 'yaxis' или 'xaxis' соответственно.

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

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

Пример: 'FrequencyResolution',1

Пределы частоты для вычисления спектра Гильберта, заданные как разделенная разделенными запятой парами, состоящая из 'FrequencyLimits'и вектор с целым числом 1 на 2. FrequencyLimits указывается в Гц.

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

Задайте FrequencyResolution в Гц. Если 'FrequencyResolution' не задан, значение (fhigh - flow )/100 выводится из FrequencyLimits. Здесь fhigh верхний предел FrequencyLimits и flow - нижний предел.

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

MinThreshold устанавливает элементы hs 0, когда соответствующие элементы 10log10(hs) меньше MinThreshold.

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

свернуть все

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

Значения частоты сигнала, возвращенные как вектор. hht использует вектор частоты f и временной вектор t создать график спектра Гильберта.

Математически, f обозначается как: f = flow: fres: fhigh, где fres - разрешение частоты.

Значения времени сигнала, возвращенные в виде вектора или a duration массив. hht использует временной вектор t и вектор частоты f создать график спектра Гильберта.

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

  • Массив, если imf задается как массив.

  • A duration массив, если imf задается как равномерно выбранное расписание.

Мгновенная частота каждого МВФ, возвращаемая в виде вектора, матрицы или расписания.

imfinsf имеет одинаковое число столбцов как imf и возвращается как:

  • Вектор, если imf задается как вектор.

  • Матрица, если imf задается как матрица.

  • Расписание, если imf задается как равномерно выбранное расписание.

Мгновенная энергия каждого МВФ, возвращаемая в виде вектора, матрицы или расписания.

imfinse имеет одинаковое число столбцов как imf и возвращается как:

  • Вектор, если imf задается как вектор.

  • Матрица, если imf задается как матрица.

  • Расписание, если imf задается как равномерно выбранное расписание.

Алгоритмы

Преобразование Гильберта-Хуанга полезно для выполнения частотно-временного анализа нестационарных и нелинейных данных. Процедура Гильберта-Хуана состоит из следующих этапов:

  1. emd разлагает набор данных x на конечное число функций внутреннего режима.

  2. Для каждой функции внутреннего режима, xi, функции hht:

    1. Использование hilbert для вычисления аналитического сигнала, zi(t)=xi(t)+jH{xi(t)}, где H {xi} - преобразование Гильберта xi.

    2. Выражает zi как zi(t)=ai(t)ejθi(t), где ai (t) является мгновенной амплитудой иθi(t) - мгновенная фаза.

    3. Вычисляет мгновенную энергию, |ai(t)|2и мгновенную частоту, ωi(t)dθi(t)/dt. Если задана частота дискретизации, hht новообращенные ωi(t) на частоту в Гц.

    4. Выводит мгновенную энергию в imfinse и мгновенную частоту в imfinsf.

  3. При вызове без выходные аргументы, hht строит график энергии сигнала как функции времени и частоты с цветом, пропорциональным амплитуде.

Ссылки

[1] Хуан, Норден Е и Сэмюэль С. П. Шен. Преобразование Гильберта-Хуанга и его приложения. 2-й ред. Том 16. Междисциплинарные математические науки. WORLD SCIENTIFIC, 2014. https://doi.org/10.1142/8804.

[2] Хуан, Норден Э., Чжаохуа У, Стивен Р. Лонг, Кеннет К. Арнольд, Сяньяо Чэнь и Карин Бланк. «НА МГНОВЕННОЙ ЧАСТОТЕ». Усовершенствования в области адаптивного анализа данных 01, № 02 (апрель 2009 года): 177-229. https://doi.org/10.1142/S1793536909000096.

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

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