Документация

fn = modalfit(frf,f,fs,mnum) оценивает естественные частоты mnum режимы системы с измеренными функциями частотной характеристики frf заданные на частотах f и для частоты дискретизации fs.

fn = modalfit(frf,f,fs,mnum,Name,Value) задает дополнительные опции, используя аргументы пары "имя-значение".

[fn,dr,ms] = modalfit(___) также возвращает коэффициенты затухания и векторы формы режима, соответствующие каждой собственной частоте в fn, с использованием любой комбинации входов из предыдущих синтаксисов.

[fn,dr,ms,ofrf] = modalfit(___) также возвращает восстановленный массив функции частотной характеристики на основе оцененных модальных параметров.

[___] = modalfit(sys,f,mnum,Name,Value) оценивает модальные параметры идентифицированной модели sys. Используйте команды оценки, такие как ssest (System Identification Toolbox) или tfest (System Identification Toolbox) для создания sys начиная с измеренной функции частотной характеристики или с входных и выходных сигналов во временной области. Этот синтаксис позволяет использовать 'DriveIndex', 'FreqRange', и 'PhysFreq' Аргументы пары "имя-значение". Обычно это требует меньше данных, чем синтаксисы, которые используют непараметрические методы. Для использования этого синтаксиса необходимо иметь лицензию System Identification Toolbox™.

Примеры

Функция частотной характеристики системы SISO

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

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

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

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

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

$\begin{array}{c} x (k + 1) = A x (k) + B u (k), \\ y (k) = C x (k) + D u (k), \end{array}$

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

$A = \exp (A_{c} Δ t), B = A_{c}^{- 1} (A - I) B_{c}, C = [\begin{array}{c} 1 & 0 \end{array}], D = 0,$

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

$A_{c} = [\begin{array}{c} 0 & 1 \\ - 1 & - b \end{array}], B_{c} = [\begin{array}{c} 0 \\ 1 \end{array}] .$

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)

Figure contains an axes. The axes contains an object of type line.

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

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])

Figure contains an axes. The axes contains 2 objects of type line.

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

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

fn = 0.1593

dr = 0.0043

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

theo = 1/(2*pi)

theo = 0.1592

Модальные параметры, использующие метод рациональной функции методом наименьших квадратов

Вычислите модальные параметры модуля Space Station, начиная с его массива функции частотной характеристики (FRF).

Загрузите структуру, содержащую массив FRF с тремя входами/тремя выходами. Отбор проб системы производится при частоте 320 Гц.

load modaldata SpaceStationFRF

frf = SpaceStationFRF.FRF; 
f = SpaceStationFRF.f; 
fs = SpaceStationFRF.Fs;

Извлеките модальные параметры самых низких 24 режимов с помощью метода рациональной функции методом наименьших квадратов.

[fn,dr,ms,ofrf] = modalfit(frf,f,fs,24,'FitMethod','lsrf');

Сравните восстановленный массив FRF с измеренной.

for ij = 1:3
 for ji = 1:3
    subplot(3,3,3*(ij-1)+ji)
    loglog(f,abs(frf(:,ji,ij)))
    hold on
    loglog(f,abs(ofrf(:,ji,ij)))
    hold off
    axis tight
    title(sprintf('In%d -> Out%d',ij,ji))
    if ij==3
        xlabel('Frequency (Hz)')
    end
 end
end

Модальные параметры двухфазного генератора

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

Идеальная одномерная колебательная система состоит из двух масс, $m_{1}$ и $m_{2}$ , ограниченный между двумя стенками. Модули такие, что $m_{1} = 1$ и $m_{2} = μ$ . Каждая масса присоединена к ближайшей стенке пружиной с упругой константой $k$ . Идентичная пружина соединяет две массы. Три демпфера препятствуют движению масс, прикладывая к ним силы, пропорциональные скорости с постоянной демпфирования $b$ . Выборка датчиков $r_{1}$ и $r_{2}$ , смещения масс, при $F_{s} = 50$ Гц.

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

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

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

$\begin{array}{c} x (k + 1) = A x (k) + B u (k), \\ y (k) = C x (k) + D u (k), \end{array}$

где $x = {[\begin{array}{c} r_{1} & v_{1} & r_{2} & v_{2} \end{array}]}^{T}$ - вектор состояния, $r_{i}$ и $v_{i}$ являются соответственно местоположением и скоростью $i$ th-я масса, $u = {[\begin{array}{c} u_{1} & u_{2} \end{array}]}^{T}$ является вектором входа движущих сил, и $y = {[\begin{array}{c} r_{1} & r_{2} \end{array}]}^{T}$ - вектор выхода. Матрицы пространства состояний

$A = \exp (A_{c} Δ t), B = A_{c}^{- 1} (A - I) B_{c}, C = [\begin{array}{c} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{array}], D = [\begin{array}{c} 0 & 0 \\ 0 & 0 \end{array}],$

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

