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

Description

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

〔従来技術及びその問題点〕

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

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

〔発明が解決しようとする問題点〕

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

〔問題点を解決するための手段及び作用〕

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

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

〔発明の実施例〕

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

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

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

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

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

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

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

デバイス３の伝達関数の測定値H_MD（_ｊ）は第５図に
より定義することができる。第５図はデバイス３の出力
が入力、伝達関数及び付加的なノイズの項に関係してい
ることを示している。従って、伝達関数H_MD（_ｊ）は
次式により陰に定義することができる。

Ｙ（_ｊ）＝H_MD（_ｊ）Ｘ（_ｊ）＋N_F（_ｊ）
（１） この式は上線を「グループ平均」と定義したとき次のよ
うに書き替えることができる。

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

ここで、 の項は無視できる程小さいと仮定し、また自己スペクト
ル及び重なりスペクトルは次のように定義している。

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

第４図のステップ107で、式（６）て定義したH
_MD（_ｊ）がメモリに記憶される。ステップ63では、測
定された伝達関数の集合したコヒーレンス（ensemble c
oherence）が測定される。この結果は極及びゼロ点の測
定におけるノイズ誤差を補正するために後で使用され
る。一般に、ランダム変数Ｘに関する分散は次式で定義
される。

σ^２＝▲▼−^２ （７） 分散は周知の多数の方法のいずれかによっても推定する
ことができるが、1957年３月グッドマン（N.R.Goodma
n）のニューヨーク大学博士論文「二次元静的ガウス過
程のスペクトル、コスペクトル及び直交スペクトルの合
同推定（On The Joint Estimation Of The Spectra,Cos
pectrum,And Quadrature Spectrum Of A Two−Dimensio
nal Stationary Gaussian Process）」の第131ページに
ある式4.92及び4.93によって厳密に求めることができ
る。これらの式は第二種ベッセル関数に書き替えること
ができるので、これを積分して、 が得られる。ここでｎは１より大きく、またｎはグルー
プ（ensemble）中の測定値の個数である。式（８）に使
用しているコヒーレンス関数は、次式で定義される。

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

Ｅ（_ｊ）＝Ｗ（_ｊ［H_E（i2π_ｊ）−H
_MD（_ｊ）］ （10） ここで、H_MDは実数の周波数点について測定され、またH
_Eは複素数ｓ上の関数として求められるが、H_MDはｓに一
般化でき、H_MDとH_Eを直接比較対照することができる。
当業者ならば、ｓとi2πの解析関数を等価であると見
なす解析接続の手続によりこの一般化を行うことができ
るであろう。

最小二乗近似解析で定義される全重み付け誤差ε′は ここで、誤差と重み付け関数は各周波数点で評価され、
合計は関心のある周波数範囲にわたって取られる。重み
付け関数Ｗ（）は最小二乗近似解析において所望の
の領域をエンファサイズ（強調）するのに使用される。
H_E（ｓ）は２つの多項式の有理式である。

Ｐ（ｓ）の係数がc₀,c₁,…,c_mでありまたＱ（ｓ）の係
数がd₀,d₁,…,d_nであると仮定すれば、最小二乗近似解
析はc₀,c₁,…,c_m及びd₀,d₁,…,d_nの各々に関するε′の
偏導関数を得、これらの偏導関数を０に等しいとおいて
同次方程式（homogeneous equations）を作り出すこと
により行われる。この最小二乗近似解析によれば、同次
方程式を連立させて解くことにより、所望の最小の重み
付き誤差ε′を見つけることができる。

