# tfestimate
传递函数估计
函数库: TySignalProcessing
提示
推荐使用 ty_tfestimate 以获取更好的性能表现和使用体验。
# 语法
txy = tfestimate(x, y)
txy = tfestimate(x, y, window)
txy = tfestimate(x, y, window, noverlap)
txy = tfestimate(x, y, window, noverlap, nfft)
txy = tfestimate(___,"mimo")
txy, w = tfestimate(___; nargout = 2)
txy, f = tfestimate(___, fs; nargout = 2)
txy, w = tfestimate(x, y, window, noverlap, w; nargout = 2)
txy, f = tfestimate(x, y, window, noverlap, f, fs; nargout = 2)
___ = tfestimate(x, y, ___, freqrange, nargout=nargout)
___ = tfestimate(___, Estimator = est, nargout=nargout)
tfestimate(___; plot=true)
# 说明
nargout 表示输出参数的个数。当只输出一个参数时可以省略。
txy = tfestimate(x, y) 找到在一组频率下评估的输入信号 x 和输出信号 y 之间的传递函数估计。
- 如果x和y都是向量，那么它们必须具有相同的长度；
- 如果其中一个信号是矩阵，另一个是向量，那么向量的长度必须等于矩阵中的行数。该函数展开向量并返回逐列传递函数估计的矩阵；
- 如果 x 和 y 是具有相同行数但不同列数的矩阵，则 txy 是组合所有输入和输出信号的多输入/多输出（MIMO）传递函数。txy 是一个三维数组。如果 x 有 m 列，y 有 n 列，那么 txy 有 n 个列和 m 个页面。有关详细信息，请参见传递函数；
- 如果 x 和 y 是大小相等的矩阵，则 tfestimate 按列运算：txy[:,n], = tfestimate(x[:,n],y[:,n])。要获得 MIMO 估计，请将 "MIMO" 附加到参数列表中。
txy = tfestimate(x, y, window) 使用 window 将 x 和 y 划分为段并执行加窗。
txy = tfestimate(x, y, window, noverlap) 使用相邻段之间重叠的 noverlap 样本。
txy = tfestimate(x, y, window, noverlap, nfft) 使用 nfft 采样点来计算离散傅立叶变换。
txy = tfestimate(___, "mimo") 计算矩阵输入的 mimo 传递函数。此语法可以包括以前语法中输入参数的任何组合。
txy, w = tfestimate(___; nargout = 2) 返回归一化频率 w 的向量，在该向量处估计传递函数。
txy, f = tfestimate(___, fs; nargout = 2) 返回频率向量 f，该频率向量用估计传递函数的采样率 fs 表示。fs 必须是 tfestimate 的第六个数字输入。要输入采样率并仍然使用前面可选参数的默认值，请将这些参数指定为空[]。
txy, w = tfestimate(x, y, window, noverlap, w; nargout = 2) 返回 w 中指定的归一化频率下的传递函数估计。
txy, f = tfestimate(x, y, window, noverlap, f, fs; nargout = 2) 返回 f 中指定频率下的传递函数估计。
__ = tfestimate(x, y, __, freqrange, nargout=nargout) 返回 freqrange 指定的频率范围内的传递函数估计。频率范围的有效选项为 "onesided", "twosided" 和 "centered"。
___ = tfestimate(___, Estimator = est, nargout=nargout) 使用估计器 est 估计传递函数。est 的有效选项为 "H1" 和 "H2"。
tfestimate(___; plotfig = rue) 在当前图形窗口中绘制传递函数估计。
# 示例
两个序列之间的传递函数
计算并绘制两个序列 x 和 y 之间的传递函数估计。序列 x 由高斯白噪声组成。y 是用归一化截止频率为 0.2 πrad/sample 的 30 阶低通滤波器对 x 进行滤波的结果。使用矩形窗口设计滤波器。为传递函数估计指定 500 Hz 的采样率和长度为 1024 的汉明窗口。
using TyPlot
using TySignalProcessing
using TyControlSystems
using TyMath
rng = MT19937ar(1234)
h = fir1(30, 0.2, rectwin(31))
x = randn(rng, 16384, 1)
y, _ = filter1(h, 1, x)
fs = 500
tfestimate(x, y, 1024, [], [], fs; plotfig=true)
通过返回变量中的传递函数估计值并以分贝为单位绘制其绝对值，可以获得相同的结果。
figure()
Txy, f = tfestimate(x, y, 1024, [], [], fs; nargout=2)
plot(f, mag2db(abs.(Txy)))
SISO系统传递函数
估计简单的单输入/单输出系统的传递函数，并将其与定义进行比较。
一维离散时间振荡系统由单位质量 m 组成，通过单位弹性常数的弹簧连接到墙上。一个传感器在
生成 2000 个时间样本。定义采样间隔
using TyPlot
using TySignalProcessing
using TyMath
using TyBase
Fs = 1
dt = 1 / Fs
N = 2000
t = dt * (0:(N - 1))
b = 0.01
系统可以用状态空间模型来描述
其中
I 是 2×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(size(A))) * Bc
C = [-1 -b]
D = 1
在一半的测量间隔内，质量由随机输入驱动。使用状态空间模型来计算系统从全零初始状态开始的时间演变。将质量的加速度绘制为时间的函数。
rng = MT19937ar(5489)
u = zeros(1, N)
u[1:Int(N / 2)] = randn(rng, 1, Int(N / 2))
y = 0
x = [0; 0]
y = Matrix{Float64}(undef, 1, 2000)
for k in 1:N
global x
y[k] = (C * x .+ D * u[k])[1]
x = A * x + B * u[k]
end
plot(t, y)
将系统的传递函数估计为频率的函数。使用 2048 个 DFT 点，并指定形状因子为 15 的 Kaiser 窗口。使用相邻线段之间重叠的默认值。
nfs = 2048
wind = kaiser(N, 15)
txy, ft = tfestimate(u, y, wind, [], nfs, Fs; nargout=2)
离散时间系统的频率响应函数可以表示为系统时域传递函数的 Z 变换，在单位圆上进行评估。验证 tfestimate 计算的估计值是否与此定义一致。
plot(ft, 20 * log10.(abs.(txy)))
使用 tfestimate 的内置功能绘制估算。
figure()
tfestimate(u, y, wind, [], nfs, Fs; plotfig=true)
MIMO系统的传递函数
估计一个简单的多输入/多输出系统的传递函数。
理想的一维振荡系统由两个质量组成，
生成 30000 个时间样本，相当于 600 秒。定义采样间隔
using TyPlot
using TyBase
using TySignalProcessing
using TyMath
Fs = 50
dt = 1 / Fs
N = 30000
t = dt * (0:(N - 1))
系统可以用状态空间模型来描述
其中
I 是 4×4 恒等式，并且连续时间状态空间矩阵是
设置 k=400, b=0, 和 μ=1/10。
k = 400
b = 0
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 = [-2*k -2*b k b; k/m b/m -2 * k/m -2 * b/m];
D = [1 0; 0 1/m]
在整个测量过程中，质量由随机输入驱动。使用状态空间模型来计算系统从全零初始状态开始的时间演变。
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
使用输入和输出数据来估计作为频率函数的系统传递函数。指定 "mimo" 选项以生成所有四个传递函数。使用 5000 个样本 Hann 窗口将信号划分为多个段。指定 2500 个相邻线段之间的重叠样本和
wind = hann(5000)
nov = 2500
q, fq = tfestimate(u', y', wind, nov, 2^14, Fs, "mimo"; nargout=2)
通过设置 b=0.1 为系统添加阻尼。计算具有相同驱动力的阻尼系统的时间演化。计算
b = 0.1
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)
B = Ac \ (A - eye(4)) * Bc
C = [-2*k -2*b k b; k/m b/m -2 * k/m -2 * b/m]
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
tfestimate(u', y', wind, nov, [], Fs, "mimo", "Estimator", "H2"; plotfig=true)
legend(["I1, O1", "I1, O2", "I2, O1", "I2, O2"]; loc="northeast")
# 输入参数
x - 输入信号向量 | 矩阵
输入信号，指定为向量或矩阵。
示例： cos.(pi/4*(0:159))+randn(160) 指定加入高斯白噪声中的正弦曲线。
复数支持： 是
y - 输出信号向量 | 矩阵
输出信号，指定为向量或矩阵。
复数支持： 是
window - 窗口整数 | 向量 | []
窗口，指定为整数或向量。使用窗口将信号划分为多个段。
- 如果 window 是一个整数，那么 tfestimate 将 x 和 y 划分为长度为 window 的分段，并用该长度的 Hammin g窗对每个分段进行开窗；
- 如果 window 是一个向量，那么 tfestimate 将 x 和 y 划分为与向量长度相同的片段，并使用 window 对每个片段进行开窗。
如果 x 和 y 的长度不能精确地划分为具有 noverlap 重叠样本的整数个分段，则相应地截断信号。
若指定 window 为空，则 tfestimate 使用 Hamming 窗口，使得 x 和 y 被划分为具有 noverlap 重叠样本的八个片段。
示例： hann(N+1) 和 (1-cos.(2pi(0:N)'/N))/2 都指定长度为 N+1 的 hann 窗口。
noverlap - 重叠样本数正整数 | []
重叠样本数，指定为正整数。
- 若 window 是标量，那么 noverlap 必须小于 window；
- 若窗口是一个向量，那么 noverlap 必须小于窗口的长度。
如果将 noverlap 指定为空，那么 tfestimate 将使用一个在段之间产生 50% 重叠的数字。
nfft - DFT点数正整数 | []
DFT点数，指定为正整数。如果您将 nfft 指定为空，则 tfestimate 将此参数设置为
fs - 采样率正标量
采样率，指定为正标量。采样率是指每单位时间的采样数。如果时间单位是秒，那么采样率的单位是赫兹。
w - 归一化频率向量
归一化频率，指定为至少包含两个元素的行或列向量。归一化频率以 rad/sample 为单位。
示例： w=[pi/4 pi/2]
f - 频率向量
频率，指定为至少包含两个元素的行或列向量。频率以单位时间的周期为单位。单位时间由采样率 fs 指定。如果 fs 以采样/秒为单位，则 f 以 Hz 为单位。
示例： fs=1000；f=[100 200]
freqrange - 传递函数估计的频率范围"onesided" | "twosided" | "centered"
传递函数估计的频率范围，指定为 "onesided"、"twosided" 或 "centered" 中的一个。对于实值信号，默认为 "onesided"，对于复值信号，默认为 "twosided"。
"onesided" - 返回两个实值输入信号 x 和 y 之间传递函数的单边估计。如果 nfft 是偶数，则 txy 有 nfft/2+1 行，并且是在区间 [0, π] rad/sample 上计算的。如果 nfft 是奇数，则 txy 有 (nfft+1)/2 行，间隔为 [0, π) rad/sample。如果指定 fs，则偶数 nfft 的对应间隔为 [0, fs/2] cycles/unit 时间，奇数 nfft 的相应间隔为 [0，fs/2) cycles/unit 时间；
"twosided" - 返回两个实值或复值输入信号 x 和 y 之间传递函数的双侧估计。在这种情况下，txy 有 nfft 行，并且是在 [0, 2π) rad/sample 的间隔上计算的。如果指定 fs，则间隔为 [0, fs) cycles/unit 时间；
"centered" - 返回两个实值或复值输入信号 x 和 y 之间传递函数的居中双侧估计。在这种情况下，txy 具有 nfft 行，并且在区间(–π, π] rad/sample（偶数 nfft）和 (–π, π) rad/sample（奇数 nfft）上计算。如果指定 fs，则偶数 nfft 的相应间隔为 (–fs/2, fs/2] cycles/unit 时间，奇数 nfft 的对应间隔为 (-fs/2, fs/2) cycles/unit 时间。
est-传递函数估计器"H1"（默认值）| "H2"
传递函数估计器，指定为"H1"或"H2"。
- 当噪声与输入信号不相关时，使用 "H1"。
- 当噪声与输出信号不相关时，使用 "H2"。在这种情况下，输入信号的数量必须等于输出信号的数量。
有关详细信息，请参见传递函数。
# 输出参数
txy - 传递函数估计向量 | 矩阵 | 三维阵列
传递函数估计，以向量、矩阵或三维数组形式返回。
w-归一化频率向量
归一化频率，作为实值向量返回。
f - 循环频率向量
循环频率，作为实值向量返回。
# 详细信息
传递函数
输入 x 和输出 y 之间的关系由线性时不变传递函数 txy 建模。在频域中，
对于单输入/单输出系统，传递函数的
估计由下式给出 其中
是 x 和 y 的互功率谱密度， 是 x 的功率谱密度。该估计假设噪声与系统输入不相关。 对于多输入/多输出（MIMO）系统，
估计器变为 对于m个输入和n个输出，其中：
是第 k 个输入和第 i 个输出的互功率谱密度； 是第 k 个和第 i 个输入的互功率谱密度。
对于两个输入和两个输出，估计器是矩阵
对于单输入/单输出系统，传递函数的
估计由下式给出 其中
是 y 的功率谱密度，并且 是 x 和 y 的互功率谱密度的复共轭。该估计假设噪声与系统输出不相关。 对于 MIMO 系统，
估计器仅对相等数量的输入和输出定义明确：n=m。估计器变为 其中：
是第 k 个和第 i 个输出的交叉功率谱密度； 是第 i 个输入和第 k 个输出的交叉功率谱密度的复共轭。
# 算法
tfestimate 使用 Welch 的平均周期图方法。详见 pwelch。
# 参考文献
[1] 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.