JP5852935B2 - 伝達関数推定装置、伝達関数推定方法、および、伝達関数推定プログラム - Google Patents

Description

本発明は、伝達関数推定装置、伝達関数推定方法、および、伝達関数推定プログラムに関する。

従来、被測定物を複数の分系に分け、分系の伝達関数から被測定物全体（以下、「全系」という。）の伝達関数を推定する伝達関数合成法が公知である。

長松昭男、大熊政明「部分構造合成法」、培風館、ｐｐ．９１〜９５（１９９１）

ところで、例えば複雑な装置を被測定物とし、ある力を入力したときの加速度に係る周波数応答特性には、共振成分や反共振成分が多い。共振成分や反共振成分が多いと、カーブフィッティング等により分系の伝達関数を算出するための演算負荷が過大であり、分系の伝達関数を算出することが実質的にできないことがある。従来の伝達関数合成法では、分系の伝達関数が算出できないと、全系の伝達関数を推定することができなかった。また、分系の伝達関数を算出すべく、一部の共振成分や反共振成分のみに着目してカーブフィッティングを行うと、着目した周波数成分以外の情報が欠落するため、分系の伝達関数を算出して全系の伝達関数を推定した場合、情報が欠落した分だけ推定の精度が低くなってしまう。
本発明は、上述の課題に鑑みてなされたものであり、その目的は、演算負荷を低減し、精度よく全系の伝達関数を演算可能な伝達関数推定装置、伝達関数推定方法、および、伝達関数推定プログラムを提供することにある。

本発明の伝達関数推定装置は、複数の分系から構成される被測定物における伝達関数を推定するものである。伝達関数推定装置は、推定式取得手段と、周波数応答特性取得手段と、四則演算判断手段と、乗除演算手段と、変換手段と、加減演算手段と、再変換手段と、を備える。

推定式取得手段は、分系の数および接続方法に応じて決定される被測定物の伝達関数を推定する推定式を取得する。周波数応答特性取得手段は、分系の伝達関数に係るゲインおよび位相を含む周波数領域のデータである周波数応答特性を取得する。
四則演算判断手段は、推定式を構成する複数の演算のうちの一の演算について、加減算か乗除算かを判断する。
四則演算判断手段により乗除算であると判断された場合、乗除演算手段は、周波数領域にて周波数応答特性の乗除算を行う。

四則演算判断手段により加減算であると判断された場合、変換手段は、加減算すべき周波数応答特性を周波数領域から時間領域に変換した時間領域のデータである時間応答特性を算出する。また、加減演算手段は、時間領域にて時間応答特性の加減算を行う。そして再変換手段は、加減演算手段により算出された加減算結果を周波数領域に再変換する。
推定式を構成する全ての演算が終了したか否かを判断し、終了していないと判断された場合、演算が終了していない残りの他の演算について、四則演算判断手段の処理を行うとともに、四則演算判断手段の判断結果に応じた乗除演算手段の処理、または、変換手段、加減算手段および再変換手段の処理を行う。

本発明では、被測定物の全系の伝達関数の推定式の演算に周波数応答特性を用いている。周波数応答特性の乗除算は、周波数領域で演算可能なので、周波数領域で演算する。一方、周波数応答特性の加減算は、周波数領域での演算が困難なので、周波数応答特性を時間応答特性に変換し、時間領域にて演算している。これにより、被測定物の分系の伝達関数を算出することなく、全系の伝達関数の推定式を演算することができる。分系の伝達関数を算出する必要がないので、演算負荷を低減することができる。また、被測定物における周波数応答特性をそのまま用いているので、例えば一部の共振成分のみに着目してカーブフィッティングする場合とは異なり、情報の欠落がないため全系の伝達関数を精度よく推定することができる。これにより、被測定物の構造が複雑である場合でも、振動や音などを精度よくシミュレーション可能となる。

なお、乗除演算手段にて演算に用いられる「周波数応答特性」および加減演算手段にて演算に用いるべく変換手段にて時間領域に変換される「周波数応答特性」には、周波数領域にて乗除算がなされた演算結果、および、時間領域にて加減算された後に周波数領域に再変換された演算結果を含むものとする。

また、上記と同様の推定式取得ステップと、周波数応答特性取得ステップと、四則演算判断ステップと、乗除演算ステップと、変換ステップと、加減演算ステップと、再変換ステップと、を備える伝達関数の推定方法としてもよい。さらにまた、上記と同様の推定式取得手段、周波数応答特性取得手段、四則演算判断手段、乗除演算手段、変換手段、加減演算手段、および、再変換手段としてコンピュータを機能させる伝達関数推定プログラムとしてもよい。
このような伝達関数の推定方法および伝達関数推定プログラムであっても、上記伝達関数推定装置と同様の効果を奏する。

本発明の第１実施形態による伝達関数推定装置の構成を示すブロック図である。 本発明の第１実施形態による被測定物に入力された力を説明する図である。 本発明の第１実施形態による被測定物から出力された加速度を説明する図である。 本発明の第１実施形態による被測定物を説明する説明図である。 本発明の第１実施形態による周波数応答特性（Ｇ₂₁）を説明する図である。 本発明の第１実施形態による周波数応答特性（Ｇ₂₂）を説明する図である。 本発明の第１実施形態による周波数応答特性（Ｇ₃₃）を説明する図である。 本発明の第１実施形態による周波数応答特性（Ｇ₄₃）を説明する図である。 本発明の第１実施形態による伝達関数推定方法を説明するフローチャートである。 本発明の第１実施形態による周波数応答特性のゲインおよび位相を説明する図である。 本発明の第１実施形態による周波数応答特性の周波数領域から時間領域への変換を説明する図である。 本発明の第１実施形態による周波数応答特性および伝達関数を説明する図である。 本発明の第１実施形態による時間領域での加算を説明する図である。 本発明の第１実施形態による時間領域における加減算結果および伝達関数同士の加減算結果を説明する図である。 本発明の第１実施形態による時間応答特性（ｇ₂₂）を説明する図である。 本発明の第１実施形態による時間応答特性（ｇ₃₃）を説明する図である。 本発明の第１実施形態による時間領域における加算結果を説明する図である。 本発明の第１実施形態による時間領域における加算結果の時間領域から周波数領域への再変換を説明する図である。 本発明の第１実施形態による周波数領域における乗算を説明する図である。 本発明の第１実施形態による周波数領域における除算を説明する図である。 本発明の第２実施形態による仮想入力波形であるステップ応答を説明する図である。 本発明の第２実施形態によるステップ応答のゲインおよび位相を説明する図である。 本発明の第２実施形態による周波数応答特性のゲインおよび位相を説明する図である。 本発明の第２実施形態による周波数応答特性の仮想入力波形による補正を説明する図である。 本発明の第２実施形態による時間応答特性を説明する図である。 本発明の第２実施形態による加減算結果を周波数領域に変換したゲインおよび位相を説明する図である。 本発明の第２実施形態による周波数領域に変換された加減算結果の仮想入力波形による再補正を説明する図である。 単位インパルス応答、ステップ応答、ランプ応答のゲインおよび位相を説明する図である。