具体的に説明すれば、式（11）において合計を取られる
各項|E（f_j）|²に式（10）を代入し、さらにこれに式
（12）のH_Eの定義を代入する。ここで、各f_jは具体的な
数値であるから、上述の代入によって計算される各項に
は記号としては多項式Ｐ（ｓ）,Q（ｓ）の係数c₀,c₁,‥
‥,c_m,d₀,d₁,・・・,d_nしか含まれておらず、またこれ
らの係数についての有理式になっている。従って、この
ようにして計算された各項を合計した値、すなわちε′
もこれらの係数の有理式になっている。ここで求めよう
としているのは、式（11）の値ε′を最小にする（つま
り現実の測定結果H_MD（f_i）をプロットした曲線に対し
て式（12）で表わされる関数が表わす曲線が最小二乗近
似の意味で最も良く当てはまるようにする）係数c₀,c₁,
‥‥,c_m,d₀,d₁,・・・,d_nである。当業者周知の通り、
最小二乗近似解析でこれらの係数を求めるには、先ず、
上述の代入・合計によって得られたこれらの係数につい
ての有理式で表現されているε′をこれらの係数の夫々
について偏微分した偏導関数を求める。次に、これらの
偏導関数を夫々０とおいて得られる式をこれらの係数に
ついての連立方程式として、これらの係数について解け
ばよい。実際に計算してみると、これらの偏導関数を０
と置いた式は全てこれらの係数の一次方程式になるの
で、この連立方程式は通常の一次連立方程式になり、簡
単にこれらの係数を求めることができる。

最小二乗近似解析の一般的な解説は、たとえばラルスト
ン（Ralston）及びラビノウイッツ（Rabinowitz）著
「数値解析の第一過程」第２版の第６章に述べられてい
る。

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

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

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

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

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

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

測定データ上の高ノイズ領域の影響を最小にするように
重み付け関数を補正することも望ましいことである。第
４図のステップ111で測定され、また式（９）で定義さ
れたコヒーレンス関数により、H_MD（_ｊ）の各測定周
波数領域データ点のノイズを推定できる。ここで、ノイ
ズ補正係数 を乗ずることにより、コヒーレンス関数が０に近い値を
取る高ノイズ領域をデエンファサイズする、つまりその
影響を弱める。係数を0.05とすれば、きれいなデータの
点をノイズの多いデータ点より20倍も強調することがで
きる。位相導関数の絶対値を含んでいる、周波数軸に沿
う各点での最終的な結果は である。

ステップ133で、ノイズを補正された重み付け関数はオ
ーバーフローしないように１に正規化され、ステップ13
5で、完全に無視されてしまうデータがないように、0.0
2のしきい値が導入される。従って、重み付け関数は、 と定義される。ただし、上式の右辺中の分数式の分母の
右下に付いている添字MAXは、今考えている範囲内で当
該添字が付けられている式が取る値の最大値を式の値と
して採用するということを表わす。このようにして得ら
れた範囲内の最大値で割ることによって、上述の１へ正
規化を行っている。

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

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

次に各データ点の振幅を関連する線分の振幅から差し引
く。この残りをしきい値と比較して、妥当なピークが検
出されたか否かを確かめる。

という代表的なしきい値を用いれば、この残りをその点
での標準偏差と比較することができ、偽のノイズピーク
を排除することができる。予め決めた個数の連続点での
上述の残りがしきい値より大きい場合には、真のピーク
値が検出されたという最終判定を下してよい。なお、４
つの連続点を用いれば、正確なピーク検出を行い得るこ
とが分かっている。

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

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

2T_m（ｘ）t_n（ｘ）＝T_n+m（ｘ）＋T_n-m（ｘ） ただし ｎｍ （21） 所望のＰ（ｓ）とＱ（ｓ）のチェビシェフ多項式は上に
引用したNBSのハンドブックの782ページに書かれている
帰納的関係（recurrence relationship）に従って容易
につくることができる。多項式Ｐ（ｓ）とＱ（ｓ）は複
素変数ｓの多項式であり、実数、虚数あるいは複素数の
係数をもっている。当業者ならは、純実数のチェビシェ
フ多項式に対して与えられる帰納的関係を使用して複素
チェビシェフ多項式を容易に作ることができる。H
_MD（_ｊ）の時間データの測定値がベースバンドの場合
のように純実数である場合には、複素チェビシェフ多項
式の使用を容易にするためエルミート対称（Hermitian
symmetry）を使用することができる。モーダル解析を行
う場合には、エルミート対称の制約を捨てて基底多項式
を修正するのが好ましい。Ｐ（ｓ）とＱ（ｓ）の基底多
項式を作るためには、２つの多項式の次数を初めに推定
することが必要である。この推定は伝達関係の極とゼロ
点の数を推定することと等価である。各多項式に対して
最初に次数を１と推定すれば、１つまたはもっと多くの
極を持つ多項式を近似することができる。解析対象の装
置についての情報をはじめにもっと多く利用できる場合
には、最初の推定次数をもっと高くしてよい。

