Генерировать квадратичную чирпу с гауссовой модуляцией. Укажите частоту дискретизации 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)

Вычислите 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.

Используйте вычисленные ЧС для построения графика гильбертова спектра квадратичного чирпа. Ограничьте диапазон частот от 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) подшипника:

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

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

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

Добавьте белый гауссов шум к сигналам. Задайте дисперсию шума $$\mathrm{}\mathrm{}\mathrm{}{\mathrm{}}^{\mathrm{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 (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])