CN102903370B - 多速率系统的高速频率响应识别法和装置 - Google Patents

Abstract

本发明提供多速率系统的高速频率响应识别法和装置，通过FIR模型表示控制对象，基于M系列信号的周期M_p和多速率比P生成数据长度M_p×P-1量的M系列信号，基于前述M系列信号和输入前述M系列信号而从控制对象得到的输出数据，计算控制对象的脉冲响应推定值，通过对脉冲响应推定值进行离散傅立叶变换，识别前述控制对象的频率响应。

Description

多速率系统的高速频率响应识别法和装置

技术领域

本发明的实施方式涉及多速率系统的高速频率响应识别法和高速频率响应识别装置，涉及例如在控制输出的采样周期为控制输入的采样周期的偶数倍的多重输入多速率系统中，高速地识别控制对象在控制输入侧的采样周期下的频率响应的算法。

背景技术

在磁盘装置的头定位控制系统等所用的多重输入多速率系统中，有时需要在控制系统设计时识别控制对象在控制输入侧的采样周期下的频率响应(传递函数)。在多重输入多速率系统中，由于控制输出(观测输出)被减少(降采样)，所以一般的频率响应的识别方法不能够直接应用。因而，需要考虑了多重输入多速率系统的特征的识别算法。

发明内容

作为本发明的一方式的多速率系统的高速频率响应识别法，识别控制输出的采样周期为控制输入的采样周期的偶数P倍的多重输入多速率系统的、控制对象在控制输入侧的采样周期下的频率响应，包括：M系列信号生成步骤、脉冲响应推定步骤、频率响应识别步骤。

前述M系列信号生成步骤，基于作为模拟白信号的一种的M系列信号的周期Mp，生成数据取得长度Mp×P-1量的M系列信号。

前述脉冲响应推定步骤，基于前述生成的M系列信号和输入前述生成的M系列信号而从前述控制对象得到的输出数据，计算前述控制对象的脉冲响应推定值。

前述频率响应识别步骤，通过对前述脉冲响应推定值进行离散傅立叶变换，识别前述控制对象的频率响应。

附图说明

图1是表示本发明的实施方式的识别算法中的数据取得长度L的图。

图2是表示本发明的实施方式的执行步骤的流程图。

图3是表示计算机模拟中的输入输出关系的框线图。

图4是表示本发明的实施方式的满足了数学式(15)的情况下的频率响应识别结果的图。

图5是表示本发明的实施方式的不满足数学式(15)的情况下的频率响应识别结果的图。

图6是本发明的实施方式所涉及的装置的功能框图。

具体实施方式

本发明的实施方式所涉及的识别算法，作为识别输入使用作为模拟白信号的一种的M系统信号，作为识别模型使用FIR(Finite ImpulseResponse，有限脉冲响应)模型。此外，以多速率比为偶数P倍、即控制输出的采样周期为控制输入的采样周期的2、4、6…倍为条件。

在此，已知多重输入多速率系统的传递函数由以下数学式(1)提供。

$Y\left(z\right)=\frac{1}{P}\underset{n=0}{\overset{P-1}{\mathrm{\Σ}}}{G}_0\left(z{\mathrm{\φ}}^n\right)\·U\left(z{\mathrm{\φ}}^n\right),\mathrm{\φ}=\exp \left(j\frac{2\mathrm{\π}}{P}\right)\text{---}\left(1\right)$

在数学式(1)中，Y(z)为输出的传递函数，U(z)为输入的传递函数，G₀(z)为控制对象的在输入侧的采样周期下的传递函数，P为多速率比。在本实施方式中应该识别的控制对象的传递函数为数学式(1)中的G₀(z)。数学式(1)，例如在多速率比P＝2时，成为下式。

$Y\left(z\right)=\frac{1}{2}{G}_0\left(z\right)U\left(z\right)+\frac{1}{2}{G}_0(-z)U(-z)---\left(2\right)$

以下，假定使用多速率比P＝2，说明本实施方式的识别算法。

首先，在本实施方式中，使用未知参数将G₀(z)的识别模型设定为以下数学式(3)的FIR模型。

