modalfit

Модальные параметры, полученные из частотных характеристик.

Библиотека

EngeeDSP

Синтаксис

Вызов функции

fn,dr,ms,ofrf = modalfit(frf,f,fs,mnum,out=:data) — оценивает собственные частоты мод mnum системы с измеренными частотными характеристиками frf, определенными на частотах f, для частоты дискретизации fs. Используйте modalfrf для построения матрицы частотных характеристик на основе измеренных данных. Предполагается, что матрица frf представлена в формате динамической гибкости (восприимчивости).

Также возвращает коэффициенты демпфирования dr, векторы форм колебаний ms, соответствующие каждой собственной частоте в fn, и массив реконструированных частотных характеристик ofrf на основе оцененных модальных параметров.

fn,dr,ms,ofrf = modalfit(___,Name,Value,out=:data) — задает дополнительные параметры с помощью аргументов типа «имя-значение».

modalfit(___,out=:plot) — строит графики частотных характеристик.

Аргументы

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

# frf — частотные характеристики
вектор | матрица | массив

Details

Частотные характеристики, заданные как вектор, матрица или массив. Аргумент frf имеет размер на на , где — количество частотных интервалов, — количество сигналов отклика, а — количество сигналов возбуждения.

Используйте modalfrf для построения матрицы частотных характеристик на основе измеренных данных.

Типы данных	`Int64`, `Float32`, `Float64`
Поддержка комплексных чисел	Да

# f — частоты
вектор

Details

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

Типы данных

Int64, Float32, Float64

# fs — частота дискретизации данных измерений
положительный скаляр

Details

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

Типы данных

Int64, Float32, Float64

# mnum — количество мод
положительное целое число

Details

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

Типы данных

Int64, Float32, Float64

Входные аргументы «имя-значение»

Укажите необязательные пары аргументов в виде Name,Value, где Name — имя аргумента, а Value — соответствующее значение. Аргументы типа «имя-значение» должны располагаться после других аргументов, но порядок пар не имеет значения. Можно указать несколько пар «имя-значение».

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

Пример: modalfit(frf,f,fs,mnum,"FitMethod","pp","FreqRange",[0 500]) или modalfit(frf,f,fs,mnum, FitMethod="pp",FreqRange=[0 500]) использует метод выделения пиков для выполнения аппроксимации и ограничивает диапазон частот от 0 до 500 Гц.

# Feedthrough — наличие пропускания в оцениваемой передаточной функции
false (по умолчанию) | true

Details

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

Типы данных

Bool

# FitMethod — метод аппроксимации
"lsce" (по умолчанию) | "lsrf" | "pp"

Details

Метод аппроксимации, заданный как "lsce", "lsrf" или "pp".

"lsce" — Метод наименьших квадратов для комплексных экспонент. Если указано значение "lsce", то fn представляет собой вектор с mnum элементами, не зависящими от размера frf.
"lsrf" — оценка рациональной функции по методу наименьших квадратов. Если указано "lsrf", то fn представляет собой вектор с mnum элементами, не зависящими от размера frf. Алгоритм описан [3]. Этот алгоритм, как правило, требует меньшего объема данных, чем непараметрические подходы, и является единственным, который работает для неравномерных f.
"pp" — Метод выделения пиков. Для frf, вычисленного на основе сигналов возбуждения и сигналов отклика, fn представляет собой массив размером mnum на на , содержащий одну оценку fn и одну оценку dr на frf.

# FreqRange — частотный диапазон
двухэлементный вектор

Details

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

Типы данных

Int64, Float32, Float64

# PhysFreq — собственные частоты физических мод
вектор

Details

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

Типы данных

Int64, Float32, Float64

# DriveIndex — индексы функции частотных характеристик в точке возбуждения
[1 1] (по умолчанию) | двухэлементный вектор из положительных целых чисел

Details

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

Пример: "DriveIndex",[2 3] указывает, что частотная характеристика в точке возбуждения равна frf(:,2,3).

Типы данных

Int64, Float32, Float64

# out — тип выходных данных
:plot (по умолчанию) | :data | :all

Details

Тип выходных данных:

:plot — функция возвращает график;
:data — функция возвращает данные;
:all — функция возвращает данные и график.

Для этого аргумента имя и значение разделяются знаком равенства (=).

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

# fn — собственные частоты
матрица | массив

Details

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

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

# dr — коэффициенты демпфирования
матрица | массив

Details

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

# ms — векторы форм колебаний
матрица

Details

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

# ofrf — реконструированные частотные характеристики
вектор | матрица | массив

Details

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

Примеры

Частотная характеристика одноканальной системы