以下、本発明による伝達関数推定装置、プログラム、および、伝達関数推定方法を図面に基づいて説明する。
（第１実施形態）
図１に示すように、本実施形態の伝達関数推定装置１０は、被測定物３０の伝達関数を推定するものであって、測定装置２０から各種データを取得可能に構成されている。本実施形態では被測定物３０の振動特性を取得する場合を例に説明を行う。
測定装置２０は、インパルスハンマ２１、加速度ピックアップ２２、チャージアンプ２３、および、ＦＦＴアナライザ２４を備える。

インパルスハンマ２１は、被測定物３０の入力点Ｐ１に打撃力を入力する。入力され打撃力（図２参照）は、ＦＦＴアナライザ２４に取得される。
加速度ピックアップ２２は、被測定物３０において、入力点Ｐ１に入力された力が伝達されて出力点Ｐ４から出力される加速度データを取得する。加速度ピックアップ２２にて取得された加速度信号は、チャージアンプ２３により増幅される。増幅された加速度データは、図３に示す如くである。増幅された加速度データは、ＦＦＴアナライザ２４にて高速フーリエ変換（Fast Fourier Transform、以下「ＦＦＴ変換」という。）により、周波数領域のデータに変換される。周波数領域のデータである周波数応答特性Ｇ_mnは、所定の周波数（例えば１Ｈｚ）毎に、周波数ｆ_i、ゲインＧ（ｆ_i）、および、位相∠Ｇ（ｆ_i）が関連付けられた行列データである。ＦＦＴアナライザ２４にて算出された周波数応答特性Ｇ_mnは、伝達関数推定装置１０に出力される。

図４に示すように、本実施形態の被測定物３０は、２つの分系３１、３２から構成される。第１分系３１は、入力点Ｐ１および出力点Ｐ２を有する。第２分系３２は、入力点Ｐ３および出力点Ｐ４を有する。第１分系３１の出力点Ｐ２と第２分系３２の入力点Ｐ３とは、剛結合により接続される。これにより、第１分系３１の入力点Ｐ１に入力された力Ｆ₁は、第１分系３１および第２分系３２に伝達され、出力点Ｐ４から加速度Ａ₄として出力される。
なお、第１分系３１の出力点Ｐ２および第２分系３２の入力点Ｐ３は、「接続点」に対応している。

ここで、被測定物３０の全系の伝達関数Ｈ_{41_sim}を推定する推定式について、図４を参照して説明する。以下、点Ｐｎに加わる力をＦ_n、点Ｐｎの加速度をＡ_n、点Ｐｎに加わる力からＰｍの加速度までの伝達関数をＨ_mnと記載する。

各点Ｐ１〜Ｐ４の加速度Ａ₁〜Ａ₄は、以下の式（１）〜（４）で表される。
Ａ₁＝Ｈ₁₁Ｆ₁＋Ｈ₁₂Ｆ₂ ・・・（１）
Ａ₂＝Ｈ₂₁Ｆ₁＋Ｈ₂₂Ｆ₂ ・・・（２）
Ａ₃＝Ｈ₃₃Ｆ₃ ・・・（３）
Ａ₄＝Ｈ₄₃Ｆ₃ ・・・（４）

また、本実施形態では、第１分系３１の出力点Ｐ２と第２分系３２の入力点Ｐ３とは剛結合にて接続されているので、以下の式（５）が成り立つ。
Ａ₂＝Ａ₃ ・・・（５）

さらにまた、内力のつり合いから、以下の式（６）が成り立つ。
Ｆ₂＋Ｆ₃＝０ ・・・（６）

上記式（１）〜（６）を連立して解くと、被測定物３０の全系の伝達関数Ｈ_{41_sim}を推定する推定式は、式（７）のように表される。

なお、全系の伝達関数を推定する推定式は、分系の数や接続方法に応じて導出可能であるので、推定したい被測定物３０に応じて、伝達関数推定装置１０に予め格納しておく。

式（７）に示すように、入力点Ｐ１に入力された力Ｆ１から出力点Ｐ４の加速度までの伝達関数Ｈ₄₁を算出するためには、分系の伝達関数として、伝達関数Ｈ₂₁、Ｈ₂₂、Ｈ₃₃、Ｈ₄₃が必要である。そこで本実施形態では、図１に示す測定装置２０を用い、伝達関数Ｈ₂₁、Ｈ₂₂、Ｈ₃₃、Ｈ₄₃に係る周波数応答特性Ｇ₂₁、Ｇ₂₂、Ｇ₃₃、Ｇ₄₃を取得する。なお、本実施形態では、測定装置２０を用いて周波数応答特性Ｇ₂₁、Ｇ₂₂、Ｇ₃₃、Ｇ₄₃を取得しているが、例えば有限要素法などを用い、内部的に取得するように構成してもよい。

図４（ｃ）に示すように、入力点Ｐ１にある力を加えたときの出力点Ｐ２の加速度Ａ₂を取得し、ＦＦＴアナライザ２４にてＦＦＴ変換した周波数応答特性Ｇ₂₁を取得する（図５参照）。また、出力点Ｐ２にある力を加えたときの出力点Ｐ２の加速度Ａ₂を取得し、ＦＦＴアナライザ２４にてＦＦＴ変換した周波数応答特性Ｇ₂₂を取得する（図６参照）。また、入力点Ｐ３にある力を加えたときの入力点Ｐ３の加速度Ａ₃を取得し、ＦＦＴアナライザ２４にてＦＦＴ変換した周波数応答特性Ｇ₃₃を取得する（図７参照）。また、入力点Ｐ３にある力を加えたときの出力点Ｐ４の加速度Ａ₄を取得し、ＦＦＴアナライザ２４にてＦＦＴ変換した周波数応答特性Ｇ₄₃を取得する（図８参照）。

