Сгенерируйте сигнал, который напоминает вибрации, вызванные реваншированием двигателя автомобиля. Сигнал дискретизируется на частоте 1 кГц в течение 30 секунд и содержит три гармонических компонентов порядков 1 , 2.4 и 3 с амплитудами 5, 4 и 0,5 соответственно. Встройте сигнал в единично-дисперсионный белый Гауссов шум и сохраните его в расписании MATLAB ®. Умножьте текущую частоту на 60, чтобы получить профиль об/мин. Постройте график профиля RPM.

fs = 1000;
t = (0:1/fs:30)';

fit = @(a,x) (t-x).^6.*exp(-(t-x)).*((t-x)>=0)*a';

fis = fit([0.4 1 0.6 1],[0 6 13 17]);
phi = 2*pi*cumtrapz(t,fis);

ol = [1 2.4 3];
amp = [5 4 0.5]';
vib = cos(phi.*ol)*amp + randn(size(t));

xt = timetable(seconds(t),vib);

plot(t,fis*60)

Используйте rpmtrack функция для вывода профиля частоты вращения из сигнала вибрации. Используйте четыре точки в 5 секундах интервалов, чтобы задать гребень, соответствующий порядку 2.4. Отображение сводных данных расписания выхода.

ndx = (5:5:20)*fs;
order = ol(2);

p = [t(ndx) order*fis(ndx)];

rpmest = rpmtrack(xt,order,p);

summary(rpmest)

RowTimes:

    tout: 30001x1 duration
        Values:
            Min           0 sec     
            Median        15 sec    
            Max           30 sec    
            TimeStep      0.001 sec 

Variables:

    rpm: 30001x1 double

        Values:

            Min       2.2204e-16
            Median        4327.2
            Max           8879.8

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

hold on
plot(seconds(rpmest.tout),rpmest.rpm,'.-')
plot(t(ndx),fis(ndx)*60,'ok')
hold off
legend('Original','Reconstructed','Ridge points','Location','northwest')

Используйте извлеченный профиль об/мин, чтобы сгенерировать карту порядка об/мин сигнала.

rpmordermap(vib,fs,rpmest.rpm)

Восстановите и постройте график форм волны во временной области, которые составляют сигнал. Изменение масштаба временного интервала, возникающего после распада переходных процессов.

xrc = orderwaveform(vib,fs,rpmest.rpm,ol);

figure
plot(t,xrc)
legend([repmat('Order = ',[3 1]) num2str(ol')])
xlim([5 20])

Оцените профиль частоты вращения блейда вентилятора, так как он замедляется после выключения.

Промышленный вентилятор на крыше вращается с частотой вращения 20 000 об/мин. Сопротивление воздуха (с незначительным вкладом от трения подшипника) заставляет ротор вентилятора остановиться приблизительно за 6 секунд. Высокоскоростная камера измеряет координату X одного из блейдов вентилятора со скоростью 1 кГц.

fs = 1000;
t = 0:1/fs:6-1/fs;

rpm0 = 20000;

Идеализируйте блейд вентилятора как точечную массу, окружающую центр ротора в радиусе 50 см. Блейд испытывает силу сопротивления, пропорциональную скорости, что приводит к следующему выражению для угла фазы:

где ${\mathit{f}}_0$ является начальной частотой и $\mathit{T}=0.75$ второе - время затухания.

a = 0.5;
f0 = rpm0/60;
T = 0.75;

phi = 2*pi*f0*T*(1-exp(-t/T));

Вычислите и постройте графики координат X и Y блейдов. Добавьте белый Гауссов шум отклонения $0.1^2$.

x = a*cos(phi) + randn(size(phi))/10;
y = a*sin(phi) + randn(size(phi))/10;

plot(t,x,t,y)

Используйте rpmtrack функция для определения профиля RPM. Напечатать

rpmtrack(x,fs)

в командной строке, чтобы открыть интерактивный рисунок.

Используйте ползунок, чтобы настроить разрешение частот частотно-временной карты приблизительно на 11 Гц. Примите, что компонент сигнала соответствует порядку 1, и установите время окончания извлечения гребня 3,0 секунды. Используйте перекрестный курсор во временной карте и кнопку «Добавить», чтобы добавить три точки, лежащие на гребне. (Также дважды щелкните курсор, чтобы добавить точки в выбранных местоположениях.) Нажмите Estimate, чтобы отследить и извлечь профиль RPM.

Проверьте, что профиль RPM затухает экспоненциально. На вкладке Экспорт нажмите Экспорт и выберите Generate MATLAB Script. Скрипт появится в редакторе.

% MATLAB Code from rpmtrack GUI

% Generated by MATLAB 9.4 and Signal Processing Toolbox 8.0

% Generated on 18-Dec-2017 19:00:35

% Set sample rate
fs = 1000.0000;
% Set order of ridge of interest
order = 1.0000;
% Set ridge points on ridge of interest
ridgePoints = [...
    0.4501 179.8246;...
    0.9944 88.5965;...
    2.4161 11.4035];
% Estimate RPM
[rpmOut,tOut] = rpmtrack(x,fs,order,ridgePoints,...
    'Method','stft',...
    'FrequencyResolution',11.1612,...
    'PowerPenalty',Inf,...
    'FrequencyPenalty',0.0000,...
    'StartTime',0.0000,...
    'EndTime',3.0000);

Запустите скрипт. Отобразите профиль RPM на семилогарифмическом графике.

semilogy(tOut,rpmOut)
ylim([500 20000])