JP2599670B2 - 伝達関数の極・ゼロ点決定方法 - Google Patents

Description

【発明の詳細な説明】

【０００１】

【産業上の利用分野】本発明は線形システムの特性を把
握するために伝達関数の極とゼロ点を解析するための方
法に関する。

【０００２】

【従来技術及びその問題点】一般に、線型システムの伝
達関数は、印加された刺激信号に応答してそのシステム
がラプラス変換（ s 平面）領域上で取る挙動を記述す
るものである。特に、電気的、機械的または電気機械的
システムの設計に際し、システムを所望の用途に対して
最適化するように s 平面内でシステムの伝達関数の各
極とゼロ点の位置と大きさを知ることは重要である。シ
ステムの安定性が特に重要な場合には、極の位置を知る
ことは極めて重要である。極とゼロ点の測定値がシステ
ム設計者に役立つためには、測定を迅速かつ正確に行う
とともに、ノイズと歪による誤差を極力少なくすること
が重要である。

【０００３】線形システムの伝達関数を測定するための
各種システムが今までに多数提案され、製作されてい
る。このような従来技術のシステムは、例えば 1976 年
8 月3 日発行された米国特許第 3,973,112 号及び 197
7 年 9 月 6 日発行の米国特許第 4,047,002 号に開示
されている。このようなシステムは、極及びゼロ点の測
定値を与えてくれるが、この測定値は、解析されるシス
テムの応答が非線形であること、及び測定データにノイ
ズが含まれていることにより、固有の誤差を含んでい
た。

【０００４】

【発明が解決しようとする問題点】本発明は上記の欠点
を解決するためになされたもので、極とゼロ点の測定を
迅速に行うとともに、ノイズによって引き起こされる誤
差を補正するものである。

【０００５】

【問題点を解決するための手段及び作用】極とゼロ点の
測定はランダムなノイズまたはガウスノイズのような所
望の刺激信号を測定対象のデバイスに印加することによ
り開始される。刺激信号と応答信号は時間領域でサンプ
ルされ、高速フーリエ変換 (FFT) により周波数領域に
変換される。これにより、刺激信号の自己 (Auto-) ス
ペクトルと、刺激信号及び応答信号の重なり (Cross-)
スペクトルが測定される。ノイズとシステム応答の非線
形性に起因する極及びゼロ点測定の誤差を最小にするた
めに、刺激信号と応答信号の測定値を組み合わせ、自己
スペクトルと重なりスペクトルをグループの平均 (ense
mble average) として測定することができる。システム
の伝達関数の測定値は重なりスペクトル及び自己スペク
トルから求められ、測定データのノイズレベルは刺激信
号と応答信号のグループの平均値から推定される。な
お、測定の時間遅れは必要に応じてこれを除去できる。

【０００６】測定された伝達関数は、分子・分母が夫々
変数 s のチェビシェフ多項式 (Chebyshev polynomial
s) である有理式 (rational fraction) である推定伝達
関数により近似 (fit) される。極・ゼロ点測定の精度
を増すために、伝達関数の一定の部分を強調する重み付
け関数 (weighting function) が求められる。次に分子
・分母の多項式の係数は、測定伝達関数データの推定伝
達関数への重み付け最小二乗近似 (weighted least squ
ares fit) として求められる。また、測定された伝達関
数に対する推定伝達関数の近似品質が求められる。もし
近似が不充分な場合には、分子・分母の多項式の次数を
変えて新しい係数を求め、再び近似度を試験する。充分
な近似が達成されると、推定伝達関数の分子・分母のチ
ェビシェフ多項式を通常の多項式に変換し、求根器 (ro
ot solver) を用いて推定伝達関数の極とゼロ点を与え
る 2 つの多項式の根を見付ける。極とゼロ点はそのデ
バイスの設計者の希望に応じてこれを表示することがで
きる。

【０００７】

【発明の実施例】図１は本発明の一実施例により構成さ
れた測定装置 1 の外観図を示す。刺激信号 x(t) は第
1 チャネルのケーブル 9 を介してデバイス 3 に加えら
れ、デバイス 3 からの応答信号 y(t) は第 2 チャネル
のケーブル 11 を経て測定装置1 に導入される。ここ
で、測定対象のデバイス 3 は任意の電気的、機械的あ
るいは電気機械的システムであり、例えばサーボモータ
でよい。なお、外部変換器を使用する場合には、デバイ
ス 3 に機械的振動を加え、それに対するデバイス3 か
らの応答信号を監視して、モーダル解析を行うことが可
能である。ユーザはキーボード 5 を通して各種の指令
を測定装置 1 に入力することができ、測定の結果は CR
T のような表示器 7 に表示することができる。

【０００８】図２は測定装置 1 の電気的ブロック図で
ある。図において、アナログ信号源27 は、キーボード
5 を介してユーザの制御の下に、ディジタル信号源 39
で駆動されて刺激信号 x(t) を発生する。この刺激信号
x(t) は前述のように第１チャネルのケーブル 9 を介
してデバイス 3 に導入される。応答信号 y(t) と刺激
信号 x(t) 入力部 21, 23 を通してそれぞれ受信され、
局部発振器 37 で設定されるサンプリングレートでアナ
ログ・ディジタル変換器 29, 31 によりディジタル化さ
れる。サンプルされた信号 x(t_j) と y(t_j) に現れるエ
イリアシング誤差を減少するには、ディジタルフイルタ
33 を使用することができる。高速フーリエ変換器 (FT
T) 43 は刺激信号と応答信号を周波数領域 (frequency
domain) の信号 X(f_j) と Y(f_j) に変換するのに使用さ
れ、これらの周波数領域の信号はバス 51 を経て RAM 4
7 に記憶される。浮動小数点プロセッサ 45 は、信号 X
(f_j)と Y(f_j) の重なりスペクトル及び自己スペクトル
を求めるため、またグループ平均を取る場合には信号 Y
(f_j) 上のノイズを測定するために使用される。プロセ
ッサ 45 は、CPU 41 の制御の下に、以下に説明するよ
うにデバイス 3 の伝達関数の極とゼロ点を測定するの
に使用される。また、前述のバス 51 は、例えば二重バ
ス系から構成されるが、データと命令を伝送するのに使
用される。

【０００９】図３はデバイス 3 の伝達関数の極とゼロ
点を測定する間に測定装置 1 が行う動作の流れ図であ
る。最初のステップ 61 で、刺激信号 x(t) がデバイス
3 に加えられ、応答信号 y(t) が記録される。2 つの
信号はディジタル化され、X(f_j) 及び Y(f_j) として周
波数領域の信号に変換される。刺激信号と応答信号の各
測定を多数回繰り返して、自己スペクトル及び重なりス
ペクトルをグループ平均として測定し、デバイス 3 の
応答中のノイズと非線形性の影響を減らすようにするこ
とができる。自己スペクトル及び重なりスペクトルが測
定され、ステップ63 で、測定された伝達関数上の測定
ノイズを推定するために測定分散が求められる。

【００１０】ステップ 67 で、帰納的関係 (recurrence
relationship) を利用して、推定伝達関数 H_E(s) の分
子・分母のチェビシェフ多項式 P(s) と Q(s) を生成す
る。ステップ 65 で重み付け関数を求める。この重み付
け関数は、測定された伝達関数を推定伝達関数で近似す
るのに用いる最小二乗解析を重み付けするために使用す
る。ステップ 69 で、分子・分母の 2 つのチェビシェ
フ多項式の係数を、測定伝達関数に対して重み付け最小
二乗近似した推定伝達関数として定める。ステップ 231
で実際の近似の程度を求め、一連の近似判定基準との
合致を解析する。近似が不充分であればこれら多項式の
次数を変え、新しい係数を求める。もしも近似が充分で
あれば、両チェビシェフ多項式をステップ 77 で通常の
多項式に変換し、これらの普通の多項式の根をステップ
79 で見出す。推定伝達関数の極とゼロ点（これらの多
項式の根）は希望に応じてステップ 81 により表示され
る。