${\hat{G}}_0\left(z\right){\hat{b}}_0+{\hat{b}}_1{z}^{-1}+{\hat{b}}_2{z}^{-2}+...+{\hat{b}}_n{z}^{-n}---\left(3\right)$

此外，如果将输入的时间序列设定为u[0]、u[1]、…、u[m]，则数学式(2)的U(z)与数学式(3)同样地，能够由以下的数学式(4)的FIR模型表示。

U(z)＝u[0]+u[1]z^-1+u[2]z^-2+...+u[m]z^-m.............................(4)

进而，成为下式。

$\begin{array}{c}{\hat{G}}_0(-z)={\hat{b}}_0-{\hat{b}}_1{z}^{-1}+{\hat{b}}_2{z}^{-2}-...+{(-1)}^n{\hat{b}}_n{z}^{-n}\\ U(-z)=u\left[0\right]-u\left[1\right]z^{-1}+u\left[2\right]z^{-2}-...+{(-1)}^{m}u\left[m\right]z^{-m}\end{array}---\left(5\right)$

在此，将数学式(2)的G₀(z)置换为的情况下的输出、即模型输出的各时刻k＝0、1、2、…的时间序列能够利用数学式(3)-(5)如以下那样通过单纯的卷积来计算。

(在k＝0时)

$\begin{array}{c}\hat{y}\left[0\right]=\frac{1}{2}\left[\begin{array}{ccccc}{\hat{b}}_0& {\hat{b}}_1& {\hat{b}}_2& ...& {\hat{b}}_n\end{array}\right]\·{\left[\begin{array}{cccccccc}u\left[0\right]& 0& 0& 0& 0& 0& ...& 0\end{array}\right]}^T\\ +\frac{1}{2}\left[\begin{array}{ccccc}{\hat{b}}_0& -{\hat{b}}_1& {\hat{b}}_2& ...& {(-1)}^n{\hat{b}}_n\end{array}\right]\·{\left[\begin{array}{cccccccc}u\left[0\right]& 0& 0& 0& 0& 0& ...& 0\end{array}\right]}^T\\ ={\hat{b}}_{0}u\left[0\right]\end{array}$

(在k＝1时)

$\begin{array}{c}\hat{y}\left[1\right]=\frac{1}{2}\left[\begin{array}{ccccc}{\hat{b}}_0& {\hat{b}}_1& {\hat{b}}_2& ...& {\hat{b}}_n\end{array}\right]\·{\left[\begin{array}{cccccccc}u\left[1\right]& u\left[0\right]& 0& 0& 0& 0& ...& 0\end{array}\right]}^T\\ +\frac{1}{2}\left[\begin{array}{ccccc}{\hat{b}}_0& {-\hat{b}}_1& {\hat{b}}_2& ...& {(-1)}^n{\hat{b}}_n\end{array}\right]\·{\left[\begin{array}{cccccccc}-u\left[1\right]& u\left[0\right]& 0& 0& 0& 0& ...& 0\end{array}\right]}^T\\ =0\end{array}$

(在k＝2时)

$\begin{array}{c}\hat{y}\left[2\right]=\frac{1}{2}\left[\begin{array}{ccccc}{\hat{b}}_0& {\hat{b}}_1& {\hat{b}}_2& ...& {\hat{b}}_n\end{array}\right]\·{\left[\begin{array}{cccccccc}u\left[2\right]& u\left[1\right]& u\left[0\right]& 0& 0& 0& ...& 0\end{array}\right]}^T\\ +\frac{1}{2}\left[\begin{array}{ccccc}{\hat{b}}_0& {-\hat{b}}_1& {\hat{b}}_2& ...& {(-1)}^n{\hat{b}}_n\end{array}\right].{\left[\begin{array}{cccccccc}u\left[2\right]& -u\left[1\right]& u\left[0\right]& 0& 0& 0& ...& 0\end{array}\right]}^T\\ ={\hat{b}}_{0}u\left[2\right]+{\hat{b}}_{1}u\left[1\right]+{\hat{b}}_{2}u\left[0\right]\end{array}$

(在k＝3时)

