特許 5854551 リアルタイム周波数解析方法

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】5854551

(24)【登録日】2015年12月18日

(45)【発行日】2016年2月9日

(54)【発明の名称】リアルタイム周波数解析方法

(51)【国際特許分類】

   G01H 13/00 20060101AFI20160120BHJP

   G01R 23/15 20060101ALI20160120BHJP

   G01H 17/00 20060101ALI20160120BHJP

【ＦＩ】

   G01H13/00

   G01R23/15 T

   G01H17/00 D

【請求項の数】6

【全頁数】18

(21)【出願番号】特願2011-191146(P2011-191146)

(22)【出願日】2011年9月1日

(65)【公開番号】特開2013-53883(P2013-53883A)

(43)【公開日】2013年3月21日

【審査請求日】2014年8月18日

(73)【特許権者】

【識別番号】000142067

【氏名又は名称】株式会社共和電業

(74)【代理人】

【識別番号】100082636

【弁理士】

【氏名又は名称】真田 修治

(72)【発明者】

【氏名】坂田 光児

【審査官】 田中 秀直

(56)【参考文献】

【文献】 特開平０８−１３６４０９（ＪＰ，Ａ）

【文献】 特開２００９−１７５１３０（ＪＰ，Ａ）

【文献】 特開２００７−２２５５６６（ＪＰ，Ａ）

【文献】 特表２００９−５２４８１２（ＪＰ，Ａ）

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ｇ０１Ｈ １／００−１７／００

Ｇ０１Ｒ ２３／１５

(57)【特許請求の範囲】

【請求項1】

加速度、荷重等の物理量を測定するセンサの出力を一定の時間間隔でサンプリングして取り出すＡ／Ｄ変換器の出力時系列データ、または該出力時系列データを一旦記憶装置に集録したデータに対し、任意の周波数成分の振動特性を抽出するリアルタイム周波数解析法であって、
任意の１つの中心周波数近傍の狭帯域の周波数成分の波形を抽出する狭帯域周波数成分波形抽出ステップと、
前記狭帯域周波数成分波形抽出ステップで抽出された波形の絶対値を求める絶対値回路の機能を有する絶対値算出ステップと、
前記絶対値算出ステップの出力絶対値波形の、任意時間毎の移動平均を求める移動平均算出ステップと、
を備え、
前記狭帯域周波数成分波形抽出ステップは、
用いる時系列データの現時点のデータおよび２サンプル前のデータと、現時点より１サンプル前および２サンプル前の当該ステップの出力から、現時点の前記当該ステップの出力を算出する機能を有し、
用いる時系列データのサンプリング時間間隔と、任意の中心周波数と、中心周波数に対する１より小さい任意の通過帯域比率とを用いて該中心周波数に対する通過帯域比率分の周波数だけ中心周波数から低い周波数側と高い周波数側にずれたカットオフ周波数と、前記任意の通過帯域比率に応じた遮断特性を有し、
前記絶対値算出ステップは、
入力された値が正値の場合はそのままの値を出力し、また、入力された値が負値の場合は負の符号を外した値を出力する機能を有し、
前記移動平均算出ステップは、
前記狭帯域周波数成分波形抽出ステップにおいて設定した任意の中心周波数に相当する周期の整数倍の時間を、前記狭帯域周波数成分振幅抽出ステップに適用するデータのサンプリング時間間隔で除したデータ数分の逆数を当該ステップに用いる時系列データの現時点のデータに乗じ、この値と、１から移動平均を取るデータ数の逆数を差し引いた値を現時点の１サンプル前の当該移動平均ステップの算出結果に乗じた値と、を足し合わせた値を当該移動平均ステップの現時点の算出結果とすることで、順次１サンプリング間隔毎の平均値を算出する機能を有し、
前記移動平均算出ステップの出力結果が前記狭帯域周波数成分波形抽出ステップにおいて設定した任意の１つの中心周波数近傍の周波数成分の時々刻々の振幅の実行値に相当する値を任意の１つの周波数成分の振動特性として抽出するに際し、掛け算の回数および一連のステップで記憶を必要とするデータの数を少なくして、演算負荷と必要記憶容量を小さく抑えるように構成したことを特徴とするリアルタイム周波数解析方法。

【請求項2】

請求項１に記載のリアルタイム周波数解析方法の全ステップを、所定の複数個の中間周波数毎に並列に処理できるステップを有し、前記複数並列の処理方法の各々の処理方法で用いる狭帯域周波数成分波形抽出ステップの中心周波数と、中心周波数に対する１より小さい任意の通過帯域比率とは、各々異なる値を設定することができ、設定した任意の複数の中心周波数近傍の周波数成分の時々刻々の振幅の実行値に相当する値を任意の複数の周波数成分の振動特性として抽出することを特徴とするリアルタイム周波数解析方法。

【請求項3】

前記解析対象への入力信号を検出する第１のセンサの出力データと前記解析対象からの出力信号を検出する前記第２のセンサの出力データに対して、請求項１に記載の前記狭帯域周波数成分波形抽出ステップと、前記絶対値算出ステップと、前記移動平均算出ステップとを適用する入出力狭帯域周波数成分波形抽出ステップと、
前記入出力狭帯域周波数成分波形抽出ステップの、前記第１のセンサと前記第２のセンサとに対応した出力の出力比を求める除算ステップと、
を備え、
前記第１のセンサと前記第２のセンサとの出力に適用する、請求項１記載の前記狭帯域周波数成分波形抽出ステップと、前記絶対値算出ステップと、前記移動平均算出ステップにおける各処理方法では、同じ中心周波数と、同じ通過帯域比率と、同じ整数周期分の移動平均データ数とを参照することで、前記解析対象に設定した任意の中心周波数に対する伝達特性応答倍率を、任意の周波数成分の振動特性として抽出することを特徴とする請求項１に記載のリアルタイム周波数解析方法。

