spa

使用频谱分析，以固定的频率分辨率估计频率响应

全页折叠

语法

G = spa(data)

G = spa(data,winSize,freq)

G = spa(data,winSize,freq,maxSize)

sys = spa(___,Name,Value)

说明

G = spa(data) 根据时域或频域数据，估计频率响应及其不确定度，并估计噪声频谱 data。

data 包含时域或频域的输入/输出数据，或时间序列数据。data 可以是 timetable 的形式，也可以是逗号分隔的数值矩阵对、iddata 对象，或者可以包含复数的 idfrd 对象。

如果 data 是一个时间序列（即不包含输入），则 spa(data) 将返回输出频谱及其不确定度。spa 使用汉窗，在 0（不包括 0）到 π 之间 128 个等间隔频率值处计算频谱。

如果 data 采用时间表的形式，软件会将最后一个变量解释为唯一的输出变量。要更改此解释，请使用 InputName 和 OutputName 名称参量。

示例

G = spa(data,winSize,freq) 使用大小为 winSize 的汉窗，对 freq 中包含的频率进行频率响应估计。

示例

G = spa(data,winSize,freq,maxSize) 将输入/输出数据划分为多个片段，每个片段包含少于 maxSize 个元素。使用此语法可提高计算性能。

sys = spa(___,Name,Value) 使用由一个或多个名称-值参量指定的附加模型选项。

例如，使用 sys = spa(data,'InputName',["u1","u3"],'OutputName',["y1","y4"]) 指定输入和输出信号变量名称。

此语法可与上述任一输入参量组合结合使用。

示例

全部折叠

估计频率响应

打开实时脚本

估计时间表 tt3 中输入/输出数据的频率响应。使用默认的固定分辨率，即在 0（不包括 0）和 $π$ 之间设置 128 个等间隔的对数频率值。

load sdata3 tt3; 
g = spa(tt3); 
bode(g)

MATLAB figure

估计指定频率下的频率响应

打开实时脚本

生成对数间隔向量 f。

f = logspace(-2,pi,128);

估计矩阵 umat3 和 ymat3 中输入/输出数据的频率响应。将窗大小指定为 []，即可获得默认的滞后窗大小。

load sdata3 umat3 ymat3;
g = spa(umat3,ymat3,[],f);

绘制波德响应曲线和频谱，并标出 3 个标准差的置信区间。

h = bodeplot(g);
showConfidence(h,3)

MATLAB figure

figure
h = spectrumplot(g);
showConfidence(h,3)

MATLAB figure

输入参数

全部折叠

`data` — 估计数据
时间表 | 数值矩阵对 | `iddata` 对象 | `idfrd` 对象 | 时间表的元胞数组 | 由矩阵组成的元胞数组

均匀采样估计数据，可指定为时间表、以逗号分隔的矩阵对或数据对象，具体如下文所述。

时间表

将 data 指定为使用等间隔时间向量的 timetable。data 包含表示输入和输出通道的变量。

在利用时间表进行估计时，您可以使用所有变量，也可以指定要使用的子集通道。

对于多试验数据，将数据指定为时间表的 Ne×1 元胞数组，其中 Ne 是试验的数量。所有试验的采样时间必须匹配。

要选择用于估计的特定输入和输出通道，请使用 InputName 和 OutputName 名称-值参量。

对于时间序列数据，可以指定 data 作为仅包含一个输出变量的单变量时间表，或者使用 InputName 和 OutputName 名称-值参量来指定仅对输出变量进行估计。

例如，sys = spa(data,'OutputName',"y2",'InputName',[]) 使用输出变量 y2（不使用任何输入变量）根据时间表 data 估计时间序列模型 sys。

以逗号分隔的矩阵对

将 data 指定为由 u、y 组成的以逗号分隔的矩阵对，其中包含均匀采样的输入和输出时域值。基于矩阵的数据不包含采样时间信息。该软件假设采样时间为一秒。不建议对连续时间系统使用基于矩阵的数据。

对于 SISO 系统，请将 u、y 指定为长度为 Ns 的列向量，其中 Ns 为采样点数。