なお、周波数応答特性Ｇ₂₁、Ｇ₄₃が「分系の入力点に入力された力が、当該分系における出力点へ伝達される物理量に係る周波数応答特性である分系内周波数応答特性」に対応し、周波数応答特性Ｇ₂₂、Ｇ₃₃が「他の分系と接続される接続点に入力された力が、当該接続点に伝達される物理量に係る接続点周波数応答特性」に対応する。なお、本実施形態における「物理量」は、加速度Ａ₂、Ａ₃、Ａ₄である。
なお、周波数応答特性Ｇ₂₁、Ｇ₂₂、Ｇ₃₃、Ｇ₄₃には、上述の通り、所定の周波数毎にゲインおよび位相が関連付けられているが、図５〜図８では位相については記載を省略し、周波数とゲインとの関係についてのみを図示している。

ところで、従来、伝達関数Ｈ_{41_sim}の推定式（式（７））を演算するためには、得られた周波数応答特性Ｇ₂₁、Ｇ₂₂、Ｇ₃₃、Ｇ₄₃から、カーブフィッティング等により伝達関数Ｈ₂₁、Ｈ₂₂、Ｈ₃₃、Ｈ₄₃を算出する必要があった。しかしながら、周波数応答特性Ｇ₂₁、Ｇ₂₂、Ｇ₃₃、Ｇ₄₃において、共振成分や反共振成分が多い場合、適切な伝達関数の次数や伝達関数の各項の初期値を設定することが難しいためにカーブフィッティングの精度が低くなったり、場合によっては演算負荷が過大であるためにカーブフィッティングができず伝達関数が算出できなかったりする。また、例えば一部の周波数における共振成分のみに着目して伝達関数を算出する場合、得られた伝達関数は、実際の周波数応答特性との乖離が大きいという問題点がある。

そこで本実施形態では、周波数応答特性Ｇ₂₁、Ｇ₂₂、Ｇ₃₃、Ｇ₄₃から伝達関数Ｈ₂₁、Ｈ₂₂、Ｈ₃₃、Ｈ₄₃を算出せず、式（７）中の伝達関数Ｈ₂₁、Ｈ₂₂、Ｈ₃₃、Ｈ₄₃に替えて、周波数応答特性Ｇ₂₁、Ｇ₂₂、Ｇ₃₃、Ｇ₄₃を用いて全系の伝達関数Ｈ_{41_sim}の推定式を演算している。本実施形態における全系の伝達関数Ｈ_{41_sim}の推定方法を図９に示すフローチャートに基づいて説明する。

最初のステップＳ１０１（以下、「ステップ」を省略し、単に記号「Ｓ」で示す。）では、予め格納されている全系の伝達関数Ｈ_{41_sim}を推定する推定式（式（７））を取得する。
Ｓ１０２では、Ｓ１０１で取得された推定式（式（７））の演算に必要な周波数応答特性Ｇ₂₁、Ｇ₂₂、Ｇ₃₃、Ｇ₄₃を、打撃試験、振動試験、または、有限要素法等により取得する。なお、周波数応答特性Ｇ₂₁、Ｇ₂₂、Ｇ₃₃、Ｇ₄₃が伝達関数推定装置１０の記憶部等に予め格納されており、演算部が内部的に取得するように構成してもよい。

Ｓ１０３では、伝達関数Ｈ_{41_sim}を推定する推定式（式（７））の各演算において、周波数応答特性の加算または減算を行うか否かを判断する。加算または減算を行わないと判断された場合（Ｓ１０３：ＮＯ）、Ｓ１０７へ移行する。加算または減算を行うと判断された場合（Ｓ１０３：ＹＥＳ）、Ｓ１０４へ移行する。

Ｓ１０４では、周波数領域においては、周波数応答特性同士の加減算は困難なので、加減算する周波数応答特性を時間応答特性に変換する。例えば式（７）では、周波数応答特性Ｇ₂₂、Ｇ₃₃を加算しているので、周波数応答特性Ｇ₂₂、Ｇ₃₃を時間応答特性ｇ₂₂、ｇ₃₃に変換する。周波数応答特性を時間応答特性に変換する変換方法については、後述する。
Ｓ１０５では、Ｓ１０４にて変換された時間応答特性ｇ₂₂、ｇ₃₃を用い、時間領域において加減算する。
Ｓ１０６では、Ｓ１０５で得られた加減算結果であるｇ₂₂＋ｇ₃₃をフーリエ変換により周波数領域に再変換し、Ｇ₂₂＋Ｇ₃₃を得る。

周波数応答特性の加算または減算を行わないと判断された場合（Ｓ１０３：ＮＯ）、すなわち、乗算または除算を行う場合に移行するＳ１０７では、周波数領域において周波数応答特の乗除算を行う。周波数領域における周波数応答特性同士の乗除算は、ゲインの加減算として取り扱うことができるので、容易に演算可能である。
Ｓ１０８では、Ｓ１０２で取得した推定式について、全ての計算が終わったか否かを判断する。全ての計算が終わっていないと判断された場合（Ｓ１０８：ＮＯ）、Ｓ１０３へ戻る。全ての計算が終わったと判断された場合（Ｓ１０８：ＹＥＳ）、本処理を終了する。

ここで、Ｓ１０４における周波数応答特性を周波数領域から時間領域に変換する変換方法について、図１０および図１１に基づいて説明する。
図１０（ａ）に示すように、ある周波数応答特性Ｇ_mnにおいて、周波数ｆ₁＝１ＨｚのときのゲインがＧ（ｆ₁）、周波数ｆ₂＝２ＨｚのときのゲインがＧ（ｆ₂）、周波数ｆ₃＝３ＨｚのときのゲインがＧ（ｆ₃）、といった具合に、１Ｈｚ毎に周波数ｆ_iとゲインＧ（ｆ_i）とが関連づけられている。
また、図１０（ｂ）に示すように、ある周波数応答特性Ｇ_mnにおいて、周波数ｆ₁＝１Ｈｚのときの位相が∠Ｇ（ｆ₁）、周波数ｆ₂＝２Ｈｚのときの位相が∠Ｇ（ｆ₂）、周波数ｆ₃＝３Ｈｚのときの位相が∠Ｇ（ｆ₃）といった具合に、１Ｈｚ毎に周波数ｆ_iと位相∠Ｇ（ｆ_i）とが関連づけられている。