【請求項4】

前記解析対象への入力信号を検出する前記第１のセンサの出力データと前記第２のセンサの出力データに対して、請求項１に記載の前記狭帯域周波数成分波形抽出ステップと、前記絶対値算出ステップと、前記移動平均算出ステップとを適用する入出力狭帯域周波数成分波形抽出ステップと、
前記入出力狭帯域周波数成分波形抽出ステップの、前記第１のセンサと前記第２のセンサとに対応した出力の出力比を求める除算ステップと、
前記除算ステップを前記第１のセンサに対応するデータに適用した結果の出力と、前記除算ステップを前記第２のセンサに対応するデータに適用した結果の出力とを乗算する伝達特性応答倍率算出ステップと、
請求項１に記載の前記狭帯域周波数成分波形抽出ステップを前記第１のセンサの出力データに適用した結果の出力と、前記狭帯域周波数成分波形抽出ステップを前記第２のセンサの出力データに適用した結果の出力とを乗算する波形乗算ステップと、
前記波形乗算ステップの出力波形に対し任意時間毎の移動平均を求める移動平均算出ステップと、
前記移動平均算出ステップの算出結果の２倍を、前記伝達特性応答倍率算出ステップの算出結果で除した結果の逆余弦を求める位相算出ステップと、
を備えたことを特徴とする請求項１に記載のリアルタイム周波数解析方法。

【請求項5】

請求項３に記載のリアルタイム周波数解析方法の全ステップを、所定の複数個の中間周波数毎に並列に処理するステップを有することを特徴とする請求項３に記載のリアルタイム周波数解析方法。

【請求項6】

請求項４に記載のリアルタイム周波数解析方法の全ステップを、所定の複数個の中間周波数毎に並列に処理するステップを有することを特徴とする請求項４に記載のリアルタイム周波数解析方法。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、リアルタイム周波数解析方法に関し、特に、計算のアルゴリズムを簡素化して処理装置や記憶装置の負担を少なくしたリアルタイム周波数解析方法に関する。

【背景技術】

【0002】

従来、時々刻々に周波数成分が変化する信号に対する周波数解析として使われてきた技術には、主に２つの技術が有ることが電気情報通信学会誌８／’９６（非特許文献１）に記載されている。以下では、この２つの技術について説明する。
（１） 短時間フーリエ変換を用いた解析方法
この方法は、基本的には、定常信号を解析するためのフーリエ変換［下記（１６）式参照］の処理を、図１０に示すように、０〜T秒間の信号をＮ分割したＴ／Ｎの長さ毎（ｉ＝１〜Ｎ）の信号に適用し、その結果を、該時間区間の周波数成分とする解析法である。

【0003】

【数1】

【0004】

この方法では、フーリエ変換が積分変換となるため、演算に多大の時間を要する。また、演算を早くする高速フーリエ変換（FFT）のアルゴリズムを使ったとしても、時々刻々の演算は複雑となり、また、ＦＦＴの制約により、周波数分解能ΔＦは、計算に用いる信号の長さの逆数、即ちΔｆ＝Ｎ／Ｔにならざるを得ないため、図１０において、周波数の時間変化を細かく見るために分割数Ｎを大きくする（個々の計算する信号の時間長さを短くする）と、周波数分解能は荒くならざるを得ず、これとは逆に、周波数分解能を細かくしようとする（個々の計算に用いる信号の長さを長くする）と、周波数成分の時間変化の分解能を荒くせざるを得ないといった欠点を有している。
（２） ウェーブレット変換［下記（１７）式参照］を用いた解析方法

【0005】

【数2】

【0006】

この方法も、前述のフーリエ変換と同様で、（１６）式でのｅの乗数である−ｊωｔの代わりに、マザーウェーブレットと呼ばれる関数を、元の信号に掛けた積分変換であるため、計算のアルゴリズムが複雑となり、計算にも多大の時間を要することになる。
このマザーウェーブレットと呼ばれる関数には、幾つかの種類の関数が提案されているが、例えば、スケールと呼ばれている数ａの逆数である１／ａが周波数に対応し、数ｂが時間のシフトに相当するような、ａ，ｂ対応の関数も提案されている。
しかしながら、この解析方法も、前述の計算時間の問題だけでは済まず、周波数分解能が周波数によって倍々に変化することになるため、周波数が高くなるに従って周波数分解能が荒くならざるを得ないことになり、また時間分解能は、低周波になるに従って荒くならざるを得ないといった欠点を有している。
この分野の特許文献としては、例えば特許文献１（特開２０００−１８６９８４号公報）に、簡単な構成および処理によって、洗掘等に起因する橋脚の支持力の変化を適切に且つ定常的に監視し、運行規制の的確な発令／解除を可能とする技術が開示されている。具体的には、ＦＦＴ処理部が、橋脚の天端に取着された加速度計による検出データに高速フーリエ変換を施してフーリエスペクトルを得る。平均化処理部は、フーリエスペクトルを所定の時間についての加算平均により時間的な平均を求め、所定の周波数幅について周波数軸上で移動平均を求めて、フーリエスペクトル波形を平滑化する。重み付け処理部は、健全時の橋脚の伝達関数に基づく１自由度系の応答倍率曲線に基づく重み付け関数で平滑化フーリエスペクトルを重み付けする。面積算定正規化処理部は、重み付けフーリエスペクトルの曲線とベースラインとで囲まれる領域の面積を求めて正規化し、正規化面積減小評価部による閾値との比較評価に供するものとしている。

