(57)【要約】【課題】演算負荷を低減し、精度よく全系の伝達関数を演算可能な伝達関数推定装置、伝達関数推定方法、および、伝達関数推定プログラムを提供する。【解決手段】被測定物の伝達関数を推定する推定式を取得し（Ｓ１０１）、周波数領域のデータである周波数応答特性を取得する（Ｓ１０２）。推定式の演算において、演算毎に加減算か乗除算かを判断し（Ｓ１０３）、加減算であると判断された場合（Ｓ１０３：ＹＥＳ）、加減算すべき周波数応答特性Ｇ22、Ｇ33を周波数領域から時間領域に変換した時間応答特性ｇ22、ｇ33を算出し（Ｓ１０４）、時間領域にて時間応答特性ｇ22、ｇ33の加減算を行う（Ｓ１０５）。そして、算出された加減算結果ｇ22＋ｇ33を周波数領域に再変換し、Ｇ22＋Ｇ33を得る（Ｓ１０６）。これにより、分系の伝達関数を算出することなく全系の伝達関数を推定する推定式を演算することができる。【選択図】  図９
2014021603.tif 000002
【特許請求の範囲】
【請求項１】
複数の分系（３１、３２）から構成される被測定物（３０）における伝達関数（Ｈ41_sim）を推定する伝達関数推定装置（１０）であって、
前記分系の数および接続方法に応じて決定される前記被測定物の伝達関数を推定する推定式（式（７））を取得する推定式取得手段（Ｓ１０１）と、
前記被測定物の入出力に係るゲインおよび位相を含む周波数領域のデータである周波数応答特性を取得する周波数応答特性取得手段（Ｓ１０２）と、
前記推定式の演算において、演算毎に加減算か乗除算かを判断する四則演算判断手段（Ｓ１０３）と、
前記四則演算判断手段により乗除算であると判断された場合、前記周波数領域にて前記周波数応答特性の乗除算を行う乗除演算手段（Ｓ１０７）と、
前記四則演算判断手段により加減算であると判断された場合、加減算すべき前記周波数応答特性を前記周波数領域から時間領域に変換した時間領域のデータである時間応答特性を算出する変換手段（Ｓ１０４）と、
前記時間領域にて前記時間応答特性の加減算を行う加減演算手段（Ｓ１０５）と、
前記加減演算手段により算出された加減算結果を前記周波数領域に再変換する再変換手段（Ｓ１０６）と、
を備えることを特徴とする伝達関数推定装置。
【請求項２】
前記変換手段は、
前記ゲインおよび前記位相に基づく三角関数である時間応答を所定の周波数毎に算出する三角関数算出手段と、
前記三角関数算出手段により算出された全ての周波数における前記時間応答を加算して前記時間応答特性を算出する加算手段と、
を有することを特徴とする請求項１に記載の伝達関数推定装置。
【請求項３】
前記周波数応答特性の前記時間領域への変換に用いる仮想入力波形における仮想入力ゲインおよび仮想入力位相を取得する仮想入力波形取得手段と、
前記仮想入力ゲインおよび前記仮想入力位相で前記周波数応答特性を補正する補正手段と、
前記仮想入力ゲインおよび前記仮想入力位相で前記再変換手段により前記周波数領域に再変換された前記加減算結果を再補正する再補正手段と、
をさらに備え、
前記三角関数算出手段は、前記補正手段により前記仮想入力ゲインで補正された補正ゲイン、および、前記補正手段により前記仮想入力位相で補正された補正位相に基づく三角関数を前記時間応答として算出することを特徴とする請求項２に記載の伝達関数推定装置。
【請求項４】
前記周波数応答特性取得手段は、
それぞれの前記分系における入力点（Ｐ１、Ｐ３）に入力された力が、当該分系における出力点（Ｐ２、Ｐ４）へ伝達される物理量に係る分系内周波数応答特性（Ｇ21、Ｇ43）と、
他の前記分系と接続される接続点（Ｐ２、Ｐ３）に入力された力が、当該接続点に伝達される物理量に係る接続点周波数応答特性（Ｇ22、Ｇ33）と、
を前記周波数応答特性として取得することを特徴とする請求項１〜３のいずれか一項に記載の伝達関数推定装置。
【請求項５】
複数の分系（３１、３２）から構成される被測定物（３０）における伝達関数（Ｈ41_sim）を推定する伝達関数推定方法であり、コンピュータの機能として実現される伝達関数の推定方法であって、
前記分系の数および接続方法に応じて決定される前記被測定物の伝達関数を推定する推定式（式（７））を取得する推定式取得ステップ（Ｓ１０２）と、
前記被測定物の入出力に係るゲインおよび位相を含む周波数領域のデータである周波数応答特性を取得する周波数応答特性取得ステップ（Ｓ１０１）と、
前記推定式の演算において、演算毎に加減算か乗除算かを判断する四則演算判断ステップ（Ｓ１０３）と、
前記四則演算判断ステップにより乗除算であると判断された場合、前記周波数領域にて前記周波数応答特性の乗除算を行う乗除演算ステップ（Ｓ１０７）と、
前記四則演算判断ステップにより加減算であると判断された場合、加減算すべき前記周波数応答特性を前記周波数領域から時間領域に変換した時間領域のデータである時間応答特性を算出する変換ステップ（Ｓ１０４）と、
前記時間領域にて前記時間応答特性の加減算を行う加減演算ステップ（Ｓ１０５）と、
前記加減演算ステップにより算出された加減算結果を前記周波数領域に再変換する再変換ステップ（Ｓ１０６）と、
を備えることを特徴とする伝達関数推定方法。
【請求項６】
複数の分系（３１、３２）から構成される被測定物（３０）における伝達関数（Ｈ41_sim）を推定する伝達関数推定プログラムであって、
前記分系の数および接続方法に応じて決定される前記被測定物の伝達関数を推定する推定式（式（７））を取得する推定式取得手段（Ｓ１０２）、
前記被測定物の入出力に係るゲインおよび位相を含む周波数領域のデータである周波数応答特性を取得する周波数応答特性取得手段（Ｓ１０１）、
前記推定式の演算において、演算毎に加減算か乗除算かを判断する四則演算判断手段（Ｓ１０３）、
前記四則演算判断手段により乗除算であると判断された場合、前記周波数領域にて前記周波数応答特性の乗除算を行う乗除演算手段（Ｓ１０７）、
前記四則演算判断手段により加減算であると判断された場合、加減算すべき前記周波数応答特性を前記周波数領域から時間領域に変換した時間領域のデータである時間応答特性を算出する変換手段（Ｓ１０４）、
前記時間領域にて前記時間応答特性の加減算を行う加減演算手段（Ｓ１０５）、
および、前記加減演算手段により算出された加減算結果を前記周波数領域に再変換する再変換手段（Ｓ１０６）、
としてコンピュータを機能させることを特徴とする伝達関数推定プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、伝達関数推定装置、伝達関数推定方法、および、伝達関数推定プログラムに関する。
【背景技術】
【０００２】
従来、被測定物を複数の分系に分け、分系の伝達関数から被測定物全体（以下、「全系」という。）の伝達関数を推定する伝達関数合成法が公知である。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】長松昭男、大熊政明「部分構造合成法」、培風館、ｐｐ．９１〜９５（１９９１）
【発明の概要】
【発明が解決しようとする課題】
【０００４】
ところで、例えば複雑な装置を被測定物とし、ある力を入力したときの加速度に係る周波数応答特性には、共振成分や反共振成分が多い。共振成分や反共振成分が多いと、カーブフィッティング等により分系の伝達関数を算出するための演算負荷が過大であり、分系の伝達関数を算出することが実質的にできないことがある。従来の伝達関数合成法では、分系の伝達関数が算出できないと、全系の伝達関数を推定することができなかった。また、分系の伝達関数を算出すべく、一部の共振成分や反共振成分のみに着目してカーブフィッティングを行うと、着目した周波数成分以外の情報が欠落するため、分系の伝達関数を算出して全系の伝達関数を推定した場合、情報が欠落した分だけ推定の精度が低くなってしまう。
本発明は、上述の課題に鑑みてなされたものであり、その目的は、演算負荷を低減し、精度よく全系の伝達関数を演算可能な伝達関数推定装置、伝達関数推定方法、および、伝達関数推定プログラムを提供することにある。
【課題を解決するための手段】
【０００５】
本発明の伝達関数推定装置は、複数の分系から構成される被測定物における伝達関数を推定するものである。伝達関数推定装置は、推定式取得手段と、周波数応答特性取得手段と、四則演算判断手段と、乗除演算手段と、変換手段と、加減演算手段と、再変換手段と、を備える。
【０００６】
推定式取得手段は、分系の数および接続方法に応じて決定される被測定物の伝達関数を推定する推定式を取得する。周波数応答特性取得手段は、被測定物の入出力に係るゲインおよび位相を含む周波数領域のデータである周波数応答特性を取得する。
四則演算判断手段は、推定式の演算において、演算毎に加減算か乗除算かを判断する。
四則演算判断手段により乗除算であると判断された場合、乗除演算手段は、周波数領域にて周波数応答特性の乗除算を行う。
【０００７】
四則演算判断手段により加減算であると判断された場合、変換手段は、加減算すべき周波数応答特性を周波数領域から時間領域に変換した時間領域のデータである時間応答特性を算出する。また、加減演算手段は、時間領域にて時間応答特性の加減算を行う。そして再変換手段は、加減演算手段により算出された加減算結果を周波数領域に再変換する。
【０００８】
本発明では、被測定物の全系の伝達関数の推定式の演算に周波数応答特性を用いている。周波数応答特性の乗除算は、周波数領域で演算可能なので、周波数領域で演算する。一方、周波数応答特性の加減算は、周波数領域での演算が困難なので、周波数応答特性を時間応答特性に変換し、時間領域にて演算している。これにより、被測定物の分系の伝達関数を算出することなく、全系の伝達関数の推定式を演算することができる。分系の伝達関数を算出する必要がないので、演算負荷を低減することができる。また、被測定物における周波数応答特性をそのまま用いているので、例えば一部の共振成分のみに着目してカーブフィッティングする場合とは異なり、情報の欠落がないため全系の伝達関数を精度よく推定することができる。これにより、被測定物の構造が複雑である場合でも、振動や音などを精度よくシミュレーション可能となる。
【０００９】
なお、乗除演算手段にて演算に用いられる「周波数応答特性」および加減演算手段にて演算に用いるべく変換手段にて時間領域に変換される「周波数応答特性」には、周波数領域にて乗除算がなされた演算結果、および、時間領域にて加減算された後に周波数領域に再変換された演算結果を含むものとする。
【００１０】
また、上記と同様の推定式取得ステップと、周波数応答特性取得ステップと、四則演算判断ステップと、乗除演算ステップと、変換ステップと、加減演算ステップと、再変換ステップと、を備える伝達関数の推定方法としてもよい。さらにまた、上記と同様の推定式取得手段、周波数応答特性取得手段、四則演算判断手段、乗除演算手段、変換手段、加減演算手段、および、再変換手段としてコンピュータを機能させる伝達関数推定プログラムとしてもよい。
このような伝達関数の推定方法および伝達関数推定プログラムであっても、上記伝達関数推定装置と同様の効果を奏する。
【図面の簡単な説明】
【００１１】
【図１】本発明の第１実施形態による伝達関数推定装置の構成を示すブロック図である。
【図２】本発明の第１実施形態による被測定物に入力された力を説明する図である。
【図３】本発明の第１実施形態による被測定物から出力された加速度を説明する図である。
【図４】本発明の第１実施形態による被測定物を説明する説明図である。
【図５】本発明の第１実施形態による周波数応答特性（Ｇ21）を説明する図である。
【図６】本発明の第１実施形態による周波数応答特性（Ｇ22）を説明する図である。
【図７】本発明の第１実施形態による周波数応答特性（Ｇ33）を説明する図である。
【図８】本発明の第１実施形態による周波数応答特性（Ｇ43）を説明する図である。
【図９】本発明の第１実施形態による伝達関数推定方法を説明するフローチャートである。
【図１０】本発明の第１実施形態による周波数応答特性のゲインおよび位相を説明する図である。
【図１１】本発明の第１実施形態による周波数応答特性の周波数領域から時間領域への変換を説明する図である。
【図１２】本発明の第１実施形態による周波数応答特性および伝達関数を説明する図である。
【図１３】本発明の第１実施形態による時間領域での加算を説明する図である。
【図１４】本発明の第１実施形態による時間領域における加減算結果および伝達関数同士の加減算結果を説明する図である。
【図１５】本発明の第１実施形態による時間応答特性（ｇ22）を説明する図である。
【図１６】本発明の第１実施形態による時間応答特性（ｇ33）を説明する図である。
【図１７】本発明の第１実施形態による時間領域における加算結果を説明する図である。
【図１８】本発明の第１実施形態による時間領域における加算結果の時間領域から周波数領域への再変換を説明する図である。
【図１９】本発明の第１実施形態による周波数領域における乗算を説明する図である。
【図２０】本発明の第１実施形態による周波数領域における除算を説明する図である。
【図２１】本発明の第２実施形態による仮想入力波形であるステップ応答を説明する図である。
【図２２】本発明の第２実施形態によるステップ応答のゲインおよび位相を説明する図である。
【図２３】本発明の第２実施形態による周波数応答特性のゲインおよび位相を説明する図である。
【図２４】本発明の第２実施形態による周波数応答特性の仮想入力波形による補正を説明する図である。
【図２５】本発明の第２実施形態による時間応答特性を説明する図である。
【図２６】本発明の第２実施形態による加減算結果を周波数領域に変換したゲインおよび位相を説明する図である。
【図２７】本発明の第２実施形態による周波数領域に変換された加減算結果の仮想入力波形による再補正を説明する図である。
【図２８】単位インパルス応答、ステップ応答、ランプ応答のゲインおよび位相を説明する図である。
【発明を実施するための形態】
【００１２】
以下、本発明による伝達関数推定装置、プログラム、および、伝達関数推定方法を図面に基づいて説明する。
（第１実施形態）
図１に示すように、本実施形態の伝達関数推定装置１０は、被測定物３０の伝達関数を推定するものであって、測定装置２０から各種データを取得可能に構成されている。本実施形態では被測定物３０の振動特性を取得する場合を例に説明を行う。
測定装置２０は、インパルスハンマ２１、加速度ピックアップ２２、チャージアンプ２３、および、ＦＦＴアナライザ２４を備える。
【００１３】
インパルスハンマ２１は、被測定物３０の入力点Ｐ１に打撃力を入力する。入力され打撃力（図２参照）は、ＦＦＴアナライザ２４に取得される。
加速度ピックアップ２２は、被測定物３０において、入力点Ｐ１に入力された力が伝達されて出力点Ｐ４から出力される加速度データを取得する。加速度ピックアップ２２にて取得された加速度信号は、チャージアンプ２３により増幅される。増幅された加速度データは、図３に示す如くである。増幅された加速度データは、ＦＦＴアナライザ２４にて高速フーリエ変換（Fast Fourier Transform、以下「ＦＦＴ変換」という。）により、周波数領域のデータに変換される。周波数領域のデータである周波数応答特性Ｇmnは、所定の周波数（例えば１Ｈｚ）毎に、周波数ｆi、ゲインＧ（ｆi）、および、位相∠Ｇ（ｆi）が関連付けられた行列データである。ＦＦＴアナライザ２４にて算出された周波数応答特性Ｇmnは、伝達関数推定装置１０に出力される。
【００１４】
図４に示すように、本実施形態の被測定物３０は、２つの分系３１、３２から構成される。第１分系３１は、入力点Ｐ１および出力点Ｐ２を有する。第２分系３２は、入力点Ｐ３および出力点Ｐ４を有する。第１分系３１の出力点Ｐ２と第２分系３２の入力点Ｐ３とは、剛結合により接続される。これにより、第１分系３１の入力点Ｐ１に入力された力Ｆ1は、第１分系３１および第２分系３２に伝達され、出力点Ｐ４から加速度Ａ4として出力される。
なお、第１分系３１の出力点Ｐ２および第２分系３２の入力点Ｐ３は、「接続点」に対応している。
【００１５】
ここで、被測定物３０の全系の伝達関数Ｈ41_simを推定する推定式について、図４を参照して説明する。以下、点Ｐｎに加わる力をＦn、点Ｐｎの加速度をＡn、点Ｐｎに加わる力からＰｍの加速度までの伝達関数をＨmnと記載する。
【００１６】
各点Ｐ１〜Ｐ４の加速度Ａ1〜Ａ4は、以下の式（１）〜（４）で表される。
Ａ1＝Ｈ11Ｆ1＋Ｈ12Ｆ2  ・・・（１）
Ａ2＝Ｈ21Ｆ1＋Ｈ22Ｆ2  ・・・（２）
Ａ3＝Ｈ33Ｆ3  ・・・（３）
Ａ4＝Ｈ43Ｆ3  ・・・（４）
【００１７】
また、本実施形態では、第１分系３１の出力点Ｐ２と第２分系３２の入力点Ｐ３とは剛結合にて接続されているので、以下の式（５）が成り立つ。
Ａ2＝Ａ3  ・・・（５）
【００１８】
さらにまた、内力のつり合いから、以下の式（６）が成り立つ。
Ｆ2＋Ｆ3＝０  ・・・（６）
【００１９】
上記式（１）〜（６）を連立して解くと、被測定物３０の全系の伝達関数Ｈ41_simを推定する推定式は、式（７）のように表される。
【数１】
2014021603.tif 000003
【００２０】
なお、全系の伝達関数を推定する推定式は、分系の数や接続方法に応じて導出可能であるので、推定したい被測定物３０に応じて、伝達関数推定装置１０に予め格納しておく。
【００２１】
式（７）に示すように、入力点Ｐ１に入力された力Ｆ１から出力点Ｐ４の加速度までの伝達関数Ｈ41を算出するためには、分系の伝達関数として、伝達関数Ｈ21、Ｈ22、Ｈ33、Ｈ43が必要である。そこで本実施形態では、図１に示す測定装置２０を用い、伝達関数Ｈ21、Ｈ22、Ｈ33、Ｈ43に係る周波数応答特性Ｇ21、Ｇ22、Ｇ33、Ｇ43を取得する。なお、本実施形態では、測定装置２０を用いて周波数応答特性Ｇ21、Ｇ22、Ｇ33、Ｇ43を取得しているが、例えば有限要素法などを用い、内部的に取得するように構成してもよい。
【００２２】
図４（ｃ）に示すように、入力点Ｐ１にある力を加えたときの出力点Ｐ２の加速度Ａ2を取得し、ＦＦＴアナライザ２４にてＦＦＴ変換した周波数応答特性Ｇ21を取得する（図５参照）。また、出力点Ｐ２にある力を加えたときの出力点Ｐ２の加速度Ａ2を取得し、ＦＦＴアナライザ２４にてＦＦＴ変換した周波数応答特性Ｇ22を取得する（図６参照）。また、入力点Ｐ３にある力を加えたときの入力点Ｐ３の加速度Ａ3を取得し、ＦＦＴアナライザ２４にてＦＦＴ変換した周波数応答特性Ｇ33を取得する（図７参照）。また、入力点Ｐ３にある力を加えたときの出力点Ｐ４の加速度Ａ4を取得し、ＦＦＴアナライザ２４にてＦＦＴ変換した周波数応答特性Ｇ43を取得する（図８参照）。
【００２３】
なお、周波数応答特性Ｇ21、Ｇ43が「分系の入力点に入力された力が、当該分系における出力点へ伝達される物理量に係る周波数応答特性である分系内周波数応答特性」に対応し、周波数応答特性Ｇ22、Ｇ33が「他の分系と接続される接続点に入力された力が、当該接続点に伝達される物理量に係る接続点周波数応答特性」に対応する。なお、本実施形態における「物理量」は、加速度Ａ2、Ａ3、Ａ4である。
なお、周波数応答特性Ｇ21、Ｇ22、Ｇ33、Ｇ43には、上述の通り、所定の周波数毎にゲインおよび位相が関連付けられているが、図５〜図８では位相については記載を省略し、周波数とゲインとの関係についてのみを図示している。
【００２４】
ところで、従来、伝達関数Ｈ41_simの推定式（式（７））を演算するためには、得られた周波数応答特性Ｇ21、Ｇ22、Ｇ33、Ｇ43から、カーブフィッティング等により伝達関数Ｈ21、Ｈ22、Ｈ33、Ｈ43を算出する必要があった。しかしながら、周波数応答特性Ｇ21、Ｇ22、Ｇ33、Ｇ43において、共振成分や反共振成分が多い場合、適切な伝達関数の次数や伝達関数の各項の初期値を設定することが難しいためにカーブフィッティングの精度が低くなったり、場合によっては演算負荷が過大であるためにカーブフィッティングができず伝達関数が算出できなかったりする。また、例えば一部の周波数における共振成分のみに着目して伝達関数を算出する場合、得られた伝達関数は、実際の周波数応答特性との乖離が大きいという問題点がある。
【００２５】
そこで本実施形態では、周波数応答特性Ｇ21、Ｇ22、Ｇ33、Ｇ43から伝達関数Ｈ21、Ｈ22、Ｈ33、Ｈ43を算出せず、式（７）中の伝達関数Ｈ21、Ｈ22、Ｈ33、Ｈ43に替えて、周波数応答特性Ｇ21、Ｇ22、Ｇ33、Ｇ43を用いて全系の伝達関数Ｈ41_simの推定式を演算している。本実施形態における全系の伝達関数Ｈ41_simの推定方法を図９に示すフローチャートに基づいて説明する。
【００２６】
最初のステップＳ１０１（以下、「ステップ」を省略し、単に記号「Ｓ」で示す。）では、予め格納されている全系の伝達関数Ｈ41_simを推定する推定式（式（７））を取得する。
Ｓ１０２では、Ｓ１０１で取得された推定式（式（７））の演算に必要な周波数応答特性Ｇ21、Ｇ22、Ｇ33、Ｇ43を、打撃試験、振動試験、または、有限要素法等により取得する。なお、周波数応答特性Ｇ21、Ｇ22、Ｇ33、Ｇ43が伝達関数推定装置１０の記憶部等に予め格納されており、演算部が内部的に取得するように構成してもよい。
【００２７】
Ｓ１０３では、伝達関数Ｈ41_simを推定する推定式（式（７））の各演算において、周波数応答特性の加算または減算を行うか否かを判断する。加算または減算を行わないと判断された場合（Ｓ１０３：ＮＯ）、Ｓ１０７へ移行する。加算または減算を行うと判断された場合（Ｓ１０３：ＹＥＳ）、Ｓ１０４へ移行する。
【００２８】
Ｓ１０４では、周波数領域においては、周波数応答特性同士の加減算は困難なので、加減算する周波数応答特性を時間応答特性に変換する。例えば式（７）では、周波数応答特性Ｇ22、Ｇ33を加算しているので、周波数応答特性Ｇ22、Ｇ33を時間応答特性ｇ22、ｇ33に変換する。周波数応答特性を時間応答特性に変換する変換方法については、後述する。
Ｓ１０５では、Ｓ１０４にて変換された時間応答特性ｇ22、ｇ33を用い、時間領域において加減算する。
Ｓ１０６では、Ｓ１０５で得られた加減算結果であるｇ22＋ｇ33をフーリエ変換により周波数領域に再変換し、Ｇ22＋Ｇ33を得る。
【００２９】
周波数応答特性の加算または減算を行わないと判断された場合（Ｓ１０３：ＮＯ）、すなわち、乗算または除算を行う場合に移行するＳ１０７では、周波数領域において周波数応答特の乗除算を行う。周波数領域における周波数応答特性同士の乗除算は、ゲインの加減算として取り扱うことができるので、容易に演算可能である。
Ｓ１０８では、Ｓ１０２で取得した推定式について、全ての計算が終わったか否かを判断する。全ての計算が終わっていないと判断された場合（Ｓ１０８：ＮＯ）、Ｓ１０３へ戻る。全ての計算が終わったと判断された場合（Ｓ１０８：ＹＥＳ）、本処理を終了する。
【００３０】
ここで、Ｓ１０４における周波数応答特性を周波数領域から時間領域に変換する変換方法について、図１０および図１１に基づいて説明する。
図１０（ａ）に示すように、ある周波数応答特性Ｇmnにおいて、周波数ｆ1＝１ＨｚのときのゲインがＧ（ｆ1）、周波数ｆ2＝２ＨｚのときのゲインがＧ（ｆ2）、周波数ｆ3＝３ＨｚのときのゲインがＧ（ｆ3）、といった具合に、１Ｈｚ毎に周波数ｆiとゲインＧ（ｆi）とが関連づけられている。
また、図１０（ｂ）に示すように、ある周波数応答特性Ｇmnにおいて、周波数ｆ1＝１Ｈｚのときの位相が∠Ｇ（ｆ1）、周波数ｆ2＝２Ｈｚのときの位相が∠Ｇ（ｆ2）、周波数ｆ3＝３Ｈｚのときの位相が∠Ｇ（ｆ3）といった具合に、１Ｈｚ毎に周波数ｆiと位相∠Ｇ（ｆi）とが関連づけられている。
【００３１】
本実施形態では、入力された入力波形が単位インパルス応答であるものと仮定し、周波数応答特性Ｇmnを時間応答特性ｇmnに変換する。すなわち、本実施形態では、仮想入力波形を単位インパルス応答としている、ということである。
仮想入力波形として単位インパルス応答を用いた場合、時間領域に変換された時間応答特性は、全周波数帯の情報を位相遅れなく均等に有しており、仮想入力波形の周波数特性に依存しない。また、単位インパルス応答の伝達関数は、Ｈ（ｓ）＝１であって、図２８中の実線Ｌ１で示すように、全周波数帯において、ゲインが０ｄＢ、かつ、位相が０ｄｅｇであるので、仮想入力波形の周波数特性によるゲインおよび位相の補正を行う必要がなく、周波数応答特性ＧmnのゲインＧ（ｆn）および位相∠Ｇ（ｆn）をそのまま演算に用いることができるので、演算の手間が少なくすむ。
【００３２】
ここで、周波数ｆ1＝１ＨｚのときのゲインＧ（ｆ1）および位相∠Ｇ（ｆ1）に基づいて時間応答ｇ（１）を算出すると、以下の式（８）のようになる。時間応答ｇ（１）は、図１１（ａ）に示す如くである。
ｇ（１）＝｜Ｇ（ｆ1）｜×ｃｏｓ（２πｆ1ｔ＋∠Ｇ（ｆ1））  ・・・（８）
【００３３】
同様に、周波数ｆ2＝２Ｈｚのとき、および、ｆ3＝３Ｈｚのときの時間応答ｇ（２）およびｇ（３）は、以下の式（９）、（１０）のように表される。また、時間応答ｇ（２）は図１１（ｂ）、時間応答ｇ（３）は図１１（ｃ）に示す如くである。
ｇ（２）＝｜Ｇ（ｆ2）｜×ｃｏｓ（２πｆ2ｔ＋∠Ｇ（ｆ2））  ・・・（９）
ｇ（３）＝｜Ｇ（ｆ3）｜×ｃｏｓ（２πｆ3ｔ＋∠Ｇ（ｆ3））  ・・・（１０）
なお、周波数ｆｉのときの時間応答ｇ（ｉ）は、式（１１）のように表される。
ｇ（ｉ）＝｜Ｇ（ｆｉ）｜×ｃｏｓ（２πｆiｔ＋∠Ｇ（ｆi））  ・・・（１１）
【００３４】
各周波数における時間応答は、上述のように三角関数で表されており、容易に加算することができる。例えば、ｇ（１）＋ｇ（２）＋ｇ（３）は、図１１（ｄ）に示す如くとなる。
このように、所定の周波数帯において、各周波数における時間応答を加算することにより、周波数領域の周波数応答特性Ｇmnを時間領域に変換した時間応答特性ｇmnを算出することができる。例えば、周波数間隔を１、周波数帯を１〜Ｎとしたとき、周波数応答特性Ｇmnを時間領域に変換した時間応答特性ｇmnは、以下の式（１２）のように表すことができる。
【００３５】
【数２】
2014021603.tif 000004
【００３６】
なお、周波数間隔は、ＦＦＴアナライザ２４の分解能等に応じて適宜設定可能である。また、周波数帯は、シミュレーションしたい周波数帯を適宜設定可能である。例えば周波数間隔を０．１Ｈｚ、シミュレーションする周波数帯の幅が数千Ｈｚであれば、加算する三角関数の数は数万のオーダーである。数万のオーダーの三角関数の加算は、現在世に普及している一般的なＰＣの能力で十分演算可能であることから、周波数応答特性Ｇ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の推定式（式（７））の演算について具体的に説明する。
式（７）の分母を見ると、Ｈ22とＨ33を加算している（図９中のＳ１０３：ＹＥＳ）。そこで、伝達関数Ｈ22およびＨ33に対応する周波数応答特性であるＧ22（図６参照）およびＧ33（図７参照）を、時間領域に変換し、時間応答特性ｇ22（図１５参照）およびｇ33（図１６参照）を得る（Ｓ１０４）。また、得られた時間応答特性ｇ22、ｇ33を用い、ｇ22＋ｇ33を算出し（Ｓ１０５、図１７参照）、算出されたｇ22＋ｇ33をＦＦＴ変換により周波数領域に再変換し、Ｇ22＋Ｇ33を得る（Ｓ１０６、図１８参照）。
【００４２】
また、式（７）の分子を見ると、Ｈ21とＨ43を乗算している（Ｓ１０３：ＮＯ）。乗除算は周波数領域で演算可能なので、伝達関数Ｈ21およびＨ43に対応する周波数応答特性であるＧ21およびＧ43は、時間領域に変換せず、そのまま周波数領域にて乗算し、Ｇ21×Ｇ43を得る（Ｓ１０７、図１９参照）。なお、周波数領域における乗算は、ゲインの加算により算出される。
【００４３】
さらに式（７）では、Ｈ２１×Ｈ４３をＨ２２＋Ｈ３３で除しているので（Ｓ１０３：ＮＯ）、周波数領域にて除算し、（Ｇ21×Ｇ43）／（Ｇ22＋Ｇ33）を得る（Ｓ１０７、図２０参照）。なお、周波数領域における除算は、ゲインの減算により算出される。
【００４４】
このように、加減算のときは周波数応答特性を時間領域に変換した時間応答特性を用いて時間領域にて演算し、乗除算のときは周波数応答特性を用いて周波数領域で演算することにより、分系の伝達関数を求めることなく、全系の伝達関数Ｈ41_simの推定式（式（７））を演算することができる。
【００４５】
以上詳述したように、伝達関数推定装置１０は、複数の分系３１、３２から構成される被測定物３０における伝達関数Ｈ41_simを推定する。伝達関数推定装置１０では、以下の処理を行う。すなわち、分系３１、３２の数および接続方法に応じて決定される被測定物３０の伝達関数Ｈ41_simを推定する推定式（式（７））を取得する（図９中のＳ１０１）。また、被測定物３０の入出力に係るゲインおよび位相を含む周波数領域のデータである周波数応答特性Ｇ21、Ｇ22、Ｇ33、Ｇ43を取得する（Ｓ１０２）。
【００４６】
推定式（式（７））の演算において、演算毎に加減算か乗除算かを判断し（Ｓ１０３）、乗除算であると判断された場合（Ｓ１０３：ＮＯ）、周波数領域にて周波数応答特性Ｇ21、Ｇ43の乗除算を行う（Ｓ１０７）。加減算であると判断された場合（Ｓ１０３：ＹＥＳ）、加減算すべき周波数応答特性Ｇ22、Ｇ33を周波数領域から時間領域に変換した時間応答特性ｇ22、ｇ33を算出し（Ｓ１０４）、時間領域にて時間応答特性ｇ22、ｇ33の加減算を行う（Ｓ１０５）。そして、算出された加減算結果ｇ22＋ｇ33を周波数領域に再変換し、Ｇ22＋Ｇ33を得る（Ｓ１０６）。
【００４７】
本実施形態では、被測定物３０の全系の伝達関数Ｈ41_simの推定式である式（７）の演算に周波数応答特性Ｇ21、Ｇ22、Ｇ33、Ｇ43を用いている。周波数応答特性の乗除算は、周波数領域で演算可能なので、周波数領域で演算する。一方、周波数応答特性の加減算は、周波数領域では演算が困難なので、周波数応答特性を時間領域に変換し、時間領域にて演算している。これにより、被測定物３０の分系の伝達関数Ｈ21、Ｈ22、Ｈ33、Ｈ43を算出することなく、全系の伝達関数Ｈ41_simの推定式である式（７）を演算することができる。また、また、被測定物３０における周波数応答特性Ｇ21、Ｇ22、Ｇ33、Ｇ43をそのまま用いているので、例えば一部の共振成分のみに着目してカーブフィッティングする場合とは異なり、全系の伝達関数Ｈ41_simを精度よく推定することができる。これにより、例えば車両の電動パワーステアリング装置のように被測定物３０の構造が複雑である場合でも、振動や音などを精度よくシミュレーション可能となる。
【００４８】
また本実施形態では、周波数応答特性Ｇmnの時間領域への変換に用いる仮想入力波形として単位インパルス応答を用い、周波数応答特性ＧmnのゲインＧ（ｆi）および位相∠Ｇ（ｆi）に基づく三角関数である時間応答ｇ（ｆi）を所定の周波数毎に算出し、算出された全ての周波数における時間応答を加算して時間応答特性ｇmnを算出する。本実施形態において仮想入力波形として用いられている単位インパルス応答は、ゲインおよび位相が全周波数帯で０である。したがって、時間領域への変換の際、周波数応答特性ＧmnのゲインＧ（ｆi）および位相∠Ｇ（ｆi）をそのまま用いることができるので、演算が簡素であり、演算負荷を低減することができる。
【００４９】
本実施形態では、それぞれの分系３１、３２に入力された力が、それぞれの分系３１、３２における出力点Ｐ２、Ｐ４へ伝達される物理量に係る分系内周波数応答特性である周波数応答特性Ｇ21、Ｇ43を取得する。また、他の分系３２と接続される接続点である出力点Ｐ２に入力された力が、当該出力点Ｐ２に伝達される物理量に係る接続点周波数応答特性である周波数応答特性Ｇ22、および、他の分系３１と接続される接続点である入力点Ｐ３に入力された力が、当該入力点Ｐ３に伝達される物理量に係る接続点周波数応答特性である周波数応答特性Ｇ33を取得する。これらの周波数応答特性Ｇ21、Ｇ43、Ｇ22、Ｇ33を用いることにより、被測定物３０の全系の伝達関数Ｈ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は、以下の式（１６）のように表すことができる（図２５参照）。
【００５９】
【数３】
2014021603.tif 000005
【００６０】
算出された時間応答特性ｇ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により再補正することにより、単位インパルス応答以外のゲインおよび位相が０ではない波形を仮想入力波形として用いることが可能となる。
【００６７】
なお、本実施形態では、仮想入力波形としてステップ応答を用いた例を説明した。仮想入力波形には、例えば、低周波領域におけるデータの精度をさらに高めるべく、ランプ応答を用いてもよい。ランプ応答は、伝達関数がＨ（ｓ）＝１／ｓ2であり、図２８に実線Ｌ３で示すように、ゲインの傾きが−４０ｄＢ／ｄｅc、位相が全周波数帯で−１８０ｄｅｇであって、ステップ応答よりもさらに高周波振動成分を励起させない。このように、シミュレーションしたい周波数帯等に適した仮想入力波形を適宜選択することができる。
また、実験データにより即したシミュレーションを行いたい場合には、インパルスハンマ２１により入力された実際の入力波形（図２参照）をＦＦＴ変換し、仮想入力波形として用いてもよい。
【００６８】
（他の実施形態）
（ア）上記実施形態では、被測定物の分系の数は２つであったが、他の実施形態では分系の数は３以上であってもよい。また、分系同士の接続形態は、直線状の接続に限らず、例えば分岐している等、どのような形態であってもよい。また、上記実施形態では、１つの分系における入力点および出力点は各１つずつであった。他の実施形態では、１つの分系に複数の入力点および出力点を有するように構成してもよい。
また、上記実施形態では、被測定物の分系同士は剛結合にて接続されていたが、他の実施形態では、剛結合に限らない。この場合、例えば上記式（５）において、所定の定数や伝達関数が係数として入ることになる。
【００６９】
（イ）上記実施形態では、周波数応答特性から伝達関数に変換することなく四則演算を行う方法について説明した。他の実施形態では、予め伝達関数がわかっている場合や、周波数応答特性から比較的容易に伝達関数を導出可能な場合には、伝達関数を用いて演算してもよい。例えば、伝達関数と周波数応答特性との加減算を行う場合は、伝達関数から周波数応答特性を算出し、上記実施形態の如く演算することができる。また、伝達関数のみからなる演算が含まれる場合、当該演算においては伝達関数にて演算を行ってもよい。これにより、演算負荷をさらに低減することができる。
また、各種演算に用いられる周波数応答特性は、分系の周波数応答特性（例えば、Ｇ21、Ｇ22、Ｇ33、Ｇ43等）に限らず、周波数領域にて乗除算がなされた演算結果、および、時間領域にて加減算された後に周波数領域に変換された演算結果であってもよい。
【００７０】
（ウ）上記実施形態では、被測定物に打撃力を加えたときの力の伝達に係る伝達関数を推定していたが、他の実施形態では、例えば熱等、他の伝達に係る伝達関数の推定に用いてもよい。
（エ）上記実施形態にて説明した伝達関数を推定する推定方法は、伝達関数推定装置１０にて実行されるものである。伝達関数推定装置１０においては、これらの処理をソフトウェアで処理する伝達関数推定プログラムを備えていてもよいし、例えば一部の処理をハードウェアにて実行するように構成されていてもよい。
以上、本発明は、上記実施形態になんら限定されるものではなく、発明の趣旨を逸脱しない範囲において種々の形態で実施可能である。
【符号の説明】
【００７１】
１０・・・伝達関数推定装置
２０・・・測定装置
３０・・・被測定物
３１・・・第１の分系（分系）
３２・・・第２の分系（分系）
Ｐ１・・・入力点（入力点）
Ｐ２・・・出力点（出力点、接続点）
Ｐ３・・・入力点（入力点、接続点）
Ｐ４・・・出力点（出力点）
【図１】
2014021603.tif 000006
【図２】
2014021603.tif 000007
【図３】
2014021603.tif 000008
【図４】
2014021603.tif 000009
【図５】
2014021603.tif 000010
【図６】
2014021603.tif 000011
【図７】
2014021603.tif 000012
【図８】
2014021603.tif 000013
【図９】
2014021603.tif 000014
【図１０】
2014021603.tif 000015
【図１１】
2014021603.tif 000016
【図１２】
2014021603.tif 000017
【図１３】
2014021603.tif 000018
【図１４】
2014021603.tif 000019
【図１５】
2014021603.tif 000020
【図１６】
2014021603.tif 000021
【図１７】
2014021603.tif 000022
【図１８】
2014021603.tif 000023
【図１９】
2014021603.tif 000024
【図２０】
2014021603.tif 000025
【図２１】
2014021603.tif 000026
【図２２】
2014021603.tif 000027
【図２３】
2014021603.tif 000028
【図２４】
2014021603.tif 000029
【図２５】
2014021603.tif 000030
【図２６】
2014021603.tif 000031
【図２７】
2014021603.tif 000032
【図２８】
2014021603.tif 000033