Exponenta Banner

hht

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

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

Синтаксис

hs = hht(imf)
hs = hht(imf,fs)
[hs,f,t] = hht(___)
[hs,f,t,imfinsf,imfinse] = hht(___)
[___] = hht(___,Name,Value)
hht(___)
hht(___,freqlocation)

Описание

пример

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 - ось по умолчанию.

Примеры

свернуть все

Скрипт Open Live Script

Сгенерируйте Модулируемый гауссовым образом квадратичный щебет. Задайте частоту дискретизации 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)

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

emd(q)

Вычислите IMFs сигнала. Используйте '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
The decomposition stopped because the number of extrema of the residual signal is less than 'MaxNumExtrema'.

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

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

Скрипт Open Live Script

Загрузите и визуализируйте неустановившийся непрерывный сигнал, состоявший из синусоидальных волн с отличным изменением в частоте. Вибрация отбойного молотка и звук фейерверков являются примерами неустановившихся непрерывных сигналов. Сигнал производится на уровне 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');

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

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

hht(imf,fs)

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

Скрипт Open Live Script

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

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

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

emd(whale,'MaxNumIMF',3)

Вычислите первые три IMFs сигнала. Используйте '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
The decomposition stopped because 'MaxNumIMF' was reached.

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

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

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

[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')

Скрипт Open Live Script

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

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

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

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

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

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

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

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

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

Скрипт 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')

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

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

Скрипт Open Live Script

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

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

imfGood = emd(yGood,'MaxNumIMF',5,'Display',1);
Current IMF  |  #Sift Iter  |  Relative Tol  |  Stop Criterion Hit  
      1      |        3     |      0.01727   |  SiftMaxRelativeTolerance
      2      |        3     |     0.059819   |  SiftMaxRelativeTolerance
      3      |        4     |      0.14846   |  SiftMaxRelativeTolerance
      4      |        1     |     0.020337   |  SiftMaxRelativeTolerance
      5      |        2     |     0.023529   |  SiftMaxRelativeTolerance
The decomposition stopped because 'MaxNumIMF' was reached.

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

emd(yGood,'MaxNumIMF',5)

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

imfBad = emd(yBad,'MaxNumIMF',5,'Display',1);
Current IMF  |  #Sift Iter  |  Relative Tol  |  Stop Criterion Hit  
      1      |        2     |     0.038341   |  SiftMaxRelativeTolerance
      2      |        3     |     0.085968   |  SiftMaxRelativeTolerance
      3      |        6     |      0.17468   |  SiftMaxRelativeTolerance
      4      |        1     |      0.01543   |  SiftMaxRelativeTolerance
      5      |        2     |     0.025733   |  SiftMaxRelativeTolerance
The decomposition stopped because 'MaxNumIMF' was reached.

emd(yBad,'MaxNumIMF',5)

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

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

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

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

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

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

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

свернуть все

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

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

Частота дискретизации в виде положительной скалярной величины. Если fs не предоставляется, нормированная частота 2π используется для расчета Гильбертов спектр. Если 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 является разрешением частоты.

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

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

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

  • duration массив, если imf задан как однородно произведенное расписание.

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

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

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

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

  • Расписание, если imf задан как однородно произведенное расписание.

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

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

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

  • Матрица A, если 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 редактора. Междисциплинарные Математические Науки. НАУЧНЫЙ МИР, 2014. https://doi.org/10.1142/8804.

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

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

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

| |

Темы

Введенный в R2018a

Документация Signal Processing Toolbox
Поддержка