本実施形態では、入力された入力波形が単位インパルス応答であるものと仮定し、周波数応答特性Ｇ_mnを時間応答特性ｇ_mnに変換する。すなわち、本実施形態では、仮想入力波形を単位インパルス応答としている、ということである。
仮想入力波形として単位インパルス応答を用いた場合、時間領域に変換された時間応答特性は、全周波数帯の情報を位相遅れなく均等に有しており、仮想入力波形の周波数特性に依存しない。また、単位インパルス応答の伝達関数は、Ｈ（ｓ）＝１であって、図２８中の実線Ｌ１で示すように、全周波数帯において、ゲインが０ｄＢ、かつ、位相が０ｄｅｇであるので、仮想入力波形の周波数特性によるゲインおよび位相の補正を行う必要がなく、周波数応答特性Ｇ_mnのゲインＧ（ｆ_n）および位相∠Ｇ（ｆ_n）をそのまま演算に用いることができるので、演算の手間が少なくすむ。

ここで、周波数ｆ₁＝１ＨｚのときのゲインＧ（ｆ₁）および位相∠Ｇ（ｆ₁）に基づいて時間応答ｇ（１）を算出すると、以下の式（８）のようになる。時間応答ｇ（１）は、図１１（ａ）に示す如くである。
ｇ（１）＝｜Ｇ（ｆ₁）｜×ｃｏｓ（２πｆ₁ｔ＋∠Ｇ（ｆ₁）） ・・・（８）

同様に、周波数ｆ₂＝２Ｈｚのとき、および、ｆ₃＝３Ｈｚのときの時間応答ｇ（２）およびｇ（３）は、以下の式（９）、（１０）のように表される。また、時間応答ｇ（２）は図１１（ｂ）、時間応答ｇ（３）は図１１（ｃ）に示す如くである。
ｇ（２）＝｜Ｇ（ｆ₂）｜×ｃｏｓ（２πｆ₂ｔ＋∠Ｇ（ｆ₂）） ・・・（９）
ｇ（３）＝｜Ｇ（ｆ₃）｜×ｃｏｓ（２πｆ₃ｔ＋∠Ｇ（ｆ₃）） ・・・（１０）
なお、周波数ｆｉのときの時間応答ｇ（ｉ）は、式（１１）のように表される。
ｇ（ｉ）＝｜Ｇ（ｆｉ）｜×ｃｏｓ（２πｆ_iｔ＋∠Ｇ（ｆ_i）） ・・・（１１）

各周波数における時間応答は、上述のように三角関数で表されており、容易に加算することができる。例えば、ｇ（１）＋ｇ（２）＋ｇ（３）は、図１１（ｄ）に示す如くとなる。
このように、所定の周波数帯において、各周波数における時間応答を加算することにより、周波数領域の周波数応答特性Ｇ_mnを時間領域に変換した時間応答特性ｇ_mnを算出することができる。例えば、周波数間隔を１、周波数帯を１〜Ｎとしたとき、周波数応答特性Ｇ_mnを時間領域に変換した時間応答特性ｇ_mnは、以下の式（１２）のように表すことができる。

なお、周波数間隔は、ＦＦＴアナライザ２４の分解能等に応じて適宜設定可能である。また、周波数帯は、シミュレーションしたい周波数帯を適宜設定可能である。例えば周波数間隔を０．１Ｈｚ、シミュレーションする周波数帯の幅が数千Ｈｚであれば、加算する三角関数の数は数万のオーダーである。数万のオーダーの三角関数の加算は、現在世に普及している一般的なＰＣの能力で十分演算可能であることから、周波数応答特性Ｇ_mnを時間応答特性ｇ_mnに変換可能である。

次に、周波数応答特性Ｇ_mnから伝達関数Ｈ_mnを算出せず、周波数応答特性Ｇ_mnを時間応答特性ｇ_mnに変換して時間領域にて加減算して周波数領域に再変換した結果と、周波数応答特性Ｇ_mnから伝達関数Ｈ_mnを算出して伝達関数同士を加減算した結果と、が一致することについて、図１２〜図１４に示す周波数応答特性Ｇ_x、Ｇ_yを用いて説明する。なお、図１２および図１４において、周波数応答特性Ｇ_x、Ｇ_yについては、ゲインのみを示し、位相については省略する。

図１２（ａ）に示すように、実線で示す周波数応答特性Ｇ_xについて、カーブフィッティングにより伝達関数Ｈ_xを算出すると、破線で示す如くとなる。また、図１２（ｂ）に示すように、実線で示す周波数応答特性Ｇ_yについて、カーブフィッティングにより伝達関数Ｈ_yを算出すると、破線で示す如くとなる。

また、周波数応答特性Ｇ_xを時間領域に変換した時間応答特性ｇ_x（図１３（ａ）参照）、および、周波数応答特性Ｇ_yを時間領域に変換した時間応答特性ｇ_y（図１３（ｂ）参照）を算出する。算出された時間応答特性ｇ_x、ｇ_yを用い、時間領域にて、ｇ_x＋ｇ_xを算出する（図１３（ｃ）参照）。時間領域における時間応答特性同士の加減算は、加速度の加減算として取り扱うことができるので、容易に演算可能である。
そして、算出されたｇ_x＋ｇ_yをＦＦＴ変換により周波数領域に変換し、周波数領域におけるＧ_x＋Ｇ_yを得る。得られたＧ_x＋Ｇ_yは、図１４中に実線で示す如くである。

また、周波数応答特性Ｇ_xの伝達関数Ｈ_xと周波数応答特性Ｇ_yの伝達関数Ｈ_yとを加算したＨ_x＋Ｈ_yは、図１４中に破線で示す如くである。図１４に示すように、時間領域にて加算して周波数領域に再変換したＧ_x＋Ｇ_yと、周波数領域にて伝達関数同士を加算したＨ_x＋Ｈ_yとは、略一致している。
そこで、本実施形態では、伝達関数同士の加減算に替えて、周波数応答特性を時間応答特性に変換して加減算を行い、時間領域における加減算結果を周波数領域に再変換している。

ここで、被測定物３０の全系の伝達関数Ｈ_{41_sim}の推定式（式（７））の演算について具体的に説明する。
式（７）の分母を見ると、Ｈ₂₂とＨ₃₃を加算している（図９中のＳ１０３：ＹＥＳ）。そこで、伝達関数Ｈ₂₂およびＨ₃₃に対応する周波数応答特性であるＧ₂₂（図６参照）およびＧ₃₃（図７参照）を、時間領域に変換し、時間応答特性ｇ₂₂（図１５参照）およびｇ₃₃（図１６参照）を得る（Ｓ１０４）。また、得られた時間応答特性ｇ₂₂、ｇ₃₃を用い、ｇ₂₂＋ｇ₃₃を算出し（Ｓ１０５、図１７参照）、算出されたｇ₂₂＋ｇ₃₃をＦＦＴ変換により周波数領域に再変換し、Ｇ₂₂＋Ｇ₃₃を得る（Ｓ１０６、図１８参照）。