【００１１】図３に示す流れ図の各ステップについて以
下個別に詳細に説明する。ステップ61 で、デバイス 3
の実際の伝達関数が測定される。図４は、ステップ 61
を構成する個々のステップを順に示している。刺激信号
x(t) はデバイス 3 に加えられ、信号 x(t) と応答信
号 y(t) はレート f_s でサンプルされる。多数の信号 x
(t) の任意のものがデバイス 3 を刺激するのに使用さ
れる。特に、ランダムノイズはデバイスの 3 の応答中
の小さな非線形性の効果が後で除かれるようになってい
るので望ましい。サーボモータの設計のような多くの用
途では、関心のある周波数範囲は直流 (DC) から 100kH
z までのように低い。このベースバンド (baseband) の
場合には、サンプル周波数 f_s はナイキスト (Nyquist)
の関係と理想的部品を用いて 200kHz まで低くするこ
とができる。現実にはフイルタ33 は理想的でないか
ら、幾らかの誤差が発生する。これはサンプル周波数 f
_sを例えば 256kHz まで上げるか、及び／またはデータ
の周波数領域における上の部分を無視することによって
補正することができる。

【００１２】ステップ 95 で、時間領域の信号 x(t_j)
と y(t_j) は高速フーリエ変換で周波数領域の信号 X
(f_j)と Y(f_j) に変換される。時間領域における各信号
から 2048サンプルを取る場合、周波数領域では 1024
の複素データ点が信号毎に発生する。エイリアシング誤
差の影響を補償するには周波数領域中で上側 200 個程
度のデータ点を無視することが必要となることもある。
それぞれの周波数データ点は複素数であるから、実数部
と虚数部とから成り、これはまた極形式で振幅と位相で
表すことができる。なお、方程式とそこにおける記法は
今後特に指示しない限り複素変数で表してある。

