(57)【特許請求の範囲】
【請求項１】
制御出力のサンプリング周期が、制御入力のそれの偶数P倍である入力多重マルチレート系の、制御入力側のサンプリング周期における制御対象の周波数応答を同定する方法であって、
擬似白色信号の一種であるM系列信号の周期Mpに基づき、データ取得長Mp×P−１分のM系列信号を生成するM系列信号生成ステップと、
前記生成されたM系列信号と、前記生成されたM系列信号を入力して前記制御対象から得られる出力データに基づき、前記制御対象のインパルス応答推定値を計算するインパルス応答推定ステップと、
前記インパルス応答推定値を離散フーリエ変換することによって、前記制御対象の周波数応答を同定する周波数応答同定ステップと
を備えた方法。
【請求項２】
前記インパルス応答推定ステップは、
u[k]を、前記M系列信号における時刻kの信号、y[k]を、u[k]を入力したときの前記制御対象の出力データ、σu2をu[k]の分散としたとき、
【数１】
0005566969.tif 000046
により、前記インパルス応答推定値を計算する、
ことを特徴とする請求項１に記載の方法。
【請求項３】
制御出力のサンプリング周期が、制御入力のそれの偶数P倍である入力多重マルチレート系の、制御入力側のサンプリング周期における制御対象の周波数応答を同定する装置であって、
擬似白色信号の一種であるM系列信号の周期Mpに基づき、データ取得長Mp×P−１分のM系列信号を生成するM系列信号生成部と、
前記生成されたM系列信号と、前記生成されたM系列信号を入力して前記制御対象から得られる出力データ基づき、前記制御対象のインパルス応答推定値を計算するインパルス応答推定部と、
前記インパルス応答推定値を離散フーリエ変換することによって、前記制御対象の周波数応答を同定する周波数応答同定部と
を備えた装置。
【請求項４】
前記インパルス応答推定部は、
u[k]を、前記M系列信号における時刻kの信号、y[k]を、u[k]を入力したときの前記制御対象の出力データ、σu2をu[k]の分散としたとき、
【数２】
0005566969.tif 000047
により、前記インパルス応答推定値を計算する、
ことを特徴とする請求項３に記載の装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明の実施形態は、マルチレート系の高速周波数応答同定法および高速周波数応答同定装置に関し、たとえば制御出力のサンプリング周期が制御入力のサンプリング周期の偶数倍となる入力多重マルチレート系において、制御入力側のサンプリング周期における制御対象の周波数応答を高速に同定するアルゴリズムに関する。
【背景技術】
【０００２】
磁気ディスク装置のヘッド位置決め制御系などに用いられている入力多重マルチレート系においては、制御系設計時に制御入力側のサンプリング周期における制御対象の周波数応答（伝達関数）を同定することが必要となる場合がある。入力多重マルチレート系では、制御出力（観測出力）が間引き（ダウンサンプリング）されるため、一般的な周波数応答の同定方法はそのままでは適用出来ない。従って、入力多重マルチレート系の特徴を考慮した同定アルゴリズムが必要となる。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】F.Ding and T.Chen, “Identification of dual-rate systems based on finite impulse response models”, International Journal of Adaptive Control and Signal Processing, 18-7, pp.589-598, (2004)
【発明の概要】
【発明が解決しようとする課題】
【０００４】
本発明の一側面は、制御出力のサンプリング周期が制御入力のそれの偶数倍である入力多重マルチレート系において、制御入力側のサンプリング周期における制御対象の周波数応答を同定することを目的とする。
【課題を解決するための手段】
【０００５】
本発明の一態様としてのマルチレート系の高速周波数応答同定法は、制御出力のサンプリング周期が、制御入力のそれの偶数P倍である入力多重マルチレート系の、制御入力側のサンプリング周期における制御対象の周波数応答を同定する方法であって、M系列信号生成ステップと、インパルス応答推定ステップと、周波数応答同定ステップとを備える。
【０００６】
前記M系列信号生成ステップは、擬似白色信号の一種であるM系列信号の周期Mpに基づき、データ取得長Mp×P−１分のM系列信号を生成する。
【０００７】
前記インパルス応答推定ステップは、前記生成されたM系列信号と、前記生成されたM系列信号を入力して前記制御対象から得られる出力データに基づき、前記制御対象のインパルス応答推定値を計算する。
【０００８】
前記周波数応答同定ステップは、前記インパルス応答推定値を離散フーリエ変換することによって、前記制御対象の周波数応答を同定する。
【図面の簡単な説明】
【０００９】
【図１】本発明の実施形態の同定アルゴリズムにおけるデータ取得長Lを示す図。
【図２】本発明の実施形態の実行手順を示すフローチャート。
【図３】計算機シミュレーションにおける入出力関係を示したブロック線図。
【図４】本発明の実施形態の数式(15)を満たした場合の周波数応答同定結果を示す図。
【図５】本発明の実施形態の数式(15)を満たさない場合の周波数応答同定結果を示す図。
【図６】本発明の実施形態に係わる装置の機能ブロック図。
【発明を実施するための形態】
【００１０】
本発明の実施形態に係る同定アルゴリズムは、同定入力として擬似白色信号の一種であるM系列信号、同定モデルとしてFIR（Finite Impulse Response）モデルを用いる。また、マルチレート比は偶数P倍、すなわち、制御出力のサンプリング周期が制御入力のそれの2, 4, 6…倍であることが条件となる。
【００１１】
ここで、入力多重マルチレート系の伝達関数は、以下数式(1)で与えられることが知られている。
【数１】
0005566969.tif 000002
数式(1)において、Y(z)は出力の伝達関数、U(z)は入力の伝達関数、G0(z)は入力側のサンプリング周期における制御対象の伝達関数、Pはマルチレート比である。本実施形態で同定すべき制御対象の伝達関数は、数式(1)におけるG0(z)となる。数式(1)は、例えばマルチレート比P = 2のとき、
【数２】
0005566969.tif 000003
となる。以下のでは、マルチレート比P = 2を用いることを仮定して、本実施形態の同定アルゴリズムを説明する。
【００１２】
まず、本実施形態ではG0(z) の同定モデル
【数３】
0005566969.tif 000004
を、未知パラメータ
【数４】
0005566969.tif 000005
を用いて、以下数式(3)のFIRモデルとする。
【数５】
0005566969.tif 000006
また、入力の時系列をu[0], u[1], …, u[m]とすれば、数式(2)のU(z)は数式(3)と同様に、以下数式(4)のFIRモデルで表すことが出来る。
【数６】
0005566969.tif 000007
さらに、
【数７】
0005566969.tif 000008
となる。
【００１３】
ここで、数式(2)のG0(z)を
【数８】
0005566969.tif 000009
に置き換えた場合の出力、すなわちモデル出力の各時刻k = 0, 1, 2,…における時系列
【数９】
0005566969.tif 000010
は、数式(3)-(5)より以下のように単純な畳み込みによって計算することが出来る。
【数１０】
0005566969.tif 000011
【数１１】
0005566969.tif 000012
【数１２】
0005566969.tif 000013
【数１３】
0005566969.tif 000014
【数１４】
0005566969.tif 000015
（以下同様）
以上より、モデル出力の時系列
【数１５】
0005566969.tif 000016
は、kが奇数のとき出力が間引き（ダウンサンプリング）された応答になることが分かる。
【００１４】
同定モデル
【数１６】
0005566969.tif 000017
の未知パラメータ
【数１７】
0005566969.tif 000018
を推定するには、実際のハードウェア（例えば磁気ディスク装置のヘッド位置決め制御系）に入力信号u[0], u[1], …, u[m]を印加した場合の出力の時系列y[k]と、モデル出力の時系列
【数１８】
0005566969.tif 000019
との誤差e[k]の二乗和を最小化する問題（最小二乗法）を解けば良い。マルチレート比P = 2の入力多重マルチレート系では、kが奇数のときは
【数１９】
0005566969.tif 000020
とも間引きされ0になるのでe[k] (k = 0, 1, 2,…) は、以下数式(6)のようになる。
【数２０】
0005566969.tif 000021
従って、解くべき最小二乗法は、
【数２１】
0005566969.tif 000022
となる。数式(7)の解は、未知パラメータ
【数２２】
0005566969.tif 000023
と入力u[0], u[1], u[2],…から成る以下数式(8)のベクトル、
【数２３】
0005566969.tif 000024
を用いて以下数式(9)で与えられる。
【数２４】
0005566969.tif 000025
【００１５】
次に、数式(9)を実際に計算することを考える。本実施形態ではu[k]に擬似白色信号の一種であるM系列を用いるので、数式(9)の総和回数Nを、M系列の周期Mpより、
【数２５】
0005566969.tif 000026
と選ぶと、数式(9)右辺の1番目の括弧内、
【数２６】
0005566969.tif 000027
は、
【数２７】
0005566969.tif 000028
となる。ここで、M系列の周期性を考慮すると数式(12)は、
【数２８】
0005566969.tif 000029
と変形出来、M系列1周期分の自己相関行列となることが分かる。従って、白色信号の自己相関の性質より数式(13)の非対角項は、以下数式(14)のように近似的に0とすることが出来る。
【数２９】
0005566969.tif 000030
【００１６】
数式(14)のσu2はu[k]の分散である。以上より、数式(9)右辺の逆行列演算は単純な除算に置き換えることが出来、計算負荷を大幅に低減することが出来る。すなわち、計算負荷の大きい逆行列演算を行う必要はない。数式(14)を成立させるために必要最小限の同定実験時の取得データ長Lは、
【数３０】
0005566969.tif 000031
となる。同定実験時の取得データ長の概念を図１に示す。図１では、マルチレート比P = 2、M系列の周期Mp=7の場合のデータ取得長が示される。
【００１７】
以上より、マルチレート比P = 2の場合の未知パラメータ推定アルゴリズムは、
【数３１】
0005566969.tif 000032
となる。一般的な偶数倍マルチレート比（P = 2, 4, 6, …）の場合は、
【数３２】
0005566969.tif 000033
となる。
【００１８】
以上の数式(15),(17)が本実施形態の基本的なアルゴリズムである。また、数式(17)によって計算される
【数３３】
0005566969.tif 000034
はG0(z)のインパルス応答推定値（FIRモデルの分子多項式係数）であるため、
【数３４】
0005566969.tif 000035
を離散フーリエ変換すればG0(z)の周波数応答を求めることが出来る。離散フーリエ変換はFFT（Fast Fourier Transform）を用いて容易に実行可能である。
【００１９】
以上示した本実施形態のアルゴリズムは、同定入力信号をM系列信号と規定しているため、一般的に正弦波掃引法などと比較して同定実験に要する時間を短縮することが出来る。また、数式(17)は単純な加算及び乗除算で構成されているため、本実施形態を適用するハードウェアに搭載されるマイクロプロセッサや外部の計算機に、簡易な記述で実装することが出来る。
【００２０】
本実施形態の実行手順を示したフローチャートは図２となる。この実行手順に従った処理は、ハードウェアに搭載されるマイクロプロセッサや、外部の計算機で実行される。
【００２１】
ステップS1において、取得したデータ長Lが、Mp×P−１に達したかを判定する。
【００２２】
Mp×P−１に達していないときは、次の時刻のM系列信号を生成し（S2）、同定対象（例：磁気ディスク装置）にM系列信号を入力する（S3）。なお、M系列信号の生成アルゴリズムは事前に与えておく。
【００２３】
制御対象（同定対象）の出力を観測し（S4）、M系列と出力データを、内部メモリまたは外部記憶装置等の記憶部に保存する（S5）。データ長Lが、Mp×P−1に達するまで、以上のステップS1〜S4を繰り返す。
【００２４】
データ長LがMp×P-1に達したら、すなわち、データ長L分の信号および出力データを取得したら、制御対象のインパルス応答推定値を計算する。具体的に、制御対象の同定モデルであるＦＩＲモデルのインパルス応答推定値
【数３５】
0005566969.tif 000036
を、数式(17) により計算する（S6）。そして、
【数３６】
0005566969.tif 000037
を離散フーリエ変換する（S7）。離散フーリエ変換の結果が、求めるべき同定対象の周波数応答となる。
【００２５】
図６に、本実施形態に係わる装置の機能ブロック図を示す。
【００２６】
M系列信号生成部１１は、擬似白色信号の一種であるM系列信号を、データ取得長L＝Mp×P−１分、生成する。
【００２７】
生成したM系列信号の各時刻の信号がそれぞれ同定対象１２に入力され、出力データが制御対象（同定対象）１２から出力される。
【００２８】
記憶部１５は、生成部１１で生成されたM系列信号と、各時刻の信号に対応して同定対象１２から得られた出力データを記憶する。
【００２９】
インパルス応答推定部１３は、記憶部１５に記憶されたM系列信号および出力データに基づき、制御対象１２のインパルス応答推定値を計算する。この計算は、数式（A）により行う。
【数３７】
0005566969.tif 000038
【００３０】
数式(A)は、前述した数式(17)と同じものである。前述した通り、u[k]は、時刻ｋ（k番目）で取得したM系列信号の成分、y[k]は、u[k]を同定対象１２の入力としたときに得られる出力データ、σu2はu[k]の分散である。
【００３１】
フーリエ変換部１４は、
【数３８】
0005566969.tif 000039
を離散フーリエ変換することによって、制御対象の周波数応答を同定する。
【００３２】
以上、本実施形態によれば、入力多重マルチレート系において、制御入力側のサンプリング周期における制御対象の周波数応答を、高速かつ簡易な演算によって同定することが可能となる。
【００３３】
以下、本発明の実施例を計算機シミュレーションによって示す。
【００３４】
まず、複数の共振モードを持つ機械システムを想定した、以下数式(18)の離散時間系伝達関数を、計算機シミュレーションにおける同定対象とする。
【数３９】
0005566969.tif 000040
【００３５】
サンプリング周期は50[μs]とする。この同定対象にサンプリング周期50[μs]でM系列信号を入力し、マルチレート比P = 2でダウンサンプリングされた出力（サンプリング周期100[μs]）を計算機シミュレーションによって得る。このときの同定対象の入出力関係を示したブロック線図を図３に示す。M系列のシフトレジスタ段数は13段（周期Mp = 8191）とした。また、数式(17)で推定するインパルス応答
【数４０】
0005566969.tif 000041
の長さは1024とした。これは、数式(8)においてn = 1023、すなわち、
【数４１】
0005566969.tif 000042
とすることに相当する。最後に、数式(17)を実行し得られたインパルス応答推定値
【数４２】
0005566969.tif 000043
に対して1024点のFFTを実行し同定対象の周波数応答を得る。ここで、入出力の取得データ長を本実施形態(15)式に従い、
【数４３】
0005566969.tif 000044
とした場合と、数式(15)式に従わず任意に設定した値、
【数４４】
0005566969.tif 000045
とした場合とで数式(17)の計算を行い、取得データ長が同定精度に与える影響を調べる。
【００３６】
以上の設定の下で計算機シミュレーションを実行した結果（同定対象の周波数応答）を図４および図５に示す。図４および図５において、実線が同定された周波数応答、点線が同定対象の真の周波数応答である。また、図中5000[Hz]に表示されている破線はサンプリング周期100[μs]時のナイキスト周波数を示している。本実施形態の同定手法と異なり、マルチレート系の特性を陽に考慮しない同定手法では、ナイキスト周波数以上（破線右側）の周波数応答を同定することは困難である。
【００３７】
図４より、本実施形態の数式(15)を満たす取得データ長Lを設定した場合は、同定された周波数応答が真の周波数応答と良く一致しており、ナイキスト周波数以上の帯域においても正確に同定出来ていることが分かる。一方、図５より、数式(15)を満たさない任意の取得データ長Lを設定した場合、周波数応答の大まかな一致は見られるものの、数式(15)を満たす場合に比べて、真の周波数応答との誤差が大きくなっていることが分かる。以上の計算機シミュレーションによる実施例により、本実施形態(15),(17)式の効果が検証された。
【００３８】
なお、磁気ディスク装置のヘッド位置決め制御系のように、ナイキスト周波数以上に共振モードを有する制御対象に対して、本実施形態のアルゴリズムによって得られる制御入力側のサンプリング周期における周波数応答は、共振モードを安定化するために用いられるマルチレートノッチフィルターの設計に利用することが出来る。
【図１】
0005566969.tif 000048
【図２】
0005566969.tif 000049
【図３】
0005566969.tif 000050
【図４】
0005566969.tif 000051
【図５】
0005566969.tif 000052
【図６】
0005566969.tif 000053