また、式（７）の分子を見ると、Ｈ₂₁とＨ₄₃を乗算している（Ｓ１０３：ＮＯ）。乗除算は周波数領域で演算可能なので、伝達関数Ｈ₂₁およびＨ₄₃に対応する周波数応答特性であるＧ₂₁およびＧ₄₃は、時間領域に変換せず、そのまま周波数領域にて乗算し、Ｇ₂₁×Ｇ₄₃を得る（Ｓ１０７、図１９参照）。なお、周波数領域における乗算は、ゲインの加算により算出される。

さらに式（７）では、Ｈ２１×Ｈ４３をＨ２２＋Ｈ３３で除しているので（Ｓ１０３：ＮＯ）、周波数領域にて除算し、（Ｇ₂₁×Ｇ₄₃）／（Ｇ₂₂＋Ｇ₃₃）を得る（Ｓ１０７、図２０参照）。なお、周波数領域における除算は、ゲインの減算により算出される。

このように、加減算のときは周波数応答特性を時間領域に変換した時間応答特性を用いて時間領域にて演算し、乗除算のときは周波数応答特性を用いて周波数領域で演算することにより、分系の伝達関数を求めることなく、全系の伝達関数Ｈ_{41_sim}の推定式（式（７））を演算することができる。

以上詳述したように、伝達関数推定装置１０は、複数の分系３１、３２から構成される被測定物３０における伝達関数Ｈ_{41_sim}を推定する。伝達関数推定装置１０では、以下の処理を行う。すなわち、分系３１、３２の数および接続方法に応じて決定される被測定物３０の伝達関数Ｈ_{41_sim}を推定する推定式（式（７））を取得する（図９中のＳ１０１）。また、被測定物３０の入出力に係るゲインおよび位相を含む周波数領域のデータである周波数応答特性Ｇ₂₁、Ｇ₂₂、Ｇ₃₃、Ｇ₄₃を取得する（Ｓ１０２）。

推定式（式（７））の演算において、演算毎に加減算か乗除算かを判断し（Ｓ１０３）、乗除算であると判断された場合（Ｓ１０３：ＮＯ）、周波数領域にて周波数応答特性Ｇ₂₁、Ｇ₄₃の乗除算を行う（Ｓ１０７）。加減算であると判断された場合（Ｓ１０３：ＹＥＳ）、加減算すべき周波数応答特性Ｇ₂₂、Ｇ₃₃を周波数領域から時間領域に変換した時間応答特性ｇ₂₂、ｇ₃₃を算出し（Ｓ１０４）、時間領域にて時間応答特性ｇ₂₂、ｇ₃₃の加減算を行う（Ｓ１０５）。そして、算出された加減算結果ｇ₂₂＋ｇ₃₃を周波数領域に再変換し、Ｇ₂₂＋Ｇ₃₃を得る（Ｓ１０６）。

本実施形態では、被測定物３０の全系の伝達関数Ｈ_{41_sim}の推定式である式（７）の演算に周波数応答特性Ｇ₂₁、Ｇ₂₂、Ｇ₃₃、Ｇ₄₃を用いている。周波数応答特性の乗除算は、周波数領域で演算可能なので、周波数領域で演算する。一方、周波数応答特性の加減算は、周波数領域では演算が困難なので、周波数応答特性を時間領域に変換し、時間領域にて演算している。これにより、被測定物３０の分系の伝達関数Ｈ₂₁、Ｈ₂₂、Ｈ₃₃、Ｈ₄₃を算出することなく、全系の伝達関数Ｈ_{41_sim}の推定式である式（７）を演算することができる。また、また、被測定物３０における周波数応答特性Ｇ₂₁、Ｇ₂₂、Ｇ₃₃、Ｇ₄₃をそのまま用いているので、例えば一部の共振成分のみに着目してカーブフィッティングする場合とは異なり、全系の伝達関数Ｈ_{41_sim}を精度よく推定することができる。これにより、例えば車両の電動パワーステアリング装置のように被測定物３０の構造が複雑である場合でも、振動や音などを精度よくシミュレーション可能となる。

また本実施形態では、周波数応答特性Ｇ_mnの時間領域への変換に用いる仮想入力波形として単位インパルス応答を用い、周波数応答特性Ｇ_mnのゲインＧ（ｆ_i）および位相∠Ｇ（ｆ_i）に基づく三角関数である時間応答ｇ（ｆ_i）を所定の周波数毎に算出し、算出された全ての周波数における時間応答を加算して時間応答特性ｇ_mnを算出する。本実施形態において仮想入力波形として用いられている単位インパルス応答は、ゲインおよび位相が全周波数帯で０である。したがって、時間領域への変換の際、周波数応答特性Ｇ_mnのゲインＧ（ｆ_i）および位相∠Ｇ（ｆ_i）をそのまま用いることができるので、演算が簡素であり、演算負荷を低減することができる。

本実施形態では、それぞれの分系３１、３２に入力された力が、それぞれの分系３１、３２における出力点Ｐ２、Ｐ４へ伝達される物理量に係る分系内周波数応答特性である周波数応答特性Ｇ₂₁、Ｇ₄₃を取得する。また、他の分系３２と接続される接続点である出力点Ｐ２に入力された力が、当該出力点Ｐ２に伝達される物理量に係る接続点周波数応答特性である周波数応答特性Ｇ₂₂、および、他の分系３１と接続される接続点である入力点Ｐ３に入力された力が、当該入力点Ｐ３に伝達される物理量に係る接続点周波数応答特性である周波数応答特性Ｇ₃₃を取得する。これらの周波数応答特性Ｇ₂₁、Ｇ₄₃、Ｇ₂₂、Ｇ₃₃を用いることにより、被測定物３０の全系の伝達関数Ｈ_{41_sim}の推定式（式（７））を適切に演算することができる。

なお、本実施形態では、図９中のＳ１０１が「推定式取得手段」の機能としての処理に相当し、Ｓ１０２が「周波数応答特性取得手段」の機能としての処理に相当し、Ｓ１０３が「四則演算判断手段」の機能としての処理に相当し、Ｓ１０７が「乗除演算手段」の機能としての処理に相当し、Ｓ１０４が「変換手段」の機能としての処理に相当し、Ｓ１０５が「加減演算手段」の機能としての処理に相当し、Ｓ１０６が「再変換手段」の機能としての処理に相当する。