对于 MIMO 系统，请将 u、y 指定为一组矩阵对，其维度如下：

u - N_s×N_u，其中 N_u 表示输入的数量
y - N_s×N_y，其中 N_y 表示输出数量

对于多试验数据，将 u、y 指定为一对 1×N_e 元胞数组，其中 N_e 为试验数。

对于时间序列系统，请指定空的 u，即 []、y。

数据对象

一个估计数据对象，指定为时域或频域的 iddata 对象，或包含均匀采样输入和输出值的 idfrd 对象。数据对象可以有一个或多个输出通道，以及零个或多个输入通道。软件默认将模型的采样时间设置为估计数据的采样时间。

有关使用估计数据类型的详细信息，请参阅Data Domains and Data Types in System Identification Toolbox。

`winSize` — 窗大小
`[]` (默认) | 标量整数

汉恩窗大小（也称为滞后大小），以标量整数形式指定。默认情况下，该函数将窗大小设置为 min(length(data)/10,30)。

`freq` — 频率
行向量

用于估计频谱响应的频率，指定为以 rad/TimeUnit 为单位的行向量，其中 TimeUnit 指 data 的 TimeUnit 属性。默认情况下，该函数将 freq 设置为一个包含 128 个值的向量，这些值位于 (0,π] 区间内，并按对数等间距分布。对于离散时间模型，应将 freq 设置在奈奎斯特频率范围内。

`maxSize` — 最大分段大小
`250e3` (默认) | 正整数

data 内各分段的最大大小，以正整数形式指定。如果省略此参量，该函数将使用 data 中的完整数据集进行估计，而不是对数据进行分段。

名称-值参数

全部折叠

以 Name1=Value1,...,NameN=ValueN 的形式指定可选参量对组，其中 Name 是参量名称，Value 是对应的值。名称-值参量必须出现在其他参量之后，但对各个参量对组的顺序没有要求。

如果使用的是 R2021a 之前的版本，请使用逗号分隔每个名称和值，并用引号将 Name 引起来。

示例: sys = spa(data,'InputName',"u2")

`InputName` — 输入通道名称
字符串 | 字符向量 | 字符串数组 | 字符向量元胞数组

输入通道名称，指定为字符串、字符向量、字符串数组或字符向量元胞数组。

如果您使用时间表作为数据源，则 InputName 中的名称必须是时间表变量的子集。

示例: sys = spa(tt,__,'InputName',["u1" "u2"]) 从时间表 u1 中选择变量 u2 和 tt 作为输入通道用于估计。

`OutputName` — 输出通道名称
字符串 | 字符向量 | 字符串数组 | 字符向量元胞数组

输出通道名称，指定为字符串、字符向量、字符串数组或字符向量单元元胞数组。

如果您使用时间表作为数据源，则 OutputName 中的名称必须是时间表变量的子集。

示例: sys = spa(tt,__,'OutputName',["y1" "y3"]) 从时间表 y1 中选择变量 y3 和 tt 作为输出通道用于估计。

输出参量

全部折叠

`G` — 频率响应与噪声频谱
`idfrd` 对象

包含不确定性和频谱的频率响应，指定为一个 idfrd 对象。对于时间序列数据，G 表示估计的频谱和标准差。

有关估计结果和所用选项的信息存储在模型的 Report 属性中。Report 具有以下字段。

报告字段描述

Status

模型状态摘要，指示模型是通过构造创建的还是通过估计获得的

Method

使用估计命令

windowSize

汉恩窗的大小。

DataUsed

用于估计的数据的属性。包含以下字段的结构：

字段	描述
`Name`	数据集名称。
`Type`	数据类型。
`Length`	数据采样的数量。
`Ts`	采样时间这相当于 `data.Ts`。
`InterSample`	输入采样间行为。以下值之一： `'zoh'` - 零阶保持在采样之间维持分段常数输入信号。 `'foh'` - 一阶保持在采样之间维持分段线性输入信号。 `'bl'` - 带限行为指定连续时间输入信号在奈奎斯特频率以上的功率为零。对于离散时间模型，`InterSample` 的值对估计结果没有影响。
`InputOffset`	在估计过程中，从时域输入数据中去除了偏移量。
`OutputOffset`	在估计过程中，从时域输出数据中去除了偏移量。

