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

emd(q)

Вычислите МВФ сигнала. Используйте '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])

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

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

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

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

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

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

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

hht(imf,fs)

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

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

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

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

emd(whale,'MaxNumIMF',3)

Вычислите первые три МВФ сигнала. Используйте '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])

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

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

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

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

Используйте вычисленные МВФ, чтобы построить график спектра многокомпонентного сигнала. Ограничьте частотную область значений [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 (Signal Processing Toolbox) функция, чтобы смоделировать влияния как периодический 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

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

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 (Signal Processing Toolbox), чтобы выполнить эмпирическое разложение исправного сигнала подшипника. Вычислите первые пять функций внутреннего режима (МВФ). Используйте '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)

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

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