$A_{c} = [\begin{array}{c} 0 & 1 & 0 & 0 \\ - 2 k & - 2 b & k & b \\ 0 & 0 & 0 & 1 \\ k / μ & b / μ & - 2 k / μ & - 2 b / μ \end{array}], B_{c} = [\begin{array}{c} 0 & 0 \\ 1 & 0 \\ 0 & 0 \\ 0 & 1 / μ \end{array}] .$

Набор $k = 400$ , $b = 0.1$ , и $μ = 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

Figure contains 4 axes. Axes 1 with title Input 1, Output 1 contains 2 objects of type line. Axes 2 with title Input 1, Output 2 contains 2 objects of type line. Axes 3 with title Input 2, Output 1 contains 2 objects of type line. Axes 4 with title Input 2, Output 2 contains 2 objects of type line.

Постройте график оценок с помощью синтаксиса 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

Модальные параметры системы MIMO

Вычислите естественные частоты, коэффициенты затухания и формы режима для системы с двумя входами/тремя выходами, возбуждаемой несколькими всплесками случайного шума. Каждый пакет длится 1 секунду, и есть 2 секунды между концом каждого пакета и началом следующего. Данные отбираются с частотой дискретизации 4 кГц.

Загрузите файл данных. Постройте график входных сигналов и выхода сигналов.

load modaldata

subplot(2,1,1)
plot(Xburst)
title('Input Signals')
subplot(2,1,2)
plot(Yburst)
title('Output Signals')

Figure contains 2 axes. Axes 1 with title Input Signals contains 2 objects of type line. Axes 2 with title Output Signals contains 3 objects of type line.

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

burstLen = 12000;
[frf,f] = modalfrf(Xburst,Yburst,fs,burstLen);

Визуализируйте диаграмму стабилизации и верните стабильные естественные частоты. Задайте максимальный порядок модели 30 режимов.

figure
modalsd(frf,f,fs,'MaxModes',30);

Figure contains an axes. The axes with title Stabilization Diagram contains 4 objects of type line. These objects represent Stable in frequency, Stable in frequency and damping, Not stable in frequency, Averaged response function.

Изменение масштаба графика. Усредненная функция отклика имеет максимумы на 373 Гц, 852 Гц и 1371 Гц, которые соответствуют физическим частотам системы. Сохраните максимумы в переменную.

phfr = [373 852 1371];

Вычислите модальные параметры с помощью алгоритма комплексной экпоненты методом наименьших квадратов (LSCE). Задайте порядок модели из 6 режимов и задайте физические частоты для 3 режимов, определенных на диаграмме стабилизации. Функция генерирует один набор собственных частот и коэффициентов затухания для каждой входной ссылки.

[fn,dr,ms,ofrf] = modalfit(frf,f,fs,6,'PhysFreq',phfr);

Постройте график восстановленных функций частотной характеристики и сравните их с исходными таковыми.