○係数の決定 第３図のステップ69で、Ｐ（ｓ）とＱ（ｓ）の係数が求
められる。第８図はステップ69を構成する各種のステッ
プを一層詳細に示している。式（10）に定義した点毎の
誤差項は次のように書き直される。

ただし、記号は便宜上省いてあり、また添字ｊは周波
数軸に沿ってサンプルしたデータを使用していることを
示す。分母の多項式Ｑ（ｓ）を両辺にかけて下式を得
る。

E_jE_j＝W_j［P_j−H_MDjQ_j］ （23） 式（６）を用い、また新しい項A_jとB_jを A_j＝G_YXj,B_j＝G_XXj （24） と定義することにより、（23）を次のように簡単化して
書き直すことができる。

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

E_jQ_jB_j＝W_j［B_jP_j−A_jQ_j］ （26） 最後に、この乗算を行わない場合、_ｊをつぎのように
定義する。_ｊ ＝E_jQ_j （27） 第８図のステップ161で、ユーザが予め指定した極とゼ
ロ点を伝達関数の推定値から除去する。この操作は必須
でないが、既知の極とゼロ点をユーザが予め指定すると
未知の極とゼロ点の解析を更に正確に行うことができ
る。次式を参照して、既知のゼロ点をＰ（ｓ）多項式か
ら除き、また既知の極をＱ（ｓ）多項式から除く。

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

W_j,U_j,A_j/B_j及びV_jの各項は推定される多項式_j,_ｊ
に適用される各点毎の重み付けの項（測定データを含
む）と見なすことができる。簡単化のため新しい項W_Pj
とW_qjを定義できる。

従って、式（32）を次のように簡単化して書換えること
ができる。_ｊ ＝W_{Pj ｊ}−W_{qj ｊ} （34） 第８図のステップ163で、測定したデータ、重み付け関
数及び推定伝達関数の多項式を、第１図及び第２図に示
す測定装置１を用いての解析を容易にするため行列記法
に変換する。この変換は当業者には容易に行うことがで
きる。係数c_kとd_kを有する２組の直交多項式φ_ｋとθ_ｋ
を上に述べた帰納的関係を用いて作り、_ｊと_ｊを表
すのに使用する。

チェビシェフ多項式はデータに依存しない（従って測定
毎に求め直す必要はない）ので、φ_ｋとθ_ｋをチェビシ
ェフ多項式として構成することが好ましい。これら２つ
の多項式は異なる次数ｍとｎをもっているが、ここでは
便宜上ｍ＝ｎでかつ両多項式ともθ_ｋと書かれると仮定
する。当業者ならば、ｍとｎが等しくない場合に一般化
することができるだろう。式（35）は次のように書き直
すことができる。

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

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

＝θＣ （38） ＝θＤ （39） また、誤差方程式（34）は次のように書き替えることが
できる。

＝W_PθＣ−W_qθＤ （40） さて、誤差の二乗の和は行列記法で次のように書くこと
ができる。

ただし、肩文字Ｔは複素共役転置を表す。

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

０＝θ^TW_P ^T（W_PθＣ−W_qθＤ）＝θ^TW_P ^T （42） ０＝θ^TW_q ^T（W_PθＣ−W_qθＤ）＝θ^TW_q ^T （43） 明確化のため、上記の式（42）と（43）を次のように書
き替える。

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

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

０＝θ^TW_q ^TW_PθＣ−θ^Ｔ（W_q ^TW_q−σ^２）θＤ （46） 第８図のステップ169で、解行列が構成され、最小二乗
近似解析の最小化のためのチェビシェフ多項式の係数が
求められる。方程式（42）と（43）は方程式（44）と
（45）のように書き直すことができる。ここでは、もし
望むならば、ノイズバイアスが上に述べたように除かれ
ていると仮定される。新しい関数F,H,H^T及びＧを次のよ
うに定義する。

Ｆ＝θ^TW_P ^TW_Pθ （47） H^T＝θ^TW_P ^TW_qθ （48） Ｈ＝θ^TW_q ^TW_Pθ （49） Ｇ＝θ^TW_q ^TW_qθ （50） このように定義することにより、次の関係式 GD−HC＝０ （式（45）より） （51） −H^TD＋FC＝０ （式（44）より） （52） が成立し、第10図に示す、４つの部分を持つ行列Ｍ及び
２つの部分を持つベクトルＶの行列乗算として示されて
いる同次方程式系を定義する。

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

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