$\begin{array}{c}\hat{y}\left[3\right]=\frac{1}{2}\left[\begin{array}{ccccc}{\hat{b}}_0& {\hat{b}}_1& {\hat{b}}_2& ...& {\hat{b}}_n\end{array}\right]\·{\left[\begin{array}{cccccccc}u\left[3\right]& u\left[2\right]& u\left[1\right]& u\left[0\right]& 0& 0& ...& 0\end{array}\right]}^T\\ +\frac{1}{2}\left[\begin{array}{ccccc}{\hat{b}}_0& {-\hat{b}}_1& {\hat{b}}_2& ...& {(-1)}^n{\hat{b}}_n\end{array}\right]\·{\left[\begin{array}{cccccccc}-u\left[3\right]& u\left[2\right]& -u\left[1\right]& u\left[0\right]& 0& 0& ...& 0\end{array}\right]}^T\\ =0\end{array}$

(在k＝4时)

$\begin{array}{c}\hat{y}\left[4\right]=\frac{1}{2}\left[\begin{array}{ccccc}{\hat{b}}_0& {\hat{b}}_1& {\hat{b}}_2& ...& {\hat{b}}_n\end{array}\right]\·{\left[\begin{array}{cccccccc}u\left[4\right]& u\left[3\right]& u\left[2\right]& u\left[1\right]& u\left[0\right]& 0& ...& 0\end{array}\right]}^T\\ +\frac{1}{2}\left[\begin{array}{ccccc}{\hat{b}}_0& {-\hat{b}}_1& {\hat{b}}_2& ...& {(-1)}^n{\hat{b}}_n\end{array}\right].{\left[\begin{array}{cccccccc}u\left[4\right]& -u\left[3\right]& u\left[2\right]& -u\left[1\right]& u\left[0\right]& 0& ...& 0\end{array}\right]}^T\\ ={\hat{b}}_{0}u\left[4\right]+{\hat{b}}_{1}u\left[3\right]+{\hat{b}}_{2}u\left[2\right]+{\hat{b}}_{3}u\left[1\right]+{\hat{b}}_{4}u\left[0\right]\end{array}$

(以下同样)。

根据以上，可以理解，模型输出的时间序列在k为奇数时成为输出减少(降采样)的响应。

为了推定识别模型的未知参数只要解决将输入信号u[0]、u[1]、…、u[m]施加于实际的硬件(例如磁盘装置的头定位控制系统)的情况下的输出的时间序列与模型输出的时间序列的误差e[k]的平方和最小化的问题(最小二乘法)即可。在多速率比P＝2的多重输入多速率系统中，在k为奇数时y[k]、都减小、成为0，所以e[k](k＝0、1、2、…)成为以下数学式(6)。

$\begin{array}{c}e\left[0\right]=y\left[0\right]-\hat{y}\left[0\right]=y\left[0\right]-{\hat{b}}_{0}u\left[0\right]\\ e\left[1\right]=y\left[1\right]-\hat{y}\left[1\right]=0\\ e\left[2\right]=y\left[2\right]-\hat{y}\left[2\right]=y\left[2\right]-\left({\hat{b}}_{0}u\right[2]+{\hat{b}}_{1}u[1]+{\hat{b}}_{0}u[0\left]\right)\\ e\left[3\right]=y\left[3\right]-\hat{y}\left[3\right]=0\\ e\left[4\right]=y\left[4\right]-\hat{y}\left[4\right]=y\left[4\right]-\left({\hat{b}}_{0}u\right[4]+{\hat{b}}_{1}u[3]+{\hat{b}}_{2}u[2]+{\hat{b}}_{3}u[1]+{\hat{b}}_{4}u[0\left]\right)\\ .\\ .\\ .\end{array}---\left(6\right)$

从而，应该解开的最小二乘法成为下式。

$\underset{l=0}{\overset{N}{\mathrm{\Σ}}}e{\left[2l\right]}^2\→\min ---\left(7\right)$

数学式(7)的解，使用包含未知参数和输入u[0]、u[1]、u[2]、…的以下数学式(8)的下述向量：

$\begin{array}{c}\widehat{\mathrm{\θ}}={\left[\begin{array}{ccccc}{\hat{b}}_0& {\hat{b}}_1& {\hat{b}}_2& ...& {\hat{b}}_n\end{array}\right]}^T\\ u\left[k\right]={\left[\begin{array}{ccccc}u\left[k\right]& u[k-1]& u[k-2]& ...& u[k-n]\end{array}\right]}^T\end{array}---\left(8\right)$