【００１３】測定が行われる物理的制約のため、あるい
は変換器その他の機器が理想的でないため、測定したデ
ータに時間遅れが存在することもあり得る。これら時間
遅れは望ましくない実在しない極及びゼロ点が解析に入
り込まないようにデータから除去すべきである。除去は
周波数領域のデータに exp(i2πft) (ここでは i は-1
の平方根、t は時間遅れの推定値であり、また exp(x)
は e の x 乗を表わす）を乗ずるというような多くの方
法のいずれかによって行うことができる。この方法の詳
細はロナルド・ブレースウエル (Ronald Bracewell) 著
「フーリエ変換とその応用」の第 2 版 104ページに述
べられている。

【００１４】デバイス 3 の伝達関数の測定値 H_MD(f_j)
は図５により定義することができる。図５はデバイス 3
の出力が入力、伝達関数及び付加的なノイズの項に関
係していることを示している。従って、伝達関数 H_MD(f
_j) は次式により陰に定義することができる。

【００１５】

【数１】

【００１６】この式は上線を「グループ平均」と定義し
たとき次のように書き替えることができる。

【００１７】

【数２】

【００１８】すなわち、

【００１９】

【数３】

【００２０】なおここで、アスタリスクを上付き添字と
して付けることにより共役複素数を表わす。デバイス 3
の伝達関数は、自己スペクトル及び重なりスペクトル
を用いて次のように定義される。

【００２１】

【数４】

【００２２】ここで、\overline{N_FX^*} (ここで \overl
ine は、その直後に置かれた中括弧 {} 内の文字列に連
続した上線を付けることを表わす記法である) の項は無
視できる程小さいと仮定し、また自己スペクトル及び重
なりスペクトルは次のように定義している。

【００２３】

【数５】

【００２４】伝達関数の測定に及ぼすノイズ項の影響を
最小にするために、自己スペクトルと重なりスペクトル
は多数の刺激−応答測定値にわたるグループ平均値とし
て求められる。すなわち、何回かの測定により得られた
各周波数スペクトルの平均としての周波数スペクトルを
得るということである。このグループ中の測定値の数が
増すにつれて、非相関ノイズ項の位相角がランダムなた
め、ノイズ項の平均は 0になる。もし刺激信号としてラ
ンダムノイズを使用すると、この平均化はデバイス 3
における応答特性の非線形性を線形化するという別の効
果を生ずる。というのは、刺激信号と非線形応答信号の
間には相関関係がないからである。10回から 1,000 回
までの測定のグループを用いると良好な結果が得られる
ことは分かっており、またこのグループ中の測定の回数
はユーザが選択して差し支えない。このようなグループ
を用いると、デバイス 3 の伝達関数は下記の商として
測定することができる。

【００２５】

【数６】

【００２６】図４のステップ 107 で、式 (6) で定義し
た H_MD(f_j) がメモリに記憶される。ステップ 63 で
は、測定された伝達関数の集合したコヒーレンス (ense
mble coherence) が測定される。この結果は極及びゼロ
点の測定におけるノイズ誤差を補正するために後で使用
される。一般に、ランダム変数 X に関する分散は次式
で定義される。

【００２７】

【数７】

【００２８】分散は周知の多数の方法のいずれかによっ
ても推定することができるが、1957年3 月グッドマン
(N. R. Goodman) のニューヨーク大学博士論文「二次元
静的ガウス過程のスペクトル、コスペクトル及び直交ス
ペクトルの合同推定 (On TheJoint Estimation Of The
Spectra, Cospectrum, And Quadrature Spectrum OfA T
wo-Dimensional Stationary Gaussian Process)」の第
131 ページにある式4.92 及び 4.93 によって厳密に求
めることができる。これらの式は第二種ベッセル関数に
書き替えることができるので、これを積分して、

【００２９】

【数８】

【００３０】が得られる。ここで n は 1 より大きく、
また n はグループ (ensemble) 中の測定値の個数であ
る。式 (8) に使用しているコヒーレンス関数は、次式
で定義される。

【００３１】

【数９】

【００３２】○重み付け関数

【００３３】ユーザは制御装置 5 によりユーザ定義重
み付け関数を入力することもできる。ユーザ定義重み付
け関数によれば、ユーザは、測定精度が高いまたは低い
ことが判っている区域を考慮して、あるいはユーザにと
って特に重要な区域を考慮して、周波数軸に沿うある所
望の領域を強調することができる。ユーザ定義重み付け
関数を入力した場合には、この関数の曲線が表示装置 7
に表示され、ステップ 65 を実行する必要はなくな
る。この関数を入力しない場合には、図３のステップ 6
5 で、最小二乗近似解析に使用する重み付け関数 W(f)
を発生する。図６と図７はそれぞれステップ 65 を構成
する各種のステップを示しており、これらのステップを
実行して 2 つの重み付け関数のうちの 1 つを作り出す
ことができる。図６に示すステップで発生した重み付け
関数は、周波数軸に添った領域で H_MD(f) が急速に変化
している部分を強調する。従って、この重み付けによる
近似は極及びゼロ点の両者の回りの領域で共に良好であ
り、この重み付け関数はサーボ解析に使用するのに最も
良く適している。最小二乗近似は伝達関数の測定値 H_MD
(f_j) と伝達関数の推定値 H_E(s) の間の各周波数データ
点で評価した誤差の二乗の和を最小にすることである。
各点で評価された重み付き誤差は、j をサンプリングイ
ンデクスとして、次のように定義できる。

【００３４】

【数１０】

【００３５】ここで、H_MD は実数の周波数点について測
定され、また H_E は複素数 s 上の関数として求められ
るが、H_MD は s に一般化でき、H_MD と H_E を直接比較
対照することができる。当業者ならば、s と i2πf の
解析関数を等価であると見なす解析接続の手続によりこ
の一般化を行うことができるであろう。

【００３６】最小二乗近似解析で定義される全重み付け
誤差 ε' は

【００３７】

【数１１】

【００３８】ここで、誤差と重み付け関数は各周波数点
で評価され、合計は関心のある周波数範囲にわたって取
られる。重み付け関数 W(f) は最小二乗近似解析におい
て所望の f の領域をエンファサイズ（強調）するのに
使用される。H_E(s) は 2 つの多項式の有理式である。

【００３９】

【数１２】

【００４０】P(s) の係数が c₀, c₁, ..., c_m でありま
た Q(s) の係数が d₀, d₁, ..., d_nであると仮定すれ
ば、最小二乗近似解析は c₀, c₁, ..., c_m 及び d₀,
d₁, ...,d_n の各々に関する ε' の偏導関数を得、これ
らの偏導関数を 0 に等しいとおいて同次方程式 (homog
eneous equations) を作り出すことにより行われる。こ
の最小二乗近似解析によれば、同次方程式を連立させて
解くことにより、所望の最小の重み付き誤差 ε' を見
つけることができる。最小二乗近似解析の一般的な解説
は、たとえばラルストン (Ralston) 及びラビノウイッ
ツ (Rabinowitz)著「数値解析の第一過程」第 2 版の第
6 章に述べられている。

【００４１】図６のステップ 121 〜 137 は、図３のス
テップ 67 で概説した重み付け関数W(f) を作り出す過
程を詳細に示している。本質的に W(f) は周波数軸に沿
う各点で評価した H_MD(f) の位相導関数である。隣接す
る 3 つの周波数データ点 (x₁, x₂, x₃) について、点
x₂ での位相導関数は ph(x₃)-ph(x₁) であると定義され
る。ただし、ph(x) はデータ点 x での位相である。位
相導関数を使用すれば周波数軸上で H_MD(f) の極及びゼ
ロ点の生起する確率が最も高い領域に最大の重みを付
け、またノイズスパイクを最小に重み付けすることがで
きる。

【００４２】図４のステップ 105 〜 107 で H_MD(f) の
複素周波数領域で測ったデータ点の実数及び虚数の成分
が測定され、記憶される。ステップ 121 で、測定した
各データ点での位相導関数が求められる。複素対数を用
いて、位相導関数は次のように表すことができる。

【００４３】

【数１３】

【００４４】

【数１４】

【００４５】ただし、位相情報だけが欲しいので、

【００４６】

【数１５】

【００４７】従って、x₃x₁ ^* の位相を評価することが必
要である。ここで小角度近似を使用するのが便利であ
り、これにより計算時間を最小限にすることができる。
ph(x₃x₁ ^*) が tan^-1(ph(x₃ x₁ ^*) ) に等しいという近似
を -2 < tan^-1(ph(x₃ x₁ ^*) ) < 2の領域で使用すること
ができる。ここで、tan^-1(ph(x₃ x₁ ^*)) はラジアンで表
現されているものとする。この領域の外では tan^-1(ph
(x₃ x₁ ^*)) は +2 または-2 と近似される。小角度近似
のその他の利点は、位相の変化の速さが制限されるため
に、H_MD(f) 上のノイズがさらに小さくなることであ
る。

【００４８】ステップ 121 で得た位相導関数はノイズ
を多く含んでいる傾向があるから、ステップ 123 で導
関数を平滑化する。平滑化は各位相導関数データ点をそ
の前後夫々 7 点ずつと平均することにより容易に行う
ことができる。従って、W(f)上の各点は実際には隣接す
る 15 個の位相導関数データ点の平均であり、その結果
ノイズスパイクが減少している。前後夫々の側に 7 つ
の点を使用して平均化したが、その数は希望に応じて 5
ないし 15 まで変えても良好な結果が得られることが
わかっている。サーボの解析では、ノイズが大きく、デ
ータの減衰特性がゆるやかなので、前後夫々の側に 13
個の点を使用した。モーダル解析では、ノイズが小さ
く、またデータの減衰が大きいので、各側に夫々 7 個
の点を取れば充分であった。

【００４９】正の重み付け関数だけが欲しいのであるか
ら、ステップ 125 で絶対値を取っている。ステップ 12
7 で平方根を取っても良い。サーボ解析では、比較的滑
らかな重み付け関数が普通望ましいので、全てのデータ
を使用するとともに、平方根もよく用いられる。逆に、
モーダル解析では、典型的には比較的 Q の高い装置を
解析するので、普通は平方根を求めることはない。

【００５０】ステップ 129 で、重み付け関数に関して
別の平滑化を行ってもよい。平滑化は上述のステップ 1
23 と同じ態様で行われる。

【００５１】測定データ上の高ノイズ領域の影響を最小
にするように重み付け関数を補正することも望ましいこ
とである。図４のステップ 111 で測定され、また式
(9) で定義されたコヒーレンス関数により、H_MD(f_j) の
各測定周波数領域データ点のノイズを推定できる。ここ
で、ノイズ補正係数

【００５２】

【数１６】

【００５３】を乗ずることにより、コヒーレンス関数が
0 に近い値を取る高ノイズ領域をデエンファサイズす
る、つまりその影響を弱める。係数を 0.05 とすれば、
きれいなデータの点をノイズの多いデータ点より 20 倍
も強調することができる。位相導関数の絶対値を含んで
いる、周波数軸に沿う各点での最終的な結果は

【００５４】

【数１７】

【００５５】である。

【００５６】ステップ 133 で、ノイズを補正された重
み付け関数はオーバーフローしないように 1 に正規化
され、ステップ 135 で、完全に無視されてしまうデー
タがないように、0.02 のしきい値が導入される。従っ
て、重み付け関数は、

【００５７】

【数１８】

【００５８】と定義される。ただし、上式の右辺中の分
数式の分母の右下に付いている添字 MAX は、今考えて
いる範囲内で当該添字が付けられている式が取る値の最
大値を式の値として採用するということを表わす。この
ようにして得られた範囲内の最大値で割ることによっ
て、上述の 1 への正規化を行っている。

【００５９】最後に、ステップ 137 で、ユーザは W(f)
を修正することもできる。W(f) は、0 から +1 まで変
わるが、これは測定装置 1 の表示器 7 に周波数の関数
として表示される。ユーザが修正を行わない場合には、
測定装置 1 は周波数の最高及び最低のデータ点から開
始し、周波数範囲の中心に向かって少なくとも 0.04の
W(f) が夫々の側に見付かるまで走査していく。こうし
て見つかった２つの端点の上下(つまり外側)にあるデー
タは無視される。また、測定時間を短縮するため、また
はある所望のデータを除外してしまわないようにするた
め、ユーザは表示器 7 上にあり制御装置 5 で制御され
る 2 つのカーソルを用いて、端点を指定することもで
きる。

【００６０】図７は、ステップ 63 を構成する詳細な各
種のステップを示している。これらのステップにより上
述の 6 図に示すステップに代えて他の重み付け関数を
発生することができる。この代わりの重み付け関数はモ
ーダル解析では好ましいものである。というのは、モー
ダル解析では多くの場合、ノイズと歪がゼロ点の近くで
しばしば発生するからである。モーダル解析値はたいて
いのサーボ解析値よりも鋭いピークを伴っていることが
しばしばある。この代わりの重み付け関数はピークの位
置を非常に正確に求める傾向があり、極とゼロ点の情報
をこれから得ることができる。ステップ 141 で、その
二乗平均値に正規化された伝達関数データ H_MD(f) の二
乗値が計算され、これが記憶される。次にステップ 145
で、H_MD(f) のピークが求められる。この判定は、短い
線分を H_MD(f) の各部分に次々に当てはめていくことに
よって行われる。800 個のデータ点に対して、典型的に
は4 種の線分長 (8, 16, 32, 64) が用いられる。これ
らの線分は 2:1 の割合でオーバーラップさせ、データ
中のどのピークも無視されないようにすることが好まし
い。

【００６１】次に各データ点の振幅を関連する線分の振
幅から差し引く。この残りをしきい値と比較して、妥当
なピークが検出されたか否かを確かめる。6^1/2σ とい
う代表的なしきい値を用いれば、この残りをその点での
標準偏差と比較することができ、偽のノイズピークを排
除することができる。予め決めた個数の連続点での上述
の残りがしきい値より大きい場合には、真のピーク値が
検出されたという最終判定を下してよい。なお、4 つの
連続点を用いれば、正確なピーク検出を行い得ることが
分かっている。

【００６２】ステップ 147 で重み付け関数の放物線が
ステップ 145 で求めた各ピークの中心の回りに構成さ
れる。放物線の中心は (S+T)/2 の位置にある。ここで
S と Tは夫々ピークがしきい値を越した最外側のデータ
点の 3 点だけ向こう側にある点である。基線からの放
物線の最大の高さは 3L/{2(T-S)} である。ただし、Lは
周波数データ点の総数である。このように、放物線は幅
が高さと逆比例するように構成されるので、広すぎたり
狭すぎたりするピークは無視される。ステップ151 で、
重み付け関数は最大ピーク値で割ることにより正規化さ
れる。ステップ 153 で、H_MD(f) の元の測定データの大
きさの二乗値をステップ 141 で測定された平均二乗値
で割ることにより正規化する。なお、分散の少ない領域
を比較的強調することが望ましいこともある。その場合
には、伝達関数の大きさの二乗値を伝達関数の大きさの
二乗値と伝達関数の分散推定値の和で割った商の k 乗
を乗算する。この関数は非常に広いピークが不注意に無
視されることがないように、W(f) に使用できる。最後
に、ステップ 155 でこれら 2 つの関数を加算し、最終
的な重み付け関数 W(f) を作る。

【００６３】○基底多項式の作成 図３のステップ 67 で、伝達関数の推定値 H_E(s) の分
子及び分母多項式がラプラス変数 s を用いて作られ
る。一般に、2 つの基底多項式 (basis polynomial) を
任意の形で作ることは可能であるが、P(s) と Q(s) に
チェビシェフ基底多項式を使用すれば、明確な利点のあ
ることが分かっている。数値的悪条件(ill-conditionin
g) は最小二乗近似解析での一般的な問題であり、有限
精度プロセッサを使用するときは、大きな数値誤差を生
ずる可能性がある。というのは、非常に大きな数と非常
に小さな数の間での加算をしなければならないからであ
る。チェビシェフ多項式には -1 から +1 までの区間に
おいて、その大きさが +1 と -1 の間で変化する性質が
あるので、加算される項はすべて相対的な大きさが同じ
であり、従って、数値的悪条件は避けられる。チェビシ
ェフ多項式を使用する上での他の固有の利点は、高次の
係数が分離する (decouple) 傾向があるということであ
る。すなわち、問題とされている周波数区間で、各種の
チェビシェフ多項式はある重み付け関数を用いること
で、互に直交しており、そのため各チェビシェフ多項式
の係数は他のチェビシェフ多項式の係数から独立に、つ
まり「分離して」、定めることができる。また、他の重
み付け関数を用いた場合でも、係数間の相互作用が小さ
いという意味で、良好に分離されている。このような分
離関係は各種の直交多項式を使っても得られる。この直
交多項式はデータ依存性があり、従って、測定値ごとに
多項式を計算しなおさなければならないという欠点があ
る。たとえば、ラルストン (Ralston) 及びラビノウイ
ッツ (Rabinowitz) 共著の「数値解析の第一過程」の第
2 版の第 6.4 章の多項式の説明を参照のこと。さらに
チェビシェフ多項式の説明を参照のこと。さらにチエシ
ビシエフ多項式と正弦・余弦関数が似ているため、チェ
ビシェフの積和関係が存在し、これによってチェビシェ
フ多項式の時間がかかる乗算をはるかに簡単な和の計算
に置き換えるように公式化することができる。例えば積
和関係は、1964 年 6 月国立標準局(NBS) から発行され
た「数学関数ハンドブック」の 782 ぺージ(式 22.7.24
)に述べられている。国立標準局の方程式には印刷上の
誤りがあり、それは次のように書き直さなければならな
い。

【００６４】

【数１９】

【００６５】所望の P(s) と Q(s) のチェビシェフ多項
式は上に引用した NBS のハンドブックの 782 ページに
書かれている帰納的関係 (recurrence relationship)
に従って容易につくることができる。多項式 P(s) と Q
(s) は複素変数 s の多項式であり、実数、虚数あるい
は複素数の係数をもっている。当業者ならは、純実数の
チェビシェフ多項式に対して与えられる帰納的関係を使
用して複素チェビシェフ多項式を容易に作ることができ
る。H_MD(f_j) の時間データの測定値がベースバンドの場
合のように純実数である場合には、複素チェビシェフ多
項式の使用を容易にするためエルミート対称 (Hermitia
n symmetry) を使用することができる。モーダル解析を
行う場合には、エルミート対称の制約を捨てて基底多項
式を修正するのが好ましい。P(s) と Q(s) の基底多項
式を作るためには、2 つの多項式の次数を初めに推定す
ることが必要である。この推定は伝達関係の極とゼロ点
の数を推定することと等価である。各多項式に対して最
初に次数を 1 と推定すれば、1 つまたはもっと多くの
極を持つ多項式を近似することができる。解析対象の装
置についての情報をはじめにもっと多く利用できる場合
には、最初の推定次数をもっと高くしてよい。

【００６６】○係数の決定 図３のステップ 69 で、P(s) と Q(s) の係数が求めら
れる。図８はステップ69 を構成する各種のステップを
一層詳細に示している。式 (10) に定義した点毎の誤差
項は次のように書き直される。

【００６７】

【数２０】

【００６８】ただし、記号 f は便宜上省いてあり、ま
た添字 j は周波数軸に沿ってサンプルしたデータを使
用していることを示す。分母の多項式 Q(s) を両辺にか
けて下式を得る。

【００６９】

【数２１】

【００７０】式 (6) を用い、また新しい項 A_j と B_j
を

【００７１】

【数２２】

【００７２】と定義することにより、式 (23) を次のよ
うに簡単化して書き直すことができる。

【００７３】

【数２３】

【００７４】B_j を両辺にかけてもよい。そのようにし
た場合は次のようになる。

【００７５】

【数２４】

【００７６】最後に、この乗算を行わない場合、\tilde
{E}_j (ここで \tilde は、その直後の中括弧 {} 内の文
字列の上に、例えば次の式 (27) の左辺のようなティル
ダを付けることを表わす記法である) をつぎのように定
義する。

【００７７】

【数２５】

【００７８】図８のステップ 161 で、ユーザが予め指
定した極とゼロ点を伝達関数の推定値から除去する。こ
の操作は必須ではないが、既知の極とゼロ点をユーザが
予め指定すると未知の極とゼロ点の解析を更に正確に行
うことができる。次式を参照して、既知のゼロ点を P
(s) 多項式から除き、また既知の極を Q(s) 多項式から
除く。

【００７９】

【数２６】

【００８０】ここで z_k は予め指定したゼロ点であり、
p_k は予め指定した極である。予め指定した除去すべき
極とゼロ点を用いて式 (27) は次のように書き換えられ
る。

【００８１】

【数２７】

【００８２】W_j, U_j, A_j/B_j 及び V_j の各項は推定され
る多項式 \tilde{P}_j, \tilde{Q}_jに適用される各点毎
の重み付けの項(測定データを含む)と見なすことができ
る。簡単化のため新しい項 W_{p_j} (p に下付き添字 j を
つけたものを W の添字として付加した記号を表わす。
以下同様である。) と W_{q_j} を定義できる。

【００８３】

【数２８】

【００８４】従って、式 (32) を次のように簡単化して
書換えることができる。

【００８５】

【数２９】

【００８６】図８のステップ 163 で、測定したデー
タ、重み付け関数及び推定伝達関数の多項式を、図１及
び図２に示す測定装置 1 を用いての解析を容易にする
ため行列記法に変換する。この変換は当業者には容易に
行うことができる。係数 c_k とd_k を有する 2 組の直交
多項式 φ_k と θ_k を上に述べた帰納的関係を用いて作
り、\tilde{P}_j と \tilde{Q}_j を表すのに使用する。

【００８７】

【数３０】

【００８８】チェビシェフ多項式はデータに依存しない
(従って測定毎に求め直す必要はない)ので、φ_k と θ_k
をチェビシェフ多項式として構成することが好まし
い。これら 2 つの多項式は異なる次数 m と n をもっ
ているが、ここでは便宜上 m=n でかつ両多項式とも θ
_k と書かれると仮定する。当業者ならば、m と n が等
しくない場合に一般化することができるだろう。式 (3
5) は次のように書き直すことができる。

【００８９】

【数３１】

【００９０】行列及びベクトルの記法は図９に示す通り
である。行列 θ を調べれば、第 1列は T₀(f) から成
る。ここで T_k(f) は f の関数である k 次のチェビシ
ェフ多項式である。従って、チェビシェフ多項式の帰納
的関係に関して上に述べたように、T₀(f) = 1, T₁(f) =
f, T₂(f) = 2f²-1, 等々となる。行列 θ の各行は、
高速フーリエ変換した H_MD(f_j) のデータをサンプルし
た周波数に対応する特定のデータの測定周波数に対応す
る。上述の通り通常の多項式の代わりにチェビシェフ多
項式を使用することにより精度が落ちることがなくな
る。

【００９１】図９に示す行列及びベクトル記号の定義を
用いれば、式 (37) と (36) は次のように書き直すこと
ができる。

【００９２】

【数３２】

【００９３】また、誤差方程式 (34) は次のように書き
替えることができる。

【００９４】

【数３３】

【００９５】さて、誤差の二乗の和は行列記法で次のよ
うに書くことができる。

【００９６】

【数３４】

【００９７】ただし、肩文字 T は複素共役転置を表
す。

【００９８】図８のステップ 165 で、上述の最小二乗
近似解析がマトリックスの形で行なわれる。近似の目標
は C 及び D のベクトルの要素に関して ε を最小にす
ることである。この最小化は C 及び D の各要素に関し
て ε を微分し、その結果を0 とおいて、これを満足す
るようなベクトル C 及び D について解けば達成され
る。得られる方程式は次のようになる。

【００９９】

【数３５】

【０１００】明確化のため、上記の式 (42) と (43) を
次のように書き替える。

【０１０１】

【数３６】

【０１０２】図８のステップ 167 で、ノイズバイアス
が最小二乗近似解析から除かれる。式 (33) から、W_{q_j}
の項は A_j を含んでおり、これは式 (24) で測定の各
周波数における G_YX データとして定義されている。図
５を参照して上に説明した通り、G_YX には付加ノイズ項
N_F(f_j) が入っている。式 (45) で、W_q の大きさの二
乗は W_q ^TW_q として得られ、その結果 G_YX 上のノイズが
二乗されている。この二乗を行うプロセスは測定データ
の各周波数点でのノイズバイアスを作り出す作用をして
いる。というのは、二乗することによって一方向のバイ
アスを生ずる効果が得られるからである。先行技術の装
置で極・ゼロ点解析にノイズバイアスが存在することを
認識したものはなく、また先行技術の装置でノイズバイ
アスをなんとかして除去しようとしたものもない。最小
二乗近似解析でノイズバイアスの影響を除かないと、伝
達関係の極とゼロ点の位置を求める際にバイアスと不確
実さが入ってくる。このことは、装置の安定性が重要な
場合には決定的に重要である。なぜなら、ノイズバイア
スがあれば支配的な極の位置を s 領域の誤った側の半
面にあると間違えて判定する可能性があるからである。

【０１０３】コヒーレンス関数と分散は式 (8) と (9)
で先に定義してあり、これは測定装置 1 で行われる 3
スペクトル測定 (tri-spectra measurement)、すなわち
式 (5) で定義された３つのスペクトル G_XX, G_YY, G_YX
の測定、から求められる。分散は二乗の大きさで表した
ノイズの電力であるから、式 (46) に示す W_q ^TW_q の項
から容易に減算することができる。このようにして、各
測定周波数点毎に分散を引き去ることにより、ノイズバ
イアスが補正される。式 (45) についてこのような補正
を行うことによって下式を得る。

【０１０４】

【数３７】

【０１０５】図８のステップ 169 で、解行列が構成さ
れ、最小二乗近似解析の最小化のためのチェビシェフ多
項式の係数が求められる。方程式 (42) と (43) は方程
式 (44) と (45) のように書き直すことができる。ここ
では、もし望むならば、ノイズバイアスが上に述べたよ
うに除かれていると仮定される。新しい関数 F, H, H^T
及び G を次のように定義する。

【０１０６】

【数３８】

【０１０７】このように定義することにより、次の関係
式

【０１０８】

【数３９】

【０１０９】が成立し、図１０に示す、4 つの部分を持
つ行列 M 及び 2 つの部分を持つベクトル V の行列乗
算として示されている同次方程式系を定義する。

【０１１０】極とゼロ点の最大可能数の解析を行ってい
る場合には (m = n = 40)、D 及びC ベクトルの各々は
40 要素の長さがあり、行列 M には 6,400 の要素が含
まれている。このような多数の要素を記憶しておくこと
は効率的でなく、必要な記憶装置の大きさ及び解行列 M
の各構成部分を測定データより式 (21) 等を用いて再
度計算して構成する (regenerate) ために必要な時間を
最小にするには、先に述べたチェビシェフの積和関係を
利用することができる。このデータ圧縮を行うと、解行
列 M を再現するのに必要な要素の数は 6,400 から 320
に減少する。更に、行列 M の第 3 象限は第 2 象限の
転置に過ぎず、従って行列の解には、どちらか 1 つだ
けを完全に再構成すればよいということが分かる。他方
は一方から容易に再構成することができる。また、加え
て、行列を解く間、一時には行列の 4 つの象限のうち
の 2 つだけしか再現されている必要がない。

【０１１１】チェビシェフの積和関係は以下の通りであ
る。

【０１１２】

【数４０】

【０１１３】従って、

【０１１４】

【数４１】

【０１１５】のようなもっと具体的な関係を上記のデー
タ圧縮に使用することができる。[なお、式 (54) では
f についての和がとられているが、これは測定データの
フーリエ変換によって得られる全ての周波数データ点
{f_i}についての和をとることを意味している（たとえ
ば、式 (11) では f のインデクスとして j を用いてい
るが、式 (54) では添字 j はすでにチェビシェフ多項
式の次数を表わすために使用されているので、ここでは
インデクスとして i を使用していることに注意された
い）。その結果、行列 M の各要素が行列内の他の 2 つ
の合計結果を足したものになる。

【０１１６】図１０から、列ベクトル V は行列 M の各
行に直交しなければならないことが分かる。従って M₁V
=0, M₂V=0, ..., M_nV=0 となる。ここで M_k は行列 M
の第k 行を表す。M が (n × n) 行列である場合に直交
解ベクトル V が存在するためには、M のすくなくとも
1 つの行が他の (n-1) 行の一次結合でなければならな
い。ノイズのない系だったならこれを求めることもでき
る。しかし、実際の極・ゼロ点解析では行列 M のすべ
ての要素は測定データがある程度ノイズで汚染されてい
る。ノイズ汚染の影響により、M のどの行も正確には他
の行の一次結合とはならず、純理論的な意味では解ベク
トルが存在しないことになる。

【０１１７】図８のステップ 171 で、行列 M の行が互
いに直交化される。図１１はステップ 171 を構成する
各ステップを詳細に示している。すなわちステップ 193
で、対角要素の大きさが最小である M の行が選択さ
れ、行列 M から削除される。こうして得られる行列 M
は (n-1) × n 行列になる。1 本の列が削除されるので
M の (n-1) 本の行に直交している長さ n のベクトル V
を得ることができる。ノイズに最も汚染されている行
を削除することにより、この削除が測定装置 1の性能の
精度に及ぼす影響を最小にすることが望ましい。基底多
項式のクラスとして直交多項式のひとつであるチェビシ
ェフ多項式が使用されているので、この行列の各行はそ
の行番号と等しい次数のチェビシェフ多項式から成る。
そこで、各行の対角要素はその行で表される特定のチェ
ビシェフ項の電力 (power) の推定値と見なすことがで
きる。対角値が最小の行を削除することによって、最小
電力を有する行、従って相対的なノイズ汚染量が最大の
行が削除される。ステップ171 中の各ステップ 195 〜
203 は行列 M の残りの (n-1) 行をグラム・シュミット
(Gram-Schmidt) 法を用いて互いに直交化することから
構成されている。グラム・シユミツト直交化法の詳細な
説明は、例えば、上に引用したラルストンとラビノウイ
ッツの文献の 256 ページに述べられている。そこで述
べられているように、各行の要素は、最小二乗の順序の
意味で連続する各行から次々と除かれる。直交化を行っ
た結果、(n-1) 行が互いに直交する。従って、長さ n
の列ベクトルがベクトル V と M の (n-1) 本の行から
成る (n × n) 行列の第 n 行と見なされる場合には、M
の (n-1) 本の行に直交するベクトル V が存在する。従
って解が存在する。

【０１１８】図８のステップ 173 で、ベクトル V のラ
ンダムな要素を発生するのに乱数発生器を使用してい
る。ステップ 175 で、上述のグラム・シユミット直交
化法を使用してベクトル V から M の直交化される (n-
1) 本の行の各々を取り除いている。ステップ 175 が終
了すると、n 本の行のすべて (M の (n-1) 本の行とベ
クトル V) が相互に直交する。ランダムに発生した初め
のベクトル V が行列 Mの直交化された (n-1) 本の行の
うちの 1 本あるいは複数本の一次結合である可能性は
極めてわずかである。ステップ 177 で、直交化された
ベクトル V の正当性が試験される。事実、もし発生さ
れた V が M の (n-1) 本の行のうちの 1 つの一次結合
であったとすれば、直交化ベクトル V は非常に小さい
であろう。この場合には、V は無効と判定され、別のラ
ンダムなベクトルを発生する。正当性を判定する経験的
に良好とされる判定基準は、直交化されたベクトル V
の大きさを評価し、これをプロセッサ 45 の丸めノイズ
レベルと比較することである。V の大きさが丸めノイズ
レベルに近付く場合には、無効であると判定される。更
に精度を高めるには、16 桁のプロセッサ 45 に対し
て、ベクトル V はその二乗した大きさが 10^-12 未満で
ある場合、無効と判定される。

【０１１９】○次数の選択 図１２Ａ及び図１２Ｂは測定装置 1 により行われ、ま
た図３のステップ 231に示してある自動次数選択の各ス
テップを詳細に示している。次数選択の目標は、H_E(f)
と H_MD(f) の近似度がある近似基準を満足するまで極と
ゼロ点の次数 (n 及び m) をインクリメントして行くこ
とである。次に、悪い近似が連続して 2つ見つかるまで
m をデクリメントして行く。その結果得られる m と n
は希望する次数である。判定基準を満たす近似がどう
しても得られない場合には、最良の近似条件が記憶され
て提示される。

【０１２０】ステップ 251 で、良好な近似か否かを定
めるための近似判定基準が決まる。実行される特定の解
析に依存するのだが、使用できる判定基準が多数存在す
る。1 つは周波数範囲を通じて近似の程度を 1 点ずつ
試験し、誤差限界外にある点の総数が合格限界点数より
少なければ良好な近似が得られたと判定することであ
る。

【０１２１】誤差限界を計算する第一歩は上に測定した
分散から各周波数点での正規化された分散を計算するこ
とである。1 点毎に正規化された分散は各点での分散を
各点での H_MD(f_j) の大きさで割ったものと定義され
る。次に T_j は、各点毎に、0.01または正規化された分
散のいずれか大きい方と定義される。各点毎の誤差限界
EL_j は

【０１２２】

【数４２】

【０１２３】である。ここで、大きさを二乗した項は、
正規化された分散または 0.01 のいずれか使用している
ものの正規化の結果をもとに戻すために使用されるもの
である。誤差限界は H_MD(f) から両方向に広がっている
窓である。

【０１２４】誤差限界は分散に基づいているから、不規
則な変動を除くために平滑化することが望ましい。13
点平滑化（各側の夫々 6 点ずつとともに各点を平均す
る）を使用することができる。元の誤差限界でのピーク
を保存するためには、各点毎に元の誤差限界と平滑化し
た誤差限界のいずれか大きい方を選択して最終的誤差限
界として使用することができる。

【０１２５】合格限界点数 AN は、周波数範囲内で誤差
限界外にはみ出している点の数であって、これだけの点
がはみ出していてもなお良好な近似を持つという、最大
点数である。この点の数は、経験的に求めるか、あるい
は、

【０１２６】

【数４３】

【０１２７】を使用することができる。F_TOT はその周
波数範囲内の周波数点の総数である。誤差限界平滑化を
行う場合には、以下に説明する近似決定において周波数
範囲内の最高及び最低の 6 点を無視することに注意し
なければならない。

【０１２８】ステップ 253 で初めに次数を m = n = 1
と推定する。ユーザは初期次数推定値を 1 以外に指定
することができることに注意しなければならない。な
お、ユーザは、図１２Ａ及び図１２Ｂのステップが全く
行われないように m と n を指定することを選択でき
る。ユーザは図１２Ａ及び図１２Ｂのステップで m と
nの最大許容値を制限する m_MAX または n_MAX を指定す
ることもできる。

【０１２９】ステップ 255 で、行列 M が m 及び n の
初期の推定次数を使用して再構成される。行列の全パラ
メータが記憶されている場合には、再構成を直接行って
もよい。データ圧縮を使用する場合には、上記の式 (5
3) 及び (54) を使用して、保持されている最大 320 の
パラメータから再構成することもできる。一旦、行列が
再構成されると、図９及び図１０に示す D 及び C ベク
トルの係数が上述のように解かれる。ステップ 257
で、伝達関数の推定値 H_E(f) が上述のようにチェビシ
ェフ係数から再構成される。H_MD(f) と H_E(f) の差が各
点毎に測定され、1点ずつ誤差限界と比較されて、限界
外の点の数を得て記憶する。この数はステップ 259 で
AN と比較され、近似が良好であるか否かが判定され
る。

【０１３０】近似が良好であると判定されれば、m と n
の現在の値が m_BEST 及び n_BESTとして記憶され、ステ
ップ 289 が次に実行される。近似が不良であれば、ス
テップ 273 で、これが今までの最良の近似であるか否
かが判定される。この判定は今回得られた誤差限界外の
点の数と以前に得られた点の数を比較して行う。以前に
得られた点の数は、ステップ 287 から前に戻って繰り
返しが行われていれば得られているはずである。これが
今までの最良の近似であるという判定がなされれば、m
と n は m_BEST 及び n_BEST として記憶され、m と n が
歩進される。これが今までの最良の近似でない場合に
は、記憶は行われず、m と n が歩進される。

【０１３１】ステップ 279 〜 285 で、m と n はユー
ザが課した最大限界値に、または測定装置 1 による最
大限界値 20 に保持される。m と n が共に最大限界値
である場合には、mと n の値がステップ 289 で記憶さ
れる。そうでない場合には、ステップ 255 が再度実行
され、歩進した m と n を用いて他の近似が評価され
る。このように、m と n は、良好な近似が得られるま
で、あるいは m と n が共に最大限界に達するまで歩進
を続ける。

【０１３２】ステップ 301 〜 303 で、ゼロ点の個数 m
が 0 より多ければ、サーボ系もモード系 (modal syst
em) も共に、典型的には、その極の数はゼロ点の個数よ
りも少なくとも 1 つ多いので、m はデクリメントされ
る。このようにゼロ点の個数を減らし、得られた近似を
再評価することによって一層正確な近似に到達すること
ができる。最良の近似は、ゼロ点の個数を 1、及び 2
減らした場合のどちらの場合も共に不適当な近似とな
る、そのようなゼロ点の数が得られたときに達成された
と考えられる。ステップ 303 で、m が 0 に等しいと分
かった場合にはステップ 333 が次に実行される。そう
でない場合には、ステップ 305 〜 307 で、行列 M が
m と n の現在の値について再構成され、チェビシェフ
係数について解かれ、H_E(f) のH_MD(f) への近似度が各
点で評価される。

【０１３３】ステップ 309 では、近似が近似判定基準
に対して評価される。このとき近似が良好ならば m と
n が記憶され、m はステップ 301 で再びデクリメント
される。近似が不良である場合、または m=0 の場合に
は、ステップ 333 が実行される。m が 0 でない場合に
は、m はステップ 315 で再び減算される。ステップ 31
7 と 319 で、新しい解が求められ、新しい近似が評価
される。近似が今度は良好であれば、m と n は記憶さ
れ、ステップ 301 で m がデクリメントされる。一方、
近似がまたも不良であれば、ステップ 333 で m は先に
記憶されている m_BEST に再設定される。この場合、今
度は m を 1、及び 2 だけデクリメントすれば、いずれ
も不良近似を生ずる値になっている。n は m を最初に
デクリメントした以前に良好な近似を与えた値になって
いる。

【０１３４】図１３Ａ及び図１３Ｂは図１２Ａ及び図１
２Ｂに示す諸ステップの代わりに図３のステップ 231
で使用することができる自動次数選択の各ステップを示
している。概観すれば、図１３Ａ及び図１３Ｂの諸ステ
ップはいくつかの直列意志決定ブロックに分けることが
できる。これらは多項式の次数を変えるべきか否か、ま
たは新しい多項式係数を求めるべきか否かを決定するも
のである。ブロック 531 では、分子の次数 m を以前に
デクリメントした際に誤差対信号比が劣化していた場合
に、分子の次数 m をインクリメントして新しい係数を
求める。ブロック 533 では、誤差対信号比がノイズ対
信号比を予め定めた量だけ超過し、これにより誤差が測
定データ上のノイズによってのみ生じたのではないこと
が示された場合に、m と n を共にインクリメントし、
新しい係数が求められる。ブロック 535 では、高次側
の係数の測定データに及ぼす影響がある最小量より少な
い場合、m を減少させて新しい係数を求める。ブロック
537 では n が最大許容次数になっても得られる近似が
不適当である場合には、新しいデータを取り、極・ゼロ
点解析を全く初めからやり直す。最後に、ブロック 539
は、近似が充分良好であることを確認する。ブロック
531 〜 537 で近似が条件付きで良好と判定され、かつ
次数 m と n が等しい場合には、近似は充分良好である
と見なされる。また、誤差がバイアスされ過ぎていない
と判定されても、その近似も充分良好であると見なされ
る。そうでなければ、m と n は共にインクリメントさ
れ、新しい係数が定められる。

【０１３５】ブロック 531 のステップ 401 で、分子 P
(s) の次数 m が図３のステップ 231 の今回の実行の間
に減らされていたか否かの判定が行われる。判定の結果
によってフラグ U0 がセットされる。m がこれまでに減
らされていれば、重み付きの誤差対信号比 V0 がステッ
プ 405 で計算され、ステップ 407 で評価される。V0は

【０１３６】

【数４４】

【０１３７】と定義される。V0 がより高次の m につい
て先に行った近似の評価から得られたV0 (V1 として記
憶されている) より 5% 以上大きい場合には m とスレ
ッショルド (S1 、以下に説明する) を共にインクリメ
ントしてステップ 69 を繰り返す。ステップ 407 で行
われる判定により、m を減らしたとき、近似の質（誤差
対信号比で測定する）が低下することが示されれば、m
をインクリメントし元の値に戻してステップ 69 を繰り
返す。

【０１３８】ブロック 533 のステップ 431 で、誤差対
信号比が 5% 増加しなかった、あるいは m がこれまで
減らされていないと仮定すれば、ノイズ対信号比 R が
計算される。R は分散の定義をノイズ電力として使用
し、上記の式 (57) の分母を使用して、次のように定義
される。

【０１３９】

【数４５】

【０１４０】ステップ 441 で、ノイズ対信号比は誤差
対信号比の倍数と比較される。倍数 S0は近似の質を測
定データ上のノイズと関係づけるのに使用される。S0
は経験的に 2 となるように選ばれてきたが、S0 は、測
定データに大きなノイズが乗っている場合に、近似の基
準を緩くするようにステップ 465 で増すことができ
る。

【０１４１】誤差対信号比 (V0) がノイズ対信号比 R
の S0 倍よりも大きくなった場合には、近似は不充分で
あると見なされる。分母の次数が最大許容数より少なけ
れば、m と n は共に 1 だけインクリメントされ、ステ
ップ 69 が繰り返される。これにより P(s) と Q(s) の
次数が共に増加して、推定伝達関数が測定された伝達関
数のデータをもっと良く近似するようになる。近似の質
は、m または nが測定された伝達関数データの真の次数
を越すまでは向上していくはずである。測定装置 1 で
は、m と n の各々の最大許容次数が 20 (合計では、組
み合わされる正負の周波数に対してそれぞれ 40) であ
る。ステップ 443 で、分母の最大次数 20 に達する
と、ブロック 535 に入る。また、近似の質がステップ
441 で条件付き良好と判定された場合にも、ブロック 5
35 に入る。

【０１４２】このブロック 535 のステップ 449 で、ス
レッショルド次数 S1 が求められる。ここで、S1 は次
のように定義される。

【０１４３】

【数４６】

【０１４４】また、S1 は式 (59) の不等式を満足する
M の最大値である。T_M は次のように定義される。

【０１４５】

【数４７】

【０１４６】ここで、C_k は分子のチェビシェフ多項式
P(s) の係数である。なお、T0 は次のように定義され
る。

【０１４７】

【数４８】

【０１４８】ここで、新しい重み付け関数 L(s) は次の
ように定義される。

【０１４９】

【数４９】

【０１５０】L(s) を表す式中、W(s) の項はピークを強
調するために設けられており、Q(s)項は分母の影響を除
き、H_MD(f_j) 項は大きさの小さい測定データの区域の影
響を最小限にしている。

【０１５１】ステップ 451 と 453 において、分子の次
数 m はスレッショルドの次数 S1-1まで減らされる。分
母の次数 n は図１３Ａ及び図１３Ｂの流れ図では決し
て減らされることがないことに注意しなければならな
い。分子の次数は、分子を高すぎる次数にしておいてし
まったせいで偽りのゼロ点が入り込むことがないように
するため、しきい値まで減らされる。S1 を決めるにあ
たっては、P(s) 中の比較的効果が小さい係数のみを除
去するようにし、これにより P(s) の次数が過度に減ら
されないように定める。式 (59) の係数 1.5 は最適数
のゼロ点を近似するため変更してよい。

【０１５２】ブロック 537 のステップ 461 で、ステッ
プ 441 のノイズ対信号比と誤差対信号比の比較が繰り
返される。誤差対信号比がノイズ対信号比の予め定めた
倍数を超え、かつ n が最大許容次数に達したか、ある
いは大きい場合には、測定データが汚染されていると判
定される。もしそうなら、データ汚染を償うため、ノイ
ズ対信号比 S0 の許容倍数を増して、極・ゼロ点解析全
体を初めからやり直す。

【０１５３】分母の次数 m が最大許容次数に到達も超
過もしなかった場合には、ブロック539 のステップ 491
を実行する。分子と分母の次数が等しくない場合に
は、近似が良好であると判定し、図１３Ａ及び図１３Ｂ
のステップを終了して、図３のステップ 77 を実行す
る。また、分子と分母の次数が等しい場合には、近似の
質の最終的チェックを行う。チェックの結果が成功であ
る場合、あるいはチェックが不成功でありかつ分母の次
数が許容最大次数と等しいかそれより大きい場合には、
近似は充分に良好であると判定される。チェックが不成
功でかつ分母の次数が最大許容次数より小さい場合に
は、分子と分母の次数を共にインクリメントし、ステッ
プ 69 で新しい係数を決める。

【０１５４】ブロック 539 のステップ 493 の最終チェ
ックはオプションであるが、ここでは H_MD(f_j) の各周
波数点での H_E(s) と H_MD(f_j) の間の誤差の極性を解析
する。もし近似が充分に良好であるなら、これらの誤差
は 0 を中心として集中しなければならず、またその極
性が頻繁に変わらなければならない。最終チェックはH_E
(s) 及び H_MD(f_j) の実数部と虚数部について別々に行
われる。チェックが成功であるとされる判定基準は、A1
個よりも多くの隣接する点の誤差が同じ極性でないと
いうことである。ただし、その大きさがプロセッサ 45
の計算ノイズの大きさより小さい誤差点は考慮されな
い。一貫して良好な近似が行えるようにするため、A1
を定めるにあたっては、近似が良好であってもそのうち
の予め定めた割合、例えば、5% がステップ 493 で排除
されるように選択される。A1 は使用する周波数データ
点の数（測定装置 1 では 800 点）と統計的に関連して
いる。なお、A1 は経験的に決めてもよい。

【０１５５】○多項式の変換 ステップ 335 で、行列 M を n, m の既に求められた最
良近似について再構成し、その後チェビシェフ係数につ
いて解く。これで図３のステップ 77 に示すような普通
の多項式への変換が行われる。図３のステップ 77 で、
H_E(s) のチェビシェフ多項式 P(s) と Q(s) が普通の多
項式に変換される。この変換は、当業者であれば、式
(10) に関して上に述べたチェビシェフ帰納的関係を利
用して容易に行うことができる。変換の結果、P(s) と
Q(s) は複素変数 s で表された普通の多項式となる。

【０１５６】○求根 ステップ 79 では、定義により H_E(s) のゼロ点と極で
ある P(s) と Q(s) の根を求めるために求根手段が用い
られている。この決定を行うために多数の既知の求根手
段を使用することができる。特に、ヒューレット・パッ
カード社の HP-9825 形卓上コンピュータに使用され、
同社の部品番号 9825--10001 として入手できる説明書
「汎用ユーティリティ・ルーチン」の 213 〜 224 ペー
ジに記されている求根手段は有効かつ有用である。求根
手段についてのこれ以上の説明は、1967 年 4 月発行の
Journal of the Association for Computing Machiner
y、第 14 巻第 2 号に掲載されたムーア (J. B. Moore)
著の「多項方程式を解く収束アルゴリズム」に述べら
れている。

【０１５７】最後に、図３のステップ 81 で、推定され
た伝達関数 H_E(s) の極とゼロ点が測定装置 1 の表示器
7 上に表示される。多数の旧来の表示手順のいずれも
有効に使用することができ、極とゼロ点は例えば、直交
座標または極座標の形でこれを表示することができる。

【０１５８】

【発明の効果】以上詳述したように、本発明によれば、
極・ゼロ点の測定が迅速かつ正確に行なわれ、また測定
値に含まれるノイズによる誤差を補正して測定対象の伝
達関数の極・ゼロを求めることができる。

【０１５９】

【図面の簡単な説明】

【図１】本発明の一実施例に基づいて構成された極・ゼ
ロ測定装置を示す外観図である。

【図２】図１に示した装置の電気的ブロック図である。

【図３】図１の装置の動作を示す流れ図である。

【図４】図３中の流れ図の一部を示す詳細な流れ図であ
る。

【図５】被測定デバイスと測定ノイズの理想モデルを示
す図である。

【図６】図３に示した流れ図の一部を示す詳細な流れ図
である。

【図７】図３に示した流れ図の一部を示す詳細な流れ図
である。

【図８】図３に示した流れ図の一部を示す詳細な流れ図
である。

【図９】本発明の実施例で用いられるベクトル及び行列
の構成を示す図である。

【図１０】本発明の実施例中に用いられる行列の区分を
示す図である。

【図１１】図３に示した流れ図の一部を示す詳細な流れ
図である。

【図１２Ａ】図３に示した流れ図の一部を示す詳細な流
れ図である。

【図１２Ｂ】図３に示した流れ図の一部を示す詳細な流
れ図である。

【図１３Ａ】図３に示した流れ図の一部を示す詳細な流
れ図である。

【図１３Ｂ】図３に示した流れ図の一部を示す詳細な流
れ図である。

【符号の説明】

1：測定装置 3：デバイス 5：キーボード 7：表示器 21,23：入力部 25：トリガ源 27：アナログ源 29,31：アナログディジタル変換器 33：ディジタルフイルタ 35：制御器 37：局部発振器 39：ディジタル信号源 41：CPU 44：FFT 45：プロセッサ 47：RAM 49：ROM 51：バス

Claims (7)

(57)【特許請求の範囲】
1. 【請求項１】以下のステップ（ａ）〜（ｉ）を設けてな
   る測定対象の推定伝達関数の極・ゼロ点決定方法： （ａ）測定対象に刺激信号を与えて応答信号を発生させ
   る； （ｂ）前記刺激信号及び応答信号を検出する； （ｃ）前記刺激信号及び応答信号から測定された伝達関
   数を決める； （ｄ）推定伝達関数についての極の個数の初期値とゼロ
   点の個数の初期値を選択する； （ｅ）前記極の個数の初期値とゼロ点の個数の初期値か
   ら前記推定伝達関数を求める； （ｆ）前記測定された伝達関数と推定伝達関数を比較
   し、前記比較から近似の程度を判定する； （ｇ）前記推定伝達関数の前記極の個数とゼロ点の個数
   を、最良の近似が得られるまで増加させていく； （ｈ）ステップ（ｇ）で求められた前記最良の推定伝達
   関数の前記ゼロ点の個数を、それに続く所定回数の相次
   ぐ前記ゼロ点の減少の結果がいずれも良い近似を与えな
   くなるまで減少させていく； （ｉ）前記ステップ（ｈ）で得られた推定伝達関数の根
   として極及びゼロ点を計算する。
2. 【請求項２】前記極の個数の初期値と前記ゼロ点の個数
   の初期値はともに１であることを特徴とする請求項１記
   載の方法。
3. 【請求項３】前記ステップ（ｆ）は下記のステップ（ｆ
   −１）〜（ｆ−４）を有することを特徴とする請求項２
   記載の方法： （ｆ−１）所望の複数の周波数点の各々における誤差限
   界を定める； （ｆ−２）前記周波数点の各々において前記測定された
   伝達関数と前記推定伝達関数の間で減算を行う； （ｆ−３）前記周波数点の各々において前記減算の結果
   を前記誤差限界と比較する； （ｆ−４）前記誤差限界を上回った前記結果の個数を数
   える。
4. 【請求項４】前記ステップ（ｆ）は下記のステップ（ｆ
   −５）及び（ｆ−６）を有することを特徴とする請求項
   ３記載の方法； （ｆ−５）前記誤差限界を上回った前記結果の個数を記
   録する； （ｆ−６）前記記録された個数が予め定められた許容個
   数よりも少ない場合に良い近似が得られたと判定する。
5. 【請求項５】前記所定回数は２回であることを特徴とす
   る請求項４記載の方法。
6. 【請求項６】以下のステップ（ａ）〜（ｉ）を設けてな
   る測定対象の推定伝達関数の極・ゼロ点決定方法： （ａ）測定対象に刺激信号を与えて応答信号を発生させ
   る； （ｂ）前記刺激信号及び応答信号を検出する； （ｃ）前記刺激信号及い応答信号から測定された伝達関
   数を決める； （ｄ）推定伝達関数についての極の個数の初期値とゼロ
   点の個数の初期値を選択する； （ｅ）前記初期値個数の極と初期値個数のゼロ点を有す
   る推定伝達関数を求める； （ｆ）前記測定された伝達関数と推定伝達関数の間の誤
   差対信号比を測定する； （ｇ）前記誤差対信号比が前記測定された伝達関数のノ
   イズ対信号比より大きい場合には、前記推定伝達関数の
   前記極の個数とゼロ点の個数を増加させる； （ｈ）前記増加した個数の極とゼロ点を持つ推定伝達関
   数を求める； （ｉ）推定伝達関数の根として極及びゼロ点を計算す
   る。
7. 【請求項７】前記ステップ（ｈ）の後に以下のステップ
   （ｊ） 〜（ｌ）を有することを特徴とする請求項６記
   載の方法： （ｊ）前記推定伝達関数の高次の係数の影響を判定す
   る； （ｋ）前記影響が予め定められた水準よりも小さい場合
   には前記推定伝達関数のゼロ点の個数を減少させる； （ｌ）ゼロ点を減少させた後の個数の極及びゼロ点を持
   つ推定伝達関数を求める。