2T_jT_k＝Ｔ_|j−k|＋Ｔ_|j＋k| （53） 従って、 のようなもっと具体的な関係を上記のデータ圧縮に使用
することができる。［なお、式（54）ではについての
和がとられているが、これは測定データのフーリエ変換
によって得られる全ての周波数データ点｛_ｉ｝につい
ての和をとることを意味している（たとえば、式（11）
ではのインデクスとしてｊを用いているが、式（54）
では添字ｊはすでにチェビシェフ多項式の次数を表わす
ために使用されているので、ここではインデクスとして
ｉを使用していることに注意されたい）。その結果、行
列Ｍの各要素が行列内の他の２つの合計結果を足したも
のになる。

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

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

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

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

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

誤差限界を計算する第一歩は上に測定した分散から各周
波数点での正規化された分散を計算することである。１
点毎に正規化された分散は各点での分散を各点でのH_MD
（_ｊ）の大きさで割ったものと定義される。次にT
_jは、各点毎に、0.01または正規化された分散のいずれ
か大きい方と定義される。各点毎の誤差限界EL_jは EL_j＝±10T_j|H_MD（_ｊ）|² （55） である。ここで、大きさを二乗した項は、正規化された
分散または0.01のいずれか使用しているものの正規化の
結果をもとに戻すために使用されるものである。誤差限
界はH_MD（）から両方向に広がっている窓である。

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

合格限界点数ANは、周波数範囲内で誤差限界外にはみ出
している点の数であって、これだけの点がはみ出してい
てもなお良好な近似を持つという、最大点数である。こ
の点の数は、経験的に求めるか、あるいは、 を使用することができる。F_TOTはその周波数範囲内の周
波数点の総数である。誤差限界平滑化を行う場合には、
以下に説明する近似決定において周波数範囲内の最高及
び最低の６点を無視することに注意しなければならな
い。

ステップ253で初めに次数をｍ＝ｎ＝１と推定する。ユ
ーザは初期次数推定値を１以外に指定することができる
ことに注意しなければならない。なお、ユーザは、第12
A図及び第12B図のステップが全く行われないようにｍと
ｎを指定することを選択できる。ユーザは第12A図及び
第12B図のステップでｍとｎの最大許容値を制限するm
_MAXまたはn_MAXを指定することもできる。

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

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

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

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

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

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

ブロック531のステップ401で、分子Ｐ（ｓ）の次数ｍが
第３図のステップ231の今回の実行の間に減らされてい
たか否かの判定が行われる。判定の結果によってフラグ
U0がセットされる。ｍがこれまでに減らされていれば、
重み付きの誤差対信号比V0がステップ405で計算され、
ステップ407で評価される。V0は と定義される。V0がより高次のｍについて先に行った近
似の評価から得られたV0（V1として記憶されている）よ
り５％以上大きい場合にはｍとスレッショルド（S1、以
下に説明する）を共にインクリメントしてステップ69を
繰り返す。ステップ407で行われる判定により、ｍを減
らしたとき、近似の質（誤差対信号比で測定する）が低
下することが示されれば、ｍをインクリメントし元の値
に戻してステップ69を繰り返す。

ブロック533のステップ431で、誤差対信号比が５％増加
しなかった、あるいはｍがこれまで減らされていないと
仮定すれば、ノイズ対信号比Ｒが計算される。Ｒは分散
の定義をノイズ電力として使用し、上記の式（57）の分
母を使用して、次のように定義される。

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

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

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

T_S1-11.5T0 （59） また、S1は式（59）の不等式を満足するＭの最大値であ
る。T_Mは次のように定義される。

ここで、C_kは分子のチェビシェフ多項式Ｐ（ｓ）の係数
である。なお、T0は次のように定義される。

ここで、新しい重み付け関数Ｌ（ｓ）は次のように定義
される。