，由以下数学式(9)提供。

$\widehat{\mathrm{\θ}}={\left(\underset{l=0}{\overset{N}{\mathrm{\Σ}}}u\right[2l\left]u{\left[2l\right]}^T\right)}^{-1}\left(\underset{l=0}{\overset{N}{\mathrm{\Σ}}}u\right[2l\left]y\right[2l\left]\right)---\left(9\right)$

接着，考虑实际计算数学式(9)。在本实施方式中，对于u[k]使用作为模拟白信号的一种的M系列，所以若通过M系列的周期M_p将数学式(9)的总和次数N选择为：

N＝M_p-1....................................................(10)

，则数学式(9)右边的第一个括弧内

$\mathrm{\Λ}=\underset{l=0}{\overset{N}{\mathrm{\Σ}}}u\left[2l\right]u{\left[2l\right]}^T---\left(11\right)$

成为下式。

$\mathrm{\Λ}=u\left[0\right]u{\left[0\right]}^T+u\left[2\right]u{\left[2\right]}^T+...+u\left[2(M_p-2)\right]u{\left[2(M_p-2)\right]}^T+u\left[2(M_p-1)\right]u{\left[2(M_p-1)\right]}^T---\left(12\right)$

在此，若考虑M系列的周期性，则数学式(12)能够变形为：

$\begin{array}{c}\mathrm{\Λ}=u\left[0\right]u{\left[0\right]}^T+u\left[2\right]u{\left[2\right]}^T+...+u[M_p-1]u{[M_p-1]}^T\\ +u\left[1\right]u{\left[1\right]}^T+u\left[3\right]u{\left[3\right]}^T+...+u[M_p-2]u{[M_p-2]}^T\\ =\underset{l=0}{\overset{M_p-1}{\mathrm{\Σ}}}u\left[l\right]u{\left[l\right]}^T\end{array}---\left(13\right)$

，可以看出成为M系列1周期量的自相关矩阵。因而，根据白信号的自相关的性质，数学式(13)的非对角项能够如以下数学式(14)那样近似地设定为0。

数学式(14)的σ_u ²是u[k]的方差。根据以上，数学式(9)右边的逆矩阵运算能够置换为单纯的除法，能够大幅地降低计算负荷。也就是说，不需要进行计算负荷大的逆矩阵运算。使数学式(14)成立所需最小限的识别实验时的取得数据长度L成为：

L＝M_pP-1....................................................(15)

。将识别实验时的取得数据长度的概念示于图1。在图1中，示出多速率比P＝2、M系列的周期Mp＝7的情况下的数据取得长度。

根据以上，多速率比P＝2的情况下的未知参数推定算法成为：

$\widehat{\mathrm{\θ}}=\frac{\left(\underset{l=0}{\overset{M_p-1}{\mathrm{\Σ}}}u\right[2l\left]y\right[2l\left]\right)}{{\mathrm{\σ}}_u^2}---\left(16\right)$

。一般的偶数倍多速率比(P＝2、4、6、…)的情况成为下式。

$\widehat{\mathrm{\θ}}=\frac{\left(\underset{l=0}{\overset{M_p-1}{\mathrm{\Σ}}}u\right[\mathrm{Pl}\left]y\right[\mathrm{Pl}\left]\right)}{{\mathrm{\σ}}_u^2}---\left(17\right)$

以上的数学式(15)、(17)是本实施方式的基本的算法。此外，根据数学式(17)计算的是G₀(z)的脉冲响应推定值(FIR模型的分子多项式系数)，所以如果对进行离散傅立叶变换，则能够求出G₀(z)的频率响应。离散傅立叶变换能够使用FFT(Fast Fourier Transform，快速傅立叶变换)容易地执行。

以上所示的本实施方式的算法，由于将识别输入信号规定为M系列信号，所以一般与正弦波扫频法比较能够缩短识别实验所需的时间。此外，由于数学式(17)由单纯的加法及乘除法构成，所以能够以简单的描述安装于应用本实施方式的硬件所搭载的微处理器和/或外部的计算机。

表示本实施方式的执行步骤的流程图为图2。按照该执行步骤的处理，由硬件所搭载的微处理器和/或外部的计算机执行。

