modalfrf

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

Библиотека

EngeeDSP

Синтаксис

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

frf,f,coh = modalfrf(x,y,fs,window,out=:data) — вычисляет матрицу частотных передаточных функций frf на основе сигналов возбуждения x и сигналов отклика y, частота дискретизации которых равна fs. Выходной аргумент frf представляет собой оценку , вычисленную с использованием метода Уэлча с окном window для оконного преобразования сигналов. Аргументы x и y должны иметь одинаковое количество строк. Если x или y представляют собой матрицы, то каждый столбец представляет собой один сигнал.

Предполагается, что отклик системы y содержит измерения ускорения. Для вычисления частотных характеристик на основе измерений перемещения или скорости используйте аргумент Sensor. Функция modalfrf всегда выводит частотные характеристики в формате динамической гибкости (восприимчивости) независимо от типа датчика.

Также функция возвращает вектор частот f, соответствующий каждой частотной характеристике, и матрицу множественной когерентности coh.

frf,f,coh = modalfrf(x,y,fs,window,noverlap,out=:data) — использует количество отсчетов перекрытия между смежными сегментами noverlap.

frf,f,coh = modalfrf(___,Name,Value,out=:data) — задает дополнительные параметры с помощью аргументов типа «имя-значение» для любого из предыдущих вариантов синтаксиса.

modalfrf(___,out=:plot) — строит графики частотных характеристик. Графики ограничены первыми четырьмя сигналами возбуждения и четырьмя сигналами отклика.

Аргументы

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

# x — сигналы возбуждения
вектор | матрица

Details

Сигналы возбуждения, заданные как вектор или матрица.

Типы данных

Int64, Float32, Float64

# y — сигналы отклика
вектор | матрица

Details

Сигналы отклика, заданные как вектор или матрица.

Типы данных

Int64, Float32, Float64

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

Details

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

Типы данных

Int64, Float32, Float64

# window — окно
целое число | вектор

Details

Окно, заданное как целое число или вектор. Используйте аргумент window для разбиения сигнала на сегменты:

Если window — целое число, то modalfrf разделяет x и y на сегменты длиной window и применяет к каждому сегменту прямоугольное окно этой длины.
Если window — вектор, то modalfrf разделяет x и y на сегменты, длина которых равна длине вектора window, и применяет к каждому сегменту окно window.

Если длину x и y невозможно точно разделить на целое число сегментов с noverlap перекрывающимися отсчетами, то сигналы соответствующим образом усекаются.