【先行技術文献】

【特許文献】

【0007】

【特許文献1】「特開２０００−１８６９８４号公報」

【非特許文献】

【0008】

【非特許文献1】電気情報通信学会誌、Vol.79,No.8,pp.820-830,1996年8月「時間−周波数解析の展望（IV）−ウェーブレットとその分類−」

【発明の概要】

【発明が解決しようとする課題】

【0009】

ところで、リアルタイム周波数解析方法は、構造物の振動を利用したヘルスモニタリングに使用されることが多いが、その普及には構造物に取り付ける振動センサ内にてある程度振動特性を評価できるところまで演算処理した結果が出力できるような安価なスマートセンサの実現が必要であると考えられる。
このためのスマートセンサ内での演算処理は、高度な処理能力や大きな記憶容量は望めないであろうことから、負荷が少なく、また、記憶容量を多く必要としない処理法が必要となる。
但し、構造物の振動データから、この振動特性の変化を見ようとした場合、少なくとも振動データの周波数成分の経時変化を求めることができる時間周波数解析的な手法が必要になる。
しかしながら、上記背景技術で述べた従来のリアルタイム周波数解析方法にあっては、前述のとおり、時間周波数解析としては短時間フーリエ変換によるランニングスペクトル、ウェーブレット解析等が通常使われているが、これらの解析手法は、演算ステップ数がリアルタイム処理を行うには多めであり、かつ、ある程度の時間データが必要なため、処理装置の処理能力や記憶容量の点で、安価なスマートセンサに組込むのは難しいと考えられる。

【0010】

また、周波数分解能に関して、必要な周波数帯域内で必要な周波数分解能を確保することが難しい。
これらの点から、従来のリアルタイム周波数解析方法は、構造物のリアルタイムヘルスモニタリングに用いるには不向きであるという問題点があった。
本発明は、上記従来の問題点に鑑みてなされたものであって、振動の物理現象と直感的に結びつく１自由度振動系を基にした複数の狭帯域バンドパスフィルタを用いる解析法を使用することを可能にして、必要な周波数分解能を確保しつつ、処理装置に対して負荷が少なく、また、記憶容量を多く必要としない、処理能力の高くない処理装置に適用可能で、特に構造物の振動特性評価のために好適なリアルタイム周波数解析方法を提供することを目的としている。

【課題を解決するための手段】

【0011】

請求項１に記載した本発明に係るリアルタイム周波数解析方法は、上述した目的を達成するために、加速度、荷重等の物理量を測定するセンサの出力を一定の時間間隔でサンプリングして取り出すＡ／Ｄ変換器の出力時系列データ、または該出力時系列データを一旦記憶装置に集録したデータに対し、任意の周波数成分の振動特性を抽出するリアルタイム周波数解析法であって、
任意の１つの中心周波数近傍の狭帯域の周波数成分の波形を抽出する狭帯域周波数成分波形抽出ステップと、
前記狭帯域周波数成分波形抽出ステップで抽出された波形の絶対値を求める絶対値回路の機能を有する絶対値算出ステップと、
前記絶対値算出ステップの出力絶対値波形の、任意時間毎の移動平均を求める移動平均算出ステップと、
を備え、
前記狭帯域周波数成分波形抽出ステップは、
用いる時系列データの現時点のデータおよび２サンプル前のデータと、現時点より１サンプル前および２サンプル前の当該ステップの出力から、現時点の前記当該ステップの出力を算出する機能を有し、
用いる時系列データのサンプリング時間間隔と、任意の中心周波数と、中心周波数に対する１より小さい任意の通過帯域比率とを用いて該中心周波数に対する通過帯域比率分の周波数だけ中心周波数から低い周波数側と高い周波数側にずれたカットオフ周波数と、前記任意の通過帯域比率に応じた遮断特性を有し、
前記絶対値算出ステップは、
入力された値が正値の場合はそのままの値を出力し、また、入力された値が負値の場合は負の符号を外した値を出力する機能を有し、
前記移動平均算出ステップは、
前記狭帯域周波数成分波形抽出ステップにおいて設定した任意の中心周波数に相当する周期の整数倍の時間を、前記狭帯域周波数成分振幅抽出ステップに適用するデータのサンプリング時間間隔で除したデータ数分の逆数を当該ステップに用いる時系列データの現時点のデータに乗じ、この値と、１から移動平均を取るデータ数の逆数を差し引いた値を現時点の１サンプル前の当該移動平均ステップの算出結果に乗じた値と、を足し合わせた値を当該移動平均ステップの現時点の算出結果とすることで、順次１サンプリング間隔毎の平均値を算出する機能を有し、
前記移動平均算出ステップの出力結果が前記狭帯域周波数成分波形抽出ステップにおいて設定した任意の１つの中心周波数近傍の周波数成分の時々刻々の振幅の実行値に相当する値を任意の１つの周波数成分の振動特性として抽出するに際し、掛け算の回数および一連のステップで記憶を必要とするデータの数を少なくして、演算負荷と必要記憶容量を小さく抑えるように構成したことを特徴としている。

【0013】