Details

Оценим частотную характеристику простой одноканальной (single-input/single-output, SISO) системы и сравним ее с определением.

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

modalfit 1

Сгенерируем 3000 временных отсчетов. Определим интервал дискретизации .

import EngeeDSP.Functions: modalfit, randn, modalfrf, hann, ss2tf, freqz

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

Систему можно описать с помощью модели пространства состояний.

где — вектор состояния, и — соответственно перемещение и скорость массы, — движущая сила, а — измеренный выходной сигнал. Матрицы пространства состояний:

где — единичная матрица , а матрицы пространства состояний в непрерывном времени имеют вид:

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

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

C = [1 0]
D = 0

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

using Random

Random.seed!(1234)
u = randn(1,N) / 2
u[2001:end] .= 0

y = zeros(N)
x = [0.0; 0.0]

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


t = dt * (0:N-1)
plot(t, y, grid=true)

modalfrf 2

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

wind = hann(Int64.(N/2));
u = vec(u)
frf,f,coh = modalfrf(u',y',Fs,wind,Sensor="dis",out=:data)

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

b,a = ss2tf(A,B,C,D);
ztf,fz = freqz(b,a,2048,Fs);

using Plots

frf_siso = vec(frf[:, 1, 1])
ztf_siso = vec(ztf[:, 1, 1])

plot(f, 20 * log10.(abs.(frf_siso) .+ 1e-12), label="FRF", grid=true)
plot!(fz * Fs, 20 * log10.(abs.(ztf_siso) .+ 1e-12), label="ZTF")
ylims!((-60, 40))

modalfit 3

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

fn,dr,ms,ofrf = modalfit(frf,f,Fs,1,FitMethod="PP",out=:data)

println("fn - ",fn[1])
println("dr - ",dr[1])

fn - 0.15912870854643102
dr - 0.0054516425536509875

Сравним собственную частоту с , что является теоретическим значением для незатухающей системы.

theo = 1/(2*pi)

0.15915494309189535

Алгоритмы

Метод наименьших квадратов для комплексных экспонент

Details

Метод наименьших квадратов для комплексных экспонент вычисляет импульсную характеристику, соответствующую каждой частотной характеристике, и аппроксимирует отклик набором комплексных затухающих синусоид с использованием метода Прони. Дискретизированная затухающая синусоида может быть представлена в виде

где

— частота дискретизации;
— частота синусоиды;
— коэффициент демпфирования;
и — амплитуда и фаза синусоиды.

Члены называются амплитудами, а — полюсами. Метод Прони выражает дискретную функцию в виде суперпозиции мод (и, следовательно, амплитуд и полюсов):

Полюсы являются корнями многочлена с коэффициентами , , …, :

Коэффициенты находят с помощью авторегрессионной модели вида отсчетов из :

Для нахождения полюсов алгоритм использует функцию roots. Как только полюсы известны, можно определить частоты и коэффициенты демпфирования, вычислив мнимую и вещественную части логарифмов полюсов. Заключительным шагом является нахождение амплитуд и реконструкция импульсной характеристики с помощью

Следующая простая реализация отражает суть процедуры:

import EngeeDSP.Functions: hankel, roots

N = 4
L = 2 * N
h = rand(L)
c = -hankel(h[1:N], h[L-N+1:L]) \ h[N+1:L]
x = roots([1; reverse(c)])
V = ComplexF64[x[j]^(i-1) for i in 1:L, j in 1:N]
hrec = V * (V \ h[1:L])
result = sum(h - hrec)

0.006229773214838086 + 9.05157097882403e-17im

Систему также можно построить так, чтобы она содержала отсчеты из нескольких частотных характеристик, и решить ее методом наименьших квадратов.

Метод выделения пиков

Details

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

где

— частотная характеристика;
— частота резонанса без демпфирования;
— относительное демпфирование;
— постоянная демпфирования;
— упругая постоянная;
— масса.

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

Литература

Allemang, Randall J., and David L. Brown. Experimental Modal Analysis and Dynamic Component Synthesis, Vol. III: Modal Parameter Estimation. Technical Report AFWAL-TR-87-3069. Air Force Wright Aeronautical Laboratories, Wright-Patterson Air Force Base, OH, December 1987.
Brandt, Anders. Noise and Vibration Analysis: Signal Analysis and Experimental Procedures. Chichester, UK: John Wiley & Sons, 2011.
Ozdemir, Ahmet Arda, and Suat Gumussoy. Transfer Function Estimation in System Identification Toolbox via Vector Fitting. Proceedings of the 20th World Congress of the International Federation of Automatic Control, Toulouse, France, July 2017.