envspectrum

Огибающая спектра для диагностики машинного оборудования

Описание

пример

es = envspectrum(x,fs) возвращает огибающую спектра x сигнала произведенный на уровне fs. Если x матрица, затем функция вычисляет огибающую спектра независимо для каждого столбца и возвращает результат в соответствующий столбец es.

пример

es = envspectrum(xt) возвращает огибающую спектра сигнала, сохраненного в расписании MATLAB® xt.

пример

es = envspectrum(___,Name,Value) задает дополнительные опции для любого из предыдущих аргументов пары "имя-значение" использования синтаксисов. Опции включают алгоритм, использовался для расчета сигнала конверта и диапазона частот, по которому можно оценить спектр.

[es,f,env,t] = envspectrum(___) возвращает f, вектор частот, в который es вычисляется; env, сигнал конверта; и t, времена, в который env вычисляется.

envspectrum(___) без выходных аргументов строит сигнал конверта и огибающую спектра в текущей фигуре.

Примеры

свернуть все

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

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

BPFO=12nf0[1-dpcosθ],

где 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

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

  • Первый канал является ослабленной синусоидой с затуханием постоянных 700 и частоты синусоиды 600 Гц.

  • Второй канал является другой ослабленной синусоидой с затуханием постоянных 800 и частоты синусоиды 500 Гц. Второй канал изолирует первый канал 5 миллисекундами.

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

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             48978  timetable              
  ttbenv      2000x1             48986  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             49200  timetable              
  btbenv      2000x2             49220  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])

Входные параметры

свернуть все

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

Пример: cos(pi/4*(0:159))+randn(1,160) одноканальный сигнал вектора-строки.

Пример: cos(pi./[4;2]*(0:159))'+randn(160,2) двухканальный сигнал.

Типы данных: single | double
Поддержка комплексного числа: Да

Частота дискретизации в виде положительного действительного скаляра.

Типы данных: single | double

Введите расписание. xt должен содержать увеличивающиеся конечные времена строки. Если xt представляет многоканальный сигнал, затем он должен иметь или одну переменную, содержащую матрицу или несколько переменных, состоящих из векторов.

Если расписание имеет пропавших без вести или дублирующиеся моменты времени, можно зафиксировать его с помощью советов в Чистом Расписании с Пропавшими без вести, Копией, или Неоднородные Времена (MATLAB).

Пример: timetable(seconds(0:4)',randn(5,2)) задает двухканальную, случайную переменную, произведенную на уровне 1 Гц в течение 4 секунд.

Типы данных: single | double

Аргументы в виде пар имя-значение

Задайте дополнительные разделенные запятой пары Name,Value аргументы. Name имя аргумента и Value соответствующее значение. Name должен появиться в кавычках. Вы можете задать несколько аргументов в виде пар имен и значений в любом порядке, например: Name1, Value1, ..., NameN, ValueN.

Пример: 'Метод', 'Гильберт', 'FilterOrder', 30, 'полоса', [0 fs/4] вычисляет огибающую спектра между 0 и половина частоты Найквиста с помощью полосового фильтра 30-го порядка и вычисляя конверт аналитического сигнала.

Алгоритм для вычисления конверта сигнализирует в виде разделенной запятой пары, состоящей из 'Method' и любой 'hilbert' или 'demod'. См. Алгоритмы для получения дополнительной информации.

Диапазон частот, чтобы вычислить огибающую спектра в виде разделенной запятой пары, состоящей из 'Band' и двухэлементный вектор строгого увеличения значений между 0 и частота Найквиста.

Типы данных: single | double
Поддержка комплексного числа: Да

Порядок фильтра FIR в виде разделенной запятой пары, состоящей из 'FilterOrder' и положительный целочисленный скаляр.

  • Если 'Method' 'hilbert', затем этот аргумент задает порядок КИХ-полосового фильтра.

  • Если 'Method' 'demod', затем этот аргумент задает порядок КИХ фильтр lowpass.

Типы данных: single | double

Выходные аргументы

свернуть все

Огибающая спектра, возвращенная как вектор или матрица.

Частоты, на которых вычисляется огибающая спектра, возвратились как вектор.

Сигнал конверта, возвращенный как вектор, матрица или расписание.

Если вход к envspectrum расписание, затем env также расписание. Временные стоимости env имейте тот же формат как временные стоимости входного расписания.

  • Если вход является расписанием с одной переменной, содержащей матрицу, то env имеет одну переменную, содержащую матрицу.

  • Если вход является расписанием с несколькими переменными, состоящими из векторов, то env имеет несколько переменных, состоящих из векторов.

Временные стоимости, в которых вычисляется сигнал конверта, возвратились как вектор.

Если вход к envspectrum расписание, затем t имеет тот же формат как временные стоимости входного расписания.

Алгоритмы

envspectrum первоначально удаляет смещение DC из входного сигнала, x, и затем вычисляет сигнал конверта.

  • Если 'Method' установлен в 'hilbert', функция:

    1. Полосовые фильтры сигнал. КИХ-фильтру задал порядок 'FilterOrder' и частоты среза в ba(1) и ba(2), где ba диапазон частот, заданный с помощью 'Band'.

    2. Вычисляет аналитический сигнал с помощью hilbert функция.

    3. Вычисляет сигнал конверта как абсолютное значение аналитического сигнала.

  • Если 'Method' установлен в 'demod', функция:

    1. Выполняет комплексную демодуляцию сигнала. Сигнал умножается на exp (j 2πf0t), где f 0 = (ba(1) + ba(2))/2.

    2. Фильтры lowpass демодулируемый сигнал вычислить аналитический сигнал. КИХ-фильтру задал порядок 'FilterOrder' и частота среза (ba(2)ba(1))/2.

    3. Вычисляет сигнал конверта как дважды абсолютное значение аналитического сигнала.

После вычисления сигнала конверта функция удаляет смещение DC из конверта и вычисляет огибающую спектра с помощью БПФ.

Ссылки

[1] Рэндалл, связь Роберта. Основанный на вибрации мониторинг состояния. Чичестер, Великобритания: John Wiley & Sons, 2011.

Расширенные возможности

Введенный в R2017b

Для просмотра документации необходимо авторизоваться на сайте