for k = 1:2
    for m = 1:3
        subplot(2,3,m+3*(k-1))
        plot(f/1000,10*log10(abs(frf(:,m,k))))
        hold on
        plot(f/1000,10*log10(abs(ofrf(:,m,k))))
        hold off
        text(1,-50,[['Output ';' Input '] num2str([m k]')])
        ylim([-100 -40])
    end
end
subplot(2,3,2)
title('Frequency-Response Functions')

Figure contains 6 axes. Axes 1 contains 3 objects of type line, text. Axes 2 with title Frequency-Response Functions contains 3 objects of type line, text. Axes 3 contains 3 objects of type line, text. Axes 4 contains 3 objects of type line, text. Axes 5 contains 3 objects of type line, text. Axes 6 contains 3 objects of type line, text.

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

`frf` - Функции частотной характеристики
вектор | матрица | трехмерный массив

Функции частотной характеристики, заданные как вектор, матрица или трехмерный массив. frf имеет размер p -by m -by - n, где p - количество интервалов частоты, m - количество ответных сигналов, а n - количество сигналов возбуждения, используемых для оценки передаточной функции.

Пример: tfestimate(randn(1,1000),sin(2*pi*(1:1000)/4)+randn(1,1000)/10) аппроксимирует частотную характеристику генератора.

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

`f` - Частоты
вектор

Частоты, заданные как вектор. Количество элементов f должно равняться количеству строк frf.

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

`fs` - Частота дискретизации данных измерений
положительная скалярная величина

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

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

`mnum` - Количество режимов
положительное целое число

Количество режимов, заданное в виде положительного целого числа.

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

`sys` - Идентифицированная система
модель с идентифицированными параметрами

Идентифицированная система, заданная как модель с идентифицированными параметрами. Используйте команды оценки, такие как ssest (System Identification Toolbox), n4sid (System Identification Toolbox), или tfest (System Identification Toolbox) для создания sys начиная с измеренной функции частотной характеристики или с входных и выходных сигналов во временной области. Смотрите Модальный анализ идентифицированных моделей для примера. Для использования этого входного параметра необходимо иметь лицензию System Identification Toolbox.

Пример: idss([0.5418 0.8373;-0.8373 0.5334],[0.4852;0.8373],[1 0],0,[0;0],[0;0],1) генерирует идентифицированную модель пространства состояний, соответствующий модуль массе, присоединенной к стенке пружиной модуля упругой константы и демпфером с постоянной 0,01. Перемещение массы отбирают с частотой дискретизации 1 Гц.

Пример: idtf([0 0.4582 0.4566],[1 -1.0752 0.99],1) генерирует идентифицированную модель передаточной функции, соответствующий модуль массе, присоединенной к стенке, пружиной модуля упругой константы и демпфером с постоянной 0,01. Перемещение массы отбирают с частотой дискретизации 1 Гц.

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

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

Пример: 'FitMethod','pp','FreqRange',[0 500] использует метод пикового захвата для выполнения подгонки и ограничивает частотную область значений от 0 до 500 Гц.

`'Feedthrough'` - Наличие прохода в расчетной передаточной функции
`false` (по умолчанию) | `true`

Наличие сквозного соединения в предполагаемой передаточной функции, заданной как разделенная разделенными запятой парами, состоящая из 'Feedthrough' и логическое значение. Этот аргумент доступен только в том случае, если 'FitMethod' задается как 'lsrf'.

Типы данных: logical

`'FitMethod'` - Алгоритм аппроксимации
`'lsce'` (по умолчанию) | `'lsrf'` | `'pp'`

Алгоритм аппроксимации, заданный как разделенная разделенными запятой парами, состоящая из 'FitMethod' и 'lsce', 'lsrf', или 'pp'.

'lsce' - Комплексный экспоненциальный метод методом наименьших квадратов. Если вы задаете 'lsce', затем fn - вектор с mnum элементы, независимо от размера frf.
'lsrf' - метод оценки рациональных функций методом наименьших квадратов. Если вы задаете 'lsrf', затем fn - вектор с mnum элементы, независимо от размера frf. Метод описан в [3]. Для получения дополнительной информации смотрите Оценку передаточной функции в непрерывном времени с использованием Частотного диапазона данных в непрерывном времени (System Identification Toolbox). Этот алгоритм обычно требует меньше данных, чем непараметрические подходы, и является единственным, который работает для неоднородных f.
'pp' - Метод пиковой комплектации. Для frf вычисляется из n сигналов возбуждения и m сигналов отклика, fn является mnum-by- m -by- n массив с одной оценкой fn и одну оценку dr в frf.

`'FreqRange'` - Частотная область значений
двухэлементный вектор увеличения положительных значений

Частотная область значений, заданный как разделенная разделенными запятой парами, состоящая из 'FreqRange' и двухэлементный вектор увеличения положительных значений, содержащийся в области значений, указанном в f.

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

`'PhysFreq'` - Естественные частоты для физических режимов
вектор

Естественные частоты для физических режимов, чтобы включить в анализ, заданные как разделенная разделенными запятой парами, состоящая из 'PhysFreq' и вектор значений частот в области значений, охватываемом f. Функция включает в анализ те режимы с естественными частотами, наиболее близкими к значениям, заданным в векторе. Если вектор содержит m значения частоты, то fn и dr иметь m строк каждая, и ms имеет m столбцов. Если вы не задаете этот аргумент, то функция использует всюсь частотную область значений в f.

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

`'DriveIndex'` - Индексы функции частотной характеристики приводной точки
`[1 1]` (по умолчанию) | двухэлементный вектор положительных целых чисел

Индексы функции частотной характеристики приводной точки, заданные как разделенная разделенными запятой парами, состоящая из 'DriveIndex' и двухэлементный вектор положительных целых чисел. Первый элемент вектора должен быть меньше или равен количеству откликов системы. Второй элемент вектора должен быть меньше или равен количеству системных возбуждений. Формы режима нормированы к единичному режиму на основе приводной точки.

Пример: 'DriveIndex',[2 3] задает, что функция частотной характеристики приводной точки frf(:,2,3).

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

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

`fn` - Естественные частоты
матрица | трехмерный массив

Естественные частоты, возвращенные в виде матрицы или трехмерные массивы. Размер fn зависит от выбора алгоритма аппроксимации, заданного в 'FitMethod':

Если вы задаете 'lsce' или 'lsrf', затем fn - вектор с mnum элементы, независимо от размера frf. Если в системе более mnum колебательные режимы, затем 'lsrf' метод возвращает первое mnum наименее демпфированные режимы сортировки в порядке увеличения собственной частоты.
Если вы задаете 'pp', затем fn - массив размера mnum-by - m -by - n с одной оценкой fn и одну оценку dr в frf.

`dr` - Коэффициенты затухания
матрица | трехмерный массив

Коэффициенты затухания для собственных частот в fn, возвращенный как матрица или трехмерный массив того же размера, что и fn.

`ms` - Векторы формы режима
матрица

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

`ofrf` - Восстановленные функции частотной характеристики
вектор | матрица | трехмерный массив

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

Алгоритмы