Syslab使用手册
⋮
2026a

# modalfrf

用于模态分析的频率响应函数

函数库: TySignalProcessing

# 语法

frf, = modalfrf(x, y, fs, window)
frf, = modalfrf(x, y, fs, window, noverlap)
frf, = modalfrf(___; Name=Value)
frf,f,coh = modalfrf(___)
___ = modalfrf(___; plotfig=true)

# 说明

frf, = modalfrf(x, y, fs, window) 从激励信号 x 和响应信号 y 中估计频率响应函数 FRF 的矩阵，所有这些函数都以速率 fs 采样。输出 frf 是使用 Welch 方法计算的 H1 估计值，其中信号经过 window 窗函数加窗。x 和 y 必须具有相同的行数。如果 x 或 y 是矩阵，则每列表示一个信号。 假设系统响应 y 包含加速度测量值。要从位移或速度测量值开始计算频率响应函数，请使用 Sensor 参数。modalFRF 始终以动态柔韧性（柔韧度）格式输出频率响应函数，而与传感器类型无关。

frf, = modalfrf(x, y, fs, window, noverlap) 指定相邻段之间重叠的 noverlap 样本。

frf, = modalfrf(___; Name=Value) 使用名称-值参数指定选项，允许使用先前语法中输入的任意组合。选项包括估算器、测量配置和测量系统响应的传感器类型。

frf, f, coh = modalfrf(___) 同时返回与每个频率响应函数对应的频率矢量，以及多重相干矩阵。

___ = modalfrf(___; plotfig=true) 在当前图窗中绘制频率响应函数。这些图仅限于前 4 个激发和 4 个响应。

# 示例

锤子激励的频率响应函数

可视化单输入/单输出锤式激励的频率响应函数。

加载包含以下内容的数据文档：

  • Xhammer — 输入激励信号，由定期传递的五次锤击组成；

  • Yhammer — 系统对输入的响应。Yhammer 以位移来衡量。

信号以 4 kHz 采样。绘制激励和输出信号。

using TySignalProcessing
using TyPlot
using TyBase

pkg_dir = pkgdir(TySignalProcessing)
load(pkg_dir * "/examples/VibrationAnalysis/modalfrf/modaldata.mat")
subplot(2, 1, 1)
plot(thammer, Xhammer[:])
ylabel("Force (N)")
subplot(2, 1, 2)
plot(thammer, Yhammer[:])
ylabel("Displacement (m)")
xlabel("Time (s)")

计算并显示频率响应函数。使用矩形窗口对信号进行加窗。指定窗口涵盖锤击之间的时间段。

winlen = size(Xhammer, 1)
modalfrf(Xhammer[:], Yhammer[:], fs, winlen; Sensor="dis", plotfig=true)
MIMO 频率响应函数

计算由随机噪声激励的双输入/双输出系统的频率响应函数。

加载一个包含 Xrand（输入激励信号）和 Yrand（系统响应）的数据文档。使用 5000 个样本的 Hann 窗口和相邻数据段之间 50% 的重叠来计算频率响应函数。指定输出测量值为位移。

using TySignalProcessing
using TyPlot
using TyBase

pkg_dir = pkgdir(TySignalProcessing)
load(pkg_dir * "/examples/VibrationAnalysis/modalfrf/modaldata.mat")
winlen = 5000

frf, f, coh = modalfrf(
    Xrand, Yrand, fs, hann(winlen), 0.5 * winlen; Sensor="dis", plotfig=true
)
tightlayout()
SISO 系统的频率响应函数

估计简单的单输入/单输出系统的频率响应函数，并将其与定义进行比较。

一维离散时间振荡系统由单位质量 m（以 kg 为单位）组成，通过弹性常数 k=1 N/m 的弹簧固定在壁上。传感器对 = 1 Hz 的质量位移进行采样。阻尼器通过对其施加与速度成正比的力来阻碍质量的运动，阻尼常数 b=0.01 kg/s。

using TySignalProcessing
using TyPlot
using TyMath
using TyBase

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

该系统可以用状态空间模型来描述:

式中 为状态向量， 和 分别为质量的位移和速度， 为驱动力， 为测量输出。状态空间矩阵为：

为 2 x 2 单位阵，连续时间状态空间矩阵是：

Ac = [0 1; -1 -b]
De, Ve = eigen(Ac * dt)
A = real.(Ve * diagm(exp.(De)) / Ve)
Bc = [0; 1]
B = Ac \ (A - eye(2)) * Bc

C = [1 0]
D = [0]

在前 2000 秒内，质量由随机输入驱动，然后返回静止状态。使用状态空间模型计算系统从全零初始状态开始的时间演变。绘制质量位移随时间的变化。

rng = MT19937ar(1234)
u = randn(rng, 1, N) / 2
u[2001:end] .= 0

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

plot(t, y)