详细信息

全部折叠

频率响应函数

频率响应函数描述了系统对正弦输入的稳态响应。对于线性系统，特定频率的正弦输入会产生一个同样频率的正弦输出，但其振幅和相位不同。频率响应函数描述了振幅变化和相位偏移与频率的关系。

为了更好地理解频率响应函数，请考虑以下对线性动态系统的描述：

$y (t) = G (q) u (t) + v (t)$

这里，u(t) 和 y(t) 分别表示输入信号和输出信号。G(q) 被称为系统的传递函数 - 它描述了将输入转换为输出的系统动态过程。符号 G(q)u(t) 表示以下运算：

$G (q) u (t) = \sum_{k = 1}^{\infty} g (k) u (t - k)$

q 是移位运算符，由以下方程定义：

$G (q) = \sum_{k = 1}^{\infty} g (k) q^{- k} q^{- 1} u (t) = u (t - 1)$

G(q) 是频率响应函数，当其在单位圆上求值时，G(q=e^iω)。

G(q=e^iω) 与输出频谱 ${\hat{Φ}}_{v} (ω)$ 共同构成了该系统的频域描述。

采用布莱克曼-图基法估计的频率响应函数由以下方程给出：

${\hat{G}}_{N} (e^{i ω}) = \frac{{\hat{Φ}}_{y u} (ω)}{{\hat{Φ}}_{u} (ω)}$

在此情况下，^ 表示近似数值。关于该方程的推导，请参阅 [1] 中关于非参数时域和频域方法的一章。

输出噪声频谱

输出噪声频谱（v(t) 的频谱）由以下方程给出：

${\hat{Φ}}_{v} (ω) = {\hat{Φ}}_{y} (ω) - \frac{{| {\hat{Φ}}_{y u} (ω) |}^{2}}{{\hat{Φ}}_{u} (ω)}$

该频谱方程的推导基于以下假设： $y (t) = G (q) u (t) + v (t)$ 之间存在线性关系；u(t) 与 v(t) 相互独立；且各频谱之间满足以下关系：

$Φ_{y} (ω) = {| G (e^{i ω}) |}^{2} Φ_{u} (ω) + Φ_{v} (ω)$

$Φ_{y u} (ω) = G (e^{i ω}) Φ_{u} (ω)$

在此，噪声频谱由以下方程给出：

$Φ_{v} (ω) \equiv \sum_{τ = - \infty}^{\infty} R_{v} (τ) e^{- i w τ}$

${\hat{Φ}}_{y u} (ω)$ 是输出-输入交叉频谱， ${\hat{Φ}}_{u} (ω)$ 是输入频谱。

此外，扰动 v(t) 也可以描述为经过滤波的白噪声：

$v (t) = H (q) e (t)$

这里，e(t) 是方差为 $λ$ 的白噪声，其噪声频谱由以下方程给出：

$Φ_{v} (ω) = λ {| H (e^{i ω}) |}^{2}$

算法

spa 通过以下步骤应用布莱克曼-图基 (Blackman-Tukey) 谱分析法：