また、請求項２に記載した本発明に係るリアルタイム周波数解析方法は、請求項１のリアルタイム周波数解析方法であって、請求項１のリアルタイム周波数解析方法の全ステップを、所定の複数個の中間周波数毎に並列に処理できるステップを有し、前記複数並列の処理方法の各々の処理方法で用いる狭帯域周波数成分波形抽出ステップの中心周波数と、中心周波数に対する１より小さい任意の通過帯域比率とは、各々異なる値を設定することができ、設定した任意の複数の中心周波数近傍の周波数成分の時々刻々の振幅の実行値に相当する値を任意の複数の周波数成分の振動特性として抽出することを特徴としている。

【0014】

また、請求項３に記載した本発明に係るリアルタイム周波数解析方法は、請求項１のリアルタイム周波数解析方法であって、前記解析対象への入力信号を検出する第１のセンサの出力データと前記解析対象からの出力信号を検出する前記第２のセンサの出力データに対して、請求項１に記載の前記狭帯域周波数成分波形抽出ステップと、前記絶対値算出ステップと、前記移動平均算出ステップとを適用する入出力狭帯域周波数成分波形抽出ステップと、
前記入出力狭帯域周波数成分波形抽出ステップの、前記第１のセンサと前記第２のセンサとに対応した出力の出力比を求める除算ステップと、
を備え、
前記第１のセンサと前記第２のセンサとの出力に適用する、請求項１記載の前記狭帯域周波数成分波形抽出ステップと、前記絶対値算出ステップと、前記移動平均算出ステップにおける各処理方法では、同じ中心周波数と、同じ通過帯域比率と、同じ整数周期分の移動平均データ数とを参照することで、前記解析対象に設定した任意の中心周波数に対する伝達特性応答倍率を、任意の周波数成分の振動特性として抽出することを特徴としている。

【0015】

また、請求項４に記載した本発明に係るリアルタイム周波数解析方法は、請求項１のリアルタイム周波数解析方法であって、前記解析対象への入力信号を検出する前記第１のセンサの出力データと前記第２のセンサの出力データに対して、請求項１に記載の前記狭帯域周波数成分波形抽出ステップと、前記絶対値算出ステップと、前記移動平均算出ステップとを適用する入出力狭帯域周波数成分波形抽出ステップと、
前記入出力狭帯域周波数成分波形抽出ステップの、前記第１のセンサと前記第２のセンサとに対応した出力の出力比を求める除算ステップと、
前記除算ステップを前記第１のセンサに対応するデータに適用した結果の出力と、前記除算ステップを前記第２のセンサに対応するデータに適用した結果の出力とを乗算する伝達特性応答倍率算出ステップと、
請求項１に記載の前記狭帯域周波数成分波形抽出ステップを前記第１のセンサの出力データに適用した結果の出力と、前記狭帯域周波数成分波形抽出ステップを前記第２のセンサの出力データに適用した結果の出力とを乗算する波形乗算ステップと、
前記波形乗算ステップの出力波形に対し任意時間毎の移動平均を求める移動平均算出ステップと、
前記移動平均算出ステップの算出結果の２倍を、前記伝達特性応答倍率算出ステップの算出結果で除した結果の逆余弦を求める位相算出ステップと、
を備えたことを特徴としている。

【0016】

また、請求項５に記載した本発明に係るリアルタイム周波数解析方法は、請求項３に記載のリアルタイム周波数解析方法の全ステップを、所定の複数個の中間周波数毎に並列に処理するステップを有することを特徴としている。
また、請求項６に記載した本発明に係るリアルタイム周波数解析方法は、請求項４に記載のリアルタイム周波数解析方法の全ステップを、所定の複数個の中間周波数毎に並列に処理するステップを有することを特徴としている。

【発明の効果】

【0017】

以上述べたように、本発明によれば、周波数分解能と時間分解能とを、時々刻々に変化する測定結果の現在値を示す関数に含まれる係数内の中心周波数と、通過帯域比率を設定することにより、任意に定めることが可能となり、計算アルゴリズムについても、１つの周波数に対して、前記関数を基本とした数式を時々刻々計算するだけであり、
特に、前記移動平均算出ステップでは、
前記狭帯域周波数成分波形抽出ステップにおいて設定した任意の中心周波数に相当する周期の整数倍の時間を、前記狭帯域周波数成分振幅抽出ステップに適用するデータのサンプリング時間間隔で除したデータ数分の逆数を当該ステップに用いる時系列データの現時点のデータに乗じ、この値と、１から移動平均を取るデータ数の逆数を差し引いた値を現時点の１サンプル前の当該移動平均ステップの算出結果に乗じた値と、を足し合わせた値を当該移動平均ステップの現時点の算出結果とすることで、順次１サンプリング間隔毎の平均値を算出する機能を有し、
前記移動平均算出ステップの出力結果が前記狭帯域周波数成分波形抽出ステップにおいて設定した任意の１つの中心周波数近傍の周波数成分の時々刻々の振幅の実行値に相当する値を任意の１つの周波数成分の振動特性として抽出するに際し、掛け算の回数および一連のステップで記憶を必要とするデータの数を少なくして、演算負荷と必要記憶容量を小さく抑えるように構成したことにより、従来よりも単純な計算処理となるので、演算負荷を少なくすることができ、処理能力の高くない処理装置に適用可能で、負荷が少なく、記憶容量を多く必要としない構造物の振動特性評価のための時間周波数解析法、伝達関数算出法を構築することができる。また、本発明によれば、構造物に取り付けるヘルスモニタリング用スマートセンサ内に構造物の振動データを用いた健全性評価のための時々刻々の処理アルゴリズムを組み込むことが容易になる。
また、請求項２〜請求項６の記載によれば、いずれも請求項１に規定するように、掛け算の回数および一連のステップで記憶を必要とするデータの数を少なくして、演算負荷と必要記憶量を小さく抑えるように構成することにより、処理能力の高くない処理装置に適用可能で、負荷が少なく、記憶容量を多く必要としない構造物の振動特性評価のために好適なリアルタイム周波数解析方法を提供することができ、延いては、構造物に取り付けるヘルスモニタリング用スマートセンサ内に構造物の振動データを用いた健全性評価のための時々刻々の処理アルゴリズムを組み込むことが容易になる。