（第２実施形態）
本発明の第２実施形態は、周波数応答特性Ｇ_mnを周波数領域から時間領域に変換する変換方法、および、時間領域から周波数領域に再変換する再変換方法のみが異なっているので、この点を中心に説明し、他の演算方法等についての説明は省略する。
第１実施形態では、周波数応答特性Ｇ_mnを周波数領域から時間領域への変換に際し、仮想入力波形を単位インパルス応答としていたが、第２実施形態では、仮想入力波形がステップ応答である点が異なっている。

仮想入力波形としてステップ応答を用いた場合、変換されたデータは、高周波振動成分が励起されない。そのため、低周波領域におけるデータの精度が必要な場合や高周波領域におけるデータの信頼度が低い場合などは、仮想入力波形としてステップ応答を用いることが好ましい。
また、仮想入力波形が単位インパルス応答以外である場合、すなわち仮想入力波形のゲインおよび位相が０でない場合、周波数領域から時間領域に変換するとき、および、時間領域から周波数領域に再変換するときに、用いる仮想入力波形自体のゲインおよび位相を考慮して、周波数応答特性を補正する必要がある。

そこで、仮想入力波形としてステップ応答を用いる場合における周波数応答特性Ｇ_mnを時間領域に変換する方法、および、周波数領域に再変換する方法について説明する。
本実施形態にて仮想入力波形として用いるステップ応答は、伝達関数がＨ（ｓ）＝１／ｓであって、波形は図２１に示す如くである。このようなステップ応答をＦＦＴ変換すると、図２２（ａ）および図２８（ａ）の実線Ｌ２に示すように、仮想入力ゲインＧ_inは傾きが−２０ｄＢ／ｄｅcであり、図２２（ｂ）および図２８（ｂ）の実線Ｌ２に示すように、仮想入力位相∠Ｇ_inは全周波数帯で−９０ｄｅｇの位相遅れとなる。

ここで、時間領域に変換する周波数応答特性Ｇ_mnのゲインをＧ（図２３（ａ）参照）、位相を∠Ｇ（図２３（ｂ）参照）とする。まず、時間領域に変換する前に、周波数領域において、周波数応答特性Ｇ_mnのゲインＧおよび位相∠Ｇを、ステップ応答の仮想入力ゲインＧ_inおよび仮想入力位相∠Ｇ_inで補正し、補正ゲインＧ_aおよび補正位相∠Ｇ_aを算出する。

具体的には、ある周波数ｆ_iにおける周波数応答特性Ｇ_mnのゲインＧをＧ（ｆ_i）、仮想入力ゲインＧ_inをＧ_in（ｆ_i）とし、式（１３）に示すように、周波数応答特性Ｇ_mnのゲインＧ（ｆ_i）にステップ応答の仮想入力ゲインＧ_in（ｆ_i）を加算して補正ゲインＧ_a（ｆ_i）を算出する（図２４（ａ）参照）。
Ｇ_a（ｆ_i）＝｜Ｇ_in（ｆ_i）＋Ｇ（ｆ_i）｜ ・・・（１３）

また、ある周波数ｆ_iにおける周波数応答特性Ｇ_mnの位相∠Ｇを∠Ｇ（ｆ_i）、仮想入力位相∠Ｇ_inを∠Ｇ_in（ｆ_i）とし、式（１４）に示すように、周波数応答特性Ｇ_mnの位相∠Ｇ（ｆ_i）にステップ応答の仮想入力位相∠Ｇ_in（ｆ_i）を加算して補正位相∠Ｇ_a（ｆ_i）を算出する（図２４（ｂ）参照）。
∠Ｇ_a（ｆ_i）＝∠Ｇ_in（ｆ_i）＋∠Ｇ（ｆ_i） ・・・（１４）

さらに、各周波数において、周波数ｆ_iのときの補正された補正ゲインＧ_a（ｆ_i）および補正位相∠Ｇ_a（ｆ_i）に基づく三角関数を時間応答ｇ（ｆ_i）として算出する。時間応答ｇ（ｆ_i）は、式（１５）のように表される
ｇ（ｆ_i）＝｜Ｇ_in（ｆ_i）＋Ｇ（ｆ_i）｜×ｃｏｓ（２πｆ_iｔ＋∠Ｇ（ｆ_i））
＝｜Ｇ_a（ｆ_i）｜×ｃｏｓ（２πｆ_iｔ＋∠Ｇ_a（ｆ_i））
・・・（１５）

そして、所定の周波数帯において、算出された各周波数における時間応答を加算することにより、仮想入力波形をステップ応答とし、周波数領域の周波数応答特性Ｇ_mnを時間領域に変換した時間応答特性ｇ_mnを算出することができる。例えば、周波数間隔を１、周波数帯を１〜Ｎとし、仮想入力波形がステップ応答である場合において、周波数応答特性Ｇ_mnを時間領域に変換した時間応答特性ｇ_mnは、以下の式（１６）のように表すことができる（図２５参照）。

算出された時間応答特性ｇ_mnは、第１実施形態で説明したように、時間領域における加減算に用いることができる。時間領域にて加減算された加減算結果ｇ_calaは、ＦＦＴ変換により周波数領域に再変換される。周波数領域に再変換された加減算結果には、再変換ゲインＧ_cala（図２６（ａ）参照）、および、再変換位相∠Ｇ_cala（図２６（ｂ）参照）が含まれる。

本実施形態では、仮想入力波形としてステップ応答を用いているので、周波数応答特性Ｇ_mnのゲインＧおよび位相∠ゲインＧにステップ応答の仮想入力ゲインＧ_inおよび仮想入力位相∠Ｇ_inが加算された補正ゲインＧ_aおよび補正位相∠Ｇ_aを用いて、時間領域に変換して加減算を行っている。そのため、加減算後に周波数領域に変換した再変換ゲインＧ_calaおよび位相∠Ｇ_calaにおいても、ステップ応答の仮想入力ゲインＧ_inおよび仮想入力位相∠Ｇ_inが加算された状態となっている。

そこで、周波数領域において、再変換ゲインＧ_calaおよび再変換位相∠Ｇ_calaを、ステップ応答の仮想入力ゲインＧ_inおよび仮想入力位相∠Ｇ_inで再補正する。

