Визуализируйте функцию частотной характеристики возбуждения молота с одним входом/одним выходом.

Загрузите файл данных, который содержит:

Дискретизация сигналов осуществляется на частоте 4 кГц. Постройте график возбуждения и выходных сигналов.

load modaldata

subplot(2,1,1)
plot(thammer,Xhammer(:))
ylabel('Force (N)')
subplot(2,1,2)
plot(thammer,Yhammer(:))
ylabel('Displacement (m)')
xlabel('Time (s)')

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

clf
winlen = size(Xhammer,1);
modalfrf(Xhammer(:),Yhammer(:),fs,winlen,'Sensor','dis')

Вычислите функции частотной характеристики для системы с двумя входами/двумя выходами, возбужденными случайным шумом.

Загрузите файл данных, содержащий Xrand, входной сигнал возбуждения и Yrand, а отклик системы. Вычислите функции частотной характеристики с помощью 5000-выборочного окна Ханна и 50% перекрытия между смежными сегментами данных. Задайте, что выходные измерения являются перемещениями.

load modaldata
winlen = 5000;

frf = modalfrf(Xrand,Yrand,fs,hann(winlen),0.5*winlen,'Sensor','dis');

Используйте функции графического изображения modalfrf чтобы визуализировать ответы.

modalfrf(Xrand,Yrand,fs,hann(winlen),0.5*winlen,'Sensor','dis')

Оцените функцию частотной характеристики для простой системы с одним входом/одним выходом и сравните ее с определением.

Одномерное дискретное время система состоит из модуля массы, $\mathit{m}$, прикрепленный к стенке пружиной с упругой константой $\mathit{k}=1$. Датчик дискретизирует перемещение массы при ${\mathit{F}}_{\mathrm{s}}=1$ Гц. Демпфер препятствует движению массы, прикладывая к ней силу, пропорциональную скорости, с постоянной демпфирования $\mathit{b}=0.01$.

Сгенерируйте 3000 временных выборок. Определите интервал дискретизации $\Delta \mathit{t}=1/{\mathit{F}}_{\mathit{s}}$.

Fs = 1;
dt = 1/Fs;
N = 3000;
t = dt*(0:N-1);
b = 0.01;

Система может быть описана моделью пространства состояний

где $\mathit{x}={\left[\begin{array}{cc}\mathit{r}& \mathit{v}\end{array}\right]}^{\mathit{T}}$ - вектор состояния, $\mathit{r}$ и $\mathit{v}$ являются соответственно перемещением и скоростью массы, $\mathit{u}$ является движущей силой, и $\mathit{y}=\mathit{r}$ - измеренный выход. Матрицы пространства состояний

$\mathit{I}$ является $2\times 2$ единичные, и матрицы пространства состояний в непрерывном времени

Ac = [0 1;-1 -b];
A = expm(Ac*dt);

Bc = [0;1];
B = Ac\(A-eye(2))*Bc;

C = [1 0];
D = 0;

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

rng default
u = randn(1,N)/2;
u(2001:end) = 0;

y = 0;
x = [0;0];
for k = 1:N
    y(k) = C*x + D*u(k);
    x = A*x + B*u(k);
end

plot(t,y)

Оцените модальную функцию частотной характеристики системы. Используйте половину окна Ханна, пока измеренные сигналы. Задайте, что выходом является перемещение массы.

wind = hann(N/2);

[frf,f] = modalfrf(u',y',Fs,wind,'Sensor','dis');

Функция частотной характеристики системы дискретного времени может быть выражена как Z-преобразование передаточной функции временной области системы, оцененное по модулю кругу. Сравните modalfrf оценка с определением.

[b,a] = ss2tf(A,B,C,D);

nfs = 2048;
fz = 0:1/nfs:1/2-1/nfs;
z = exp(2j*pi*fz);
ztf = polyval(b,z)./polyval(a,z);

plot(f,20*log10(abs(frf)))
hold on
plot(fz*Fs,20*log10(abs(ztf)))
hold off
grid
ylim([-60 40])

Оцените естественную частоту и коэффициент затухания для режима вибрации.

[fn,dr] = modalfit(frf,f,Fs,1,'FitMethod','PP')

fn = 0.1593

dr = 0.0043

Сравните естественную частоту с $1/2\pi$, которое является теоретическим значением для незаполненной системы.

theo = 1/(2*pi)

theo = 0.1592

Оцените функцию частотной характеристики и модальные параметры простой системы с несколькими входами/несколькими выходами.

Идеальная одномерная колебательная система состоит из двух масс, ${\mathit{m}}_1$ и ${\mathit{m}}_2$, ограниченный между двумя стенками. Модули такие, что ${\mathit{m}}_1=1$ и ${\mathit{m}}_2=\mu$. Каждая масса присоединена к ближайшей стенке пружиной с упругой константой $\mathit{k}$. Идентичная пружина соединяет две массы. Три демпфера препятствуют движению масс, прикладывая к ним силы, пропорциональные скорости с постоянной демпфирования $\mathit{b}$. Выборка датчиков ${\mathit{r}}_1$ и ${\mathit{r}}_2$, смещения масс, при ${\mathit{F}}_{\mathrm{s}}=50$ Гц.

