Загрузка и визуализация нестационарного непрерывного сигнала, состоящего из синусоидальных волн с отчетливым изменением частоты. Вибрация отбойного молотка и звук фейерверка являются примерами нестационарных непрерывных сигналов. Сигнал дискретизируется со скоростью 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{5}{2}\mathrm{cos2}\pi {\mathit{f}}_1\mathit{t}+\frac{1}{4}(\mathrm{cos2}\pi \left({\mathit{f}}_1+{\mathit{f}}_2\right)t+\mathrm{cos2}\pi \left({\mathit{f}}_1-{\mathit{f}}_2\right)t)=(2+{\cos }^2\pi {\mathit{f}}_2\mathit{t})\mathrm{cos2}\pi {\mathit{f}}_1\mathit{t}$.

Сгенерируйте две синусоиды, 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) подшипника:

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

Следующим шагом в анализе является вычисление спектра Гильберта извлеченных МВФ. Для получения дополнительной информации смотрите пример Compute Hilbert Spectrum of Vibration Signal (Signal Processing Toolbox).

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

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

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

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

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