在步骤S1，判定取得的数据长度L是否达到了Mp×P-1。

在未达到Mp×P-1时，生成下一时刻的M系列信号(S2)，将M系列信号输入于识别对象(例如磁盘装置)(S3)。另外，M系列信号的生成算法预先提供。

观测控制对象(识别对象)的输出(S4)，将M系列和输出数据保存于内部存储器或外部存储装置等存储部(S5)。反复以上的步骤S1～S4，直至数据长度L达到Mp×P-1为止。

若数据长度L达到了Mp×P-1，即若取得了数据长度L量的M系列信号及输出数据，则计算控制对象的脉冲响应推定值。具体地，根据数学式(17)计算控制对象的识别模型、即FIR模型的脉冲响应推定值(S6)。并且，对进行离散傅立叶变换(S7)。离散傅立叶变换的结果，为应该求取的识别对象的频率响应。

在图6中示出本实施方式所涉及的装置的功能框图。

M系列信号生成部11按数据取得长度L＝Mp×P-1的量，生成作为模拟白信号的一种的M系列信号。

将所生成的M系列信号的各时刻的信号分别输入于识别对象12，并从控制对象(识别对象)12输出输出数据。

存储部15存储由生成部11生成的M系列信号和与各时刻的信号对应地从识别对象12得到的输出数据。

脉冲响应推定部13基于存储于存储部15的M系列信号及输出数据，计算控制对象12的脉冲响应推定值。该计算根据数学式(A)进行。

脉冲响应推定值 $\widehat{\mathrm{\θ}}=\frac{\left(\underset{l=0}{\overset{M_p-1}{\mathrm{\Σ}}}u\right[\mathrm{Pl}\left]y\right[\mathrm{Pl}\left]\right)}{{\mathrm{\σ}}_u^2}---\left(A\right)$

数学式(A)与前述的数学式(17)相同。如前所述，u[k]是在时刻k(第k个)取得的M系列信号的分量，y[k]是在将u[k]设定为识别对象12的输入时得到的输出数据，是u[k]的方差。

傅立叶变换部14通过对进行离散傅立叶变换，识别控制对象的频率响应。

以上，根据本实施方式，在多重输入多速率系统中，能够通过高速且简单的运算识别控制对象在控制输入侧的采样周期下的频率响应。

以下，通过计算机模拟表示本发明的实施例。

首先，将假定了具有多个谐振模式的机械系统的、以下数学式(18)的离散时间序列传递函数设定为计算机模拟中的识别对象。

$P\left(z\right)=\frac{\begin{array}{c}-0.009539z^9+0.08884z^8+0.1452z^7+0.3556z^6+0.4401z^5\\ +0.4295z^4+0.3348z^3+0.1316z^2+0.07566z-0.005744\end{array}}{\begin{array}{c}{z}^{10}-0.4028z^9+1.818z^8-2.063z^7+1.624z^6-2.951z^5\\ +1.453z^4-1.915z^3+1.352z^2-0.3971z+0.6743\end{array}}---\left(18\right)$

采样周期设定为50[μs]。对该识别对象以采样周期50[μs]输入M系列信号，通过计算机模拟得到以多速率比P＝2降采样后的输出(采样周期100[μs])。在图3中示出表示此时的识别对象的输入输出关系的框线图。M系列的移位寄存器阶数设定为13阶(周期M_p＝8191)。此外，由数学式(17)推定的脉冲响应的长度设定为1024。这相当于在数学式(8)中设定为n＝1023，即： $\begin{array}{c}\widehat{\mathrm{\θ}}={\left[\begin{array}{ccccc}{\hat{b}}_0& {\hat{b}}_1& {\hat{b}}_2& ...& {\hat{b}}_{1023}\end{array}\right]}^T\\ u\left[k\right]={\left[\begin{array}{ccccc}u\left[k\right]& u[k-1]& u[k-2]& ...& u[k-1023]\end{array}\right]}^T\end{array}---\left(19\right)$

。最后，对执行数学式(17)而得到的脉冲响应推定值执行1024点的FFT，得到识别对象的频率响应。在此，在根据本实施方式的数学式(15)将输入输出的取得数据长度设定为

