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

Средний диаметр подшипника 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

Вычислите валлийские спектры сигналов. Задайте разрешение частоты 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 введите, как временные стоимости входного расписания. Выходное расписание имеет тот же размер как входное расписание.

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

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])