Ｌ（_ｊ）＝Ｗ（_ｊ）・|Q（i2π_ｊ）|²・|H_MD（
_ｊ）|⁴ （62） Ｌ（ｓ）を表す式中、Ｗ（ｓ）の項はピークを強調する
ために設けられており、Ｑ（ｓ）項は分母の影響を除
き、H_MD（_ｊ）項は大きさの小さい測定データの区域
の影響を最小限にしている。

ステップ451と453において、分子の次数ｍはスレッショ
ルドの次数S1−１まで減らされる。分母の次数ｎは第13
A図及び第13B図の流れ図では決して減らされることがな
いことに注意しなければならない。分子の次数は、分子
を高すぎる次数にしておいてしまったせいで偽りのゼロ
点が入り込むことがないようにするため、しきい値まで
減らされる。S1を決めるにあたっては、Ｐ（ｓ）中の比
較的効果が小さい係数のみを除去するようにし、これに
よりＰ（ｓ）の次数が過度に減らされないように定め
る。式（59）の係数1.5は最適数のゼロ点を近似するた
め変更してよい。

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

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

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

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

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

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

〔発明の効果〕

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

【図面の簡単な説明】

第１図は本発明の一実施例に基づいて構成された極・ゼ
ロ測定装置を示す外観図である。 第２図は第１図に示した装置の電気的ブロック図であ
る。 第３図は第１図の装置の動作を示す流れ図である。 第４図は第３図中の流れ図の一部を示す詳細な流れ図で
ある。 第５図は被測定デバイスと測定ノイズの理想モデルを示
す図である。 第６図ないし第８図は第３図に示した流れ図の一部を示
す詳細な流れ図である。 第９図は本発明の実施例で用いられるベクトル及び行列
の構成を示す図である。 第10図は本発明の実施例中に用いられる行列の区分を示
す図である。 第11図は第３図に示した流れ図の一部を示す詳細な流れ
図である。 第12A図、第12B図、第13A図及び第13B図はそれぞれ第３
図に示した流れ図の一部を示す詳細な流れ図である。 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 (32)

【特許請求の範囲】
1. 【請求項１】以下のステップ（ａ）〜（ｊ）を設けてな
   る測定対象の推定伝達関数の極・ゼロ点決定方法： （ａ）測定対象に刺激信号を与えて応答信号を発生させ
   る； （ｂ）前記刺激信号及び応答信号を検出する； （ｃ）前記刺激信号及び応答信号から複数の周波数点に
   おける自己スペクトル及び重なりスペクトルを計算す
   る； （ｄ）前記周波数点の各々における前記重なりスペクト
   ル上のノイズレベルを測定する； （ｅ）前記周波数点の各々における推定伝達関数と測定
   された伝達関数の間の誤差を表わす関数を求める； （ｆ）前記周波数点の各々において前記誤差を表わす関
   数を前記推定された伝達関数の係数について微分し、こ
   れにより前記周波数点の各々において前記重なりスペク
   トルを二乗する； （ｇ）前記周波数点の各々において前記重なりスペクト
   ルの二乗から前記測定されたノイズレベルを減算する； （ｈ）前記周波数点の各々において前記微分された誤差
   を０と置く； （ｉ）前記推定伝達関数の係数を定める； （ｊ）前記推定伝達関数の極及びゼロ点を前記推定伝達
   関数の根として求める。
2. 【請求項２】前記ステップ（ｂ）は前記検出された刺激
   信号及び応答信号を周波数領域に変換することを特徴と
   する特許請求の範囲第１項記載の方法。
3. 【請求項３】前記ステップ（ｄ）は下記のステップ（ｄ
   −１）及び（ｄ−２）を有することを特徴とする特許請
   求の範囲第２項記載の方法： （ｄ−１）前記周波数点の各々において前記自己スペク
   トル及び重なりスペクトルから前記測定対象のコヒーレ
   ンス関数を計算する； （ｄ−２）前記周波数点の各々において前記コヒーレン
   ス関数から前記推定伝達関数についての分散を計算す
   る。
4. 【請求項４】前記ステップ（ｇ）は前記周波数点の各々
   において前記重なりスペクトルの大きさの二乗から前記
   分散を減算することを特徴とする特許請求の範囲第３項
   記載の方法。
5. 【請求項５】前記刺激信号は前記周波数点上のランダム
   ノイズであることを特徴とする特許請求の範囲第４項記
   載の方法。