Пример: hann(N+1) или (1-cos(2*pi*(0:N)'/N))/2 задают окно Ханна длиной N+1.

Типы данных

Int64, Float32, Float64

# noverlap — количество перекрывающихся отсчетов
0 (по умолчанию) | положительное целое число

Details

Количество перекрывающихся отсчетов, заданное как положительное целое число.

Если window — скаляр, то noverlap должно быть меньше window.
Если window — вектор, то noverlap должно быть меньше длины window.

Типы данных

Int64, Float32, Float64

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

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

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

Пример: modalfrf(x,y,fs,window, "Sensor","vel","Estimator","H1") или modalfrf(x,y,fs,window, Sensor="vel", Estimator="H1") указывает, что входной сигнал состоит из измерений скорости и в качестве оценки выбрано "H1".

# Estimator — оценка
"H1" (по умолчанию) | "H2" | "Hv"

Details

Оценка, заданная как "H1", "H2" или "Hv":

Используйте "H1", если шум не коррелирует с сигналами возбуждения.
Используйте "H2", если шум не коррелирует с сигналами отклика. В этом случае количество сигналов возбуждения должно равняться количеству сигналов отклика.
Используйте "Hv" для минимизации расхождения между смоделированными и оцененными данными отклика путем минимизации следа матрицы ошибок. Величина — это геометрическое среднее значений и : . Измерение должно быть одноканальным (single-input/single-output, SISO).

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

Details

Наличие пропускания в модели пространства состояний, заданное как логическое значение.

Типы данных

Bool

# Measurement — конфигурация измерений
"fixed" (по умолчанию) | "rovinginput" | "rovingoutput"

Details

Конфигурация измерений для равного числа каналов возбуждения и отклика, заданная как "fixed", "rovinginput" или "rovingoutput":

Используйте "fixed", если источники возбуждения и датчики расположены в фиксированных точках системы. Каждое возбуждение влияет на каждый отклик.
Используйте "rovinginput", если измерения получены в результате испытания с подвижным возбуждением (или подвижным молотком). Один датчик находится в фиксированном месте системы. Один источник возбуждения размещается в нескольких точках и генерирует один сигнал от датчика для каждой точки. Выход функции frf(:,:,i) = modalfrf(x(:,i),y(:,i)).
Используйте "rovingoutput", если измерения получены в результате испытания с подвижным датчиком. Один источник возбуждения находится в фиксированном месте системы. Один датчик размещается в нескольких точках и реагирует на одно возбуждение в каждой точке. Выход функции frf(:,i) = modalfrf(x(:,i),y(:,i)).

# Order — порядок модели пространства состояний
1:10 (по умолчанию) | целое число | вектор-строка целых чисел

Details

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

Типы данных

Int64, Float32, Float64

# Sensor — тип датчика
"acc" (по умолчанию) | "dis" | "vel"

Details

Тип датчика, заданный как "acc", "vel" или "dis".

"acc" — указывает, что сигнал отклика системы пропорционален ускорению;
"vel" — указывает, что сигнал отклика системы пропорционален скорости;
"dis" — указывает, что сигнал отклика системы пропорционален перемещению.

Функция modalfrf всегда выводит частотную характеристику в формате динамической гибкости (восприимчивости) независимо от типа датчика.

Пример использования этого аргумента см. в разделе Незатухающий гармонический осциллятор.

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

Details

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

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

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

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

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

Details

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

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

Details

Частоты, возвращаемые в виде вектора.

# coh — матрица множественной когерентности
матрица

Details

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

Примеры

Незатухающий гармонический осциллятор

Details

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

где числитель зависит от измеряемой величины:

перемещение: ;
скорость: ;
ускорение: ;

Вычислим частотную характеристику для трех возможных типов датчиков отклика системы. Используем частоту дискретизации 2 Гц и 30 000 отсчетов белого шума в качестве входных данных.

import EngeeDSP.Functions: modalfrf, filter,hann
fs = 2
dt = 1/fs
N = 30000
u = randn(N,1)
ydis = filter((1-cos(dt))*[0 1 1],[1 -2*cos(dt) 1],u)
w=hann(Int64.(N/2))

frfd,fd,coh = modalfrf(u,ydis,fs,w,Sensor="dis",out=:data)

yvel = filter(sin(dt)*[0 1 -1],[1 -2*cos(dt) 1],u)
frfv,fv,coh = modalfrf(u,yvel,fs,hann(Int64.(N/2)),Sensor="vel",out=:data)

yacc = filter([1 -(1+cos(dt)) cos(dt)],[1 -2*cos(dt) 1],u)
frfa,fa,coh = modalfrf(u,yacc,fs,hann(Int64.(N/2)),Sensor="acc",out=:data)

frfd_vec = frfd[:]
frfv_vec = frfv[:]
frfa_vec = frfa[:]

plot(fd, abs.(frfd_vec), xscale=:log10, yscale=:log10, label="dis", grid=true)
plot!(fv, abs.(frfv_vec), xscale=:log10, yscale=:log10, label="vel", grid=true)
plot!(fa, abs.(frfa_vec), xscale=:log10, yscale=:log10, label="acc", grid=true)

modalfrf

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

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

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

Литература

Dynamic Stiffness, Compliance, Mobility, and more Siemens, last modified 2019, https://community.sw.siemens.com/s/article/dynamic-stiffness-compliance-mobility-and-more.
Brandt, Anders. Noise and Vibration Analysis: Signal Analysis and Experimental Procedures. Chichester, UK: John Wiley & Sons, 2011.
Irvine, Tom. An Introduction to Frequency Response Functions Vibrationdata, 2000, https://vibrationdata.com/tutorials2/frf.pdf.
Vold, Havard, John Crowley, and G. Thomas Rocklin. New Ways of Estimating Frequency Response Functions. Sound and Vibration. Vol. 18, November 1984, pp. 34–38.