【図面の簡単な説明】

【0018】

【図1】本発明の実施の形態１に係るリアルタイム周波数解析方法の処理手順を示すフローチャート図である。

【図2】本発明の実施の形態２に係るリアルタイム周波数解析方法の処理手順を示すフローチャート図である。

【図3】本発明の実施の形態３に係るリアルタイム周波数解析方法の処理手順を示すフローチャート図である。

【図4】本発明の実施の形態４に係るリアルタイム周波数解析方法の処理手順を示すフローチャート図である。

【図5】本発明の実施の形態５に係るリアルタイム周波数解析方法の処理手順を示すフローチャート図である。

【図6】本発明の実施の形態６に係るリアルタイム周波数解析方法の処理手順を示すフローチャート図である。

【図7】本発明の実施の形態２に係るリアルタイム周波数解析方法での処理方法を示す説明図である。

【図8】本発明の実施の形態５に係るリアルタイム周波数解析方法での処理方法を示す説明図である。

【図9】本発明の実施の形態６に係るリアルタイム周波数解析方法での処理方法を示す説明図である。

【図10】従来のリアルタイム周波数解析方法での短時間フーリエ変換を用いた解析方法を示す説明図である。

【発明を実施するための形態】

【0019】

本発明のリアルタイム周波数解析方法は、加速度、荷重等の物理量を測定するセンサの出力を一定の時間間隔でサンプリングして取り出すＡ／Ｄ変換器（図示は省略）を使用し、該Ａ／Ｄ変換器の出力時系列データ、または該出力時系列データを、一旦記憶装置に集録したデータに対し、任意の周波数成分の振動特性を抽出する方法である。
よって、本発明のリアルタイム周波数解析方法は、構造物のヘルスモニタリングのためのリアルタイム周波数解析に適した方法である。
以下、本発明のリアルタイム周波数解析方法の実施の形態について、〔実施の形態１〕〜〔実施の形態７〕の順に、図面を参照して詳細に説明する。

【0020】

〔実施の形態１〕
本実施の形態１に係るリアルタイム周波数解析方法は、設定した任意の１つの中心周波数近傍の周波数成分の時々刻々の振幅の実行値に相当する値を任意の１つの周波数成分の振動特性として抽出する方法である。
即ち、本方法は、特定の周波数成分の時間変化を計算する方法であり、例えば、回転機器周辺の振動において、回転機器の回転数と同じ周波数成分の変化を得ようとする場合、中心周波数を、時々刻々の回転機器の回転数に相当する値に変えることで、上記の目的を達成することができる。
具体的な処理ステップとしては、任意の１つの中心周波数近傍の狭帯域の周波数成分の波形を抽出する狭帯域周波数成分波形抽出ステップと、前記狭帯域周波数成分波形抽出ステップで抽出された波形の絶対値を求める絶対値回路の機能を有する絶対値算出のステップと、前記絶対値算出ステップの出力絶対値波形に対し任意時間毎の移動平均を求める移動平均算出ステップと、を有する。

【0021】

前記狭帯域周波数成分波形抽出ステップは、時系列データの現時点のデータおよび２サンプル前のデータと、現時点より１サンプル前および２サンプル前の当該ステップの出力とから、現時点の出力を算出し、その狭帯域周波数成分波形抽出の周波数特性は、時系列データのサンプリング時間間隔と、任意の中心周波数と、中心周波数に対する１より小さい任意の通過帯域比率と、を用いることにより、該中心周波数に対する通過帯域比率分の周波数だけ中心周波数から低い周波数側と高い周波数側にずれたカットオフ周波数と、前記任意の通過帯域比率に応じた遮断特性とを有する。
また、絶対値算出ステップは、入力された値が正の値の場合はそのままの値を出力し、入力された値が負の値の場合は該負を示す符号を外した値を出力する。
さらに、移動平均算出ステップは、前記狭帯域周波数成分波形抽出ステップにおいて設定した任意の中心周波数に相当する周期の整数倍の時間を、前記狭帯域周波数成分振幅抽出ステップに適用するデータのサンプリング時間間隔で除したデータ数分の逆数を当該ステップに用いる時系列データの現時点のデータに乗じ、この値と、１から移動平均を取るデータ数の逆数を差し引いた値を現時点の１サンプル前の当該移動平均ステップの算出結果に乗じた値と、を足し合わせた値を当該移動平均ステップの現時点の算出結果とすることで、順次１サンプリング間隔毎の平均値を算出する機能を有し、
前記移動平均算出ステップの出力結果が前記狭帯域周波数成分波形抽出ステップにおいて設定した任意の１つの中心周波数近傍の周波数成分の時々刻々の振幅の実行値に相当する値を任意の１つの周波数成分の振動特性として抽出するに際し、掛け算の回数および一連のステップで記憶を必要とするデータの数を少なくして、演算負荷と必要記憶容量を小さく抑えるように構成したものである。
この移動平均ステップは、信号の移動平均により、振幅実行値や平均値を求めるための移動平均計算を簡略化する方法である。この方法によれば、計算アルゴリズムの簡略化と計算速度の向上を達成することができる。
この（３）式で示す指数移動平均ステップは、移動平均（Ｚｎ）を取るデータ数の逆数を当該方法に用いる時系列データの現時点のデータに乗じ、この値と、１から移動平均を取るデータ数の逆数を差し引いた値を現時点の１サンプル前の前記指数移動平均ステップの算出結果に乗じた値とを足し合わせた値を、現時点の算出結果とするものである。
但し、（３）式の符号の意味を、下記のとおりとする。
ｚ：ｖの移動平均値；本ステップの出力
ｋ：移動平均を取るデータ数の逆数