估计系统的模态频率响应函数。使用 Hann 窗口的长度是测量信号的一半。指定输出为质量的位移。

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.(2im * pi * fz)
ztf = polyval(b, z) ./ polyval(a, z)
figure()
plot(f, 20 * log10.(abs.(frf)))
hold("on")
plot(fz * Fs, 20 * log10.(abs.(ztf)))
hold("off")
grid()
ylim([-60 40])
双体振荡器的频率响应函数

估计一个简单的多输入/多输出(MIMO)系统的频率响应函数和模态参数。

理想的一维振荡系统由两个质量 和 组成，它们被限制在两面壁之间。单位为 和 。每个质量都通过弹性常数为 的弹簧连接到最近的壁上。相同的另一个弹簧连接两个质量。三个阻尼器通过对质量施加与速度成正比的力来阻碍质量的运动，阻尼常数为 。传感器采样 和 ，质量的位移，采样频率 。

生成 30000 个时间样本，相当于 600 秒。定义采样间隔 。

using TySignalProcessing
using TyPlot
using TyMath
using TyBase

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

该系统可以用状态空间模型来描述：

其中 是状态矢量， 和 分别是第 个质量的位置和速度， 是输入驱动力的矢量， 是输出矢量。状态空间矩阵是：

是 4 × 4 单位阵，连续时间状态空间矩阵是：

其中 。

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]
De, Ve = eigen(Ac * dt)
A = real.(Ve * diagm(exp.(De)) / Ve)

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, 2)

在整个测量过程中，质量由随机输入驱动。使用状态空间模型计算系统从全零初始状态开始的时间演变。

rng = MT19937ar(1234)
u = randn(rng, 2, N)

x = [0; 0; 0; 0];
y = Array{Float64}(undef, 2, 30000)
for kk in 1:N
    global x
    y[:, kk] = C * x + D * u[:, kk]
    x = A * x + B * u[:, kk]
end

使用输入和输出数据来估计系统的传递函数与频率的函数关系。使用 15000 个样本的 Hann 窗口，相邻段之间有 9000 个重叠样本。指定测量的输出为位移。

wind = hann(5000)
nove = 2500
FRF, f, = modalfrf(u', y', Fs, wind, nove; Sensor="dis")

将理论传递函数计算为时域传递函数的 Z 变换，在单位圆处计算。

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

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