6. 【請求項６】前記ステップ（ｆ）、（ｈ）及び（ｉ）は
   行列演算を用いて行われることを特徴とする特許請求の
   範囲第５項記載の方法。
7. 【請求項７】前記ステップ（ｆ）、（ｈ）及び（ｉ）は
   下記のステップ（ｋ）〜（ｎ）を行うことを特徴とする
   特許請求の範囲第６項記載の方法： （ｋ）推定伝達関数の係数を要素とするベクトルを発生
   する； （ｌ）行列への前記ベクトルの乗算を０と置く； （ｍ）前記行列の行を直交化する； （ｎ）前記乗算を行う。
8. 【請求項８】前記ステップ（ｍ）は下記のステップ（ｍ
   −１）〜（ｍ−３）を含むことを特徴とする特許請求の
   範囲第７項記載の方法： （ｍ−１）他の全ての行の対角要素よりも小さな対角要
   素を持つ最小行を削除する； （ｍ−２）一番上にある行をそれよりも下にある全ての
   行から減算する； （ｍ−３）下にある行の各々をそれよりも下にある全て
   の行から減算する。
9. 【請求項９】前記ステップ（ｍ−２）及び（ｍ−３）は
   最小二乗近似技術に基づいてことを特徴とする特許請求
   の範囲第８項記載の方法。
10. 【請求項１０】前記ステップ（ｍ−２）は下にある行か
    ら上にある行に順に行われることを特徴とする特許請求
    の範囲第９項記載の方法。
11. 【請求項１１】前記推定伝達関数の極及びゼロ点を表示
    するステップを含むことを特徴とする特許請求の範囲第
    ９項記載の方法。
12. 【請求項１２】以下の（ａ）〜（ｌ）を設けてなる測定
    対象の推定伝達関数の極・ゼロ点決定方法： （ａ）測定対象に刺激信号を与えて応答信号を発生させ
    る； （ｂ）前記刺激信号及び応答信号を検出する； （ｃ）前記刺激信号及び応答信号から複数の周波数点に
    おける自己スペクトル及び重なりスペクトルを計算す
    る； （ｄ）前記周波数点の各々における前記重なりスペクト
    ル上のノイズレベルを測定する； （ｅ）推定伝達関数を生成する； （ｆ）前記推定伝達関数の所望の領域を強調するための
    重み付け関数を生成する； （ｇ）前記推定伝達関数に前記重み付け関数を乗算す
    る； （ｈ）前記周波数点の各々における測定された伝達関数
    と前記重み付け関数を乗算された推定伝達関数の間の誤
    差を求める； （ｉ）前記周波数点の各々において前記誤差を前記推定
    された伝達関数について微分する； （ｊ）前記周波数点の各々において前記微分された誤差
    を０と置く； （ｋ）前記推定伝達関数の係数を定める； （ｌ）前記推定伝達関数の極及びゼロ点を前記推定伝達
    関数の根として求める。
13. 【請求項１３】前記ステップ（ｂ）は前記検出された刺
    激信号及び応答信号を周波数領域に変換することを特徴
    とする特許請求の範囲第12項記載の方法。
14. 【請求項１４】前記刺激信号は前記周波数点上のランダ
    ムノイズであることを特徴とする特許請求の範囲第13項
    記載の方法。
15. 【請求項１５】ステップ（ｉ）、（ｊ）及び（ｋ）を行
    うに当って行列演算を用いることを特徴とする特許請求
    の範囲第14項記載の方法。
16. 【請求項１６】前記ステップ（ｆ）は下記のステップ
    （ｆ−１）〜（ｆ−４）を有することを特徴とする特許
    請求の範囲第15項記載の方法： （ｆ−１）前記周波数点の各々において前記測定された
    伝達関数の位相を周波数について微分する； （ｆ−２）前記周波数点の各々において前記微分された
    位相を絶対値に変換する； （ｆ−３）前記重なりスペクトル上の測定されたノイズ
    の大きい周波数点をデエンファサイズする； （ｆ−４）ステップ（ｆ−３）の結果を１に正規化す
    る。
17. 【請求項１７】前記ステップ（ｆ−１）の後に、得られ
    た微分係数を平滑化することを特徴とする特許請求の範
    囲第16項記載の方法。