根据 u(t) 和 y(t) 计算协方差和交叉协方差：
$\begin{array}{l} {\hat{R}}_{y} (τ) = \frac{1}{N} \sum_{t = 1}^{N} y (t + τ) y (t) \\ {\hat{R}}_{u} (τ) = \frac{1}{N} \sum_{t = 1}^{N} u (t + τ) u (t) \\ {\hat{R}}_{y u} (τ) = \frac{1}{N} \sum_{t = 1}^{N} y (t + τ) u (t) \end{array}$
计算协方差和交叉协方差的傅里叶变换：
$\begin{array}{l} {\hat{Φ}}_{y} (ω) = \sum_{τ = - M}^{M} {\hat{R}}_{y} (τ) W_{M} (τ) e^{- i ω τ} \\ {\hat{Φ}}_{u} (ω) = \sum_{τ = - M}^{M} {\hat{R}}_{u} (τ) W_{M} (τ) e^{- i ω τ} \\ {\hat{Φ}}_{y u} (ω) = \sum_{τ = - M}^{M} {\hat{R}}_{y u} (τ) W_{M} (τ) e^{- i ω τ} \end{array}$
其中 $W_{M} (τ)$ 是宽度（滞后量）为 M 的汉窗。您可以指定 M 来控制估计值的频率分辨率，该分辨率约等于 2π/M 弧度/采样时间。
默认情况下，该操作使用 0（不包括 0）到 π 之间 128 个等间隔的频率值，其中 w = [1:128]/128*pi/Ts，且 Ts 是该数据集的采样时间。汉窗的默认滞后长度为 M = min(length(data)/10,30)。对于默认频率，该操作采用快速傅里叶变换 (FFT)，其效率高于用户自定义频率的情况。
注意
M = γ 见 [1] 的表 6.1。标准差见 [1] 第 184 页和第 188 页。
计算频率响应函数 ${\hat{G}}_{N} (e^{i ω})$ 和输出噪声频谱 ${\hat{Φ}}_{v} (ω)$ 。
${\hat{G}}_{N} (e^{i ω}) = \frac{{\hat{Φ}}_{y u} (ω)}{{\hat{Φ}}_{u} (ω)}$
$Φ_{v} (ω) \equiv \sum_{τ = - \infty}^{\infty} R_{v} (τ) e^{- i w τ}$

spectrum 是输出通道和输入通道的频谱矩阵。也就是说，如果 z = [data.OutputData、data.InputData]、spectrum 包含 z 的矩阵值频谱作为频谱数据。

$S = \sum_{m = - M}^{M} E z (t + m) z {(t)}^{'} W_{M} (T_{s}) \exp (- i ω m)$

这里，' 表示复共轭转置。

参考

[1] Ljung, Lennart. System Identification: Theory for the User. 2nd ed. Prentice Hall Information and System Sciences Series. Upper Saddle River, NJ: Prentice Hall PTR, 1999.

版本历史记录

在 R2006a 之前推出

全部展开

R2022b: 接受时间表和矩阵形式的时域估计数据

现在，大多数估计、验证、分析和实用工具函数都接受时域输入/输出数据，其形式为包含输入和输出数据的单个时间表或分别包含输入和输出数据的一对矩阵。这些函数仍然接受 iddata 对象作为数据源，无论是时域数据还是频域数据。

另请参阅

spa

语法

说明

示例

估计频率响应

估计指定频率下的频率响应

输入参数

data — 估计数据 时间表 | 数值矩阵对 | iddata 对象 | idfrd 对象 | 时间表的元胞数组 | 由矩阵组成的元胞数组

时间表

以逗号分隔的矩阵对

数据对象

winSize — 窗大小 [] (默认) | 标量整数

freq — 频率 行向量

maxSize — 最大分段大小 250e3 (默认) | 正整数

名称-值参数

InputName — 输入通道名称 字符串 | 字符向量 | 字符串数组 | 字符向量元胞数组

OutputName — 输出通道名称 字符串 | 字符向量 | 字符串数组 | 字符向量元胞数组

输出参量

G — 频率响应与噪声频谱 idfrd 对象

详细信息

频率响应函数

输出噪声频谱

算法

参考

版本历史记录

R2022b: 接受时间表和矩阵形式的时域估计数据

另请参阅

主题

`data` — 估计数据
时间表 | 数值矩阵对 | `iddata` 对象 | `idfrd` 对象 | 时间表的元胞数组 | 由矩阵组成的元胞数组

`winSize` — 窗大小
`[]` (默认) | 标量整数

`freq` — 频率
行向量

`maxSize` — 最大分段大小
`250e3` (默认) | 正整数

`InputName` — 输入通道名称
字符串 | 字符向量 | 字符串数组 | 字符向量元胞数组

`OutputName` — 输出通道名称
字符串 | 字符向量 | 字符串数组 | 字符向量元胞数组

`G` — 频率响应与噪声频谱
`idfrd` 对象