frf = Array{ComplexF64}(undef, 2, Int(nfs / 2), 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 in 1:2
    for kj in 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("Input $jk, Output $kj")
    end
end
tightlayout()

使用不带 output 参数的 modalfrf 语法绘制估计值。

modalfrf(u', y', Fs, wind, nove; Sensor="dis", plotfig=true)
tightlayout()

# 输入参数

x - 激励信号
向量 | 矩阵

激励信号，指定为向量或矩阵。

y - 响应信号
向量 | 矩阵

响应信号，指定为矢量或矩阵。

fs - 采样率
正标量

采样率，指定为以赫兹表示的正标量。

window - 窗函数
正整数 | 向量

Window，指定为整数或向量。使用 window 将信号划分为多个段：

  • 如果 window 是一个整数，则 modalfrf 将 x 和 y 划分为长度为 window 的段，并且 windows 每个段都有一个该长度的矩形窗口；

  • 如果 window 是一个矢量，则 modalfrf 将 x 和 y 划分为与矢量长度相同的段，并使用 window 对每个段进行窗口；

  • 如果 'Estimator' 指定为 'subspace'，则 modalfrf 忽略窗口的形状，并使用其长度来确定返回的频率响应函数中的频率点数。

如果 x 和 y 的长度无法精确划分为具有 noverlap 重叠样本的整数个段，则信号会相应地截断。

示例： hann(N+1) 和 (1-cos(2pi(0:N)'/N))/2 都能指定一个长度为 N+1 的 Hann 窗。

noverlap - 重叠样本数
0（默认）| 正整数

重叠样本数，指定为正整数。

  • 如果 window 是标量，则 noverlap 必须小于 window；

  • 如果 window 是一个向量，则 noverlap 必须小于 window 的长度。

# 名称-值对参数

以 Name1=Value1,...,NameN=ValueN 的形式指定可选的参数对，其中 Name 是参数名，Value 是相应的值。名称-值参数必须出现在其他参数之后，但参数对的顺序并不重要。

示例： Sensor="vel"，Est="H1" 指定输入信号由速度测量组成，选择的估计器是 H1。

Estimator - 估计器
"H1"（默认）| "H2" | "Hv" | "subspace"

Estimator，指定为 "H1"、"H2"、"Hv" 或 "subspace"。有关 H1 和 H2 估计器的更多信息，请参阅传递函数。

  • 当噪声与激励信号无关时，使用 "H1"；

  • 当噪声与响应信号无关时，使用 "H2"。在这种情况下，激励信号的数量必须等于响应信号的数量；

  • 使用 "Hv" 可以通过最小化误差矩阵的跟踪来最小化建模响应数据和估计响应数据之间的差异。Hv 是 和 的几何平均值：。测量必须是单输入/单输出 （SISO）；

  • 使用 "subspace" 通过状态空间模型计算频率响应函数。在这种情况下，将忽略 noverlap 参数。与非参数方法相比，此方法通常需要更少的数据。

Feedthrough - 在状态空间模型中存在馈通
false（默认） | true

状态空间模型中存在馈通，指定为逻辑值。仅当 Estimator 指定为 "subspace" 时，此参数才可用。

Measurement - 测量配置
"fixed"（默认） | "rovinginput" | "rovingoutput"

相同数量的激励和响应信道的测量配置，指定为 "fixed"、"rovinginput" 或 "rovingoutput"。

  • 当系统的固定位置有激发源和传感器时，使用"fixed"。每一次激励都有助于每一次响应;

  • 当测量结果来自流动激励（或流动锤）测试时，请使用 "rovinginput"。单个传感器保存在系统的固定位置。将单个激励源放置在多个位置，并为每个位置产生一个传感器响应。函数输出 frf[：，：，i] = modalfrf(x[：，i]，y[：，i])；

  • 当测量结果来自粗纱传感器测试时，请使用 "rovingoutput"。单个激发源保存在系统的固定位置。单个传感器放置在多个位置，并响应每个位置的一次激励。函数输出 frf[：，i] = modalfrf(x[：，i]，y[：，i]]。

Order - 状态空间模型阶数
1:10（默认） | 整数 | 整数向量

状态空间模型阶数，指定为整数或整数向量。如果您指定整数向量，则函数将从指定范围内选择最佳阶数值。仅当 'Estimator' 指定为 "subspace" 时，此参数才可用。

Sensor - 传感器类型
"acc"（默认） | "dis" | "vel"

传感器类型，指定为 "acc"、"vel" 或 "dis"。

  • "acc" — 指定系统的响应信号与加速度成正比；

  • "vel" — 指定系统的响应信号与速度成正比；

  • "dis" — 指定系统的响应信号与位移成正比。

modalfrf 始终以动态柔韧性（接受度）格式输出频率响应函数，而与传感器类型无关。

示例：无阻尼谐振子

以 的速率采样的单位质量和弹性常数的简单无阻尼谐振子的运动由传递函数描述

其中分子取决于被测量的量级：

  • 位移：

  • 速度：

  • 加速度：

计算三种可能的系统响应传感器类型的频率响应函数。使用 2 Hz 的采样率和 30,000 个白噪声样本作为输入。

using TySignalProcessing
using TyPlot
using TyMath

fs = 2
dt = 1/fs
N = 30000

rng = MT19937ar(1234)
u = randn(rng,N,1)
ydis, = filter1((1-cos(dt))*[0 1 1],[1 -2*cos(dt) 1],u)
frfd,fd, = modalfrf(u,ydis,fs,hann(N/2),Sensor="dis")

yvel, = filter1(sin(dt)*[0 1 -1],[1 -2*cos(dt) 1],u)
frfv,fv, = modalfrf(u,yvel,fs,hann(N/2),Sensor="vel")

yacc, = filter1([1 -(1+cos(dt)) cos(dt)],[1 -2*cos(dt) 1],u)
frfa,fa, = modalfrf(u,yacc,fs,hann(N/2),Sensor="acc")

loglog(fd,abs.(frfd[:]),fv,abs.(frfv[:]),fa,abs.(frfa[:]))
grid("on")
legend(["dis" "vel" "acc"],loc="best")

在所有情况下，生成的频响函数都采用与位移对应的格式。速度和加速度测量分别是位移测量的第一次和第二次导数。频率响应函数在系统固有频率附近是等效的。远离固有频率，频率响应函数也不同。

# 输出参数

frf - 频率响应函数
向量 | 矩阵 | 3维数组

频率响应函数，以向量、矩阵或3维数组的形式返回。frf 的大小为 P×M×N，其中 P 是频率区间的数量，M 是响应的数量，N 是激励信号的数量。

modalfrf 始终以动态柔韧性（柔韧度）格式输出频率响应函数，而与传感器类型无关。

f - 频率
向量

频率，以向量形式返回。

coh - 多重相干矩阵
向量

多重相干矩阵，以矩阵形式返回。coh 为每个响应信号提供一列。

# 参考文献

[1] "Dynamic Stiffness, Compliance, Mobility, and more..." Siemens, last modified 2019.

[2] Brandt, Anders. Noise and Vibration Analysis: Signal Analysis and Experimental Procedures. Chichester, UK: John Wiley & Sons, 2011.

[3] Irvine, Tom. "An Introduction to Frequency Response Functions," Vibrationdata, 2000.

[4] Vold, Håvard, John Crowley, and G. Thomas Rocklin. "New Ways of Estimating Frequency Response Functions." Sound and Vibration. Vol. 18, November 1984, pp. 34–38.