$\begin{array}{c}L=M_{p}P-1\\ =8191\×2-1\\ =16381\end{array}---\left(20\right)$

的情况下和在不根据数学式(15)而将输入输出的取得数据长度设定为任意设定的值

L＝13000.....................................................(21)

的情况下进行数学式(17)的计算，检查取得数据长度对识别精度产生的影响。

将在以上的设定下执行计算机模拟的结果(识别对象的频率响应)示于图4及图5。在图4及图5中，实线是所识别的频率响应，点线是识别对象的真实的频率响应。此外，图中在5000[Hz]处示出的虚线表示采样周期100[Hz]时的奈奎斯特频率。与本实施方式的识别方法不同，在不正面考虑多速率系统的特性的识别方法中，难以识别奈奎斯特频率以上(虚线右侧)的频率响应。

从图4可以看出，在设定了满足本实施方式的数学式(15)的取得数据长度L的情况下，所识别的频率响应与真实的频率响应高度一致，在奈奎斯特频率以上的频带也能够正确地识别。另一方面，从图5可以看出，在设定了不满足数学式(15)的任意的取得数据长度L的情况下，虽然可以看出频率响应的大致的一致，但是与满足数学式(15)的情况相比，与真实的频率响应的误差较大。通过以上的计算机模拟所实现的实施例，验证了本实施方式(15)、(17)的效果。

另外，对于如磁盘装置的头定位控制系统那样在奈奎斯特频率以上具有谐振模式的控制对象，通过本实施方式的算法得到的控制输入侧的采样周期下的频率响应能够用于为了将谐振模式稳定化而使用的多速率陷波滤波器的设计。

Claims (4)

1.一种识别控制对象的频率响应的方法，其识别控制输出的采样周期为控制输入的采样周期的偶数P倍的多重输入多速率系统的、控制对象在控制输入侧的采样周期下的频率响应，包括：

M系列信号生成步骤，基于作为模拟白信号的一种的M系列信号的周期Mp，生成数据取得长度Mp×P-1量的M系列信号；

脉冲响应推定步骤，基于前述生成的M系列信号和输入前述生成的M系列信号而从前述控制对象得到的输出数据，计算前述控制对象的脉冲响应推定值；以及

频率响应识别步骤，通过对前述脉冲响应推定值进行离散傅立叶变换，识别前述控制对象的频率响应。

2.根据权利要求1所述的方法，其中，

前述脉冲响应推定步骤，在将u[k]设定为前述M系列信号中的时刻k的信号、将y[k]设定为输入了u[k]时的前述控制对象的输出数据、将σ_u ²设定为u[k]的方差时，根据下式计算前述脉冲响应推定值：

$\widehat{\mathrm{\θ}}=\frac{\left(\underset{l=0}{\overset{M_p-1}{\mathrm{\Σ}}}u\right[\mathrm{Pl}\left]y\right[\mathrm{Pl}\left]\right)}{{\mathrm{\σ}}_u^2}.$

3.一种识别控制对象的频率响应的装置，其识别控制输出的采样周期为控制输入的采样周期的偶数P倍的多重输入多速率系统的、控制对象在控制输入侧的采样周期下的频率响应，具备：

M系列信号生成部，其基于作为模拟白信号的一种的M系列信号的周期Mp，生成数据取得长度Mp×P-1量的M系列信号；

脉冲响应推定部，其基于前述生成的M系列信号和输入前述生成的M系列信号而从前述控制对象得到的输出数据，计算前述控制对象的脉冲响应推定值；以及

频率响应识别部，其通过对前述脉冲响应推定值进行离散傅立叶变换，识别前述控制对象的频率响应。

4.根据权利要求3所述的装置，其中，

前述脉冲响应推定部，在将u[k]设定为前述M系列信号中的时刻k的信号、将y[k]设定为输入了u[k]时的前述控制对象的输出数据、将σ_u ²设定为u[k]的方差时，根据下式计算前述脉冲响应推定值：

$\widehat{\mathrm{\θ}}=\frac{\left(\underset{l=0}{\overset{M_p-1}{\mathrm{\Σ}}}u\right[\mathrm{Pl}\left]y\right[\mathrm{Pl}\left]\right)}{{\mathrm{\σ}}_u^2}.$