(57)【特許請求の範囲】
【請求項１】
予め異常診断対象から測定した時系列データを定Ｑ変換して多変量サンプルを作成し、前記作成された多変量サンプルを基に診断用パラメータを作成する診断用パラメータ作成手段と、
新たに異常診断対象から測定した時系列データを定Ｑ変換して多変量サンプルを作成し、前記診断用パラメータ作成手段により作成された診断用パラメータを用いて異常を診断する診断手段と、を備え、
前記診断用パラメータ作成手段は、
時系列データを定Ｑ変換して多変量サンプルを作成する第１の定Ｑ変換装置と、
前記作成された多変量サンプルを主成分分析して変換行列を得るとともに各サンプルの主成分得点を計算する主成分分析部と、
前記主成分得点を基にサンプルの指標となる統計量を得る第１の統計値計算部と、
前記統計量を正規化して異常度とする第１の正規化部と、を備え、
前記診断手段は、
時系列データを定Ｑ変換して多変量サンプルを作成する第２の定Ｑ変換装置と、
前記作成された多変量サンプルから、前記主成分分析部で得られた変換行列を使って主成分得点を得る主成分計算部と、
前記主成分得点を基に統計量を得る第２の統計値計算部と、
前記統計量を、前記第１の正規化部で使用した正規化係数によって正規化して異常度を計算する第２の正規化部と、を備え、前記異常度に基づいて異常診断対象の異常を診断し、
前記診断用パラメータ作成手段および診断手段の第１、第２の定Ｑ変換装置は、ともに、
設定した周波数成分ＫM以下の低周波数帯域において、時系列データの離散フーリエ変換と定Ｑ変換における周波数軸係数との積和計算により求めた定Ｑ変換成分の第１の演算量と、前記ＫMを超える高周波帯域において、時系列データと定Ｑ変換における時間軸係数との積和計算により求めた定Ｑ変換成分の第２の演算量とを含む合計演算量を算出し、
前記ＫMを全周波数成分Ｋに各々設定して各ＫM毎の前記合計演算量を求め、
前記各ＫM毎の合計演算量を評価し、該演算量が最小となるＫMを、計算法境界を示す周波数成分閾値として決定する計算法境界決定手段と、
定Ｑ変換における周波数成分Ｋが前記決定された周波数成分閾値ＫM以下の低周波帯域では、時系列データの離散フーリエ変換と前記周波数軸係数との積和計算である第１の計算法を選択し、前記周波数成分Ｋが前記決定された周波数成分閾値ＫMを超える高周波帯域では、時系列データと前記時間軸係数との積和計算である第２の計算法を選択する計算法選択手段と、
時系列データに対して、前記計算法選択手段により選択された第１の計算法、第２の計算法を実施して定Ｑ変換の成分を求める定Ｑ変換手段と、
を備えたことを特徴とする異常診断装置。
【請求項２】
前記各定Ｑ変換装置は、
前記定Ｑ変換における時間軸係数を、
【数１】
0006627639.tif 000025
によって算出する時間軸係数算出部と、
前記定Ｑ変換における周波数軸係数を、時間軸係数Ｔn,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
【数２】
0006627639.tif 000026
から（１／Ｎ）Ｓn,kとして算出する周波数軸係数算出部と、を備え、
前記定Ｑ変換手段は、
前記周波数軸係数算出部で算出された周波数軸係数と時系列データの離散フーリエ変換とを積和計算して前記第１の計算法を実施する低周波帯域計算部と、
前記時間軸係数算出部で算出された時間軸係数と時系列データとを積和計算して前記第２の計算法を実施する高周波帯域計算部と、
前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求める合成部と、
を備えていることを特徴とする請求項１に記載の異常診断装置。
【請求項３】
前記各定Ｑ変換装置は、
【数１】
0006627639.tif 000027
を、予め設定した時間方向の複数の解析位置ｐm（ｍ＝１，…，Ｍ）毎に算出して時間軸係数Ｔn,k,mを求め、該Ｔn,k,mを前記定Ｑ変換における時間軸係数とする時間軸係数算出部と、
時間軸係数Ｔn,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
【数２】
0006627639.tif 000028
から得た（１／Ｎ）Ｓn,kを、前記解析位置ｐm毎に算出して周波数軸係数（１／Ｎ）Ｓn,k,mを求め、該（１／Ｎ）Ｓn,k,mを前記定Ｑ変換における周波数軸係数とする周波数軸係数算出部と、を備え、
前記定Ｑ変換手段は、
時系列データｘn全体に対して離散フーリエ変換を行ってフーリエ係数Ｘnを求め、前記解析位置ｐm毎に、前記周波数軸係数算出部で算出された周波数軸係数（１／Ｎ）Ｓn,k,mと前記フーリエ係数Ｘnとを積和計算して前記第１の計算法を実施する低周波帯域計算部と、
前記解析位置ｐm毎に、前記時間軸係数算出部で算出された時間軸係数Ｔn,k,mと時系列データｘnとを積和計算して前記第２の計算法を実施する高周波帯域計算部と、
前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求める合成部と、
を備えていることを特徴とする請求項１に記載の異常診断装置。
【請求項４】
前記各定Ｑ変換装置は、
【数１】
0006627639.tif 000029
を、予め設定した時間方向の複数の解析位置ｐm（ｍ＝１，…，Ｍ）のうち１つの解析位置ｐ1において算出して時間軸係数Ｔn,k,1を求め、該Ｔn,k,1を前記定Ｑ変換における時間軸係数とする時間軸係数算出部と、
時間軸係数Ｔn,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
【数２】
0006627639.tif 000030
から得た（１／Ｎ）Ｓn,kを、前記１つの解析位置ｐ1において算出して周波数軸係数（１／Ｎ）Ｓn,k,1を求め、該（１／Ｎ）Ｓn,k,1を前記定Ｑ変換における周波数軸係数とする周波数軸係数算出部と、を備え、
前記定Ｑ変換手段は、
時系列データｘn全体に対して離散フーリエ変換を行ってフーリエ係数Ｘnを求め、  前記解析位置ｐm毎に、
前記周波数軸係数算出部で求められた周波数軸係数（１／Ｎ）Ｓn,k,1と、
前記解析位置ｐm毎の周波数軸係数（１／Ｎ）Ｓn,k,mは前記１つの解析位置ｐ1における周波数軸係数（１／Ｎ）Ｓn,k,1に帰結できる関係にあることから導かれる、複素指数関数のべき乗
【数３】
0006627639.tif 000031
と、
前記フーリエ係数Ｘnとを積和計算して前記第１の計算法を実施する低周波帯域計算部と、
前記解析位置ｐm毎に範囲をずらして、前記時間軸係数算出部で算出された時間軸係数Ｔn,k,1と時系列データｘn-p1+pmとを積和計算して前記第２の計算法を実施する高周波帯域計算部と、
前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求める合成部と、
を備えていることを特徴とする請求項１に記載の異常診断装置。
【請求項５】
診断用パラメータ作成手段が、予め異常診断対象から測定した時系列データを定Ｑ変換して多変量サンプルを作成し、前記作成された多変量サンプルを基に診断用パラメータを作成する診断用パラメータ作成ステップと、
診断手段が、新たに異常診断対象から測定した時系列データを定Ｑ変換して多変量サンプルを作成し、前記診断用パラメータ作成手段により作成された診断用パラメータを用いて異常を診断する診断ステップと、を備え、
前記診断用パラメータ作成ステップは、
第１の定Ｑ変換装置が、時系列データを定Ｑ変換して多変量サンプルを作成するステップと、
主成分分析部が、前記作成された多変量サンプルを主成分分析して変換行列を得るとともに各サンプルの主成分得点を計算するステップと、
第１の統計値計算部が、前記主成分得点を基にサンプルの指標となる統計量を得るステップと、
第１の正規化部が、前記統計量を正規化して異常度とするステップと、を備え、
前記診断ステップは、
第２の定Ｑ変換装置が、時系列データを定Ｑ変換して多変量サンプルを作成するステップと、
主成分計算部が、前記作成された多変量サンプルから、前記主成分分析部で得られた変換行列を使って主成分得点を得るステップと、
第２の統計値計算部が、前記主成分得点を基に統計量を得るステップと、
第２の正規化部が、前記統計量を、前記第１の正規化部で使用した正規化係数によって正規化して異常度を計算するステップと、を備え、前記異常度に基づいて異常診断対象の異常を診断し、
前記第１の定Ｑ変換装置が、時系列データを定Ｑ変換して多変量サンプルを作成するステップは、
設定した周波数成分ＫM以下の低周波数帯域において、時系列データの離散フーリエ変換と定Ｑ変換における周波数軸係数との積和計算により求めた定Ｑ変換成分の第１の演算量と、前記ＫMを超える高周波帯域において、時系列データと定Ｑ変換における時間軸係数との積和計算により求めた定Ｑ変換成分の第２の演算量とを含む合計演算量を算出し、
前記ＫMを全周波数成分Ｋに各々設定して各ＫM毎の前記合計演算量を求め、
前記各ＫM毎の合計演算量を評価し、該演算量が最小となるＫMを、計算法境界を示す周波数成分閾値として決定するステップを備え、
前記第２の定Ｑ変換装置が、時系列データを定Ｑ変換して多変量サンプルを作成するステップは、
定Ｑ変換における周波数成分Ｋが前記決定された周波数成分閾値ＫM以下の低周波帯域では、時系列データの離散フーリエ変換と前記周波数軸係数との積和計算である第１の計算法を選択し、前記周波数成分Ｋが前記決定された周波数成分閾値ＫMを超える高周波帯域では、時系列データと前記時間軸係数との積和計算である第２の計算法を選択するステップと、
時系列データに対して、前記選択された第１の計算法、第２の計算法を実施して定Ｑ変換の成分を求めるステップと、を備えたことを特徴とする異常診断方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、異常診断対象機器の動作時の振動・音響などの波形データ（時系列データ）を周波数解析して異常診断を行う異常診断装置および方法に関する。
【背景技術】
【０００２】
異常診断対象機器の異常を診断する方法は、例えば特許文献１に記載のものが提案されている。
【０００３】
この特許文献１の方法では、時系列の波形データを、短い時間毎に分割して、フーリエ変換して周波数成分の多変量データのサンプルを多数作り、機器が正常と考えられる期間に予め計測した多数の多変量サンプルを主成分分析して主成分得点を求めるための固有ベクトル（ローテーション行列ともいう）を用意する。
【０００４】
診断の際には、計測データを前記同様にフーリエ変換して周波数成分の多変量データ・サンプルを作り、これを用意しておいたローテーション行列で変換して主成分得点を求め、それを正常時の主成分得点と比較することで異常診断を行っている。要するに診断毎に独立に主成分分析するのではなく、予め主成分分析しておいた結果を利用して、診断時には同じ基準で変換することにより比較を容易にするのが要点である。
【０００５】
尚、本発明で利用する定Ｑ変換の成分演算方法は、例えば非特許文献１〜４に記載されている。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特許第３３８２２４０号公報
【非特許文献】
【０００７】
【非特許文献１】Ｊｕｄｉｔｈ  Ｃ．Ｂｒｏｗｎ：Ｃａｌｃｕｌａｔｉｏｎ  ｏｆ  ａ  ｃｏｎｓｔａｎｔ  Ｑ  ｓｐｅｃｔｒａｌ  ｔｒａｎｓｆｏｒｍ，Ｊ．Ａｃｏｕｓｔ．Ｓｏｃ．Ａｍ．８９（１）：４２５−４３４，１９９１
【非特許文献２】Ｊｕｄｉｔｈ  Ｃ．Ｂｒｏｗｎ  ａｎｄ  Ｍｉｌｌｅｒ  Ｓ．Ｐｕｃｋｅｔｔｅ：Ａｎ  ｅｆｆｉｃｉｅｎｔ  ａｌｇｏｒｉｔｈｍ  ｆｏｒ  ｔｈｅ  ｃａｌｃｕｌａｔｉｏｎ  ｏｆ  ａ  ｃｏｎｓｔａｎｔ  Ｑ  ｔｒａｎｓｆｏｒｍ，Ｊ．Ａｃｏｕｓｔ．Ｓｏｃ．Ａｍ．９２（５）：２６９８−２７０１，１９９２
【非特許文献３】ｙｕｋａｒａ＿１３：[Ｐｙｔｈｏｎ]Ｃｏｎｓｔａｎｔ−Ｑ変換（対数周波数スペクトログラム）、音楽プログラミングの超入門（仮）  ｉｎ  Ｈａｔｅｎａ  Ｂｌｏｇ，２０１３−１２−０１，２０１３、インターネット＜ＵＲＬ：ｈｔｔｐ：／／ｙｕｋａｒａ−１３．ｈａｔｅｎａｂｌｏｇ．ｃｏｍ／ｅｎｔｒｙ／２０１３／１２／０１／２２２７４２＞．[２０１６−０２−２２  アクセス]
【非特許文献４】ｙｕｋａｒａ＿１３：[Ｐｙｔｈｏｎ]  高速なＣｏｎｓｔａｎｔ−Ｑ変換（ｗｉｔｈ  ＦＦＴ）、音楽プログラミングの超入門（仮）  ｉｎ  Ｈａｔｅｎａ  Ｂｌｏｇ，２０１４−０１−０５，２０１４、インターネット＜ＵＲＬ：ｈｔｔｐ：／／ｙｕｋａｒａ−１３．ｈａｔｅｎａｂｌｏｇ．ｃｏｍ／ｅｎｔｒｙ／２０１４／０１／０５／０６２４１４＞．[２０１６−０２−２２  アクセス]
【発明の概要】
【発明が解決しようとする課題】
【０００８】
フーリエ変換を基にした多変量サンプルで予め作成したローテーション行列によって長期にわたって診断しようとすると、事前の解析には多くのデータを必要とする。例えば回転機の診断においては、少なくとも低い周波数領域では１Ｈｚ程度の周波数解像度が必要とされる。
【０００９】
このためには各サンプルは少なくとも１秒程度の測定が必要であり、時間を重複してサンプル数を増やしても本質的には類似サンプルとなって多様性が不足するため、それを基に診断すると、長期的には異常と単なる時間経過による変化の区別がつけられなくなる。
【００１０】
以下に、その事例を述べる。この事例は、ある回転機の振動を長期間（約４ヶ月）にわたって測定して分析したものである。１時間に一回、約１００ｋＨｚのサンプリング周波数で５秒間の測定を行い、これから時間を重複させながら１秒強の分割データを１測定から１００サンプル作成して、各サンプルをフーリエ変換して１／６オクターブバンドの約８０系列の多変量サンプルにした。この１週間分（ただし稼動している時刻のデータに絞ったため都合８千サンプルほど）を主成分分析してローテーション行列を計算し、これを基にホテリングＴ２／Ｑ統計量をベースとした異常度を図６のように算出している。
【００１１】
対象回転機は特に問題のないものであったが、ホテリングＴ２統計量をベースとしたＩｎｄｉｃａｔｏｒ１、Ｑ統計量をベースとしたＩｎｄｉｃａｔｏｒ２とも、図６の「線形」の傾きから分かるように、時間の経過とともに増加していることが伺える。この増加の割合は、学習期間（７日）の１０倍（７０日）もあれば平均値が異常になるほどであり、長期的な診断にはとても使えない。
【００１２】
この問題は異常度の定義に依存するものではなく、全期間データを主成分分析した図７の第一主成分×第二主成分の散布状況をみても、白色の丸印で示す学習期間のデータが、黒丸印で示す全データの範囲をカバーできていないことは明白である。
【００１３】
尚、図７において、横軸ＰＣ１は第一主成分、縦軸ＰＣ２は第二主成分を各々示している。
【００１４】
この問題は、時系列データから多変量データを作成するのにフーリエ変換を使っているためにデータの多様性が確保できないためと考えられる。
【００１５】
本発明は上記課題を解決するものであり、その目的は、サンプルの多様性を確保し、長期間にわたって信頼性の高い異常診断を実施することができる異常診断装置および方法を提供することにある。
【課題を解決するための手段】
【００１６】
上記課題を解決するための請求項１に記載の異常診断装置は、予め異常診断対象から測定した時系列データを定Ｑ変換して多変量サンプルを作成し、前記作成された多変量サンプルを基に診断用パラメータを作成する診断用パラメータ作成手段と、
新たに異常診断対象から測定した時系列データを定Ｑ変換して多変量サンプルを作成し、前記診断用パラメータ作成手段により作成された診断用パラメータを用いて異常を診断する診断手段と、を備え、
前記診断用パラメータ作成手段は、
時系列データを定Ｑ変換して多変量サンプルを作成する第１の定Ｑ変換装置と、
前記作成された多変量サンプルを主成分分析して変換行列を得るとともに各サンプルの主成分得点を計算する主成分分析部と、
前記主成分得点を基にサンプルの指標となる統計量を得る第１の統計値計算部と、
前記統計量を正規化して異常度とする第１の正規化部と、を備え、
前記診断手段は、
時系列データを定Ｑ変換して多変量サンプルを作成する第２の定Ｑ変換装置と、
前記作成された多変量サンプルから、前記主成分分析部で得られた変換行列を使って主成分得点を得る主成分計算部と、
前記主成分得点を基に統計量を得る第２の統計値計算部と、
前記統計量を、前記第１の正規化部で使用した正規化係数によって正規化して異常度を計算する第２の正規化部と、を備え、前記異常度に基づいて異常診断対象の異常を診断し、
前記診断用パラメータ作成手段および診断手段の第１、第２の定Ｑ変換装置は、ともに、
設定した周波数成分ＫM以下の低周波数帯域において、時系列データの離散フーリエ変換と定Ｑ変換における周波数軸係数との積和計算により求めた定Ｑ変換成分の第１の演算量と、前記ＫMを超える高周波帯域において、時系列データと定Ｑ変換における時間軸係数との積和計算により求めた定Ｑ変換成分の第２の演算量とを含む合計演算量を算出し、
前記ＫMを全周波数成分Ｋに各々設定して各ＫM毎の前記合計演算量を求め、
前記各ＫM毎の合計演算量を評価し、該演算量が最小となるＫMを、計算法境界を示す周波数成分閾値として決定する計算法境界決定手段と、
定Ｑ変換における周波数成分Ｋが前記決定された周波数成分閾値ＫM以下の低周波帯域では、時系列データの離散フーリエ変換と前記周波数軸係数との積和計算である第１の計算法を選択し、前記周波数成分Ｋが前記決定された周波数成分閾値ＫMを超える高周波帯域では、時系列データと前記時間軸係数との積和計算である第２の計算法を選択する計算法選択手段と、
時系列データに対して、前記計算法選択手段により選択された第１の計算法、第２の計算法を実施して定Ｑ変換の成分を求める定Ｑ変換手段と、
を備えたことを特徴とする。
【００１８】
また、請求項５に記載の異常診断方法は、診断用パラメータ作成手段が、予め異常診断対象から測定した時系列データを定Ｑ変換して多変量サンプルを作成し、前記作成された多変量サンプルを基に診断用パラメータを作成する診断用パラメータ作成ステップと、
診断手段が、新たに異常診断対象から測定した時系列データを定Ｑ変換して多変量サンプルを作成し、前記診断用パラメータ作成手段により作成された診断用パラメータを用いて異常を診断する診断ステップと、を備え、
前記診断用パラメータ作成ステップは、
第１の定Ｑ変換装置が、時系列データを定Ｑ変換して多変量サンプルを作成するステップと、
主成分分析部が、前記作成された多変量サンプルを主成分分析して変換行列を得るとともに各サンプルの主成分得点を計算するステップと、
第１の統計値計算部が、前記主成分得点を基にサンプルの指標となる統計量を得るステップと、
第１の正規化部が、前記統計量を正規化して異常度とするステップと、を備え、
前記診断ステップは、
第２の定Ｑ変換装置が、時系列データを定Ｑ変換して多変量サンプルを作成するステップと、
主成分計算部が、前記作成された多変量サンプルから、前記主成分分析部で得られた変換行列を使って主成分得点を得るステップと、
第２の統計値計算部が、前記主成分得点を基に統計量を得るステップと、
第２の正規化部が、前記統計量を、前記第１の正規化部で使用した正規化係数によって正規化して異常度を計算するステップと、を備え、前記異常度に基づいて異常診断対象の異常を診断し、
前記第１の定Ｑ変換装置が、時系列データを定Ｑ変換して多変量サンプルを作成するステップは、
設定した周波数成分ＫM以下の低周波数帯域において、時系列データの離散フーリエ変換と定Ｑ変換における周波数軸係数との積和計算により求めた定Ｑ変換成分の第１の演算量と、前記ＫMを超える高周波帯域において、時系列データと定Ｑ変換における時間軸係数との積和計算により求めた定Ｑ変換成分の第２の演算量とを含む合計演算量を算出し、
前記ＫMを全周波数成分Ｋに各々設定して各ＫM毎の前記合計演算量を求め、
前記各ＫM毎の合計演算量を評価し、該演算量が最小となるＫMを、計算法境界を示す周波数成分閾値として決定するステップを備え、
前記第２の定Ｑ変換装置が、時系列データを定Ｑ変換して多変量サンプルを作成するステップは、
定Ｑ変換における周波数成分Ｋが前記決定された周波数成分閾値ＫM以下の低周波帯域では、時系列データの離散フーリエ変換と前記周波数軸係数との積和計算である第１の計算法を選択し、前記周波数成分Ｋが前記決定された周波数成分閾値ＫMを超える高周波帯域では、時系列データと前記時間軸係数との積和計算である第２の計算法を選択するステップと、
時系列データに対して、前記選択された第１の計算法、第２の計算法を実施して定Ｑ変換の成分を求めるステップと、を備えたことを特徴とする。
【００２０】
上記構成によれば、時系列データから周波数成分の多変量データを作成する際に、フーリエ変換でなく定Ｑ変換を使ったことにより、診断用パラメータ作成手段で使用するデータが比較的少なくてもサンプルの多様性を確保できるようになり、長期間にわたって同じ基準で異常診断が可能になる。これによって異常診断の信頼性が向上する。
【００２２】
また、請求項２に記載の異常診断装置は、請求項１において、前記各定Ｑ変換装置は、前記定Ｑ変換における時間軸係数を、
【００２３】
【数１】
0006627639.tif 000002
【００２４】
によって算出する時間軸係数算出部と、
前記定Ｑ変換における周波数軸係数を、時間軸係数Ｔn,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
【００２５】
【数２】
0006627639.tif 000003
【００２６】
から（１／Ｎ）Ｓn,kとして算出する周波数軸係数算出部と、を備え、
前記定Ｑ変換手段は、
前記周波数軸係数算出部で算出された周波数軸係数と時系列データの離散フーリエ変換とを積和計算して前記第１の計算法を実施する低周波帯域計算部と、
前記時間軸係数算出部で算出された時間軸係数と時系列データとを積和計算して前記第２の計算法を実施する高周波帯域計算部と、
前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求める合成部と、
を備えていることを特徴とする。
【００２７】
上記構成によれば、低周波帯域、高周波帯域のどちらの場合でも演算量の少ない計算法で定Ｑ変換の成分を演算することができ、演算の高速化が実現される。
【００２８】
また、請求項３に記載の異常診断装置は、請求項１において、前記各定Ｑ変換装置は、
【００２９】
【数１】
0006627639.tif 000004
【００３０】
を、予め設定した時間方向の複数の解析位置ｐm（ｍ＝１，…，Ｍ）毎に算出して時間軸係数Ｔn,k,mを求め、該Ｔn,k,mを前記定Ｑ変換における時間軸係数とする時間軸係数算出部と、
時間軸係数Ｔn,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
【００３１】
【数２】
0006627639.tif 000005
【００３２】
から得た（１／Ｎ）Ｓn,kを、前記解析位置ｐm毎に算出して周波数軸係数（１／Ｎ）Ｓn,k,mを求め、該（１／Ｎ）Ｓn,k,mを前記定Ｑ変換における周波数軸係数とする周波数軸係数算出部と、を備え、
前記定Ｑ変換手段は、
時系列データｘn全体に対して離散フーリエ変換を行ってフーリエ係数Ｘnを求め、前記解析位置ｐm毎に、前記周波数軸係数算出部で算出された周波数軸係数（１／Ｎ）Ｓn,k,mと前記フーリエ係数Ｘnとを積和計算して前記第１の計算法を実施する低周波帯域計算部と、
前記解析位置ｐm毎に、前記時間軸係数算出部で算出された時間軸係数Ｔn,k,mと時系列データｘnとを積和計算して前記第２の計算法を実施する高周波帯域計算部と、
前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求める合成部と、
を備えていることを特徴とする。
【００３３】
上記構成によれば、複数の時刻で各々定Ｑ変換を行う際に、低周波帯域における周波数軸での計算は、時系列データの離散フーリエ変換を、複数の時刻毎に行わなくても、ｘn全体に対して一括して離散フーリエ変換したフーリエ係数Ｘnと各解析位置での周波数軸係数との積和計算でよいため、全体としての計算が高速化される。
【００３４】
また、請求項４に記載の異常診断装置は、請求項１において、前記各定Ｑ変換装置は、
【００３５】
【数１】
0006627639.tif 000006
【００３６】
を、予め設定した時間方向の複数の解析位置ｐm（ｍ＝１，…，Ｍ）のうち１つの解析位置ｐ1において算出して時間軸係数Ｔn,k,1を求め、該Ｔn,k,1を前記定Ｑ変換における時間軸係数とする時間軸係数算出部と、
時間軸係数Ｔn,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
【００３７】
【数２】
0006627639.tif 000007
【００３８】
から得た（１／Ｎ）Ｓn,kを、前記１つの解析位置ｐ1において算出して周波数軸係数（１／Ｎ）Ｓn,k,1を求め、該（１／Ｎ）Ｓn,k,1を前記定Ｑ変換における周波数軸係数とする周波数軸係数算出部と、を備え、
前記定Ｑ変換手段は、
時系列データｘn全体に対して離散フーリエ変換を行ってフーリエ係数Ｘnを求め、  前記解析位置ｐm毎に、
前記周波数軸係数算出部で求められた周波数軸係数（１／Ｎ）Ｓn,k,1と、
前記解析位置ｐm毎の周波数軸係数（１／Ｎ）Ｓn,k,mは前記１つの解析位置ｐ1における周波数軸係数（１／Ｎ）Ｓn,k,1に帰結できる関係にあることから導かれる、複素指数関数のべき乗
【００３９】
【数３】
0006627639.tif 000008
【００４０】
と、
前記フーリエ係数Ｘnとを積和計算して前記第１の計算法を実施する低周波帯域計算部と、
前記解析位置ｐm毎に範囲をずらして、前記時間軸係数算出部で算出された時間軸係数Ｔn,k,1と時系列データｘn-p1+pmとを積和計算して前記第２の計算法を実施する高周波帯域計算部と、
前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求める合成部と、
を備えていることを特徴とする。
【００４１】
上記構成によれば、複数の時刻で各々定Ｑ変換を行う際に、１つの解析位置ｐ1における時間軸係数および周波数軸係数を利用して各時刻での計算が可能となるので、複数の時刻毎（複数の解析位置毎）に時間軸係数および周波数軸係数を用意する必要はなく、データ量が少なくなってメモリを節約することができる。
【発明の効果】
【００４２】
（１）請求項１〜５に記載の発明によれば、時系列データから周波数成分の多変量データを作成する際に、フーリエ変換でなく定Ｑ変換を使ったことにより、診断用パラメータ作成手段で使用するデータが比較的少なくてもサンプルの多様性を確保できるようになり、長期間にわたって同じ基準で異常診断が可能になる。これによって異常診断の信頼性が向上する。
（２）請求項１、２に記載の発明によれば、低周波帯域、高周波帯域のどちらの場合でも演算量の少ない計算法で定Ｑ変換の成分を演算することができ、演算の高速化が実現される。
（３）請求項３に記載の発明によれば、複数の時刻で各々定Ｑ変換を行う際に、低周波帯域における周波数軸での計算は、時系列データの離散フーリエ変換を、複数の時刻毎に行わなくても、ｘn全体に対して一括して離散フーリエ変換したフーリエ係数Ｘnと各解析位置での周波数軸係数との積和計算でよいため、全体としての計算が高速化される。
（４）請求項４に記載の発明によれば、複数の時刻で各々定Ｑ変換を行う際に、１つの解析位置ｐ1における時間軸係数および周波数軸係数を利用して各時刻での計算が可能となるので、複数の時刻毎（複数の解析位置毎）に時間軸係数および周波数軸係数を用意する必要はなく、データ量が少なくなってメモリを節約することができる。
【図面の簡単な説明】
【００４３】
【図１】本発明の実施形態例による異常診断装置の構成を示すブロック図。
【図２】本発明の実施形態例による全期間データの第一主成分×第二主成分散布図。
【図３】本発明の実施形態例による異常度の期間変化を示す説明図。
【図４】本発明の実施形態例による異常診断例を示す説明図。
【図５】本発明の実施形態例における定Ｑ変換装置の一例を示す構成図。
【図６】従来手法をベースとした異常度の期間変化を示す説明図。
【図７】従来のフーリエ変換による全期間データの第一主成分×第二主成分散布図。
【発明を実施するための形態】
【００４４】
以下、図面を参照しながら本発明の実施の形態を説明するが、本発明は下記の実施形態例に限定されるものではない。本発明では、時系列データ（波形データ）から周波数成分の多変量サンプルを生成する際に、従来のフーリエ変換の代わりに定Ｑ変換を使う。
【００４５】
定Ｑ変換（Ｃｏｎｓｔａｎｔ  Ｑ  ｔｒａｎｓｆｏｒｍ）は、フーリエ変換のように全ての周波数帯域で同じ期間のデータを解析するのではなく、全ての周波数帯域で同じ周期数になるように、周波数毎に参照するデータ数を変えて解析する。この方法では高周波帯域ほど短い期間のデータで解析するため、低周波帯域での周波数解像度を高く（例えば１Ｈｚ程度）するために長い期間（１秒以上）の測定を必要とする際に、期間を重複させ多数のサンプルをとっても高周波帯域では期間が重複しないため、多変量サンプルとしては多様性が確保される。
【実施例１】
【００４６】
具体的な診断手法の構成を図１に示す。図１において、本実施例１の手法は、予め測定して集めた時系列データから診断のためのパラメータを計算する準備フェーズと、新たに測定された時系列データを診断する診断フェーズからなる。
【００４７】
準備フェーズは、予め異常診断対象から測定して集めた時系列データを定Ｑ変換して多変量サンプルを作成し、前記作成された多変量サンプルを基に診断用パラメータを作成する診断用パラメータ作成手段１００として構成している。
【００４８】
診断フェーズは、新たに異常診断対象から測定した時系列データを定Ｑ変換して多変量サンプルを作成し、前記診断用パラメータ作成手段により作成された診断用パラメータを用いて異常を診断する診断手段２００として構成している。
【００４９】
診断用パラメータ作成手段１００は、予め異常診断対象から測定して集めた時系列データを定Ｑ変換して多変量サンプルを作成する第１の定Ｑ変換装置１０１と、前記作成された多変量サンプルを主成分分析して変換行列を得るとともに各サンプルの主成分得点を計算する主成分分析部１０２と、前記主成分得点を基にサンプルの指標となる統計量を得る第１の統計値計算部１０３と、前記統計量を正規化して異常度とする第１の正規化部１０４と、を備えている。
【００５０】
前記診断手段２００は、新たに測定された時系列データを定Ｑ変換して多変量サンプルを作成する第２の定Ｑ変換装置２０１と、前記作成された多変量サンプルから、前記主成分分析部１０２で得られた変換行列を使って主成分得点を得る主成分計算部２０２と、前記主成分得点を基に統計量を得る第２の統計値計算部２０３と、前記統計量を、前記第１の正規化部１０４で使用した正規化係数によって正規化して異常度を計算する第２の正規化部２０４と、を備え、前記異常度に基づいて異常診断対象の異常を診断するものである。
【００５１】
図１の異常診断装置は、例えばコンピュータにより構成され、通常のコンピュータのハードウェアリソース、例えばＲＯＭ、ＲＡＭ、ＣＰＵ、入力装置、出力装置、通信インターフェース、ハードディスク、記録媒体およびその駆動装置を備えている。
【００５２】
このハードウェアリソースとソフトウェアリソース（ＯＳ、アプリケーションなど）との協働の結果、異常診断装置は、図１のように第１の定Ｑ変換装置１０１、主成分分析部１０２、第１の統計値計算部１０３、第１の正規化部１０４、第２の定Ｑ変換装置２０１、主成分計算部２０２、第２の統計値計算部２０３、第２の正規化部２０４を実装する。
【００５３】
前記多変量サンプル、主成分得点、統計量、異常度、変換行列、正規化係数などのデータはハードディスクあるいはＲＡＭなどの保存手段・記憶手段に格納されるものとする。
【００５４】
次に、図１の装置の各部の動作を説明する。準備フェーズ（診断用パラメータ作成手段１００）では、予め集めた時系列データを、第１の定Ｑ変換装置１０１において順次定Ｑ変換して多数の多変量サンプルを作成する。定Ｑ変換装置１０１は、例えば非特許文献１〜４に開示されている方式を実行するものでもよく、また後述する実施例２〜４の高速演算が可能な定Ｑ変換装置を採用しても良い。尚、実施例２〜４に示す方法で高速に求める場合には、時系列データのデータ長はそろっていることが望ましい。
【００５５】
例えば、１０２．４ｋＨｚサンプリングで５秒間のデータが１００件など。定Ｑ変換の高周波領域での参照データ数の少なさにより、１回の計測から回転周期との位相や電源周期との位相がそれぞれ異なる１００〜４００程度の多様なサンプルを得ることができる。
【００５６】
次に、主成分分析部１０２において、多変量サンプルを主成分分析して変換行列を得るとともに、各サンプルの主成分得点を計算する。そして統計値計算部１０３が、主成分得点を基にホテリングＴ２／Ｑ統計量のようなサンプルの指標となる統計量を得る。
【００５７】
最後に、算出した統計量が準備フェーズの時系列データで平均０になるように、第１の正規化部１０４において統計量を正規化して異常度とする。この際に、正規化に使用した正規化係数を保存して診断フェーズ（診断手段２００）で使う。尚、異常度の取り得る数値範囲を限定するために、異常度を統計量の対数値にすることが望ましい。
【００５８】
診断フェーズ（診断手段２００）では、新たに測定した時系列データを第２の定Ｑ変換装置２０１において定Ｑ変換して多変量サンプルを作成する。
【００５９】
定Ｑ変換装置２０１は、定Ｑ変換装置１０１と同様に、例えば非特許文献１〜４に開示されている方式を実行するものでもよく、また後述する実施例２〜４の高速演算が可能な定Ｑ変換装置を採用しても良い。
【００６０】
次に、多変量サンプルに準備フェーズの主成分分析部１０２で作成した変換行列を使って、主成分計算部２０２で主成分得点を得る。そしてこれを基に、第２の統計値計算部２０３が準備フェーズと同様の計算で統計量を得る。最後に準備フェーズの第１の正規化部１０４で作成した正規化係数を使って、第２の正規化部２０４が異常度を計算し、その値により機器が正常か異常かを診断する。
【００６１】
異常診断に使う異常値の閾値を決めるのは一般的に困難であるが、異常値を統計値の対数で取る場合には、簡易の診断としては閾値をｌｏｇ３、ｌｏｇ５のような値としても良い。これは振動実効値での異常判断基準が通常時の３倍、５倍などとするのと同様である。
【００６２】
本発明の効果を示すため、図７の第一主成分×第二主成分散布図と同じデータについて、定Ｑ変換による多変量化データを主成分分析した場合の第一主成分×第二主成分散布図を図２に示す。図２によれば、図７の場合と異なり、白色の丸印で示す学習期間のデータが、黒丸印で示す全データの分布範囲を概ねカバーできていることが分かる。
【００６３】
また、図６に示したものと同様の異常度を提案手法により計算した結果を図３に示す。図３によれば、図６の場合と異なり、時間がたっても異常度が上昇傾向にないことが分かる。線形回帰による傾きはほとんどゼロであり、正の傾きを取るＩｎｄｉｃａｔｏｒ１でも、異常度の平均が閾値に達するのに１０年ほど掛かる。
【００６４】
次に、提案手法による異常診断例を図４に示す。図４は、図６や図３とは別の回転機の計算例である。こちらも同様に１時間に一回、約１００ｋＨｚのサンプリング周波数で５秒間の測定を行ったが、測定開始から２５日目に故障により緊急停止した。２５日目は稼動直後から異常度が非常に大きくなっているが、良く見ると、１３日目にＩｎｄｉｃａｔｏｒ１が注意ラインを越えていることが分かる。  
【００６５】
尚、この診断結果はフーリエ変換による手法でもほぼ同様の結果が得られる。
【実施例２】
【００６６】
本実施例２は、図１の第１、第２の定Ｑ変換装置１０１、２０１を、より精度の高い定Ｑ変換を、少ない演算量で高速に行うことができる図５の定Ｑ変換装置で構成したものである。
【００６７】
まず定Ｑ変換の成分について説明する。  
【００６８】
非特許文献１によると、定Ｑ変換は次で定式化される（表記は少し変更している）。
【００６９】
【数４】
0006627639.tif 000009
【００７０】
０以外の値は各ｋに対して、ｎ＝（Ｎ−Ｎk）／２，…，（Ｎ＋Ｎk）／２−１の範囲でのみ取り、この範囲の窓関数の形状として具体的にはハミング窓、ハニング窓、矩形窓など、多くの種類が提案されている
Ｎk：窓幅，Ｎk＝Ｑ（ｆs／ｆk）となるように取る
Ｎ≧Ｎ1≧…≧Ｎk…≧ＮK≧２Ｑとなる
Ｑ：窓の周期数，定Ｑ変換ではこれを全ての周波数で共通にする
非特許文献１、２ではＱ＝（２1/B−１）-1程度に取るようにしている
Ｂはビン数で、１オクターブを分割する数のこと
ｆs：サンプリング周波数
ｆk：第ｋ成分の周波数，ｆk＝（２1/B）k-1ｆmin
ｆmin＝ｆ1は定Ｑ変換で解析する最小の周波数で、解析する範囲から決める
ｆkは解析する最大の周波数で、ｆk≦ｆs／２となるようにする。
【００７１】
定Ｑ変換は特に低周波帯域の成分で窓幅Ｎkが非常に大きくなることにより演算量が多い。
【００７２】
これに対して非特許文献２では、高速フーリエ変換（ＦＦＴ；Ｆａｓｔ  Ｆｏｕｒｉｅｒ  ｔｒａｎｓｆｏｒｍ）を利用した演算の高速化方法が提案されている。
【００７３】
この方法ではパーセバルの公式により、定Ｑ変換を次の数式（２）のように時間軸での計算から周波数軸での計算に置き換える。
【００７４】
【数５】
0006627639.tif 000010
【００７５】
尚、明細書中の以下の文章では、前記数式（１）、（２）の各左辺の定Ｑ変換の成分を、「Ｘkcq」と表記することもある。
【００７６】
ここでＳn,kは（小さなｋに対して特に）ほとんどの要素が０に近いので、一定値以下の項を０として疎行列表現すると（実質非ゼロ項の積和計算だけになるため）非常に高速に計算できるようになる。
【００７７】
さらに（１／Ｎ）Ｓn,kは、計算するデータ長と解析する周波数範囲が決まっていれば、入力データ列ｘnに寄らず事前に計算できるためこれを予め用意しておけば、入力データの離散フーリエ変換はＦＦＴで実行するとして、定Ｑ変換を入力データに対する１回のＦＦＴと、ＦＦＴ結果ベクトルと予め用意した疎行列との積演算により計算できることになる。
【００７８】
【数６】
0006627639.tif 000011
【００７９】
離散でないフーリエ変換でも同等の公式が成り立ち、これらはパーセバルの定理、プランシュレルの定理など様々な呼び方がされる。
【００８０】
非特許文献２による高速化のアイデアは周波数軸での係数Ｓn,kが疎行列とみなせることに依拠しているが、実際にはＳn,kが疎行列というのは粗い評価である。確かにその値は一部を除いて小さいが、精度を高めると必ずしも疎にならない。特に高周波帯域でこの傾向が顕著である。実際、非特許文献２、４で係数Ｓn,kが十分小さいとしている閾値はそれぞれ０．１５、０．００５４であり、通常の評価には問題ないが有効数字７桁以上の精度を得るには十分とは言えない。このため、要求精度によってはこの高速化手法がかえって遅くなる可能性がある。
【００８１】
実際に周波数軸での計算でどの程度の演算量を必要とするかを調べるため、非特許文献４に準じた設定により時間軸での計算と周波数軸での計算でそれぞれ計算対象となる非ゼロの項数を評価する。
【００８２】
設定パラメータは、本明細書の記号法で表記してＢ＝２４、Ｑ＝２８、ｆs＝１６０００、ｆmin＝６０、Ｋ＝１６０とし、定Ｑ変換をそれぞれ矩形窓、ハミング窓、ハニング窓で計算する場合において、時間軸での計算に必要な非ゼロ項数と、周波数軸での計算で所定の精度の係数まで計算する際に必要な項数を表１に示す。
【００８３】
【表１】
0006627639.tif 000012
【００８４】
なお、非特許文献１、２に合わせるとＢ＝２４のときＱ＝３４となるが、非特許文献３、４では追加の係数ｆｒａｔｉｏを導入してＱの値を変えているので、その設定値Ｑ＝２８を使っている。また解析する最大周波数６０００ＨｚからＫ＝１６０を導いた。ｆsの値は非特許文献３，４に明記されていなかったので、最大周波数６０００Ｈｚまで解析できるサンプリング周波数１２０００Ｈｚ以上のＷＡＶ音声データとしてｆs＝１６０００Ｈｚを選択した。
【００８５】
表１中の周波数軸右の数値は係数をゼロと評価する係数の大きさの閾値であり、各数値右のパーセンテージは矩形窓の時間軸計算の項数を基準としている。
【００８６】
表１から、窓形状により程度が違うが、精度が高いほど（係数をゼロと評価する係数Ｓn,kの大きさの閾値が小さいほど）周波数軸での計算に必要な項数が増えることが分かる。特に矩形窓やハミング窓では精度が高くなると時間軸での計算より周波数軸での計算の方が積和すべき項数が大きくなり、周波数軸での計算で高速にならなくなる。
【００８７】
尚、周波数軸での計算にはＦＦＴ計算と疎行列の非ゼロ係数抽出のオーバーヘッドもあるので、必要項数に差がなければ時間軸計算が有利である。
【００８８】
より精度の高い定Ｑ変換を高速に行うには工夫が必要になる。  
【００８９】
そこで本実施例２では、定Ｑ変換を高精度でも演算量を少なく計算するために、定Ｑ変換の成分Ｘkcqを数式（２）の時間軸での計算と周波数軸での計算のどちらの方法で計算するか、予め周波数成分ｋ毎に演算量の小さい方を選択しておくことで実現する。
【００９０】
定Ｑ変換の成分Ｘkcqの計算はおおよそ低周波帯域では周波数軸での計算の演算量が小さく、高周波帯域では時間軸での計算の演算量が小さいから、ｋの境界値ＫMを決めて、ｋ≦ＫMでは周波数軸での計算、ｋ＞ＫMでは時間軸での計算、というようにすれば良い。
【００９１】
図５において、１０は予め定Ｑ変換を行うためのパラメータを決定して定Ｑ変換計算の準備をする準備部であり、２０は個別の時系列データに対する定Ｑ変換を計算する個別計算部である。
【００９２】
準備部１０は、
設定した周波数成分ＫM以下の低周波数帯域において、時系列データの離散フーリエ変換と定Ｑ変換における周波数軸係数との積和計算により求めた定Ｑ変換成分の第１の演算量と、前記ＫMを超える高周波帯域において、時系列データと定Ｑ変換における時間軸係数との積和計算により求めた定Ｑ変換成分の第２の演算量とを含む合計演算量を算出し、前記ＫMを全周波数成分ｋに各々設定して各ＫM毎の前記合計演算量を求め、前記各ＫM毎の合計演算量を評価し、該演算量が最小となるＫMを、計算法境界を示す周波数成分閾値として決定する計算法境界決定手段の機能と、
定Ｑ変換における周波数成分ｋが前記決定された周波数成分閾値ＫM以下の低周波帯域では、時系列データの離散フーリエ変換と前記周波数軸係数との積和計算である第１の計算法を選択し、前記周波数成分ｋが前記決定された周波数成分閾値ＫMを超える高周波帯域では、時系列データと前記時間軸係数との積和計算である第２の計算法を選択する計算法選択手段の機能とを具備している。
【００９３】
前記計算法境界決定手段および計算法選択手段の各機能を実行する準備部１０は、具体的には、定Ｑ変換のパラメータを決定するパラメータ決定部１１、時間軸係数Ｔn,kを計算する時間軸係数算出部１２、周波数軸係数（１／Ｎ）Ｓn,kを計算する周波数軸係数算出部１３、時間軸での計算と周波数軸での計算との境界のｋ値ＫMを決定する計算法境界決定部１４からなる。
【００９４】
個別計算部２０は、時系列データに対して、前記計算法選択手段により選択された第１の計算法、第２の計算法を実施して定Ｑ変換の成分を求める定Ｑ変換手段の機能を具備している。
【００９５】
前記定Ｑ変換手段の機能を実行する個別計算部２０は、具体的には、時系列データの長さを整える前処理部２１、高周波帯域の定Ｑ変換を計算する高周波帯域計算部２２、低周波帯域の定Ｑ変換を計算する低周波帯域計算部２３、高周波帯域での結果と低周波帯域での結果をまとめる定Ｑ変換合成部２４からなる。
【００９６】
図５の定Ｑ変換の成分演算装置は、例えばコンピュータにより構成され、通常のコンピュータのハードウェアリソース、例えばＲＯＭ、ＲＡＭ、ＣＰＵ、入力装置、出力装置、通信インターフェース、ハードディスク、記録媒体およびその駆動装置を備えている。
【００９７】
このハードウェアリソースとソフトウェアリソース（ＯＳ、アプリケーションなど）との協働の結果、定Ｑ変換の成分演算装置は、図５に示すように、パラメータ決定部１１、時間軸係数算出部１２、周波数軸係数算出部１３、計算法境界決定部１４、前処理部２１、高周波帯域計算部２２、低周波帯域計算部２３、定Ｑ変換合成部２４を実装する。
【００９８】
また各部での処理結果は、ハードディスク或いはＲＡＭなどの保存手段・記憶手段に格納され、随時利用されるものとする。
【００９９】
まず、準備部１０のパラメータ決定部１１では、予め定Ｑ変換で解析する時系列データのサンプリング周波数ｆs、解析する周波数範囲（ｆmin，ｆmax）と１オクターブの周波数を区切るビン数Ｂ、解析対象データの周期数Ｑ、窓関数を決定する。このときｆmax≦ｆs／２が必要で、Ｋ＝ｆｌｏｏｒ（Ｂ＊ｌｏｇ2（ｆmax／ｆmin））＋１となる。これを基に計算する入力データの長さをＮ＝Ｎ1＝Ｑ（ｆs／ｆmin）とする。
【０１００】
また周波数軸係数をゼロとみなすべき係数の大きさの閾値Ｔｈもここで決定しておく。
【０１０１】
時間軸係数算出部１２では、パラメータ決定部１１で決定したパラメータを基に、時間軸係数Ｔn,kを次の定義式に従って計算する。
【０１０２】
【数１】
0006627639.tif 000013
【０１０３】
周波数軸係数算出部１３では、パラメータ決定部１１で決定したパラメータを基に、周波数軸係数（１／Ｎ）Ｓn,kを、
【０１０４】
【数２】
0006627639.tif 000014
【０１０５】
の定義式に従って計算する。この際、｜Ｓn,k｜≦Ｔｈとなる項は周波数軸係数（１／Ｎ）Ｓn,k＝０とする。
【０１０６】
計算法境界決定部１４では、周波数軸での計算を行う最大のｋ値ＫM（周波数成分閾値）を求める。ＫMは０，…，Ｋの何れかの値を持ち、ＫM＝０は周波数軸での計算を行わずに非特許文献１、２と同様に全て時間軸での計算で求める場合を意味する。
【０１０７】
周波数成分閾値ＫMの決定方法としては、ＫM＝０，…，Ｋの各場合について、定Ｑ変換の成分Ｘkcqを求めるのに必要な演算量ＣX（ＫM）を評価してそれが最小となるようにＫMを決定する。
【０１０８】
定Ｑ変換の成分Ｘkcqの演算量ＣX（ＫM）（合計演算量）は、長さＮのＦＦＴの演算量ＣFFT（Ｎ）、非ゼロの時間軸係数毎の演算量ＣT、非ゼロの時間軸係数の数ＮT（ｋ）、非ゼロの周波数軸係数毎の演算量ＣF、ｋで非ゼロとみなす周波数軸係数の数ＮF（ｋ）を使っておおよそ
【０１０９】
【数７】
0006627639.tif 000015
【０１１０】
と表せる。
【０１１１】
ここでＣFFTは計算を実行する計算機やそれぞれの算法の詳細により変わり得るが、実測からの一例としては、ＣFFT（Ｎ）＝（１／４）Ｎｌｏｇ2（Ｎ）、ＣT＝１、ＣF＝１としても良い。
【０１１２】
上記演算量ＣX（ＫM）を示す数式（３）において、右辺第１項は、ＫMが非ゼロのときＣFFTとし、そうでないときは０とすることを表し、右辺第２項は、周波数成分ｋが周波数成分閾値ＫM以下の低周波帯域における周波数軸係数毎の演算量を表し、右辺第３項は、周波数成分ｋが周波数成分閾値ＫMを超える高周波帯域における時間軸係数毎の演算量を表している。
【０１１３】
このためＫMが非ゼロのときは、数式（３）の右辺第１項のＦＦＴの演算量および右辺第２項の演算量（周波数軸係数による演算量）からなる第１の演算量と、数式（３）の右辺第３項の時間軸係数による第２の演算量との合計が演算量ＣX（ＫM）となる。
【０１１４】
またＫMがゼロのときは、数式（３）の右辺第１項および右辺第２項がともに０となって右辺第３項の時間軸係数による第２の演算量のみが演算量ＣX（ＫM）となる。
【０１１５】
この演算量ＣX（ＫM）は、非ゼロの時間軸係数の数ＮT（ｋ）がｋに関して指数的に減少し、非ゼロの周波数軸係数の数ＮF（ｋ）がｋに関しておおよそ指数的に増加するため、ＫM＝１，…，Ｋの範囲ではおおよそ下に凸の特性カーブとなる。ただし、ＫM＝０ではＣFFTがなくなるためＣX（０）が小さくなる。結果として演算量ＣX（ＫM）はＫM＝０あるいはＫM＝１，…，Ｋでの最小点の何れかで最小値を取ることになる。
【０１１６】
個別計算部２０では、入力された時系列データの定Ｑ変換を次のようにして求める。  
【０１１７】
前処理部２１では、入力された時系列データの長さがＮより短ければ前後に０を補完して、Ｎより長ければ前後の値を除去して長さがＮの時系列データｘnになるようにする。このとき入力された時系列データの中央部分が時系列データｘnの中央部分に一致するようにする。
【０１１８】
次に定Ｑ変換の成分Ｘkcqを、前記計算法境界決定部１４で決定した周波数成分閾値ＫMにより、ｋ≦ＫMについては低周波帯域計算部２３で計算し、ｋ＞ＫMについては高周波帯域計算部２２で計算する。
【０１１９】
高周波帯域計算部２２は、前記時間軸係数算出部１２で算出された時間軸係数Ｔn,kと時系列データｘnを積和計算して定Ｑ変換の成分Ｘkcqを求める。各ｋについて非ゼロの時間軸係数Ｔn,kは高々Ｎk個なので、積和計算を非ゼロの範囲に限定することで効率的に計算できる。
【０１２０】
低周波帯域計算部２３は、まず時系列データｘnをＦＦＴしてフーリエ係数Ｘnを求める。次に前記周波数軸係数算出部１３で算出された周波数軸係数（１／Ｎ）Ｓn,kとフーリエ係数Ｘnを積和計算して定Ｑ変換の成分Ｘkcqを求める。各ｋについて非ゼロの周波数軸係数（１／Ｎ）Ｓn,kは少数なので、積和計算を非ゼロの範囲に限定することで効率的に計算できる。
【０１２１】
尚、周波数軸係数は多くの場合、両端（周波数ゼロ近辺とサンプリング周波数近辺）部分では非ゼロ項とゼロ項が入り混じるが中間部分ではゼロ項が続くので、非ゼロ項の範囲を少し広げて、両端部分は全て非ゼロ項として計算しても良い。この場合、非ゼロ項が少し増えるがその数はそれほど多くはなく、また非ゼロ判定ロジックが簡略化できるので演算時間は大差なくなる。
【０１２２】
最後に定Ｑ変換合成部２４で、高周波帯域計算部２２と低周波帯域計算部２３とで計算した定Ｑ変換の成分Ｘkcqをまとめて、定Ｑ変換を完了する。
【０１２３】
上記のように本実施例２によれば、従来の高速化技術の問題点である、高い精度の定Ｑ変換を計算しようとすると演算量が膨れ上り通常の時間軸での計算よりもむしろ遅くなるという問題が改善される。
【０１２４】
その効果を示すため、定Ｑ変換の計算方法（時間軸計算、周波数軸計算、提案手法）毎に、必要な演算量を表２のように試算した。試算は、時系列データのパラメータを表１と同じにして、演算量パラメータをＣFFT（Ｎ）＝（１／４）Ｎｌｏｇ2（Ｎ）、ＣT＝１、ＣF＝１として計算した定Ｑ変換の成分Ｘkcqの演算量ＣX（ＫM）を表２に示す。
【０１２５】
【表２】
0006627639.tif 000016
【０１２６】
表２中の周波数軸右の数値は係数をゼロと評価する係数の大きさの閾値であり、各数値右のパーセンテージは矩形窓の時間軸計算の演算量を基準としている。
【０１２７】
表２では、ＦＦＴ演算の演算量が原因で高速化の効果が小さめに出ているが、ハニング窓では閾値１．０ｅ−５でもおおよそ１／５（２２．５％）の演算量になっている。またＦＦＴ演算の演算量が重く周波数軸での計算で効果が出ない場合でも、自動的に時間軸での計算になるので反って遅くなるという問題は発生しない。
【実施例３】
【０１２８】
実施例２は、ある一時刻における定Ｑ変換の各成分値を求める際には効率的であるが、実際の定Ｑ変換では、時系列データの周波数成分の時間変化を捉えるために少しずつ時刻をずらしながら繰り返し定Ｑ成分を計算する。このときずらす時刻は、最短ではサンプリング周波数で、時系列データの情報を余すことなく計算するには最長でも最高周波数での窓幅にしたい。もちろん、時系列データの特徴サンプル抽出などでは、これよりも長い間隔で時刻をずらして定Ｑ変換する場合もあるが、いずれにせよ時系列データ全体に対しては多くの繰り返し計算をすることになり、まだまだ負荷が高い。
【０１２９】
多数の時刻で同時に定Ｑ変換を計算する際には更なる効率化が必要である。  
【０１３０】
そこで本実施例３では、数式（２）の周波数軸での計算で定Ｑ変換を行う場合に、定Ｑ変換を計算する時刻ごとに計算に使用する時系列データの範囲でＤＦＴ（あるいはＦＦＴ）する代わりに、複数時刻での定Ｑ変換に必要な範囲全体で一括してＤＦＴ（あるいはＦＦＴ）を行うことで高速化を図った。
【０１３１】
数式（１）で定義する定Ｑ変換はデータ列ｘn，ｎ＝０，…，Ｎ−１の中央ｘN/2での成分Ｘkcqを計算するように定式化しているが、窓関数の非ゼロ範囲を時刻方向に平行移動することで他の時刻での成分を計算するように調整できる。すなわちＴn,kとｘnとの積和のインデックスをずらすことで実現できる。
【０１３２】
数式（２）による周波数軸での計算に関しては、初めに定Ｑ変換する時刻全体で必要な全データに対して一括してＤＦＴ（あるいはＦＦＴ）を行い、予め各時刻用に用意した（１／Ｎ）Ｓn,kと順次積和することで実現する。
【０１３３】
本実施例３による定Ｑ変換の成分演算装置の具体的な計算の構成は、図５と同様であるが、各部で実行される機能が以下に示すように異なる。
【０１３４】
まず、準備部１０のパラメータ決定部１１では、予め定Ｑ変換で解析する時系列データのサンプリング周波数ｆｓ、解析する周波数範囲（ｆmin，ｆmax）と１オクターブの周波数を区切るビン数Ｂ、解析対象データの周期数Ｑ、窓関数を決定する。このときｆmax≦ｆs／２が必要で、Ｋ＝ｆｌｏｏｒ（Ｂ＊ｌｏｇ2（ｆmax／ｆmin））＋１となる。
【０１３５】
さらに入力データｘnの全長Ｎと、解析位置ｐm，ｍ＝１，…，Ｍ（時間方向の解析位置；時刻）を決定しておく。このとき、各位置ｐmでの定Ｑ変換はｘpmを中心とする長さＮ1＝Ｑ（ｆs／ｆmin）のデータ列をもとに計算し、（Ｎ1／２）≦ｐ1≦…≦ｐM≦Ｎ−（Ｎ1／２）であるものとする。すなわち、全ての解析位置ｐm，ｍ＝１，…，Ｍで定Ｑ変換の計算に必要なデータがｘ0，…，ｘN-1に含まれるものとする。尚、入力データの全点で定Ｑ変換を計算したいような場合は予め入力データを前後にゼロ拡張して上記条件を満たすように変換しておけばよい。
【０１３６】
また、周波数軸係数をゼロとみなすべき係数の大きさの閾値Ｔｈもここで決定しておく。
【０１３７】
時間軸係数算出部１２ではパラメータ決定部１１で決定したパラメータを基に、時間軸係数Ｔn,kを各解析位置ｐmについて次の定義式に従って計算する。
【０１３８】
【数１】
0006627639.tif 000017
【０１３９】
このときｐmに対するＴn,kをＴn,k,mと書くと、非ゼロ項に関してＴn,k,m＝Ｔn-pm+p1,k,1である。
【０１４０】
周波数軸係数算出部１３ではパラメータ決定部１１で決定したパラメータを基に、各解析位置ｐmについて周波数軸係数（１／Ｎ）Ｓn,k,mを定義式に従って計算する。この際、｜Ｓn,k,m｜≦Ｔｈとなる項は周波数軸係数（１／Ｎ）Ｓn,k,m＝０とする。
【０１４１】
計算法境界決定部１４では、周波数軸での計算を行う最大のｋ値ＫM（周波数成分閾値）を求める。ＫMは０，…，Ｋの何れかの値を持ち、ＫM＝０は周波数軸での計算を行わずに非特許文献１、２と同様に全て時間軸での計算で求める場合を意味する。
【０１４２】
周波数成分閾値ＫMの決定方法としては、ＫM＝０，…，Ｋの各場合について、定Ｑ変換の成分Ｘkcqを求めるのに必要な演算量ＣX（ＫM）を評価してそれが最小となるようにＫMを決定する。尚、この条件は複数のｐmに寄らないので、ｐ1（１つの解析位置）について求めておけば十分である。
【０１４３】
定Ｑ変換の成分Ｘkcqの演算量ＣX（ＫM）（合計演算量）は、長さＮのＦＦＴの演算量ＣFFT（Ｎ）、非ゼロの時間軸係数毎の演算量ＣT、非ゼロの時間軸係数の数ＮT（ｋ）、非ゼロの周波数軸係数毎の演算量ＣF、ｋで非ゼロとみなす周波数軸係数の数ＮF（ｋ）を使っておおよそ、前記実施例２で述べた数式（３）
【０１４４】
【数７】
0006627639.tif 000018
【０１４５】
と表せる。
【０１４６】
ここでＣFFTは計算を実行する計算機やそれぞれの算法の詳細により変わり得るが、実測からの一例としては、ＣFFT（Ｎ）＝（１／４）Ｎｌｏｇ2（Ｎ）、ＣT＝１、ＣF＝１としても良い。
【０１４７】
この演算量ＣX（ＫM）は、時間軸係数および周波数軸係数の、ハニング窓での項数と周波数（周波数成分ｋ；対数周波数）の関係が、非ゼロの時間軸係数の数ＮT（ｋ）がｋに関して指数的に減少し、非ゼロの周波数軸係数の数ＮF（ｋ）がｋに関しておおよそ指数的に増加するため、ＫM＝１，…，Ｋの範囲ではおおよそ下に凸の特性カーブとなる。ただし、ＫM＝０ではＣFFTがなくなるためＣX（０）が小さくなる。結果として演算量ＣX（ＫM）はＫM＝０あるいはＫM＝１，…，Ｋでの最小点の何れかで最小値を取ることになる。
【０１４８】
個別計算部２０では入力された時系列の定Ｑ変換を次のようにして求める。  
【０１４９】
前処理部２１では、入力された時系列データの長さがＮより短ければ前後にゼロを補完して、Ｎより長ければ前後の値を除去して長さがＮの時系列データｘnになるようにする。ゼロを補完する場合には、入力された時系列データの中央部分が時系列データｘnの中央部分に一致するようにする。データを除去する場合は、解析したい範囲に寄るが、単純に後ろのデータを除去するのでも良い。
【０１５０】
次に定Ｑ変換の成分Ｘkcqを、計算法境界決定部１４で決定したＫMにより、ｋ≦ＫMについては低周波帯域計算部２３で、ｋ＞ＫMについては高周波帯域計算部２２で計算する。
【０１５１】
【数８】
0006627639.tif 000019
【０１５２】
尚、前記定Ｑ変換の成分は以下、「Ｘk,mcq」と表記することもある。
【０１５３】
各ｋについて、前記時間軸係数算出部１２で算出された非ゼロの時間軸係数Ｔn,k,mは高々Ｎk個なので、積和計算を非ゼロの範囲に限定することで効率的に計算できる。
【０１５４】
低周波帯域計算部２３は、まず時系列データｘnを全体に対して一度ＦＦＴしてフーリエ係数Ｘnを求める。次にｐm毎に、前記周波数軸係数算出部１３で算出された周波数軸係数（１／Ｎ）Ｓn,k,mとフーリエ係数Ｘnを積和計算して定Ｑ変換の成分Ｘk,mcqを求める。
【０１５５】
各ｋについて非ゼロの周波数軸係数（１／Ｎ）Ｓn,k,mは少数なので、積和計算を非ゼロの範囲に限定することで効率的に計算できる。
【０１５６】
尚、周波数軸係数は多くの場合、両端（周波数ゼロ近辺とサンプリング周波数近辺）部分では非ゼロ項とゼロ項が入り混じるが中間部分ではゼロ項が続くので、非ゼロ項の範囲を少し広げて、両端部分は全て非ゼロ項として計算しても良い。この場合、非ゼロ項が少し増えるがその数はそれほど多くはなく、また非ゼロ判定ロジックが簡略化できるので演算時間は大差なくなる。
【０１５７】
最後に定Ｑ変換合成部２４で、高周波帯域計算部２２と低周波帯域計算部２３とで計算した定Ｑ変換の成分Ｘk,mcqをまとめて、定Ｑ変換を完了する。
【実施例４】
【０１５８】
前記実施例３では、時間軸係数Ｔn,k,m、周波数軸係数（１／Ｎ）Ｓn,k,mは解析位置ｐm毎に用意するので、事前に用意して保管しておくデータ量が大きくなる。本実施例４では、係数データのメモリ使用量を大幅に減らすように以下のような対策を講じた。
【０１５９】
【数９】
0006627639.tif 000020
【０１６０】
尚、数式（４）中の
【０１６１】
【数３】
0006627639.tif 000021
【０１６２】
は、複素指数関数の時刻および周波数毎に異なるべき乗を表している。
【０１６３】
本実施例４による定Ｑ変換の成分演算装置の具体的な構成は図５と同様であり、以下に図５と異なる部分のみを説明する。
【０１６４】
時間軸係数算出部１２では、パラメータ決定部１１で決定したパラメータを基に、時間軸係数Ｔn,k,1を定義式に従って求める。すなわち、
【０１６５】
【数１】
0006627639.tif 000022
【０１６６】
を、予め設定した解析位置ｐmのうち１つの解析位置（例えば１番目の解析位置）ｐ1において算出し、時間軸係数Ｔn,k,1を求める。
【０１６７】
周波数軸係数算出部１３では、パラメータ決定部１１で決定したパラメータを基に、周波数軸係数（１／Ｎ）Ｓn,k,1を定義式に従って計算する。
【０１６８】
すなわち、
【０１６９】
【数２】
0006627639.tif 000023
【０１７０】
から得た（１／Ｎ）Ｓn,kを、前記１つの解析位置ｐ1において算出して周波数軸係数（１／Ｎ）Ｓn,k,1を求める。
【０１７１】
この際、｜Ｓn,k,1｜≦Ｔｈとなる項は周波数軸係数（１／Ｎ）Ｓn,k,1＝０とする。
【０１７２】
計算法境界決定部１４では、周波数軸計算の際に複素べき乗の積和が追加されるために非ゼロの周波数軸係数毎の演算量ＣFを実施例３の場合より大きく、例えば５（複素べき乗の演算に＋３、複素べき乗の積和に＋１）とする。
【０１７３】
高周波帯域計算部２２は、複数の解析位置ｐm毎に範囲をずらして、前記時間軸係数算出部１２で算出された時間軸係数Ｔn,k,1と時系列データｘn-p1+pmを積和計算して定Ｑ変換の成分Ｘk,mcqを求める。
【０１７４】
低周波帯域計算部２３は、まず時系列データｘnを全体に対して一度ＦＦＴしてフーリエ係数Ｘnを求める。
【０１７５】
【数１０】
0006627639.tif 000024
【０１７６】
前記実施例３、４の手法によれば、定Ｑ変換を高速に実行するにあたって、先行の高速化技術（実施例２）が時間軸での計算でも周波数軸での計算でもそれらを併用する計算方法でも、一時刻毎に定Ｑ変換を計算するのに対して、複数時刻での定Ｑ変換を同時に計算する際に、ＤＦＴ（あるいはＦＦＴ）を一時刻の計算範囲毎ではなく全体に対して一度だけ実施することで、全体としての計算が高速化される。
【０１７７】
実施例４では、時間軸と周波数軸との係数データを計算時刻毎に持つ代わりに、それら複素数のべき乗計算に置き換えて（数式（４））、演算量は少し増える代わりに係数データのメモリ使用量を大幅に減らすことができる。
【符号の説明】
【０１７８】
１０…準備部
１１…パラメータ決定部
１２…時間軸係数算出部
１３…周波数軸係数算出部
１４…計算法境界決定部
２０…個別計算部
２１…前処理部
２２…高周波帯域計算部
２３…低周波帯域計算部
２４…定Ｑ変換合成部
１００…診断用パラメータ作成手段
１０１…第１の定Ｑ変換装置
１０２…主成分分析部
１０３…第１の統計値計算部
１０４…第１の正規化部
２００…診断手段
２０１…第２の定Ｑ変換装置
２０２…主成分計算部
２０３…第２の統計値計算部
２０４…第２の正規化部
【図１】
0006627639.tif 000032
【図２】
0006627639.tif 000033
【図３】
0006627639.tif 000034
【図４】
0006627639.tif 000035
【図５】
0006627639.tif 000036
【図６】
0006627639.tif 000037
【図７】
0006627639.tif 000038