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

Подшипник тангажа диаметром 12 см имеет восемь элементы качения. Каждый элемент качения имеет диаметр 2 см. Внешнее кольцо остается стационарным, внутренне кольцо вращается со скоростью 25 оборотов в секунду. Акселерометр производит измерения колебаний подшипника с частотой дискретизации 10 кГц.

fs = 10000;
f0 = 25;
n = 8;
d = 0.02;
p = 0.12;

Сигнал вибрации от исправного подшипника включает в себя несколько порядков ведущей частоты. Постройте график 0,1 секунды данных.

t = 0:1/fs:1-1/fs;
z = [1 0.5 0.2 0.1 0.05]*sin(2*pi*f0*[1 2 3 4 5]'.*t)/5;

plot(t,z)
xlim([0.4 0.5])

Дефект внешнего кольца подшипника вызывает ряд влияний 5 миллисекунд на подшипник. В конечном счете эти влияния приводят к износу подшипника. Влияния происходят на частоте мяча внешнего кольца (BPFO) подшипника,

где $$f_0$$ - скорость привода, n - количество тел качения, d - диаметр тел качения, p - диаметр тангажа подшипника, Примите контактный угол нуля и вычислите BPFO.

ca = 0;
bpfo = n*f0/2*(1-d/p*cos(ca))

bpfo = 83.3333

Моделируйте каждое влияние как синусоиду с частотой 3 кГц, оконную с плоским верхним окном. Сделать влияние периодическим путем свертки его с помощью функции гребня. Постройте график 0,1 секунды данных.

fImpact = 3000;
tImpact = 0:1/fs:5e-3-1/fs;
xImpact = sin(2*pi*fImpact*tImpact).*flattopwin(length(tImpact))'/10;

xComb = zeros(size(t));
xComb(1:fs/bpfo:end) = 1;

x = conv(xComb,xImpact,'same')/3;

plot(t,x+z)
xlim([0.4 0.5])

Добавьте белый Гауссов шум к сигналам. Задайте отклонение шума 1/30 ². Постройте график 0,1 секунды данных.

yGood = z + randn(size(z))/30;
yBad = x+z + randn(size(z))/30;
plot(t,yGood,t,yBad)
xlim([0.4 0.5])
legend('Healthy','Damaged')

Вычислите и постройте график огибающих сигналов и спектров.

envspectrum([yGood' yBad'],fs)
xlim([0 10*bpfo]/1000)

Сравните пиковые местоположения с частотами гармоник BPFO. Гармоники BPFO в огибающей спектра являются признаком износа подшипника.

harmImpact = (1:10)*bpfo;
[X,Y] = meshgrid(harmImpact,ylim);

hold on
plot(X/1000,Y,':k')
legend('Healthy','Damaged','BPFO harmonics')
hold off

Вычислите спектры сигналов Welch. Задайте разрешение частоты 5 Гц.

figure
pspectrum([yGood' yBad'],fs,'FrequencyResolution',5)
legend('Healthy','Damaged')

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

xlim([0 10*bpfo]/1000)

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

xlim((bpfo*[-10 10]+fImpact)/1000)

Сгенерируйте двухканальный сигнал, который напоминает сигналы вибрации от подшипника, который завершает вращение каждые 10 миллисекунд. Сигнал дискретизируется на частоте 10 кГц в течение 0,2 секунд, что соответствует 20 вращениям подшипника.

fs = 10000;
tmax = 20;
mlt = 0.01;
t = 0:1/fs:mlt-1/fs;

В течение каждого 10-миллисекундного интервала:

Постройте график сигнала.

y1 = sin(2*pi*600*t).*exp(-700*t);
y2 = sin(2*pi*500*t).*exp(-800*t);
y2 = [y2(51:100) y2(1:50)];

T = (0:1/fs:mlt*tmax-1/fs)';
Y = repmat([y1;y2],1,tmax)';

plot(T,Y)

Создайте массив длительности с использованием временного интервала T. Создайте расписание с массивом длительности и двухканальным сигналом.

dt = seconds(T);
ttb = timetable(dt,Y);

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

envspectrum(ttb,'Band',[100 4900])

Огибающие спектров сигналов имеют peaks с целочисленными кратными частотой повторения 1/0,01 = 0,1 кГц. Это так же, как и ожидалось. envspectrum удаляет высокочастотные синусоидальные компоненты и фокусируется на более низкочастотном повторении. Вот почему огибающая спектра является полезным инструментом для анализа вращательного машинного оборудования.

Вычислите огибающий сигнал и время, в которое он вычисляется. Проверяйте типы выходных переменных.

[~,~,ttbenv,ttbt] = envspectrum(ttb,'Band',[100 4900]);
whos ttb*

  Name           Size            Bytes  Class        Attributes

  ttb         2000x1             48977  timetable              
  ttbenv      2000x1             48985  timetable              
  ttbt        2000x1             16002  duration               

Система временного вектора duration введите, как и временные значения входа timetable. Расписание выхода имеет тот же размер, что и расписание входа.

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

btb = timetable(dt,Y(:,1),Y(:,2));

[~,~,btbenv,btbt] = envspectrum(btb,'Band',[100 4900]);
whos btb*

  Name           Size            Bytes  Class        Attributes

  btb         2000x2             49199  timetable              
  btbenv      2000x2             49219  timetable              
  btbt        2000x1             16002  duration               

Расписание выхода имеет тот же размер, что и расписание входа.

Сгенерируйте сигнал, дискретизированный в 1 кГц в течение 5 секунд. Сигнал состоит из 0,01-секундных прямоугольных импульсов, которые повторяются каждые T = 0,25 секунды. Амплитуда модулирует сигнал на синусоиду несущей частоты 150 Гц.

fs = 1e3;
tmax = 5;

t = 0:1/fs:tmax;
y = pulstran(t,0:0.25:tmax,'rectpuls',0.01);

fc = 150;
z = modulate(y,fc,fs);

Постройте график исходных и модулированных сигналов. Показывать только первые несколько циклов.

plot(t,y,t,z,'-')
grid on
axis([0 1 -1.1 1.1])

Вычислите огибающую и огибающий спектра сигнала. Определите огибающую сигнала с помощью комплексной демодуляции. Вычислите огибающую спектра на интервале 20 Гц с центром на несущей частоте.

[q,f,e,te] = envspectrum(z,fs,'Method','demod','Band',[fc-10 fc+10]);

Постройте график огибающего сигнала и огибающего спектра. Изменение масштаба интервала от 0 до 50 Гц.

subplot(2,1,1)
plot(te,e)
xlabel('Time')
title('Envelope')

subplot(2,1,2)
plot(f,q)
xlim([0 50])
xlabel('Frequency')
title('Envelope Spectrum')

Огибающий сигнал имеет тот же период времени, T = 0,25 секунды, что и исходный сигнал. Огибающая спектра имеет импульсы 1/T = 4 Гц.

Повторите расчеты, но теперь используйте hilbert функция для вычисления огибающей. Полосу пропускания - фильтрация сигнала с помощью фильтра 10-го порядка с конечной импульсной характеристикой (КИХ). Постройте график огибающего сигнала и огибающей спектра с помощью встроенной функциональности envspectrum.

envspectrum(z,fs,'Method','hilbert','FilterOrder',10)

Встройте сигнал в белый Гауссов шум отклонения 1/3. Постройте график результата.

zn = z + randn(size(z))/3;

plot(t,zn,'-')
grid on
axis([0 1 -1.1 1.1])

Вычислите и отобразите огибающий сигнал и огибающую спектра. Вычислите огибающую спектра с помощью комплексной демодуляции на интервале 10 Гц с центром на несущей частоте. Изменение масштаба интервала от 0 до 50 Гц.

envspectrum(zn,fs,'Band',[fc-5 fc+5])
xlim([0 50])