【0022】

図１は、本発明の実施の形態１に係るリアルタイム周波数解析方法の処理手順を示すフローチャート図である。
但し、ここで使用される符号の意味は、下記のとおりとする。
ｘ：使用する時系列データ
ｙ：周波数成分の波形出力
入出力の添え字は、ｎ：現在の値、ｎ−ｊ：現在よりｊ：サンプル前の値
ｖ：ｙの絶対値
ｚ：ｖの移動平均値：本ステップの出力
ｋ：移動平均を取るデータ数の逆数
ａ、ｂ、ｃ：係数
（処理手順）
まず、ステップＳ１の数式〔（１）式〕を計算し、次にステップＳ２の数式〔（２）式〕を計算し、最後にステップＳ３の数式〔（３）式〕を計算する。

【0023】

【数3】

【0024】

【数4】

【0025】

【数5】

【0026】

〔実施の形態２〕
本発明の実施の形態２に係るリアルタイム周波数解析方法は、設定した任意の複数の中心周波数近傍の周波数成分の時々刻々の振幅の実行値に相当する値を任意の複数の周波数成分の振動特性として抽出する方法、即ち、前記中心周波数についての周波数成分の時間変化を得るための方法であるが、従来法に比べ、出力関数の係数ａ、ｂ、ｃ内の中心周波数と通過帯域比率の設定を任意に行うことが可能なため、周波数毎に所望の周波数分解能と時間分解能を設定できる利点を有する。このことは、例えば構造物の損傷により構造物の固有振動数が変化するような現象においては、比較的長い時間間隔における固有振動数の微小な変化を検知する必要があるため、時間分解能は荒くし、周波数分解能を細かくするといった操作を可能にするものである。
具体的な処理ステップとしては、前述の実施の形態１に係るリアルタイム周波数解析方法の全ステップを複数並列に処理するように配置した処理方法を有し、前記複数並列に処理するように配置された前記の処理方法で用いる狭帯域周波数成分波形抽出ステップの中心周波数と、該中心周波数に対する１より小さい任意の通過帯域比率とは、各々、異なる値を設定することができるようにし、これにより、設定した任意の複数の中心周波数近傍の周波数成分の時々刻々の振幅の実行値に相当する値を任意の複数の周波数成分の振動特性として抽出できるようにしている。

【0027】

図２は、本発明の実施の形態２に係るリアルタイム周波数解析方法処理手順を示すフローチャート図である。
但し、ここでは、測定環境を示す一般的な符号、及び図２に記載の符号の意味を下記のとおりとする。
ｘ：使用する時系列データ
ｙ：周波数成分の波形出力
入出力の添え字は、ｎ：現在の値、ｎ−ｊ：現在よりｊ：サンプル前の値
ｖ：ｙの絶対値
ｚ：ｖの移動平均値（ｐ個を出力）
ｋ：移動平均を取るデータ数の逆数
（処理手順）
図２におけるステップＳ１〜Ｓ３の処理を、ｉ＝１〜ｐ個並べて並列に処理する。

【0028】

〔実施の形態３〕
本発明の実施の形態３に係るリアルタイム周波数解析方法は、任意の中心周波数に対する伝達特性応答倍率を任意の周波数成分の振動特性として抽出できるようにしたものであり、特定の周波数の入力信号に対する出力信号の応答倍率を計算できることから、例えば、回転機器取付け部と周辺の振動において、回転数の変化に対し、それに相当する周波数成分の応答倍率の変化を求めることなどの処理が可能となり、振動対策が必要な場合などの解析に有用となる。
具体的な処理ステップとしては、前記実施の形態１に係るリアルタイム周波数解析方法の全処理ステップと同じ処理ステップを、所定の複数個の中間周波数毎に配して並列に処理できるステップを有し、かつ前記複数並列の処理方法の各々の処理方法で用いる狭帯域周波数成分波形抽出ステップの中心周波数と、中心周波数に対する１より小さい任意の通過帯域比率とは、各々異なる値を設定することを可能とする。これにより、設定した任意の複数の中心周波数近傍の周波数成分の時々刻々の振幅の実行値に相当する値を任意の複数の周波数成分の振動特性として抽出することができる。
図３は、本発明の実施の形態３に係るリアルタイム周波数解析方法の処理手順を示すフローチャート図である。
但し、ここでは、測定環境を示す一般的な符号、及び図３に記載の符号の意味を下記のとおりとする。

【0029】

