Создайте моделируемый сигнал с частотой 600 Гц в течение 5 секунд. Испытываемая система увеличивает скорость вращения с 10 до 40 оборотов в секунду в течение периода наблюдения.

Создайте показания тахометра.

fs = 600;
t1 = 5;
t = 0:1/fs:t1;

f0 = 10;
f1 = 40;
rpm = 60*linspace(f0,f1,length(t));

Сигнал состоит из четырех гармонически связанных чирп с порядками 1, 0,5, 4 и 6. Чирп порядка 4 имеет вдвое большую амплитуду, чем остальные. Чтобы создать чирпы, используйте трапециевидное правило для выражения фазы как интеграла скорости вращения.

o1 = 1;
o2 = 0.5;
o3 = 4;
o4 = 6;

ph = 2*pi*cumtrapz(rpm/60)/fs;

x = [1 1 2 1]*cos([o1 o2 o3 o4]'*ph);

Визуализируйте карту частота-частота вращения сигнала.

rpmfreqmap(x,fs,rpm)

Проанализируйте смоделированные данные с акселерометра, размещенного в кабине вертолета.

Загрузите данные вертолета. Вибрационные измерения, vib, отбирают пробы со скоростью 500 Гц в течение 10 секунд. Проверка данных показывает, что она имеет линейный тренд. Устраните тенденцию, чтобы предотвратить ухудшение качества оценки частоты.

load('helidata.mat')

vib = detrend(vib);

Постройте график нелинейного профиля RPM. Ротор поднимается вверх до тех пор, пока не достигнет максимальной скорости вращения около 27600 оборотов в минуту, а затем опускается вниз.

plot(t,rpm)
xlabel('Time (s)')
ylabel('RPM')

Вычислите карту частота-частота вращения. Укажите разрешающую способность 2,5 Гц.

[map,freq,rpmOut,time] = rpmfreqmap(vib,fs,rpm,2.5);

Визуализируйте карту.

imagesc(time,freq,map)
ax = gca;
ax.YDir = 'normal';
xlabel('Time (s)')
ylabel('Frequency (Hz)')

Повторите вычисления с использованием полосы пропускания с более высоким разрешением. Постройте карту, используя встроенные функциональные возможности rpmfreqmap. Усиление частотного разрешения происходит за счет временного разрешения.

rpmfreqmap(vib,fs,rpm,1.5);

Создайте сигнал, который состоит из двух линейных чирп и квадратичной чирп, все дискретизированные при 600 Гц в течение 15 секунд. Система, которая выдает сигнал, увеличивает скорость его вращения с 10 до 40 оборотов в секунду в течение периода испытаний.

Создайте показания тахометра.

fs = 600;
t1 = 15;
t = 0:1/fs:t1;

f0 = 10;
f1 = 40;
rpm = 60*linspace(f0,f1,length(t));

Линейные чирки имеют порядки 1 и 2,5. Компонент порядка 1 имеет половину амплитуды другого. Квадратичная чирп начинается с порядка 6 и возвращается к этому порядку в конце измерения. Его амплитуда составляет 0,8. Создайте сигнал, используя эту информацию.

o1 = 1;
o2 = 2.5;
o6 = 6;

x = 0.5*chirp(t,o1*f0,t1,o1*f1)+chirp(t,o2*f0,t1,o2*f1) + ...
    0.8*chirp(t,o6*f0,t1,o6*f1,'quadratic');

Вычислите карту частота-частота вращения сигнала. Используйте амплитуду пика в каждой измерительной ячейке. Задайте разрешение 6 Гц. Окно данных с плоским верхним окном.

[map,fr,rp] = rpmfreqmap(x,fs,rpm,6, ...
    'Amplitude','peak','Window','flattopwin');

Нарисуйте карту частоты вращения в виде графика водопада.

[FR,RP] = meshgrid(fr,rp);
waterfall(FR,RP,map')

view(-6,60)
xlabel('Frequency (Hz)')
ylabel('RPM')
zlabel('Amplitude')

Постройте интерактивную карту частоты вращения с помощью вызова rpmfreqmap без выходных аргументов.

load helidata.mat
rpmfreqmap(detrend(vib),fs,rpm,2)