Сгенерируйте Модулируемый гауссовым образом квадратичный щебет. Задайте частоту дискретизации 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
Decomposition stopped because the number of extrema in the residual signal is less than the 'MaxNumExtrema' value.

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

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

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

Загрузите файл, который содержит аудиоданные от Тихоокеанского голубого кита, произведенного на уровне 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
Decomposition stopped because maximum number of intrinsic mode functions was extracted.

Используйте вычисленный 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')

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

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

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

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

Средний диаметр подшипника 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) подшипника:

где $$f_0$$ пропускная способность, $$n$$ количество прокручивающихся элементов, $$d$$ диаметр прокручивающихся элементов, $$p$$ диаметр тангажа подшипника, и $$\theta$$ угол контакта подшипника. Примите угол контакта 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/150^2$$.

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

emd(yGood,'MaxNumIMF',5)

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

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
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])