具体的には、ある周波数ｆ_iにおける再変換ゲインＧ_calaをＧ_cala（ｆ_i）、仮想入力ゲインをＧ_in（ｆ_i）とし、式（１７）に示すように、再変換ゲインＧ_cala（ｆ_i）からステップ応答の仮想入力ゲインＧ_in（ｆ_i）を減算して再補正ゲインＧ_cal（ｆ_i）を算出する（図２７（ａ）参照）。
Ｇ_cal（ｆ_i）＝Ｇ_cala（ｆ_i）−∠Ｇ_in（ｆ_i） ・・・（１７）

また、ある周波数ｆ_iにおける再変換位相∠Ｇ_calaを∠Ｇ_cala（ｆ_i）、仮想入力位相∠Ｇ_inを∠Ｇ_in（ｆ_i）とし、式（１８）に示すように、再変換位相∠Ｇ_cala（ｆ_i）からステップ応答の仮想入力位相∠Ｇ_in（ｆ_i）を減算して再補正位相∠Ｇ_cal（ｆ_i）を算出する（図２７（ｂ）参照）。
∠Ｇ_cal（ｆ_i）＝∠Ｇ_cala（ｆ_i）−∠Ｇ_in（ｆ_i） ・・・（１８）

本実施形態では、周波数応答特性Ｇ_mnの時間領域への変換に用いる仮想入力波形における仮想入力ゲインＧ_inおよび仮想入力位相∠Ｇ_inを取得し、仮想入力ゲインＧ_inおよび仮想入力位相∠Ｇ_inで周波数応答特性Ｇ_mnを補正する。また、仮想入力ゲインＧ_inで補正された補正ゲインＧ_a、および、仮想入力位相∠Ｇ_aで補正された補正位相∠Ｇ_aに基づく三角関数を時間応答ｇ（ｆ_i）として算出し、算出された各周波数における時間応答ｇ（ｆ_i）を加算することにより時間応答特性ｇ_mnを得ることができる。さらに、仮想入力ゲインＧ_inおよび仮想入力位相∠Ｇ_inで周波数領域に再変換された加減算結果である再変換ゲインＧ_calaおよび再変換位相∠Ｇ_calaを再補正する。

これにより、上記実施形態と同様の効果を奏する。また、仮想入力ゲインＧ_inおよび仮想入力位相∠Ｇ_inにより周波数応答特性Ｇ_mnを補正し、周波数領域に再変換された加減算結果を仮想入力ゲインＧ_inおよび仮想入力位相∠Ｇ_inにより再補正することにより、単位インパルス応答以外のゲインおよび位相が０ではない波形を仮想入力波形として用いることが可能となる。

なお、本実施形態では、仮想入力波形としてステップ応答を用いた例を説明した。仮想入力波形には、例えば、低周波領域におけるデータの精度をさらに高めるべく、ランプ応答を用いてもよい。ランプ応答は、伝達関数がＨ（ｓ）＝１／ｓ²であり、図２８に実線Ｌ３で示すように、ゲインの傾きが−４０ｄＢ／ｄｅc、位相が全周波数帯で−１８０ｄｅｇであって、ステップ応答よりもさらに高周波振動成分を励起させない。このように、シミュレーションしたい周波数帯等に適した仮想入力波形を適宜選択することができる。
また、実験データにより即したシミュレーションを行いたい場合には、インパルスハンマ２１により入力された実際の入力波形（図２参照）をＦＦＴ変換し、仮想入力波形として用いてもよい。

（他の実施形態）
（ア）上記実施形態では、被測定物の分系の数は２つであったが、他の実施形態では分系の数は３以上であってもよい。また、分系同士の接続形態は、直線状の接続に限らず、例えば分岐している等、どのような形態であってもよい。また、上記実施形態では、１つの分系における入力点および出力点は各１つずつであった。他の実施形態では、１つの分系に複数の入力点および出力点を有するように構成してもよい。
また、上記実施形態では、被測定物の分系同士は剛結合にて接続されていたが、他の実施形態では、剛結合に限らない。この場合、例えば上記式（５）において、所定の定数や伝達関数が係数として入ることになる。

（イ）上記実施形態では、周波数応答特性から伝達関数に変換することなく四則演算を行う方法について説明した。他の実施形態では、予め伝達関数がわかっている場合や、周波数応答特性から比較的容易に伝達関数を導出可能な場合には、伝達関数を用いて演算してもよい。例えば、伝達関数と周波数応答特性との加減算を行う場合は、伝達関数から周波数応答特性を算出し、上記実施形態の如く演算することができる。また、伝達関数のみからなる演算が含まれる場合、当該演算においては伝達関数にて演算を行ってもよい。これにより、演算負荷をさらに低減することができる。
また、各種演算に用いられる周波数応答特性は、分系の周波数応答特性（例えば、Ｇ₂₁、Ｇ₂₂、Ｇ₃₃、Ｇ₄₃等）に限らず、周波数領域にて乗除算がなされた演算結果、および、時間領域にて加減算された後に周波数領域に変換された演算結果であってもよい。

（ウ）上記実施形態では、被測定物に打撃力を加えたときの力の伝達に係る伝達関数を推定していたが、他の実施形態では、例えば熱等、他の伝達に係る伝達関数の推定に用いてもよい。
（エ）上記実施形態にて説明した伝達関数を推定する推定方法は、伝達関数推定装置１０にて実行されるものである。伝達関数推定装置１０においては、これらの処理をソフトウェアで処理する伝達関数推定プログラムを備えていてもよいし、例えば一部の処理をハードウェアにて実行するように構成されていてもよい。
以上、本発明は、上記実施形態になんら限定されるものではなく、発明の趣旨を逸脱しない範囲において種々の形態で実施可能である。

１０・・・伝達関数推定装置
２０・・・測定装置
３０・・・被測定物
３１・・・第１の分系（分系）
３２・・・第２の分系（分系）
Ｐ１・・・入力点（入力点）
Ｐ２・・・出力点（出力点、接続点）
Ｐ３・・・入力点（入力点、接続点）
Ｐ４・・・出力点（出力点）

Claims (6)

