(57)【要約】【課題】時系列データの非ガウス性又は非線形性の検出に当たり、観測データの大きさ、歪度に依存しない検定量を得る。【解決手段】時系列データ中の複数個のデータブロックをフーリエ変換して、各データブロックについて複数の周波数成分を計算し（１４）、計算された周波数成分のうちの２つの周波数（ｋ１、ｋ２）からなる周波数ペアの各々のバイコヒーレンス（γ［ｋ１，ｋ２］）を計算し（２０）、周波数ペアの各々のバイコヒーレンス（γ［ｋ１，ｋ２］）に基づいて、検定量（ｚγ）を計算する（３２）。周波数ペアの各々のバイスペクトルの正規化に当たり、周波数ペアを構成する周波数の周波数成分（Ｘｉ［ｋ１］、Ｘｉ［ｋ２］の積の絶対値の自乗の、複数のデータブロックについての平均（Ｃａ［ｋ１，ｋ２］）を用いる。【選択図】図５
2015046116.tif 000002
【特許請求の範囲】
【請求項１】
観測値を表す時系列データ中の、各々相連続するサンプル値で構成される複数個のデータブロックをフーリエ変換して、各データブロックについて複数の周波数成分を計算するフーリエ変換部と、
前記フーリエ変換部で計算された各データブロックについての周波数成分のうちの２つの周波数からなる周波数ペアの各々のバイコヒーレンスを計算するバイコヒーレンス計算部と、
前記周波数ペアの各々についての前記バイコヒーレンスに基づいて、前記時系列データについての検定量を計算する検定量計算部と、
前記検定量計算部で計算された前記検定量が予め定められた閾値よりも大きいか否かを判断する比較部とを有し、
前記比較部による判断結果が、解析の結果として出力される
データ解析装置。
【請求項２】
前記バイコヒーレンス計算部は、
前記各データブロックについての前記周波数ペアの各々のバイスペクトルの、前記複数のデータブロックについての平均を計算し、前記時系列データについてのバイスペクトルとして出力するバイスペクトル計算部と、
前記各データブロックについての前記周波数ペアの各々を構成する周波数の和に等しい周波数の成分のパワースペクトルの、前記複数のデータブロックについての平均を計算して、前記時系列データについてのパワースペクトルとして出力するパワースペクトル計算部と、
前記各データブロックについての前記周波数ペアの各々を構成する周波数の成分の積の絶対値の自乗の、前記複数のデータブロックについての平均を計算し、前記時系列データについての積自乗平均として出力する積自乗平均計算部と、
前記バイスペクトル計算部から出力された前記時系列データについての前記バイスペクトル、前記パワースペクトル計算部から出力された前記時系列データについての前記パワースペクトル、及び前記積自乗平均計算部から出力された前記時系列データについての前記積自乗平均に基づいて、前記バイコヒーレンスを計算する規格化部とを有する
ことを特徴とする請求項１に記載のデータ解析装置。
【請求項３】
前記バイスペクトル計算部は、前記各データブロックについての前記バイスペクトルをＢｉ［ｋ１，ｋ２］で表すとき、
【数２２】
2015046116.tif 000024
（但し、
ｋ１、ｋ２はそれぞれ前記周波数ペアを構成する周波数、
Ｘｉ[ｋ１]、Ｘｉ[ｋ２]、Ｘｉ[ｋ１＋ｋ２]はそれぞれ、ｉ番目のデータブロック（ｉは１乃至Ｐのいずれか）の、周波数ｋ１、ｋ２、ｋ１＋ｋ２の周波数成分）
により、前記各データブロックについての前記バイスペクトルＢｉ［ｋ１，ｋ２］を計算し、
前記時系列データについての前記バイスペクトルをＢａ［ｋ１，ｋ２］で表すとき、
【数２３】
2015046116.tif 000025
（但し、Ｐは前記時系列データを構成するデータブロックの数）
により、前記時系列データについての前記バイスペクトルＢａ［ｋ１，ｋ２］を計算する
ことを特徴とする請求項２に記載のデータ解析装置。
【請求項４】
前記積自乗平均計算部は、前記時系列データについての積自乗平均をＣａ［ｋ１，ｋ２］で表すとき、
【数２４】
2015046116.tif 000026
（但し、
ｋ１、ｋ２はそれぞれ前記周波数ペアを構成する周波数、
Ｘｉ[ｋ１]、Ｘｉ[ｋ２]はそれぞれ、ｉ番目のデータブロック（ｉは１乃至Ｐのいずれか）の、周波数ｋ１、ｋ２の周波数成分、
Ｐは前記時系列データを構成するデータブロックの数）
により、前記時系列データについての前記積自乗平均Ｃａ［ｋ１，ｋ２］を計算する
ことを特徴とする請求項２又は３に記載のデータ解析装置。
【請求項５】
前記パワースペクトル計算部は、前記時系列データについてのパワースペクトルをＳａ［ｋ１＋ｋ２］で表すとき、
【数２５】
2015046116.tif 000027
（但し、
ｋ１、ｋ２はそれぞれ前記周波数ペアを構成する周波数、
Ｘｉ［ｋ１＋ｋ２］はｉ番目のデータブロック（ｉは１乃至Ｐのいずれか）の、周波数ｋ１＋ｋ２の周波数成分であり、
Ｘｉ＊［ｋ１＋ｋ２］は、Ｘｉ［ｋ１＋ｋ２］の複素共役であり、
Ｐはデータブロックの数である。）
により、前記時系列データについての前記パワースペクトルＳａ［ｋ１＋ｋ２］を計算する
ことを特徴とする請求項２、３又は４に記載のデータ解析装置。
【請求項６】
前記規格化部は、前記バイコヒーレンスをγ［ｋ１，ｋ２］で表すとき、
【数２６】
2015046116.tif 000028
により、前記バイコヒーレンスγ［ｋ１，ｋ２］を計算する
ことを特徴とする請求項２乃至５のいずれかに記載のデータ解析装置。
【請求項７】
前記検定量計算部は、前記検定量の平均値及び分散が、各データブロックを構成するサンプル値の数、前記時系列データを構成するデータブロックの数、及び
前記データブロックのオーバーラップ率によらずに一定に保たれるようにするための調整量を用いて、前記検定量の計算を行うことを特徴とする請求項１に記載のデータ解析装置。
【請求項８】
前記検定量計算部は、前記検定量をｚγで表すとき、
【数２７】
2015046116.tif 000029
（但し、
Ｑは前記周波数ペアの数、
γｑは、前記周波数ペアの各々について、前記バイコヒーレンス計算部で計算されたバイコヒーレンス、
ｂｉａｓは、前記検定量の平均を０とするための調整値、
σ２は、前記検定量の分散を補正するための調整値）
により前記検定量ｚγを計算することを特徴とする請求項１に記載のデータ解析装置。
【請求項９】
前記調整値ｂｉａｓ及びσ２を、
【数２８】
2015046116.tif 000030
（但し、
Ｐは、前記時系列データを構成するデータブロックの数、
Ｎは各データブロックを構成するサンプル値の数、
α及びβはパラメータ）
により計算する調整値計算部をさらに有する
ことを特徴とする請求項８に記載のデータ解析装置。
【請求項１０】
白色雑音を表す時系列データを、前記観測値を表す時系列データの代わりに用いて、
前記フーリエ変換及び前記バイコヒーレンスの計算を行い、
該計算により求められたバイコヒーレンスに基づいて試験的検定量を計算し、
複数回の試行についての前記試験的検定量の平均ｚｗａ及び分散ｚｗｖを計算し、
該計算により求められた平均ｚｗａ及び分散ｚｗｖに基づいて、
【数２９】
2015046116.tif 000031
により、前記パラメータα及びβを生成するパラメータ生成装置と、
前記パラメータ生成装置で生成された前記パラメータα及びβを記憶するパラメータメモリと
をさらに有することを特徴とする請求項９に記載のデータ解析装置。
【請求項１１】
観測値を表す時系列データ中の、各々相連続するサンプル値で構成される複数個のデータブロックをフーリエ変換して、各データブロックについて複数の周波数成分を計算するフーリエ変換ステップと、
前記フーリエ変換ステップで計算された各データブロックについての周波数成分のうちの２つの周波数からなる周波数ペアの各々のバイコヒーレンスを計算するバイコヒーレンス計算ステップと、
前記周波数ペアの各々についての前記バイコヒーレンスに基づいて、前記時系列データについての検定量を計算する検定量計算ステップと、
前記検定量計算ステップで計算された前記検定量が予め定められた閾値よりも大きいか否かを判断する比較ステップとを有し、
前記比較ステップによる判断結果が、解析の結果として出力される
データ解析方法。
【請求項１２】
前記バイコヒーレンス計算ステップは、
前記各データブロックについての前記周波数ペアの各々のバイスペクトルの、前記複数のデータブロックについての平均を計算し、前記時系列データについてのバイスペクトルとして出力するバイスペクトル計算ステップと、
前記各データブロックについての前記周波数ペアの各々を構成する周波数の和に等しい周波数の成分のパワースペクトルの、前記複数のデータブロックについての平均を計算して、前記時系列データについてのパワースペクトルとして出力するパワースペクトル計算ステップと、
前記各データブロックについての前記周波数ペアの各々を構成する周波数の成分の積の絶対値の自乗の、前記複数のデータブロックについての平均を計算し、前記時系列データについての積自乗平均として出力する積自乗平均計算ステップと、
前記バイスペクトル計算ステップから出力された前記時系列データについての前記バイスペクトル、前記パワースペクトル計算ステップから出力された前記時系列データについての前記パワースペクトル、及び前記積自乗平均計算ステップから出力された前記時系列データについての前記積自乗平均に基づいて、前記バイコヒーレンスを計算する規格化ステップとを有する
ことを特徴とする請求項１１に記載のデータ解析方法。
【請求項１３】
前記バイスペクトル計算ステップは、前記各データブロックについての前記バイスペクトルをＢｉ［ｋ１，ｋ２］で表すとき、
【数３０】
2015046116.tif 000032
（但し、
ｋ１、ｋ２はそれぞれ前記周波数ペアを構成する周波数、
Ｘｉ[ｋ１]、Ｘｉ[ｋ２]、Ｘｉ[ｋ１＋ｋ２]はそれぞれ、ｉ番目のデータブロック（ｉは１乃至Ｐのいずれか）の、周波数ｋ１、ｋ２、ｋ１＋ｋ２の周波数成分）
により、前記各データブロックについての前記バイスペクトルＢｉ［ｋ１，ｋ２］を計算し、
前記時系列データについての前記バイスペクトルをＢａ［ｋ１，ｋ２］で表すとき、
【数３１】
2015046116.tif 000033
（但し、Ｐは前記時系列データを構成するデータブロックの数）
により、前記時系列データについての前記バイスペクトルＢａ［ｋ１，ｋ２］を計算する
ことを特徴とする請求項１２に記載のデータ解析方法。
【請求項１４】
前記積自乗平均計算ステップは、前記時系列データについての積自乗平均をＣａ［ｋ１，ｋ２］で表すとき、
【数３２】
2015046116.tif 000034
（但し、
ｋ１、ｋ２はそれぞれ前記周波数ペアを構成する周波数、
Ｘｉ[ｋ１]、Ｘｉ[ｋ２]はそれぞれ、ｉ番目のデータブロック（ｉは１乃至Ｐのいずれか）の、周波数ｋ１、ｋ２の周波数成分、
Ｐは前記時系列データを構成するデータブロックの数）
により、前記時系列データについての前記積自乗平均Ｃａ［ｋ１，ｋ２］を計算する
ことを特徴とする請求項１２又は１３に記載のデータ解析方法。
【請求項１５】
前記パワースペクトル計算ステップは、前記時系列データについてのパワースペクトルをＳａ［ｋ１＋ｋ２］で表すとき、
【数３３】
2015046116.tif 000035
（但し、
ｋ１、ｋ２はそれぞれ前記周波数ペアを構成する周波数、
Ｘｉ［ｋ１＋ｋ２］はｉ番目のデータブロック（ｉは１乃至Ｐのいずれか）の、周波数ｋ１＋ｋ２の周波数成分であり、
Ｘｉ＊［ｋ１＋ｋ２］は、Ｘｉ［ｋ１＋ｋ２］の複素共役であり、
Ｐはデータブロックの数である。）
により、前記時系列データについての前記パワースペクトルＳａ［ｋ１＋ｋ２］を計算する
ことを特徴とする請求項１２、１３又は１４に記載のデータ解析方法。
【請求項１６】
前記規格化ステップは、前記バイコヒーレンスをγ［ｋ１，ｋ２］で表すとき、
【数３４】
2015046116.tif 000036
により、前記バイコヒーレンスγ［ｋ１，ｋ２］を計算する
ことを特徴とする請求項１２乃至１５のいずれかに記載のデータ解析方法。
【請求項１７】
前記検定量計算ステップは、前記検定量の平均値及び分散が、各データブロックを構成するサンプル値の数、前記時系列データを構成するデータブロックの数、及び
前記データブロックのオーバーラップ率によらずに一定に保たれるようにするための調整量を用いて、前記検定量の計算を行うことを特徴とする請求項１１に記載のデータ解析方法。
【請求項１８】
前記検定量計算ステップは、前記検定量をｚγで表すとき、
【数３５】
2015046116.tif 000037
（但し、
Ｑは前記周波数ペアの数、
γｑは、前記周波数ペアの各々について、前記バイコヒーレンス計算ステップで計算されたバイコヒーレンス、
ｂｉａｓは、前記検定量の平均を０とするための調整値、
σ２は、前記検定量の分散を補正するための調整値）
により前記検定量ｚγを計算することを特徴とする請求項１１に記載のデータ解析方法。
【請求項１９】
前記調整値ｂｉａｓ及びσ２を、
【数３６】
2015046116.tif 000038
（但し、
Ｐは、前記時系列データを構成するデータブロックの数、
Ｎは各データブロックを構成するサンプル値の数、
α及びβはパラメータ）
により計算する調整値計算ステップをさらに有する
ことを特徴とする請求項１８に記載のデータ解析方法。
【請求項２０】
白色雑音を表す時系列データを、前記観測値を表す時系列データの代わりに用いて、
前記フーリエ変換及び前記バイコヒーレンスの計算を行い、
該計算により求められたバイコヒーレンスに基づいて試験的検定量を計算し、
複数回の試行についての前記試験的検定量の平均ｚｗａ及び分散ｚｗｖを計算し、
該計算により求められた平均ｚｗａ及び分散ｚｗｖに基づいて、
【数３７】
2015046116.tif 000039
により、前記パラメータα及びβを生成するパラメータ生成ステップと、
前記パラメータ生成ステップで生成された前記パラメータα及びβをパラメータメモリに記憶させる記憶ステップと
をさらに有することを特徴とする請求項１９に記載のデータ解析方法。
【請求項２１】
請求項１１乃至２０のいずれかに記載のデータ解析方法の各処理をコンピュータに実行させるためのプログラム。
【請求項２２】
請求項２１に記載のプログラムを記録した、コンピュータで読み取り可能な記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、データ解析装置及び方法、並びにプログラム及び記録媒体に関し、特に時系列データの非ガウス性又は非線形性の検出に関する。
【背景技術】
【０００２】
時系列データの非ガウス性とは、時系列データの振幅分布が、平均値を中心として左右対称でないことを意味する。一方時系列データの非線形性とは、時系列データの生成過程に、非線形項が含まれていることを意味する。非線形性のデータは、非ガウス性を有する。
取得したデータ（観測データ）の非ガウス性及び非線形性を確認するためのデータ解析方法として、データの歪度に関係する規格化バイスペクトルを用いることが提案されている（非特許文献１〜３）。
【０００３】
従来の方法の問題点は、観測データがガウス性であっても有限のデータを抽出すると、抽出の仕方によってデータが非ガウス性又は非線形性を有するように見え、誤判定が生じる可能性があることである。
【０００４】
以下、従来の方法について図１を参照して説明する。図１において、観測の結果得られた時系列データをｘ（ｔ）と表し、このｘ（ｔ）の離散化されたものをｘ［ｎ］と表す。ここで、ｎは時刻を表すインデックスである。
時系列データｘ［ｎ］から図２に示すように、複数個即ちＰ個のデータブロックＤＢ１〜ＤＢＰを生成する。データブロックの各々ＤＢｉ（ｉは１乃至Ｐのいずれか）はＮ個のサンプル値を含む。
各データブロックＤＢｉに対して、フーリエ変換部１４でＮサンプル離散フーリエ変換を行って、各周波数の成分Ｘｉ［ｋ］を求める。フーリエ変換は、下記の式（１）で表される。
【０００５】
【数１】
2015046116.tif 000003
ここで、ｋは周波数を表すインデックスである。
【０００６】
バイスペクトル計算部２２及びパワースペクトル計算部２４には、複数個即ちＱ個の周波数ペアを構成する周波数についての周波数成分Ｘｉ［ｋ］が供給される。例えば、Ｑは４９であり、４９個の周波数ペアの２次元周波数平面上における配置は図３に「○」印で示される通りである。
【０００７】
バイスペクトル計算部２２では、データブロックＤＢ１〜ＤＢＰの各々について周波数成分Ｘｉ［ｋ］を用いて、バイスペクトルＢａ［ｋ１，ｋ２］を計算する。
そのために、まず各データブロック（ｉ番目のデータブロック）ＤＢｉについてのバイスペクトルＢｉ［ｋ１，ｋ２］を、次式（２）によって計算する。
【数２】
2015046116.tif 000004
但し、Ｘｉ＊［ｋ１＋ｋ２］はＸｉ［ｋ１＋ｋ２］の複素共役を表す。
【０００８】
式（２）で得られるバイスペクトルＢｉ［ｋ１，ｋ２］は、観測データｘ［ｎ］の周波数間の相関の程度を示す指標である。
【０００９】
バイスペクトル計算部２２は、式（２）の計算を、Ｐ個のデータブロックＤＢ１〜ＤＢＰの各々に対して行い、計算の結果得られたバイスペクトルＢｉ［ｋ１，ｋ２］を平均することによって、Ｐ個のデータブロックから成る時系列データＤＳについてのバイスペクトルＢａ［ｋ１，ｋ２］を得る。バイスペクトルＢａ［ｋ１，ｋ２］を求めるための計算は、次式（３）で表される。
【００１０】
【数３】
2015046116.tif 000005
【００１１】
観測データに白色雑音しか含まれていなければ、観測データの周波数間に相関関係はないので、バイスペクトルは小さな値となる。一方、観測データに非定常成分が含まれていれば、観測データの周波数間に相関関係が現れ、式（３）で示されるバイスペクトルが大きな値を持つ。
【００１２】
式（３）で示されたバイスペクトルは、観測データに含まれる各周波数の成分の大きさ（パワー）に依存した値を持つ。そのため、このバイスペクトルそのものでは、周波数間の相関関係の大きさが判断しにくい。そこで通常は、バイスペクトルの大きさが観測データの各周波数の成分の大きさに依存しない、規格化した値が用いられる。
【００１３】
規格化のため、パワースペクトル計算部２４では、観測データのパワースペクトルＳａ［ｋ］を、次のように計算する。
【数４】
2015046116.tif 000006
ここで、
Ｘｉ＊［ｋ］は、Ｘｉ［ｋ］の複素共役である。
式（４）の計算は、周波数ｋ１、ｋ２、ｋ３（＝ｋ１＋ｋ２）の各々について行われ、従って、パワースペクトル計算部２４からは、Ｓａ［ｋ１］、Ｓａ［ｋ２］、Ｓａ［ｋ１＋ｋ２］が出力される。
【００１４】
規格化部２７は、バイスペクトル計算部２２で計算されたバイスペクトルＢａ［ｋ１，ｋ２］と、パワースペクトル計算部２４で計算されたパワースペクトルＳａ［ｋ１］、Ｓａ［ｋ２］及びＳａ［ｋ１＋ｋ２］を用いて、規格化バイスペクトルＢｎ［ｋ１，ｋ２］を、次式（５）により、計算する。
【数５】
2015046116.tif 000007
【００１５】
分母にパワースペクトルの項があることによって、観測データの各周波数の成分の大きさが規格化される。従って、式（５）の規格化バイスペクトルＢｎ［ｋ１，ｋ２］は観測データの大きさ（パワー）に依存せず、周波数間の相関に応じた値を示すものとなる。
但し、実際の計算では、式（５）のバイスペクトルＢｎ［ｋ１，ｋ２］の分散を１にするために、次式（６）のように補正することで得られる、補正された規格化バイスペクトルＢｎｃ［ｋ１，ｋ２］を用いる。
【数６】
2015046116.tif 000008
【００１６】
式（６）で表される規格化バイスペクトルＢｎｃ［ｋ１，ｋ２］は、特定の周波数ペア（ｋ１，ｋ２）について計算された値であり、Ｑ個の周波数ペアについて上記した処理が繰り返されて、Ｑ個の周波数ペアの各々について式（６）で表される規格化バイスペクトルＢｎｃ［ｋ１，ｋ２］が求められる。
【００１７】
検定量計算部２９は、上記のようにしてＱ個の周波数ペアについて求められた規格化バイスペクトルＢｎｃ［ｋ１，ｋ２］を用いて、観測データの非ガウス性及び非線形性を検出する。
【００１８】
観測データが白色雑音であれば、それぞれＱ個の周波数ペア（ｋ１，ｋ２）について求められた規格化バイスペクトルＢｎｃ［ｋ１，ｋ２］は、平均０、分散１の正規分布に従う。この性質を利用して、観測データの非定常性を検出するために、次に示す検定量ｚを計算する。
【００１９】
【数７】
2015046116.tif 000009
【００２０】
式（７）でＢｎｃｑ（但しｑは１乃至Ｑのいずれか）は、Ｑ個の周波数ペア（ｋ１，ｋ２）の各々についてのＢｎｃ［ｋ１，ｋ２］と同じものを表す。
式（７）の分子はＱ個の規格化バイスペクトルの値を加算することを表している。式（７）の分母は、検定量ｚの分散を規格化するための項である。
観測データに白色雑音のみが含まれており、Ｑが十分に大きい値であれば、検定量ｚは平均０、分散１の正規分布に従う。もし、観測データに非ガウス性又は非線形性の成分が含まれていれば、規格化バイスペクトルの値が大きくなり、その結果、検定量ｚは平均０、分散１の正規分布には従わずに大きな値となる。このことを利用して、観測データの非ガウス性又は非線形性が検出できる。
以上が規格化バイスペクトルを用いて観測データから非ガウス性及又は非線形性を検出するための従来の方法である。
【先行技術文献】
【非特許文献】
【００２１】
【非特許文献１】Ｍｅｌｖｉｎ  Ｊ．  Ｈｉｎｉｃｈ，  “Ｂｉｓｐｅｃｔｒｕｍ  ｏｆ  ｓｈｉｐ−ｒａｄｄｉａｔｅｄ  ｎｏｉｓｅ”，  Ｊ．  Ａｃｏｕｓｔ．  Ｓｏｃ．  Ａｍ．  ８５（４），  Ａｐｒｉｌ  １９８９
【非特許文献２】Ｍｅｌｖｉｎ  Ｊ．  Ｈｉｎｉｃｈ，  “Ｄｅｔｅｃｔｉｎｇ  ａ  Ｔｒａｎｓｉｅｎｔ  Ｓｉｇｎａｌ  ｂｙ  Ｂｉｓｐｅｃｔｒａｌ  Ａｎａｌｙｓｉｓ”，  ＩＥＥＥ  Ｔｒａｎｓ．  Ｏｎ  Ａｃｏｕｓｔｉｃｓ，  Ｓｐｅｅｃｈ  ａｎｄ  Ｓｉｇｎａｌ  Ｐｒｏｃｅｓｓｉｎｇ，  Ｖｏｌ．３６，  Ｎｏ．７，  Ｊｕｌｙ  １９９０
【非特許文献３】Ｍｅｌｖｉｎ  Ｊ．  Ｈｉｎｉｃｈ  ｅｔ．  ａｌ，  “Ｏｎ  ｔｈｅ  Ｐｒｉｎｃｉｐａｌ  Ｄｏｍａｉｎ  ｏｆ  ｔｈｅ  Ｄｉｓｃｒｅｔｅ  Ｂｉｓｐｅｃｔｒｕｍ  ｏｆ    ａ  Ｓｔａｔｉｏｎａｒｙ  Ｓｉｇｎａｌ”，  ＩＥＥＥ  Ｔｒａｎｓ．  Ｏｎ  Ｓｉｇｎａｌ  Ｐｒｏｃｅｓｓｉｎｇ，  Ｖｏｌ．４３，  Ｎｏ．９，  Ｓｅｐｔｅｍｂｅｒ  １９９５
【発明の概要】
【発明が解決しようとする課題】
【００２２】
上記した従来の方法は、観測データが非ガウス性又は非線形性を有すれば、観測データがある程度の歪度を有しているということに着目した処理である。しかしながら、観測データにガウス性又は線形性があっても、有限のデータを抽出すると、抽出の仕方によっては非ガウス性又は非線形性を有するように見え、そのため判定に誤りが生じる可能性がある。
【００２３】
本発明は、上記の課題の解決のためのものであり、その目的は、観測データの大きさ、歪度の影響が抑制された処理結果を得ることを可能にすることにある。
【課題を解決するための手段】
【００２４】
本発明のデータ解析装置は、
観測値を表す時系列データ中の、各々相連続するサンプル値で構成される複数個のデータブロックをフーリエ変換して、各データブロックについて複数の周波数成分を計算するフーリエ変換部と、
前記フーリエ変換部で計算された各データブロックについての周波数成分のうちの２つの周波数からなる周波数ペアの各々のバイコヒーレンスを計算するバイコヒーレンス計算部と、
前記周波数ペアの各々についての前記バイコヒーレンスに基づいて、前記時系列データについての検定量を計算する検定量計算部と、
前記検定量計算部で計算された前記検定量が予め定められた閾値よりも大きいか否かを判断する比較部とを有し、
前記比較部による判断結果が、解析の結果として出力されるものである。
【発明の効果】
【００２５】
本発明によれば、観測データの周波数間の位相的な相関関係を表す指標であるバイコヒーレンスを用い、バイコヒーレンスに基づいた検定量を用いて、観測データの非ガウス性及び非線形性を検出するので、観測データの大きさ、歪度の影響を抑制しつつ、非ガウス性又は非線形性の判定を正確に行うことができる。
【図面の簡単な説明】
【００２６】
【図１】従来のデータ解析方法を実施するための装置の概略を示すブロック図である。
【図２】時系列データからのデータブロックの生成方法を示す図である。
【図３】データ解析に用いられる周波数ペアの２次元周波数平面上での配置の一例を示す図である。
【図４】本発明のデータ解析装置の概略を示すブロック図である。
【図５】図４のデータ解析装置の検出処理部の構成例を示すブロック図である。
【図６】図５の検出処理部のコヒーレンス計算部の構成例を示すブロック図である。
【図７】図５の検出処理部で用いられるパラメータを生成する装置の構成例を示すブロック図である。
【図８】（ａ）及び（ｂ）は、模擬した過渡的信号、及びその模擬信号に白色雑音を付加した波形を示す図である。
【図９】（ａ）及び（ｂ）は、図８（ｂ）に示される信号に対して、従来の方法及び本発明の方法を適用した結果を示す図である。
【発明を実施するための形態】
【００２７】
図４は、本発明の時系列データ解析装置を示す。以下では、一例として、可聴域の音響信号を受信して、解析を行う場合について説明する。図示の時系列データ解析装置は、入力信号処理部２と、時系列データメモリ４と、検出処理部６とを有する。
【００２８】
入力信号処理部２は、Ａ／Ｄ変換器３を有し、音響信号を電気信号に変換するマイクロフォン１に接続され、該マイクロフォン１から与えられるアナログ信号ｘ（ｔ）をサンプリングして、ディジタル化して、時系列データ（観測値を表す時系列データ）ｘ［ｎ］を生成する。Ａ／Ｄ変換器３のサンプリング周波数は、例えば５０ｋＨｚであり、Ａ／Ｄ変換器３から出力される時系列データはｘ［ｎ］で表される。ｎは、時刻を示すインデックスである。
【００２９】
時系列データメモリ４は、入力信号処理部２から出力される時系列データｘ［ｎ］を記憶する。
【００３０】
検出処理部６は、図５に示すように、データブロック生成部１２と、フーリエ変換部１４と、周波数成分メモリ１６と、周波数成分選択部１８と、バイコヒーレンス計算部２０と、検定量計算部３２と、比較部３４と、パラメータメモリ３６と、調整値計算部３８とを有する。図１と同じ符号は同様のブロックを示す。データブロック生成部１２、周波数成分メモリ１６及び周波数成分選択部１８は図１には示されていないが、従来の方法でも実際には同様のものが用いられる。
【００３１】
データブロック生成部１２は、時系列データメモリ４に記憶された時系列データｘ［ｎ］から、各々複数のサンプル値で構成される、複数個即ちＰ個のデータブロックＤＢｉ（ｉ＝１〜Ｐ）を生成して、フーリエ変換部１４に供給する。
データブロックＤＢｉは、例えば、図３に示されるように、各々相連続するＮ個のサンプル値から成り、オーバーラップ率Ｌが５０％である。即ち、各データブロック（第ｉ番目のデータブロック）の後半のＮ／２個のサンプル値と、次のデータブロック（第ｉ＋１番目のデータブロック）の前半のＮ／２個のサンプル値とは同じものである。
【００３２】
フーリエ変換部１４は、Ｐ個のデータブロックＤＢ１〜ＤＢＰの各々のＮ個のサンプル値（観測データ）ｘ［０］〜ｘ［Ｎ−１］に対して離散フーリエ変換を行い、変換の結果として得られる周波数成分Ｘ１［ｋ］〜ＸＰ［ｋ］を出力する。ｋは周波数を表すインデックスである。各データブロック（ｉ番目のデータブロック）ＤＢｉ（ｉは、１乃至Ｐのいずれか）についてのフーリエ変換の結果（周波数ｋの成分）をＸｉ［ｋ］で表す。
【００３３】
フーリエ変換の結果得られる周波数成分Ｘｉ［ｋ］は、下記の式（１）で表される。式（１）は従来の方法に関して先に記載した式（１）と同じであるが、再度記載する。以下他の数式についても同様である。
【００３４】
【数８】
2015046116.tif 000010
【００３５】
ここで、ｋはＭ個の周波数の値ｆ１、ｆ２、…ｆＭのいずれかを表すインデックスである。ｆ１〜ｆＭを一般化してｆｍ（ｍは１〜Ｍ）で表す。
周波数ｆｍは、以下の値を有する。
ｍ＝１のとき、ｆｍ＝０                    （８Ａ）
２≦ｍ＜Ｍのとき、ｆｍ＝ｆ（ｍ−１）＋ΔＦ    （８Ｂ）
但し、ΔＦ＝１５６２．５Ｈｚであり、これは５０ｋＨｚ（サンプリング周波数）をＮ＝３２（各データブロック中のサンプル値の数）で割った値に等しい。
Ｍは例えば２１である。
【００３６】
フーリエ変換部１４で算出された、複数のデータブロックＤＢｉ（ｉ＝１〜Ｐ）の各々についての周波数成分Ｘｉ［ｋ］（ｋ＝ｆ１〜ｆＭ）は、周波数成分メモリ１６に記憶される。
【００３７】
周波数成分選択部１８は、周波数成分メモリ１６に記憶された複数のデータブロックＤＢｉ（ｉ＝１〜Ｐ）の各々についての周波数成分Ｘｉ［ｋ］（ｉ＝１〜Ｐ、ｋ＝ｆ１〜ｆＭ）のうち、バイコヒーレンスの計算に用いられるものを逐次読み出して、バイコヒーレンス計算部２０に供給する。
【００３８】
バイコヒーレンス計算部２０におけるバイコヒーレンスの計算は、複数個即ちＱ個の周波数ペアについて順に行われる。Ｑは例えば４９であり、４９個の周波数ペアの２次元周波数平面上における配置は例えば図３に「○」印で示される如くである。
【００３９】
周波数ペア（ｋ１，ｋ２）についてバイコヒーレンスを計算する際、複数個即ちＰ個のデータブロックの各々ＤＢｉについての、周波数ｋ１、ｋ２の周波数成分Ｘｉ［ｋ１］、Ｘｉ［ｋ２］のみならず、ｋ１とｋ２の和に等しい周波数ｋ３の周波数成分Ｘｉ［ｋ３］も用いられる。そこで周波数成分選択部１８は、各データブロックＤＢｉについての周波数ｋ１、ｋ２、ｋ３の周波数成分Ｘｉ［ｋ１］、Ｘｉ［ｋ２］、Ｘｉ［ｋ３］をバイコヒーレンス計算部２０に供給する。
上記のように、
ｋ３＝ｋ１＋ｋ２
の関係があるので、以下では、「ｋ３」を、「ｋ１＋ｋ２」と表記する。
一方、任意の周波数についての説明、或いはｋ１、ｋ２、ｋ３のすべてに当てはまる説明の際は、符号「ｋ」を用いる。
【００４０】
バイコヒーレンス計算部２０は、各周波数ペア（ｋ１，ｋ２）について周波数成分選択部１８から供給された周波数成分Ｘｉ［ｋ１］、Ｘｉ［ｋ２］、Ｘｉ［ｋ１＋ｋ２］を用いて、当該周波数ペア（ｋ１，ｋ２）のバイコヒーレンスγ［ｋ１，ｋ２］を計算する。以下に詳しく述べるように、Ｐ個のデータブロックから成る時系列データＤＳについて、Ｑ個のバイコヒーレンスの値が求められる。検定量計算部３２では、Ｑ個のバイコヒーレンスの値に基づいて時系列データＤＳについての検定量ｚγを求める。比較部３４では、求められた検定量ｚγを予め定められた閾値ｚｔｈと比較する。
【００４１】
バイコヒーレンス計算部２０は、例えば図６に示すように、バイスペクトル計算部２２、パワースペクトル計算部２４、積自乗平均計算部２６、及び規格化部２８を有する。
【００４２】
バイスペクトル計算部２２、パワースペクトル計算部２４、及び積自乗平均計算部２６は、周波数成分選択部１８から供給された周波数成分を受けて、それぞれバイスペクトル、パワースペクトル、及び積自乗平均を計算する。
【００４３】
周波数ペア（ｋ１，ｋ２）のバイコヒーレンスを計算する際は、バイスペクトル計算部２２には、周波数成分Ｘｉ［ｋ１］、Ｘｉ［ｋ２］、Ｘｉ［ｋ１＋ｋ２］（ｉ＝１〜Ｐ）が供給され、パワースペクトル計算部２４には、周波数成分Ｘｉ［ｋ１＋ｋ２］（ｉ＝１〜Ｐ）が供給され、積自乗平均計算部２６には、周波数成分Ｘｉ［ｋ１］、Ｘｉ［ｋ２］（ｉ＝１〜Ｐ）が供給される。
【００４４】
バイスペクトル計算部２２は、Ｐ個のデータブロックの各々についての、周波数ペア（ｋ１，ｋ２）のバイスペクトルを計算し、さらに該バイスペクトルの、Ｐ個のデータブロックについての平均を求め、該平均を、Ｐ個のデータブロックから成る時系列データＤＳについてのバイスペクトルとして出力する。即ち、まず各データブロック（ｉ番目のデータブロック）ＤＢｉについてのバイスペクトルＢｉ［ｋ１，ｋ２］を、次式（２）によって計算する。
【数９】
2015046116.tif 000011
ここで、
Ｘｉ［ｋ１］はｉ番目のデータブロックの周波数ｋ１の周波数成分であり、
Ｘｉ［ｋ２］はｉ番目のデータブロックの周波数ｋ２の周波数成分であり、
Ｘｉ＊［ｋ１＋ｋ２］はＸｉ［ｋ１＋ｋ２］の複素共役を表し、
Ｘｉ［ｋ１＋ｋ２］はｉ番目のデータブロックの周波数ｋ１＋ｋ２の周波数成分である。
【００４５】
式（２）で得られるバイスペクトルＢｉ［ｋ１，ｋ２］は、観測データｘ［ｎ］の周波数間の相関の程度を示す指標である。
【００４６】
バイスペクトル計算部２２は、式（２）の計算をＰ個のデータブロックＤＢｉ（ｉ＝１〜Ｐ）の各々に対して行い、計算の結果得られたＰ個のバイスペクトルＢｉ［ｋ１，ｋ２］の値を平均することによって、Ｐ個のデータブロックから成る時系列データＤＳについてのバイスペクトルＢａ［ｋ１，ｋ２］を得る。バイスペクトルＢａ［ｋ１，ｋ２］を求めるための計算式は、次式（３）で表される。
【数１０】
2015046116.tif 000012
【００４７】
式（２）と式（３）を組合せて書き直すと、次の式（９）となる。
【数１１】
2015046116.tif 000013
【００４８】
パワースペクトル計算部２４は、Ｐ個のデータブロックの各々についての、各周波数ペア（ｋ１，ｋ２）を構成する２つの周波数ｋ１、ｋ２の和に等しい周波数（ｋ１＋ｋ２）の周波数成分Ｘｉ［ｋ１＋ｋ２］のパワースペクトルを計算し、さらに該パワースペクトルの、Ｐ個のデータブロックについての平均を求めることによって、Ｐ個のデータブロックから成る時系列データＤＳについてのパワースペクトルＳａ［ｋ１＋ｋ２］を得る。即ち、時系列データＤＳについてのパワースペクトルＳａ［ｋ１＋ｋ２］を求めるための計算は、次の式（１０）で表される。
【数１２】
2015046116.tif 000014
ここで、
Ｘｉ［ｋ１＋ｋ２］はｉ番目のデータブロックの、周波数ｋ１＋ｋ２の周波数成分であり、
Ｘｉ＊［ｋ１＋ｋ２］は、Ｘｉ［ｋ１＋ｋ２］の複素共役である。
【００４９】
積自乗平均計算部２６は、Ｐ個のデータブロックの各々についての、各周波数ペア（ｋ１，ｋ２）を構成する２つの周波数ｋ１、ｋ２の周波数成分Ｘｉ［ｋ１］、Ｘｉ［ｋ２］の積の絶対値の自乗を求め、該自乗の、Ｐ個のデータブロックについての平均を求め、該平均をＰ個のデータブロックから成る時系列データＤＳについての積自乗平均Ｃａ［ｋ１，ｋ２］として出力する。積自乗平均Ｃａ［ｋ１，ｋ２］を求めるための計算は、次の式（１１）で表される。
【数１３】
2015046116.tif 000015
ここで、
Ｘｉ［ｋ１］はｉ番目のデータブロックの周波数ｋ１の周波数成分であり、
Ｘｉ［ｋ２］はｉ番目のデータブロックの周波数ｋ２の周波数成分である。
【００５０】
規格化部２８は、バイスペクトル計算部２２で計算されたバイスペクトルＢａ［ｋ１，ｋ２］と、パワースペクトル計算部２４で計算されたパワースペクトルＳａ［ｋ１＋ｋ２］と、積自乗平均計算部２６で計算された積自乗平均Ｃａ［ｋ１，ｋ２］に基づいて、バイコヒーレンスγ［ｋ１，ｋ２］を計算する。その計算は次の式（１２）で表される。
【数１４】
2015046116.tif 000016
この計算は、バイスペクトルＢａ［ｋ１，ｋ２］を、積自乗平均Ｃａ［ｋ１，ｋ２］とパワースペクトルＳａ［ｋ１＋ｋ２］の積の平方根で割ることで規格化することを意味する。
規格化部２８で算出されたバイコヒーレンスγ［ｋ１，ｋ２］は、バイコヒーレンス計算部２０の出力となる。
式（１２）で求められるバイコヒーレンスγ［ｋ１，ｋ２］は、周波数ペアごとに計算され、周波数ペアを構成する２つの周波数ｋ１及びｋ２の成分間の位相的な相関関係を表す複素数であり、その振幅が０〜１の間になるように規格化されており、周波数ｋ１とｋ２の関数である。式（５）または（６）で求められる規格化バイスペクトルは、振幅の大きさが観測データの歪度に依存している。歪度は、観測データの抽出の仕方によって大きさが変化する場合があるので、結果的に検定量の大きさに影響を及ぼす。
一方、式（１２）で求められるバイコヒーレンスは、周波数ｋ１及びｋ２の成分間の位相的な相関関係のみを抽出し、観測データの抽出の仕方によって大きさが変化する可能性のある歪度の影響を抑制している。
【００５１】
検定量計算部３２は、Ｑ個の周波数ペア（ｋ１，ｋ２）の各々についてバイコヒーレンス計算部２０で計算される、Ｐ個のデータブロックから成る時系列データＤＳについてのバイコヒーレンスγ［ｋ１，ｋ２］を、Ｑ個の周波数ペアのすべてについて積算することで、検定量ｚγを計算する。この計算は、算出される検定量ｚγの平均値及び分散が、各データブロックを構成するサンプル値の数Ｎ、時系列データＤＳを構成するデータブロックの数Ｐ、及びデータブロックのオーバーラップ率Ｌによらずに一定に保たれるようにするための調整量（下記のｂｉａｓ、σ２）を用いて行われるものであり、次の式（１３）で表される。
【数１５】
2015046116.tif 000017
【００５２】
式（１３）で、γｑ（但しｑは１乃至Ｑのいずれか）は、Ｑ個の周波数ペア（ｋ１，ｋ２）の各々についてのγ［ｋ１，ｋ２］と同じものを表す。
式（１３）の分子のｂｉａｓ項は分子の平均値を０とする。従って、ｚγの平均値を０とするための調整項、分母のσ２は検定量ｚγの分散を補正するための調整項であり、調整値計算部３８で、それぞれ以下の式（１４）、式（１５）によって求められる。
【００５３】
【数１６】
2015046116.tif 000018
【００５４】
式（１４）、式（１５）で、α及びβは定数であり、パラメータメモリ３６に記憶されている。
例えば、オーバーラップ率Ｌを５０％として、Ｐ＝３２個のデータブロックのＮ＝３２サンプル離散フーリエ変換を計算する場合には、α＝２、β＝６と設定すれば、観測データがガウス性を有し、且つ線形性を有する場合に、検定量ｚγは、平均０、分散１の正規分布に従う。
これによって検定量ｚγに対してある閾値を設定し、検定量ｚγが閾値よりも大きくなったときに、観測データが非ガウス性を有する、又は観測データが非線形性を有すると判断することが可能となる。
【００５５】
比較部３４は、検定量計算部３２で計算された検定量ｚγが予め定められた閾値ｚｔｈよりも大きいか否かを判定し、判定結果ＤＪを出力する。即ち閾値ｚｔｈよりも大きければ、非ガウス性又は非線形性があるとの判断結果を出力し、閾値ｚｔｈ以下であれば、非ガウス性も非線形性もないとの判断結果を出力する。判定結果ＤＪは、データ解析装置による解析の結果として出力され、例えば図４に示される表示部８に解析の結果が表示される。
【００５６】
上記のように、式（１３）の分子のｂｉａｓ項は分子の平均値を０とするための調整項、分母のσ２は検定量ｚγの分散を補正するための調整項であり、この２つの調整項は、Ｎ（時間方向のサンプル数）及びＰ（データブロック数）を変えた場合でも、ガウス性及び線形性の信号に対する検定量ｚγの大きさを一定に保つためのものであり、上記の式（１４）、式（１５）によって求められる。
即ち、パラメータα及びβは、Ｐ、Ｎ、Ｌに応じて定める必要がある。
【００５７】
式（１４）、式（１５）中のα、βは、以下のようにして定められる。即ち、パラメータα、βは、疑似的な白色雑音を用いることによって、実験的に求め、パラメータメモリ３６に記憶させておく。
【００５８】
具体的には、疑似白色雑音に対して式（１２）でバイコヒーレンスγｗ［ｋ１，ｋ２］を計算し、これを用いて式（１３）と同様の計算式で計算した検定量（試験的検定量）ｚｗの平均値ｚｗａと分散ｚｗｖからパラメータα、βを求める。
【００５９】
パラメータα、βの決定には、図７に示すパラメータ生成装置を用いる。このパラメータ生成装置のうち、図４、図５と同じ符号は、同じ機能を有するものであり、図４、図５の装置と共用することができる。
【００６０】
（ａ）  まず、白色雑音生成部４２により、疑似的な白色雑音を表す時系列データ（観測データと同じ符号ｘ［ｎ］で表す）を生成する。白色雑音の生成は、例えばソフトウェアにより行うことができる。生成された白色雑音を表す時系列データｘ［ｎ］を、図４に関連して説明した観測データ（観測値を表す時系列データ）の代わりに、時系列データメモリ４に記憶する。
【００６１】
（ｂ）  データブロック生成部１２により、時系列データメモリ４に記憶されている時系列データｘ［ｎ］を、各々Ｎサンプル値から成るＰ個のデータブロックを、オーバーラップ率Ｌで読み出し、フーリエ変換部１４に供給する。
フーリエ変換部１４、及びバイコヒーレンス計算部２０は、時系列データメモリ４から読み出された白色雑音を表すデータｘ［ｎ］に対して、図４を参照して説明した観測データに対する処理と同様の処理を行い、バイコヒーレンス計算部２０からは、計算されたバイコヒーレンスが出力される。出力されるバイコヒーレンスを、観測データに対する処理の場合と異なる記号γｗ［ｋ１，ｋ２］で表す。
【００６２】
（ｃ）  検定量計算部３２は、バイコヒーレンスγｗ［ｋ１，ｋ２］に基づいて検定量（試験的検定量）ｚｗを計算する。この計算は次の式（１６）で表される。
【数１７】
2015046116.tif 000019
式（１６）でγｗｑの添え字ｑは、それぞれ異なる周波数ペア（ｋ１，ｋ２）に対応するものであり、Ｑ個の周波数ペアを用いる場合、ｑは１からＱまでの値を取る。例えば、図３に示されるのと同じ配置の周波数ペアを用いる場合、Ｑは例えば４９である。
【００６３】
（ｄ）  上記（ａ）〜（ｃ）の処理を複数回繰り返し、複数回、即ちＲ回の試行で得られた検定量ｚｗｒ（ｒは１〜Ｒのいずれか）の平均及び分散を平均分散計算部４４に供給する。平均分散計算部４４は、複数回の試行で得られた検定量ｚｗｒの平均ｚｗａ、分散ｚｗｖを計算する。この計算は、次の式（１７）、（１８）で表される。
【００６４】
【数１８】
2015046116.tif 000020
【００６５】
平均分散計算部４４は、このようにして求めた平均値ｚｗａ、分散ｚｗｖをパラメータ計算部４６に供給する。パラメータ計算部４６では、平均値ｚｗａ、分散ｚｗｖに基づいてパラメータα、βを計算する。この計算は下記の式（１９）、（２０）で表される。
【数１９】
2015046116.tif 000021
【００６６】
ここで、Ｎ、Ｐはそれぞれ白色雑音を表す時系列データを読み出して検出処理部６に供給する際のデータブロック数及びサンプル数である。
【００６７】
（ｅ）  手順（ａ）〜（ｄ）までの計算を、データブロック数Ｐ、サンプル数Ｎ及びオーバーラップ率Ｌの異なる値について行って、それぞれα、βの値を求め、パラメータメモリ３６に記憶する。例えば、Ｐ、Ｎ、Ｌに対応付けて、例えばルックアップテーブルの形式で記憶する。
【００６８】
観測データに対する解析を行う際には、その時のＰ、Ｎ、Ｌに対応するα、βを読み出して、式（１４）、式（１５）の計算に用いる。
【００６９】
なお、上記のように、それぞれの設定条件（Ｐ、Ｎ、Ｌ）に対して、α及びβをテーブル形式で記憶しておく代わりに、設定条件Ｐ、Ｎ、Ｌごとにｂｉａｓ、σ２を計算してメモリに記憶させておくこととしても良い。その場合には、図５の調整値計算部３８は不要となる。
【００７０】
以上のように、本発明では、従来の方法と同様に、観測データからＮサンプル×Ｐデータブロックのデータを抽出し、Ｐ個のＮサンプル離散フーリエ変換を計算する。その後、Ｐ個のデータブロックの各々のフーリエ変換の出力を用いて、パワースペクトル、バイスペクトルを計算するところまでは、従来の方法と同じである。
【００７１】
バイコヒーレンスγ［ｋ１，ｋ２］（式（１２））と、従来の規格化バイスペクトルＢｎ［ｋ１，ｋ２］（式（６））との違いは、式中の分母にある。従来の規格化バイスペクトルＢｎ［ｋ１，ｋ２］の分母が
【数２０】
2015046116.tif 000022
であるのに対し、バイコヒーレンスγ［ｋ１，ｋ２］の分母は
【数２１】
2015046116.tif 000023
であり、両者の違いは、Ｓａ［ｋ１］Ｓａ［ｋ２］と、Ｃａ［ｋ１，ｋ２］の部分にある。この違いによって、バイコヒーレンスγ［ｋ１，ｋ２］は、式の上では信号Ｘ［ｋ１］Ｘ［ｋ２］と信号Ｘ［ｋ１＋ｋ２］間のコヒーレンスを求めているが、結果的には観測データの周波数ｋ１及びｋ２の成分間の位相的な相関関係を評価していることになる。
【００７２】
なお、バイスペクトル自体も、周波数間の位相関係を表すが、バイスペクトルは、これとともに、観測データの振幅にも依存する。これに対して、本発明で用いるバイコヒーレンスは、式（１２）に示すように規格化を行うことで、振幅の大きさを０〜１の間になるようにして、２つの周波数成分間の位相に関する情報のみを抽出するようにしたものである。
【００７３】
観測データが非ガウス性及び非線形性を有していれば、周波数ｋ１とｋ２の信号間に位相的なつながりがあるので、バイコヒーレンスは高い値を有することとなる。以上から、バイコヒーレンスによって観測データが非ガウス性を有し、かつ非線形性を有するか否かを確認することが可能となる。
【００７４】
本発明の方法の効果を確認するため、模擬的に生成した過渡的信号を用いて実験を行った。過渡的な信号は、非ガウス性及び非線形性を有している。図８（ａ）及び（ｂ）に、模擬した過渡的信号Ｓｔｒ、及び該過渡的信号Ｓｔｒに白色雑音Ｎｗｔを付加した波形を示す。図８（ａ）は生成した過渡的信号であり、過渡的信号はある起点からの経過時間５、１０、１５、２０秒の時刻に現れている。図８（ｂ）は、生成した過渡的信号Ｓｔｒに白色雑音Ｎｗｔを付加した波形であり、この波形に対して従来の方法及び本発明の方法を適用した。
【００７５】
図９（ａ）及び（ｂ）は、生成した信号（図８（ｂ）に示す信号）に対して、従来の方法及び本発明の方法を適用した結果を示す図である。図９（ａ）は、従来の方法を適用した結果を示し、図９（ｂ）は、本発明の方法を適用した結果を示す。なお、計算パラメータは以下のように設定した。
・１データブロックあたりのサンプル数Ｎ＝３２
・データブロック数Ｐ＝３２
・オーバーラップ率Ｌ＝５０％
・周波数ペア数Ｑ＝４９
・検定量ｚγの計算に用いられる調整値ｂｉａｓ＝０．０６２５
・検定量ｚγの計算に用いられる調整値σ２＝０．１８７５
なお、サンプリング周波数は５０ｋＨｚであり、
周波数ペアの２次元平面上の配置は、図３に示す如くであり、図３に示される周波数ｆｍ（ｍ＝１〜１４）はそれぞれ以下の値を有する。
ｍ＝１のとき、ｆｍ＝０
１＜ｍ＜１４のとき、ｆｍ＝ｆ（ｍ−１）＋ΔＦ
但し、ΔＦ＝１５６２．５Ｈｚである。
【００７６】
図９（ａ）及び（ｂ）の２つの処理の結果を比べると、白色雑音のみを観測しているときのレベルは同等である。一方、過渡的信号を観測しているときの応答レベルを見ると、本発明の方法は、従来の方法の概ね２倍以上の大きさである。したがって、本発明の方法は、従来の方法よりも出力ＳＮＲが高く、観測データの非ガウス性及び非線形性の検出能力が高いことが確認できる。
【００７７】
以上本発明を、データ解析装置として説明したが、データ解析装置で実施されるデータ解析方法もまた本発明の一部を成す。さらに、データ解析装置における各部の処理、又はデータ解析方法における各処理をコンピュータに実行させるためのプログラム、及びそのようなプログラムを記録したコンピュータで読み取り可能な記録媒体もまた本発明の一部を成す。
【産業上の利用可能性】
【００７８】
以上観測データが可聴域の音響信号である場合について説明したが、本発明はこれに限定されず、例えば、超音波信号の解析、機械振動の解析、脳波の解析にも適用できる。
【符号の説明】
【００７９】
２  入力信号処理部、  ４  時系列データメモリ、  ６  検出処理部、  １２  データブロック生成部、  １４  フーリエ変換部、  １６  周波数成分メモリ、  １８  周波数成分選択部、  ２０  バイコヒーレンス計算部、  ２２  バイスペクトル計算部、  ２４  パワースペクトル計算部、  ２６  積自乗平均計算部、  ２８  規格化部、  ３２  検定量計算部、  ３４  比較部、  ３６  パラメータメモリ、  ３８  調整値計算部、  ４２  白色雑音生成部、  ４４  平均分散計算部、  ４６  パラメータ計算部。
【図１】
2015046116.tif 000040
【図２】
2015046116.tif 000041
【図３】
2015046116.tif 000042
【図４】
2015046116.tif 000043
【図５】
2015046116.tif 000044
【図６】
2015046116.tif 000045
【図７】
2015046116.tif 000046
【図８】
2015046116.tif 000047
【図９】
2015046116.tif 000048