Сгенерируйте 30 000 временных выборок, эквивалентных 600 секундам. Определите интервал дискретизации $\Delta \mathit{t}=1/{\mathit{F}}_{\mathrm{s}}$.

Fs = 50;
dt = 1/Fs;
N = 30000;
t = dt*(0:N-1);

Система может быть описана моделью пространства состояний

где $$x={\left[\begin{array}{cccc}{r}_1& v_1& r_2& v_2\end{array}\right]}^T$$ - вектор состояния, $$r_i$$ и $$v_i$$ являются соответственно местоположением и скоростью $$i$$th-я масса, $$u={\left[\begin{array}{cc}{u}_1& u_2\end{array}\right]}^T$$ является вектором входа движущих сил, и $$y={\left[\begin{array}{cc}{r}_1& r_2\end{array}\right]}^T$$ - вектор выхода. Матрицы пространства состояний

$\mathit{I}$ является $4\times 4$ единичные, и матрицы пространства состояний в непрерывном времени

Набор $\mathit{k}=400$, $\mathit{b}=0.1$, и $\mu =1/10$.

k = 400;
b = 0.1;
m = 1/10;

Ac = [0 1 0 0;-2*k -2*b k b;0 0 0 1;k/m b/m -2*k/m -2*b/m];
A = expm(Ac*dt);
Bc = [0 0;1 0;0 0;0 1/m];
B = Ac\(A-eye(4))*Bc;
C = [1 0 0 0;0 0 1 0];
D = zeros(2);

Массы управляются случайным входом на протяжении всего измерения. Используйте модель пространства состояний, чтобы вычислить временную эволюцию системы, начиная с полностью нулевого начального состояния.

rng default
u = randn(2,N);

y = [0;0];
x = [0;0;0;0];
for kk = 1:N
    y(:,kk) = C*x + D*u(:,kk);
    x = A*x + B*u(:,kk);
end

Используйте входные и выходные данные, чтобы оценить передаточную функцию системы как функцию от частоты. Используйте 15000-образцевое окно Ханна с 9000 выборками перекрытия между смежными сегментами. Задайте, что измеренные выходы являются перемещениями.

wind = hann(15000);
nove = 9000;
[FRF,f] = modalfrf(u',y',Fs,wind,nove,'Sensor','dis');

Вычислите теоретическую передаточную функцию как Z-преобразование передаточной функции во временной области, оцененное в модуль круге.

nfs = 2048;
fz = 0:1/nfs:1/2-1/nfs;
z = exp(2j*pi*fz);

[b1,a1] = ss2tf(A,B,C,D,1);
[b2,a2] = ss2tf(A,B,C,D,2);

frf(1,:,1) = polyval(b1(1,:),z)./polyval(a1,z);
frf(1,:,2) = polyval(b1(2,:),z)./polyval(a1,z);
frf(2,:,1) = polyval(b2(1,:),z)./polyval(a2,z);
frf(2,:,2) = polyval(b2(2,:),z)./polyval(a2,z);

Постройте графики оценок и наложите теоретические предсказания.

for jk = 1:2
    for kj = 1:2
        subplot(2,2,2*(jk-1)+kj)
        plot(f,20*log10(abs(FRF(:,jk,kj))))
        hold on
        plot(fz*Fs,20*log10(abs(frf(jk,:,kj))))
        hold off
        axis([0 Fs/2 -100 0])
        title(sprintf('Input %d, Output %d',jk,kj))
    end
end

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

figure
modalfrf(u',y',Fs,wind,nove,'Sensor','dis')

Оцените естественные частоты, коэффициенты затухания и формы режима системы. Используйте метод пикового комплектования для вычисления.

[fn,dr,ms] = modalfit(FRF,f,Fs,2,'FitMethod','pp');
fn

fn = 
fn(:,:,1) =

    3.8466    3.8466
    3.8495    3.8495


fn(:,:,2) =

    3.8492    3.8490
    3.8552   14.4684

Сравните естественные частоты с теоретическими предсказаниями для незаполненной системы.

undamped = sqrt(eig([2*k -k;-k/m 2*k/m]))/2/pi

undamped = 2×1

    3.8470
   14.4259

Вычислите функцию частотной характеристики набора данных с двумя входами/шестью выходами, соответствующего стальной системе координат.

Загрузите структуру, содержащую входные возбуждения и измерения выходного акселерометра. Система отбирается с частотой дискретизации 1024 Гц в течение около 3,9 секунд.

load modaldata SteelFrame
X = SteelFrame.Input;
Y = SteelFrame.Output;
fs = SteelFrame.Fs;

Используйте метод subspace, чтобы вычислить функции частотной характеристики. Разделите входной и выходной сигналы на неперекрывающиеся, 1000-выборочные сегменты. Окно каждого сегмента с помощью прямоугольного окна. Задайте порядок модели 36.

[frf,f] = modalfrf(X,Y,fs,1000,'Estimator','subspace','Order',36);

Визуализация диаграммы стабилизации системы. Идентифицируйте до 15 физических режимов.

modalsd(frf,f,fs,'MaxModes',15)