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

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

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

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

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{}}_{\mathrm{f0}}$ - начальная частота, а $\mathit{T}\mathrm{=}\mathrm{}$0,75 секунды - время затухания.

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

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

Вычислите и постройте график координат x и y лезвия. Добавьте белый гауссов шум дисперсии $\mathrm{}{\mathrm{}}^{0.12}$.

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. На вкладке Экспорт (Export) щелкните Экспорт (Export) и выберите 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])