ｘ：使用する時系列データ
ｙ：周波数成分の波形出力
入出力の添え字は、ｎ：現在の値、ｎ−ｊ：現在よりｊ：サンプル前の値、ｉｎ：構造物への入力、ｏｔ：構造物からの出力
ｖ：ｙの絶対値
ｚ：ｖの移動平均値
ｋ：移動平均を取るデータ数の逆数
Ａａｍｐ＿ｉ ：任意の中心周波数に対する伝達特性応答倍率とする。
（処理手順）
まず、ステップＳ４の数式〔（４）式および（５）式〕を計算し、次にステップＳ５の数式〔（６）式および（７）式〕を計算し、次に、ステップＳ６の数式〔（８）式および（９）式〕を計算し、最後にステップＳ７の数式〔（１０）式〕を計算する。

【0030】

【数6】

【0031】

【数7】

【0032】

【数8】

【0033】

【数9】

【0034】

【数10】

【0035】

【数11】

【0036】

【数12】

【0037】

〔実施の形態４〕
本発明の実施の形態４に係るリアルタイム周波数解析方法は、設定した任意の中心周波数に対する伝達特性位相差を任意の周波数成分の振動特性として抽出する方法であり、即ち、特定の周波数の入力信号に対する出力信号の位相差を計算するものであり、前述の実施の形態３と同様の用途への適用が可能であり、複数の出力点と入力点の位相関係から振動の伝達経路の推定等に役立たせることができる。
具体的な処理ステップとしては、解析対象への入力信号を検出する第１のセンサの出力データと前記第２のセンサの出力データに対して、前記実施の形態１と同じ、狭帯域周波数成分波形抽出ステップと、絶対値算出ステップと、移動平均算出ステップとを適用する入出力狭帯域周波数成分波形抽出ステップと、前記入出力狭帯域周波数成分波形抽出ステップの、前記第１のセンサと前記第２のセンサとに対応した出力の出力比を求める除算ステップと、前記除算ステップを前記第１のセンサに対応するデータに適用した結果の出力と、前記除算ステップを前記第２のセンサに対応するデータに適用した結果の出力とを乗算する伝達特性応答倍率算出ステップと、前記実施の形態１と同じ狭帯域周波数成分波形抽出ステップを前記第１のセンサの出力データに適用した結果の出力と、前記狭帯域周波数成分波形抽出ステップを前記第２のセンサの出力データに適用した結果の出力とを乗算する波形乗算ステップと、前記波形乗算ステップの出力波形に対し任意時間毎の移動平均を求める移動平均算出ステップと、前記移動平均算出ステップの算出結果の２倍を、前記伝達特性応答倍率算出ステップの算出結果で除した結果の逆余弦を求める位相算出ステップと、を有する。

【0038】

このように構成して、前記各ステップには、前記実施の形態１と同じ中心周波数と、同じ通過帯域比率と、を用い、前記移動平均算出ステップには前記実施の形態１と同じ整数周期分の移動平均データ数を用いることで、前記構造物の設定した任意の中心周波数に対する伝達特性位相差を任意の周波数成分の振動特性として抽出することができるようにしている。
図４は、本発明の実施の形態４に係るリアルタイム周波数解析方法の処理手順を示すフローチャート図である。
但し、ここでは、測定環境を示す一般的な符号、及び図４に記載の符号の意味を下記のとおりとする。

【0039】

ｘ：使用する時系列データ
ｙ：周波数成分の波形出力
入出力の添え字は、ｎ：現在の値、ｎ−ｊ：現在よりｊ：サンプル前の値、ｉｎ：構造物への入力、ｏｔ：構造物からの出力
ｖ：ｙの絶対値
ｚ：ｖの移動平均値
ｋ：移動平均を取るデータ数の逆数
Aamp_inot：任意の中心周波数に対する構造物の入出力の振幅成分の積
Ainot ：任意の中心周波数に対する構造物の入出力の波形の積
cos φ：中心周波数に対する伝達特性位相差の余弦
φ：中心周波数に対する伝達特性位相差；
（処理手順）
まず、ステップＳ４の数式〔（４）式と（５）式〕を計算し、次に、ステップＳ５，６の数式〔（６）式、（７）式、（８）式、（９）式〕とステップＳ８〔（１１）式〕の数式を計算し、また、この計算と並列にステップＳ９〔（１２）式〕の数式とＳ１０〔（１３）式〕の数式を計算し、最後にステップＳ１１〔（１４）式〕の数式を計算する。

【0040】

【数13】

【0041】

【数14】

【0042】

【数15】

【0043】

【数16】

【0044】

〔実施の形態５〕
本発明の実施の形態５に係るリアルタイム周波数解析方法は、狭帯域伝達特性応答倍率スペクトル抽出方法であり、即ち、複数の中心周波数についての入力信号に対する出力信号の応答倍率を計算する方法である。
この方法により、解析対象の時々刻々に伝達関数の内の応答倍率の変化を見ることが可能となるので、この方法は、例えば構造物の損傷により構造物の固有振動数の変化を検知することに適用することができる。
具体的な処理ステップとしては、前記実施の形態３に係るリアルタイム周波数解析方法の全処理ステップと同じ処理ステップを複数配して並列処理できるステップを有し、前記複数並列に処理するように配置された処理方法の各々が用いる狭帯域周波数成分波形抽出ステップの中心周波数と、該中心周波数に対する１より小さい任意の通過帯域比率とは、各々異なる値を設定することができるようにしている。これにより、設定した任意の複数の任意の中心周波数に対する伝達特性応答倍率を任意の周波数成分の振動特性として抽出することができる。

