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

Терпение диаметра подачи 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             49034  timetable              
  ttbenv      2000x1             49042  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             49272  timetable              
  btbenv      2000x2             49292  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])