18. 【請求項１８】前記ステップ（ｆ−２）の後に、前記絶
    対値を平滑化することを特徴とする特許請求の範囲第17
    項記載の方法。
19. 【請求項１９】前記ステップ（ｆ−２）と前記微分係数
    の平滑化の間に前記周波数点の各々において前記絶対値
    の平方根を形成することを特徴とする特許請求の範囲第
    18項記載の方法。
20. 【請求項２０】前記重なりスペクトル上のノイズは前記
    周波数点の各々において複数の測定についてのコヒーレ
    ンスとして測定されることを特徴とする特許請求の範囲
    第19項記載の方法。
21. 【請求項２１】前記ステップ（ｆ）は下記のステップ
    （ｆ−１′）〜（ｆ−３′）を有することを特徴とする
    特許請求の範囲15項記載の方法： （ｆ−１′）測定された伝達関数中のピークを検出す
    る； （ｆ−２′）前記検出されたピークの各々に輪郭を形成
    する； （ｆ−３′）前記周波数点に渡って前記輪郭を足し合わ
    せる。
22. 【請求項２２】前記足し合わせられた輪郭を最大の輪郭
    値に正規化することを特徴とする特許請求の範囲第21項
    記載の方法。
23. 【請求項２３】前記輪郭の各々は放物線であることを特
    徴とする特許請求の範囲第22項記載の方法。
24. 【請求項２４】前記放物線の各々は前記測定された伝達
    関数の検出されたピークに中心を合わせられていること
    を特徴とする特許請求の範囲第23項記載の方法。
25. 【請求項２５】前記放物線の各々はほぼ同じ面積を有し
    ていることを特徴とする特許請求の範囲第24項記載の方
    法。
26. 【請求項２６】前記放物線の各々は当該放物線に対応す
    るピークの幅に比例する幅を有することを特徴とする特
    許請求の範囲第25項記載の方法。
27. 【請求項２７】前記放物線の各々はその幅に反比例する
    高さを有することを特徴とする特許請求の範囲第26項記
    載の方法。
28. 【請求項２８】下記のステップ（ｍ）〜（ｐ）を有する
    ことを特徴とする特許請求の範囲第27項記載の方法： （ｍ）測定された伝達関数の二乗平均値を求める； （ｎ）前記周波数点の各々において測定された伝達関数
    の絶対値の二乗を計算する； （ｏ）前記周波数点の各々において前記絶対値の二乗を
    前記二乗平均値で除算する； （ｐ）前記除算された絶対値の二乗を正規化された放物
    線の和と組み合わせる。
29. 【請求項２９】以下の（ａ）〜（ｊ）を設けてなる測定
    対象の推定伝達関数の極・ゼロ点決定方法： （ａ）測定対象に刺激信号を与えて応答信号を発生させ
    る； （ｂ）前記刺激信号及び応答信号を検出する； （ｃ）前記刺激信号及び応答信号から複数の周波数点に
    おける自己スペクトル及び重なりスペクトルを計算す
    る； （ｄ）推定伝達関数を２つのチェビシェフ多項式の有理
    式として生成する； （ｅ）前記周波数点の各々における推定伝達関数と測定
    された伝達関数の間の誤差を求める； （ｆ）前記周波数点の各々において前記誤差を前記推定
    された伝達関数について微分する； （ｇ）前記周波数点の各々において前記微分された誤差
    を０と置く； （ｈ）前記推定伝達関数の係数を定める； （ｉ）前記推定伝達関数を構成するチェビシェフ多項式
    を通常の多項式に変換する； （ｊ）前記推定伝達関数の極及びゼロ点を前記推定伝達
    関数の根として求める。
30. 【請求項３０】前記ステップ（ｂ）は前記検出された刺
    激信号及び応答信号を周波数領域に変換することを特徴
    とする特許請求の範囲第29項記載の方法。
31. 【請求項３１】前記刺激信号は前記周波数点上のランダ
    ムノイズであることを特徴とする特許請求の範囲第30項
    記載の方法。
32. 【請求項３２】前記ステップ（ｆ）、（ｇ）及び（ｈ）
    を行うに当って行列演算を用いることを特徴とする特許
    請求の範囲第31項記載の方法。