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

Подшипник с диаметром шага 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) подшипника,

где $${}_{\mathrm{f0}}$$ - скорость приведения в движение, 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])

Огибающие спектры сигналов имеют пики при целых кратных частоте повторения 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 введите, как значения времени входного расписания. Расписание вывода имеет тот же размер, что и расписание ввода.

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

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