【0045】

図５は、本発明の実施の形態５に係るリアルタイム周波数解析方法の処理手順を示すフローチャート図である。
但し、ここでは、測定環境を示す一般的な符号、及び図５に記載の符号の意味を下記のとおりとする。
ｘ：使用する時系列データ
ｙ：周波数成分の波形出力
入出力の添え字は、ｎ：現在の値、ｎ−ｊ：現在よりｊ：サンプル前の値、ｉｎ：構造物への入力、ｏｔ：構造物からの出力
ｖ：ｙの絶対値
ｚ：ｖの移動平均値
Aamp＿i ：任意の中心周波数に対する伝達特性応答倍率
（処理手順）
図５のステップＳ４〜Ｓ７の処理を、ｉ＝１〜ｐ個並べて並列に処理する。

【0046】

〔実施の形態６〕
本発明の実施の形態６に係るリアルタイム周波数解析方法は、狭帯域周波数成分位相差スペクトル抽出方法であり、即ち、複数の中心周波数についての入力信号に対する出力信号の位相特性を計算するものである。
この方法を使用すると、解析対象の時々刻々に伝達関数の内の応答倍率の変化を見ることが可能となり、よって、例えば構造物の損傷により構造物の固有振動数の変化を位相特性の変化により検知することに適用することができる。
具体的な処理ステップとしては、前記実施の形態４の全ステップと同じ処理ステッブを複数配して並列処理できるようにしたステップを有し、かつ前記複数並列の処理方法の各々の処理方法で用いる狭帯域周波数成分波形抽出ステップの中心周波数と、中心周波数に対する１より小さい任意の通過帯域比率とは、各々異なる値を設定することを可能とする。これにより、設定した任意の複数の任意の中心周波数に対する伝達特性位相差を任意の周波数成分の振動特性として抽出することができる。
図６は、本発明の実施の形態６に係るリアルタイム周波数解析方法の処理手順を示すフローチャート図である。
但し、ここでは、測定環境を示す一般的な符号および図６に記載の符号の意味を下記のとおりとする。

【0047】

ｘ：使用する時系列データ
ｙ：周波数成分の波形出力
入出力の添え字は、ｎ：現在の値、ｎ−ｊ：現在よりｊ：サンプル前の値、ｉｎ：構造物への入力、ｏｔ：構造物からの出力
ｖ：ｙの絶対値
ｚ：ｖの移動平均値
ｋ：移動平均を取るデータ数の逆数
Aamp＿inot：任意の中心周波数に対する構造物の入出力の振幅成分の積
Ainot ：任意の中心周波数に対する構造物の入出力の波形の積
cos φ：中心周波数に対する伝達特性位相差の余弦
φ：中心周波数に対する伝達特性位相差（ｐ個出力）
（処理手順）
図６におけるステップＳ４〜Ｓ６およびステップＳ８〜Ｓ１１の処理を、ｉ＝１〜ｐ個並べて並列に処理する。

【0050】

図７は、本発明の実施の形態２に係るリアルタイム周波数解析方法での処理を示す説明図である。
同図において、本発明の実施の形態１に係るリアルタイム周波数解析方法での処理は、一点鎖線内のブロックに纏めて示している。
図７に示す処理は、前述のとおり、中心周波数の周波数成分の時々刻々の振幅値を求める処理である。
図８は、本発明の実施の形態５に係るリアルタイム周波数解析方法での処理を示す説明図である。
同図において、本発明の実施の形態３に係るリアルタイム周波数解析方法での処理は、一点鎖線内のブロックに纏めて示している。
図８に示す処理は、前述のとおり、各々異なる中心周波数での本発明の実施の形態１に係るリアルタイム周波数解析方法の処理を纏めたものであり、各々の中心周波数の周波数成分の時々刻々の振幅値を求める処理である。

【0051】

図９は、本発明の実施の形態６に係るリアルタイム周波数解析方法での処理を示す説明図である。
図９において、本発明の実施の形態４に係るリアルタイム周波数解析方法での処理は一点鎖線（但し、画像変換時の画質の劣化により一点鎖線でなく、破線として見える場合もある）内のブロックに纏めて示している。
図９に示す処理は、前述のとおり、各々異なる中心周波数での本発明の実施の形態４に係るリアルタイム周波数解析方法の処理を纏めたものであり、各周波数毎の入出力の位相差を求める処理である。
なお、本発明に係るリアルタイム周波数解析方法を実現する処理の少なくとも一部をコンピュータ制御により実行するものとし、かつ、上記処理を、図１〜６のフローチャートで示した手順および（１５）式によりコンピュータに実行せしめるプログラムは、半導体メモリを始め、ＣＤ−ＲＯＭや磁気テープなどのコンピュータ読み取り可能な記録媒体に格納して配付してもよい。そして、少なくともマイクロコンピュータ，パーソナルコンピュータ，汎用コンピュータを範疇に含むコンピュータが、上記の記録媒体から上記プログラムを読み出して、実行するものとしてもよい。

【符号の説明】

【0052】

Ｔ：サンプリング時間間隔；データ集録時設定
ωｉ：中心周波数；任意設定
ζｉ（ツェータ）；通過帯域比率；任意設定
ｘ：使用する時系列データ
ｙ：周波数成分の波形出力
入出力の添え字は、ｎ：現在の値、ｎ−ｊ：現在よりｊ：サンプル前の値
ｖ：ｙの絶対値
ｚ：ｖの移動平均値

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】