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

Эта тригонометрическая идентичность представляет два различных вида одного и того же физического сигнала:

$\frac{\mathrm{}}{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}\mathit{52cos2āf1t}\frac{\mathrm{+}}{\mathrm{}}\mathrm{14\; (}{\mathit{}}_{\mathrm{cos2\delta }}{\mathit{(}}_{\mathrm{}}f1+\mathrm{f2)}t{\mathit{}}_{\mathrm{+}}{\mathit{}}_{\mathrm{}}cos2\lambda (\mathrm{}{\mathrm{f1-f2}}^{\mathrm{)}}{\mathit{t}}_{\mathrm{)}}\mathit{}=\mathrm{(2}+{\mathit{}}_{\mathrm{}}\mathit{}$cos2āf2t) cos2āf1t.

Создать две синусоиды, s и z, такой, что s - сумма трех синусоидальных волн и z - одиночная синусоидальная волна с модулированной амплитудой. Убедитесь, что два сигнала равны, вычисляя бесконечную норму их разности.

t = 0:1e-3:10;
omega1 = 2*pi*100;
omega2 = 2*pi*20;
s = 0.25*cos((omega1-omega2)*t) + 2.5*cos(omega1*t) + 0.25*cos((omega1+omega2)*t);
z = (2+cos(omega2/2*t).^2).*cos(omega1*t);

norm(s-z,Inf) 

ans = 3.2729e-13

Постройте график синусоид и выберите 1-секундный интервал, начинающийся с 2 секунд.

plot(t,[s' z'])
xlim([2 3])
xlabel('Time (s)')
ylabel('Signal')

Получить спектрограмму сигнала. Спектрограмма показывает три различных синусоидальных компонента. Анализ Фурье рассматривает сигналы как наложение синусоидальных волн.

pspectrum(s,1000,'spectrogram','TimeResolution',4)

Использовать emd вычислять функции внутреннего режима сигнала и дополнительную диагностическую информацию. Функция по умолчанию выводит таблицу, которая указывает количество итераций отсечения, относительный допуск и критерий остановки отсева для каждого МВФ. Разложение в эмпирическом режиме воспринимает сигнал как z.

[imf,~,info] = emd(s);

Количество нулевых переходов и локальных экстремумов различается максимум на единицу. Это удовлетворяет необходимому условию для того, чтобы сигнал был МВФ.

info.NumZerocrossing - info.NumExtrema

ans = 1

Постройте график МВФ и выберите 0,5-секундный интервал, начинающийся с 2 секунд. МВФ - это сигнал AM, потому что emd воспринимает сигнал как амплитудно модулированный.

plot(t,imf)
xlim([2 2.5])
xlabel('Time (s)')
ylabel('IMF')

Моделирование сигнала вибрации от поврежденного подшипника. Выполните декомпозицию в эмпирическом режиме для визуализации Т сигнала и поиска дефектов.

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

Следующим шагом в анализе является вычисление гильбертова спектра извлеченных ПН. Дополнительные сведения см. в примере «Вычислить спектр вибросигнала Гильберта» (панель инструментов обработки сигналов).

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

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

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

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

emd(X,'Interpolation','pchip','Display',1)

Current IMF  |  #Sift Iter  |  Relative Tol  |  Stop Criterion Hit  
      1      |        2     |     0.026352   |  SiftMaxRelativeTolerance
      2      |        2     |    0.0039573   |  SiftMaxRelativeTolerance
      3      |        1     |     0.024838   |  SiftMaxRelativeTolerance
      4      |        2     |      0.05929   |  SiftMaxRelativeTolerance
      5      |        2     |      0.11317   |  SiftMaxRelativeTolerance
      6      |        2     |      0.12599   |  SiftMaxRelativeTolerance
      7      |        2     |      0.13802   |  SiftMaxRelativeTolerance
      8      |        3     |      0.15937   |  SiftMaxRelativeTolerance
      9      |        2     |      0.15923   |  SiftMaxRelativeTolerance
Decomposition stopped because the number of extrema in the residual signal is less than the 'MaxNumExtrema' value.

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

Щелкните правой кнопкой мыши по белому пространству на графике, чтобы открыть окно селектора МВФ. Используйте селектор МВФ для выборочного просмотра генерируемых ИС, исходного сигнала и остатка.

Выберите IMFs, который будет показан из списка. Выберите необходимость отображения исходного сигнала и остатка на графике.

Выбранные ПС теперь отображаются на графике.

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