1. 複数の分系（３１、３２）から構成される被測定物（３０）における伝達関数（Ｈ_{41_sim}）を推定する伝達関数推定装置（１０）であって、
   前記分系の数および接続方法に応じて決定される前記被測定物の伝達関数を推定する推定式（式（７））を取得する推定式取得手段（Ｓ１０１）と、
   前記分系の伝達関数に係るゲインおよび位相を含む周波数領域のデータである周波数応答特性を取得する周波数応答特性取得手段（Ｓ１０２）と、
   前記推定式を構成する複数の演算のうちの一の演算について、加減算か乗除算かを判断する四則演算判断手段（Ｓ１０３）と、
   前記四則演算判断手段により乗除算であると判断された場合、前記周波数領域にて前記周波数応答特性の乗除算を行う乗除演算手段（Ｓ１０７）と、
   前記四則演算判断手段により加減算であると判断された場合、加減算すべき前記周波数応答特性を前記周波数領域から時間領域に変換した時間領域のデータである時間応答特性を算出する変換手段（Ｓ１０４）と、
   前記時間領域にて前記時間応答特性の加減算を行う加減演算手段（Ｓ１０５）と、
   前記加減演算手段により算出された加減算結果を前記周波数領域に再変換する再変換手段（Ｓ１０６）と、
   を備え、
   前記推定式を構成する全ての演算が終了したか否かを判断し、終了していないと判断された場合、演算が終了していない残りの他の演算について、前記四則演算判断手段の処理を行うとともに、前記四則演算判断手段の判断結果に応じた前記乗除演算手段の処理、または、前記変換手段、前記加減算手段および前記再変換手段の処理を行うことを特徴とする伝達関数推定装置。
2. 前記変換手段は、
   前記ゲインおよび前記位相に基づく三角関数である時間応答を所定の周波数毎に算出する三角関数算出手段と、
   前記三角関数算出手段により算出された全ての周波数における前記時間応答を加算して前記時間応答特性を算出する加算手段と、
   を有することを特徴とする請求項１に記載の伝達関数推定装置。
3. 前記周波数応答特性の前記時間領域への変換に用いる仮想入力波形における仮想入力ゲインおよび仮想入力位相を取得する仮想入力波形取得手段と、
   前記仮想入力ゲインおよび前記仮想入力位相で前記周波数応答特性を補正する補正手段と、
   前記仮想入力ゲインおよび前記仮想入力位相で前記再変換手段により前記周波数領域に再変換された前記加減算結果を再補正する再補正手段と、
   をさらに備え、
   前記三角関数算出手段は、前記補正手段により前記仮想入力ゲインで補正された補正ゲイン、および、前記補正手段により前記仮想入力位相で補正された補正位相に基づく三角関数を前記時間応答として算出することを特徴とする請求項２に記載の伝達関数推定装置。
4. 前記周波数応答特性取得手段は、
   それぞれの前記分系における入力点（Ｐ１、Ｐ３）に入力された力が、当該分系における出力点（Ｐ２、Ｐ４）へ伝達される物理量に係る分系内周波数応答特性（Ｇ₂₁、Ｇ₄₃）と、
   他の前記分系と接続される接続点（Ｐ２、Ｐ３）に入力された力が、当該接続点に伝達される物理量に係る接続点周波数応答特性（Ｇ₂₂、Ｇ₃₃）と、
   を前記周波数応答特性として取得することを特徴とする請求項１〜３のいずれか一項に記載の伝達関数推定装置。
5. 複数の分系（３１、３２）から構成される被測定物（３０）における伝達関数（Ｈ_{41_sim}）を推定する伝達関数推定方法であり、コンピュータの機能として実現される伝達関数の推定方法であって、
   前記分系の数および接続方法に応じて決定される前記被測定物の伝達関数を推定する推定式（式（７））を取得する推定式取得ステップ（Ｓ１０１）と、
   前記分系の伝達関数に係るゲインおよび位相を含む周波数領域のデータである周波数応答特性を取得する周波数応答特性取得ステップ（Ｓ１０２）と、
   前記推定式を構成する複数の演算のうちの一の演算について、加減算か乗除算かを判断する四則演算判断ステップ（Ｓ１０３）と、
   前記四則演算判断ステップにより乗除算であると判断された場合、前記周波数領域にて前記周波数応答特性の乗除算を行う乗除演算ステップ（Ｓ１０７）と、
   前記四則演算判断ステップにより加減算であると判断された場合、加減算すべき前記周波数応答特性を前記周波数領域から時間領域に変換した時間領域のデータである時間応答特性を算出する変換ステップ（Ｓ１０４）と、
   前記時間領域にて前記時間応答特性の加減算を行う加減演算ステップ（Ｓ１０５）と、
   前記加減演算ステップにより算出された加減算結果を前記周波数領域に再変換する再変換ステップ（Ｓ１０６）と、
   を備え、
   前記推定式を構成する全ての演算が終了したか否かを判断し、終了していないと判断された場合、演算が終了していない残りの他の演算について、前記四則演算判断ステップの処理を行うとともに、前記四則演算判断ステップの判断結果に応じた前記乗除演算ステップの処理、または、前記変換ステップ、前記加減算ステップおよび前記再変換ステップの処理を行うことを特徴とする伝達関数推定方法。
6. 複数の分系（３１、３２）から構成される被測定物（３０）における伝達関数（Ｈ_{41_sim}）を推定する伝達関数推定プログラムであって、
   前記分系の数および接続方法に応じて決定される前記被測定物の伝達関数を推定する推定式（式（７））を取得する推定式取得手段（Ｓ１０１）、
   前記分系の伝達関数に係るゲインおよび位相を含む周波数領域のデータである周波数応答特性を取得する周波数応答特性取得手段（Ｓ１０２）、
   前記推定式を構成する複数の演算のうちの一の演算について、加減算か乗除算かを判断する四則演算判断手段（Ｓ１０３）、
   前記四則演算判断手段により乗除算であると判断された場合、前記周波数領域にて前記周波数応答特性の乗除算を行う乗除演算手段（Ｓ１０７）、
   前記四則演算判断手段により加減算であると判断された場合、加減算すべき前記周波数応答特性を前記周波数領域から時間領域に変換した時間領域のデータである時間応答特性を算出する変換手段（Ｓ１０４）、
   前記時間領域にて前記時間応答特性の加減算を行う加減演算手段（Ｓ１０５）、
   および、前記加減演算手段により算出された加減算結果を前記周波数領域に再変換する再変換手段（Ｓ１０６）、
   としてコンピュータを機能させるプログラムであって、
   前記推定式を構成する全ての演算が終了したか否かを判断し、終了していないと判断された場合、演算が終了していない残りの他の演算について、前記四則演算判断手段の処理を行うとともに、前記四則演算判断手段の判断結果に応じた前記乗除演算手段の処理、または、前記変換手段、前記加減算手段および前記再変換手段の処理を行うことを特徴とする伝達関数推定プログラム。

Images