特許 6183067 データ解析装置及び方法、並びにプログラム及び記録媒体

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

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

(11)【特許番号】6183067

(24)【登録日】2017年8月4日

(45)【発行日】2017年8月23日

(54)【発明の名称】データ解析装置及び方法、並びにプログラム及び記録媒体

(51)【国際特許分類】

   G06F 17/10 20060101AFI20170814BHJP

   G01R 23/16 20060101ALI20170814BHJP

【ＦＩ】

   G06F17/10 Z

   G01R23/16 D

   G01R23/16 Z

【請求項の数】20

【全頁数】23

(21)【出願番号】特願2013-178142(P2013-178142)

(22)【出願日】2013年8月29日

(65)【公開番号】特開2015-46116(P2015-46116A)

(43)【公開日】2015年3月12日

【審査請求日】2016年5月17日

(73)【特許権者】

【識別番号】000000295

【氏名又は名称】沖電気工業株式会社

(74)【代理人】

【識別番号】100083840

【弁理士】

【氏名又は名称】前田 実

(74)【代理人】

【識別番号】100116964

【弁理士】

【氏名又は名称】山形 洋一

(74)【代理人】

【識別番号】100135921

【弁理士】

【氏名又は名称】篠原 昌彦

(72)【発明者】

【氏名】高嶋 昭一

【審査官】 小林 哲雄

(56)【参考文献】

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

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

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

【文献】 米国特許出願公開第２００７／０２８８１０３（ＵＳ，Ａ１）

【文献】 米国特許出願公開第２００９／００８２９７６（ＵＳ，Ａ１）

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

Ｇ０６Ｆ １７／１０

Ｇ０１Ｒ ２３／１６

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

【請求項1】

観測値を表す時系列データ中の、各々相連続するサンプル値で構成される複数個のデータブロックをフーリエ変換して、各データブロックについて複数の周波数成分を計算するフーリエ変換部と、
前記フーリエ変換部で計算された各データブロックについての周波数成分のうちの２つの周波数からなる周波数ペアの各々のバイコヒーレンスを計算するバイコヒーレンス計算部と、
前記周波数ペアの各々についての前記バイコヒーレンスに基づいて、前記時系列データについての検定量を計算する検定量計算部と、
前記検定量計算部で計算された前記検定量が予め定められた閾値よりも大きいか否かを判断する比較部とを有し、
前記比較部による判断結果が、解析の結果として出力され、
前記検定量計算部は、前記検定量の平均値及び分散が、各データブロックを構成するサンプル値の数、前記時系列データを構成するデータブロックの数、及び
前記データブロックのオーバーラップ率によらずに一定に保たれるようにするための調整量を用いて、前記検定量の計算を行う
データ解析装置。

【請求項2】

観測値を表す時系列データ中の、各々相連続するサンプル値で構成される複数個のデータブロックをフーリエ変換して、各データブロックについて複数の周波数成分を計算するフーリエ変換部と、
前記フーリエ変換部で計算された各データブロックについての周波数成分のうちの２つの周波数からなる周波数ペアの各々のバイコヒーレンスを計算するバイコヒーレンス計算部と、
前記周波数ペアの各々についての前記バイコヒーレンスに基づいて、前記時系列データについての検定量を計算する検定量計算部と、
前記検定量計算部で計算された前記検定量が予め定められた閾値よりも大きいか否かを判断する比較部とを有し、
前記比較部による判断結果が、解析の結果として出力され、
前記検定量計算部は、前記検定量をｚ_γで表すとき、

【数27】

（但し、
Ｑは前記周波数ペアの数、
γ_ｑは、前記周波数ペアの各々について、前記バイコヒーレンス計算部で計算されたバイコヒーレンス、
ｂｉａｓは、前記検定量の平均を０とするための調整量、
σ^２は、前記検定量の分散を補正するための調整量）
により前記検定量ｚ_γを計算する
データ解析装置。

【請求項3】

前記調整量ｂｉａｓ及びσ^２を、

【数28】

（但し、
Ｐは、前記時系列データを構成するデータブロックの数、
Ｎは各データブロックを構成するサンプル値の数、
α及びβはパラメータ）
により計算する調整量計算部をさらに有する
ことを特徴とする請求項２に記載のデータ解析装置。

【請求項4】

白色雑音を表す時系列データを、前記観測値を表す時系列データの代わりに用いて、
前記フーリエ変換及び前記バイコヒーレンスの計算を行い、
該計算により求められたバイコヒーレンスに基づいて試験的検定量を計算し、
複数回の試行についての前記試験的検定量の平均ｚｗａ及び分散ｚｗｖを計算し、
該計算により求められた平均ｚｗａ及び分散ｚｗｖに基づいて、

【数29】

により、前記パラメータα及びβを生成するパラメータ生成装置と、
前記パラメータ生成装置で生成された前記パラメータα及びβを記憶するパラメータメモリと
をさらに有することを特徴とする請求項３に記載のデータ解析装置。

【請求項5】

前記バイコヒーレンス計算部は、
前記各データブロックについての前記周波数ペアの各々のバイスペクトルの、前記複数のデータブロックについての平均を計算し、前記時系列データについてのバイスペクトルとして出力するバイスペクトル計算部と、
前記各データブロックについての前記周波数ペアの各々を構成する周波数の和に等しい周波数の成分のパワースペクトルの、前記複数のデータブロックについての平均を計算して、前記時系列データについてのパワースペクトルとして出力するパワースペクトル計算部と、
前記各データブロックについての前記周波数ペアの各々を構成する周波数の成分の積の絶対値の自乗の、前記複数のデータブロックについての平均を計算し、前記時系列データについての積自乗平均として出力する積自乗平均計算部と、
前記バイスペクトル計算部から出力された前記時系列データについての前記バイスペクトル、前記パワースペクトル計算部から出力された前記時系列データについての前記パワースペクトル、及び前記積自乗平均計算部から出力された前記時系列データについての前記積自乗平均に基づいて、前記バイコヒーレンスを計算する規格化部とを有する
ことを特徴とする請求項１乃至４のいずれかに記載のデータ解析装置。

【請求項6】

前記バイスペクトル計算部は、前記各データブロックについての前記バイスペクトルをＢ_ｉ［ｋ_１，ｋ_２］で表すとき、

【数22】

（但し、
ｋ_１、ｋ_２はそれぞれ前記周波数ペアを構成する周波数、
Ｘ_ｉ[ｋ_１]、Ｘ_ｉ[ｋ_２]、Ｘ_ｉ[ｋ_１＋ｋ_２]はそれぞれ、ｉ番目のデータブロック（ｉは１乃至Ｐのいずれか）の、周波数ｋ_１、ｋ_２、ｋ_１＋ｋ_２の周波数成分）
により、前記各データブロックについての前記バイスペクトルＢ_ｉ［ｋ_１，ｋ_２］を計算し、
前記時系列データについての前記バイスペクトルをＢａ［ｋ_１，ｋ_２］で表すとき、

【数23】

（但し、Ｐは前記時系列データを構成するデータブロックの数）
により、前記時系列データについての前記バイスペクトルＢａ［ｋ_１，ｋ_２］を計算する
ことを特徴とする請求項５に記載のデータ解析装置。

【請求項7】

前記積自乗平均計算部は、前記時系列データについての積自乗平均をＣａ［ｋ_１，ｋ_２］で表すとき、

【数24】

（但し、
ｋ_１、ｋ_２はそれぞれ前記周波数ペアを構成する周波数、
Ｘ_ｉ[ｋ_１]、Ｘ_ｉ[ｋ_２]はそれぞれ、ｉ番目のデータブロック（ｉは１乃至Ｐのいずれか）の、周波数ｋ_１、ｋ_２の周波数成分、
Ｐは前記時系列データを構成するデータブロックの数）
により、前記時系列データについての前記積自乗平均Ｃａ［ｋ_１，ｋ_２］を計算する
ことを特徴とする請求項５又は６に記載のデータ解析装置。

【請求項8】

前記パワースペクトル計算部は、前記時系列データについてのパワースペクトルをＳａ［ｋ_１＋ｋ_２］で表すとき、

【数25】

（但し、
ｋ_１、ｋ_２はそれぞれ前記周波数ペアを構成する周波数、
Ｘ_ｉ［ｋ_１＋ｋ_２］はｉ番目のデータブロック（ｉは１乃至Ｐのいずれか）の、周波数ｋ_１＋ｋ_２の周波数成分であり、
Ｘ_ｉ^＊［ｋ_１＋ｋ_２］は、Ｘ_ｉ［ｋ_１＋ｋ_２］の複素共役であり、
Ｐはデータブロックの数である。）
により、前記時系列データについての前記パワースペクトルＳａ［ｋ_１＋ｋ_２］を計算する
ことを特徴とする請求項５、６又は７に記載のデータ解析装置。

【請求項9】

前記規格化部は、前記バイコヒーレンスをγ［ｋ_１，ｋ_２］で表すとき、

【数26】

により、前記バイコヒーレンスγ［ｋ_１，ｋ_２］を計算する
ことを特徴とする請求項５乃至８のいずれかに記載のデータ解析装置。

【請求項10】

観測値を表す時系列データ中の、各々相連続するサンプル値で構成される複数個のデータブロックをフーリエ変換して、各データブロックについて複数の周波数成分を計算するフーリエ変換ステップと、
前記フーリエ変換ステップで計算された各データブロックについての周波数成分のうちの２つの周波数からなる周波数ペアの各々のバイコヒーレンスを計算するバイコヒーレンス計算ステップと、
前記周波数ペアの各々についての前記バイコヒーレンスに基づいて、前記時系列データについての検定量を計算する検定量計算ステップと、
前記検定量計算ステップで計算された前記検定量が予め定められた閾値よりも大きいか否かを判断する比較ステップとを有し、
前記比較ステップによる判断結果が、解析の結果として出力され、
前記検定量計算ステップは、前記検定量の平均値及び分散が、各データブロックを構成するサンプル値の数、前記時系列データを構成するデータブロックの数、及び
前記データブロックのオーバーラップ率によらずに一定に保たれるようにするための調整量を用いて、前記検定量の計算を行う
データ解析方法。

【請求項11】

観測値を表す時系列データ中の、各々相連続するサンプル値で構成される複数個のデータブロックをフーリエ変換して、各データブロックについて複数の周波数成分を計算するフーリエ変換ステップと、
前記フーリエ変換ステップで計算された各データブロックについての周波数成分のうちの２つの周波数からなる周波数ペアの各々のバイコヒーレンスを計算するバイコヒーレンス計算ステップと、
前記周波数ペアの各々についての前記バイコヒーレンスに基づいて、前記時系列データについての検定量を計算する検定量計算ステップと、
前記検定量計算ステップで計算された前記検定量が予め定められた閾値よりも大きいか否かを判断する比較ステップとを有し、
前記比較ステップによる判断結果が、解析の結果として出力され、
前記検定量計算ステップは、前記検定量をｚ_γで表すとき、

【数35】

（但し、
Ｑは前記周波数ペアの数、
γ_ｑは、前記周波数ペアの各々について、前記バイコヒーレンス計算ステップで計算されたバイコヒーレンス、
ｂｉａｓは、前記検定量の平均を０とするための調整量、
σ^２は、前記検定量の分散を補正するための調整量）
により前記検定量ｚ_γを計算する
データ解析方法。

【請求項12】

前記調整量ｂｉａｓ及びσ^２を、

【数36】

（但し、
Ｐは、前記時系列データを構成するデータブロックの数、
Ｎは各データブロックを構成するサンプル値の数、
α及びβはパラメータ）
により計算する調整量計算ステップをさらに有する
ことを特徴とする請求項１１に記載のデータ解析方法。

【請求項13】

白色雑音を表す時系列データを、前記観測値を表す時系列データの代わりに用いて、
前記フーリエ変換及び前記バイコヒーレンスの計算を行い、
該計算により求められたバイコヒーレンスに基づいて試験的検定量を計算し、
複数回の試行についての前記試験的検定量の平均ｚｗａ及び分散ｚｗｖを計算し、
該計算により求められた平均ｚｗａ及び分散ｚｗｖに基づいて、

【数37】

により、前記パラメータα及びβを生成するパラメータ生成ステップと、
前記パラメータ生成ステップで生成された前記パラメータα及びβをパラメータメモリに記憶させる記憶ステップと
をさらに有することを特徴とする請求項１２に記載のデータ解析方法。

【請求項14】

前記バイコヒーレンス計算ステップは、
前記各データブロックについての前記周波数ペアの各々のバイスペクトルの、前記複数のデータブロックについての平均を計算し、前記時系列データについてのバイスペクトルとして出力するバイスペクトル計算ステップと、
前記各データブロックについての前記周波数ペアの各々を構成する周波数の和に等しい周波数の成分のパワースペクトルの、前記複数のデータブロックについての平均を計算して、前記時系列データについてのパワースペクトルとして出力するパワースペクトル計算ステップと、
前記各データブロックについての前記周波数ペアの各々を構成する周波数の成分の積の絶対値の自乗の、前記複数のデータブロックについての平均を計算し、前記時系列データについての積自乗平均として出力する積自乗平均計算ステップと、
前記バイスペクトル計算ステップから出力された前記時系列データについての前記バイスペクトル、前記パワースペクトル計算ステップから出力された前記時系列データについての前記パワースペクトル、及び前記積自乗平均計算ステップから出力された前記時系列データについての前記積自乗平均に基づいて、前記バイコヒーレンスを計算する規格化ステップとを有する
ことを特徴とする請求項１０乃至１３のいずれかに記載のデータ解析方法。

【請求項15】

前記バイスペクトル計算ステップは、前記各データブロックについての前記バイスペクトルをＢ_ｉ［ｋ_１，ｋ_２］で表すとき、

【数30】

（但し、
ｋ_１、ｋ_２はそれぞれ前記周波数ペアを構成する周波数、
Ｘ_ｉ[ｋ_１]、Ｘ_ｉ[ｋ_２]、Ｘ_ｉ[ｋ_１＋ｋ_２]はそれぞれ、ｉ番目のデータブロック（ｉは１乃至Ｐのいずれか）の、周波数ｋ_１、ｋ_２、ｋ_１＋ｋ_２の周波数成分）
により、前記各データブロックについての前記バイスペクトルＢ_ｉ［ｋ_１，ｋ_２］を計算し、
前記時系列データについての前記バイスペクトルをＢａ［ｋ_１，ｋ_２］で表すとき、

【数31】

（但し、Ｐは前記時系列データを構成するデータブロックの数）
により、前記時系列データについての前記バイスペクトルＢａ［ｋ_１，ｋ_２］を計算する
ことを特徴とする請求項１４に記載のデータ解析方法。

【請求項16】

前記積自乗平均計算ステップは、前記時系列データについての積自乗平均をＣａ［ｋ_１，ｋ_２］で表すとき、

【数32】

（但し、
ｋ_１、ｋ_２はそれぞれ前記周波数ペアを構成する周波数、
Ｘ_ｉ[ｋ_１]、Ｘ_ｉ[ｋ_２]はそれぞれ、ｉ番目のデータブロック（ｉは１乃至Ｐのいずれか）の、周波数ｋ_１、ｋ_２の周波数成分、
Ｐは前記時系列データを構成するデータブロックの数）
により、前記時系列データについての前記積自乗平均Ｃａ［ｋ_１，ｋ_２］を計算する
ことを特徴とする請求項１４又は１５に記載のデータ解析方法。

【請求項17】

前記パワースペクトル計算ステップは、前記時系列データについてのパワースペクトルをＳａ［ｋ_１＋ｋ_２］で表すとき、

【数33】

（但し、
ｋ_１、ｋ_２はそれぞれ前記周波数ペアを構成する周波数、
Ｘ_ｉ［ｋ_１＋ｋ_２］はｉ番目のデータブロック（ｉは１乃至Ｐのいずれか）の、周波数ｋ_１＋ｋ_２の周波数成分であり、
Ｘ_ｉ^＊［ｋ_１＋ｋ_２］は、Ｘ_ｉ［ｋ_１＋ｋ_２］の複素共役であり、
Ｐはデータブロックの数である。）
により、前記時系列データについての前記パワースペクトルＳａ［ｋ_１＋ｋ_２］を計算する
ことを特徴とする請求項１４、１５又は１６に記載のデータ解析方法。

【請求項18】

前記規格化ステップは、前記バイコヒーレンスをγ［ｋ_１，ｋ_２］で表すとき、

【数34】

により、前記バイコヒーレンスγ［ｋ_１，ｋ_２］を計算する
ことを特徴とする請求項１４乃至１７のいずれかに記載のデータ解析方法。

【請求項19】

請求項１０乃至１８のいずれかに記載のデータ解析方法の各処理をコンピュータに実行させるためのプログラム。

【請求項20】

請求項１９に記載のプログラムを記録した、コンピュータで読み取り可能な記録媒体。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、データ解析装置及び方法、並びにプログラム及び記録媒体に関し、特に時系列データの非ガウス性又は非線形性の検出に関する。

【背景技術】

【0002】

時系列データの非ガウス性とは、時系列データの振幅分布が、平均値を中心として左右対称でないことを意味する。一方時系列データの非線形性とは、時系列データの生成過程に、非線形項が含まれていることを意味する。非線形性のデータは、非ガウス性を有する。
取得したデータ（観測データ）の非ガウス性及び非線形性を確認するためのデータ解析方法として、データの歪度に関係する規格化バイスペクトルを用いることが提案されている（非特許文献１〜３）。

【0003】

従来の方法の問題点は、観測データがガウス性であっても有限のデータを抽出すると、抽出の仕方によってデータが非ガウス性又は非線形性を有するように見え、誤判定が生じる可能性があることである。

【0004】

以下、従来の方法について図１を参照して説明する。図１において、観測の結果得られた時系列データをｘ（ｔ）と表し、このｘ（ｔ）の離散化されたものをｘ［ｎ］と表す。ここで、ｎは時刻を表すインデックスである。
時系列データｘ［ｎ］から図２に示すように、複数個即ちＰ個のデータブロックＤＢ_１〜ＤＢ_Ｐを生成する。データブロックの各々ＤＢ_ｉ（ｉは１乃至Ｐのいずれか）はＮ個のサンプル値を含む。
各データブロックＤＢ_ｉに対して、フーリエ変換部１４でＮサンプル離散フーリエ変換を行って、各周波数の成分Ｘ_ｉ［ｋ］を求める。フーリエ変換は、下記の式（１）で表される。

【0005】

【数1】

ここで、ｋは周波数を表すインデックスである。

【0006】

バイスペクトル計算部２２及びパワースペクトル計算部２４には、複数個即ちＱ個の周波数ペアを構成する周波数についての周波数成分Ｘ_ｉ［ｋ］が供給される。例えば、Ｑは４９であり、４９個の周波数ペアの２次元周波数平面上における配置は図３に「○」印で示される通りである。

【0007】

バイスペクトル計算部２２では、データブロックＤＢ_１〜ＤＢ_Ｐの各々について周波数成分Ｘ_ｉ［ｋ］を用いて、バイスペクトルＢａ［ｋ_１，ｋ_２］を計算する。
そのために、まず各データブロック（ｉ番目のデータブロック）ＤＢ_ｉについてのバイスペクトルＢ_ｉ［ｋ_１，ｋ_２］を、次式（２）によって計算する。

【数2】

但し、Ｘ_ｉ^＊［ｋ_１＋ｋ_２］はＸ_ｉ［ｋ_１＋ｋ_２］の複素共役を表す。

【0008】

式（２）で得られるバイスペクトルＢ_ｉ［ｋ_１，ｋ_２］は、観測データｘ［ｎ］の周波数間の相関の程度を示す指標である。

【0009】

バイスペクトル計算部２２は、式（２）の計算を、Ｐ個のデータブロックＤＢ_１〜ＤＢ_Ｐの各々に対して行い、計算の結果得られたバイスペクトルＢ_ｉ［ｋ_１，ｋ_２］を平均することによって、Ｐ個のデータブロックから成る時系列データＤＳについてのバイスペクトルＢａ［ｋ_１，ｋ_２］を得る。バイスペクトルＢａ［ｋ_１，ｋ_２］を求めるための計算は、次式（３）で表される。

【0010】

【数3】

【0011】

観測データに白色雑音しか含まれていなければ、観測データの周波数間に相関関係はないので、バイスペクトルは小さな値となる。一方、観測データに非定常成分が含まれていれば、観測データの周波数間に相関関係が現れ、式（３）で示されるバイスペクトルが大きな値を持つ。

【0012】

式（３）で示されたバイスペクトルは、観測データに含まれる各周波数の成分の大きさ（パワー）に依存した値を持つ。そのため、このバイスペクトルそのものでは、周波数間の相関関係の大きさが判断しにくい。そこで通常は、バイスペクトルの大きさが観測データの各周波数の成分の大きさに依存しない、規格化した値が用いられる。

【0013】

規格化のため、パワースペクトル計算部２４では、観測データのパワースペクトルＳａ［ｋ］を、次のように計算する。

【数4】

ここで、
Ｘ_ｉ^＊［ｋ］は、Ｘ_ｉ［ｋ］の複素共役である。
式（４）の計算は、周波数ｋ_１、ｋ_２、ｋ_３（＝ｋ_１＋ｋ_２）の各々について行われ、従って、パワースペクトル計算部２４からは、Ｓａ［ｋ_１］、Ｓａ［ｋ_２］、Ｓａ［ｋ_１＋ｋ_２］が出力される。

【0014】

規格化部２７は、バイスペクトル計算部２２で計算されたバイスペクトルＢａ［ｋ_１，ｋ_２］と、パワースペクトル計算部２４で計算されたパワースペクトルＳａ［ｋ_１］、Ｓａ［ｋ_２］及びＳａ［ｋ_１＋ｋ_２］を用いて、規格化バイスペクトルＢｎ［ｋ_１，ｋ_２］を、次式（５）により、計算する。

【数5】

【0015】

分母にパワースペクトルの項があることによって、観測データの各周波数の成分の大きさが規格化される。従って、式（５）の規格化バイスペクトルＢｎ［ｋ_１，ｋ_２］は観測データの大きさ（パワー）に依存せず、周波数間の相関に応じた値を示すものとなる。
但し、実際の計算では、式（５）のバイスペクトルＢｎ［ｋ_１，ｋ_２］の分散を１にするために、次式（６）のように補正することで得られる、補正された規格化バイスペクトルＢｎｃ［ｋ_１，ｋ_２］を用いる。

【数6】

【0016】

式（６）で表される規格化バイスペクトルＢｎｃ［ｋ_１，ｋ_２］は、特定の周波数ペア（ｋ_１，ｋ_２）について計算された値であり、Ｑ個の周波数ペアについて上記した処理が繰り返されて、Ｑ個の周波数ペアの各々について式（６）で表される規格化バイスペクトルＢｎｃ［ｋ_１，ｋ_２］が求められる。

【0017】

検定量計算部２９は、上記のようにしてＱ個の周波数ペアについて求められた規格化バイスペクトルＢｎｃ［ｋ_１，ｋ_２］を用いて、観測データの非ガウス性及び非線形性を検出する。

【0018】

観測データが白色雑音であれば、それぞれＱ個の周波数ペア（ｋ_１，ｋ_２）について求められた規格化バイスペクトルＢｎｃ［ｋ_１，ｋ_２］は、平均０、分散１の正規分布に従う。この性質を利用して、観測データの非定常性を検出するために、次に示す検定量ｚを計算する。

【0019】

【数7】

【0020】

式（７）でＢｎｃ_ｑ（但しｑは１乃至Ｑのいずれか）は、Ｑ個の周波数ペア（ｋ_１，ｋ_２）の各々についてのＢｎｃ［ｋ_１，ｋ_２］と同じものを表す。
式（７）の分子はＱ個の規格化バイスペクトルの値を加算することを表している。式（７）の分母は、検定量ｚの分散を規格化するための項である。
観測データに白色雑音のみが含まれており、Ｑが十分に大きい値であれば、検定量ｚは平均０、分散１の正規分布に従う。もし、観測データに非ガウス性又は非線形性の成分が含まれていれば、規格化バイスペクトルの値が大きくなり、その結果、検定量ｚは平均０、分散１の正規分布には従わずに大きな値となる。このことを利用して、観測データの非ガウス性又は非線形性が検出できる。
以上が規格化バイスペクトルを用いて観測データから非ガウス性及又は非線形性を検出するための従来の方法である。

【先行技術文献】

【非特許文献】

【0021】

【非特許文献1】Ｍｅｌｖｉｎ Ｊ． Ｈｉｎｉｃｈ， “Ｂｉｓｐｅｃｔｒｕｍ ｏｆ ｓｈｉｐ−ｒａｄｄｉａｔｅｄ ｎｏｉｓｅ”， Ｊ． Ａｃｏｕｓｔ． Ｓｏｃ． Ａｍ． ８５（４）， Ａｐｒｉｌ １９８９

【非特許文献2】Ｍｅｌｖｉｎ Ｊ． Ｈｉｎｉｃｈ， “Ｄｅｔｅｃｔｉｎｇ ａ Ｔｒａｎｓｉｅｎｔ Ｓｉｇｎａｌ ｂｙ Ｂｉｓｐｅｃｔｒａｌ Ａｎａｌｙｓｉｓ”， ＩＥＥＥ Ｔｒａｎｓ． Ｏｎ Ａｃｏｕｓｔｉｃｓ， Ｓｐｅｅｃｈ ａｎｄ Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓｉｎｇ， Ｖｏｌ．３６， Ｎｏ．７， Ｊｕｌｙ １９９０

【非特許文献3】Ｍｅｌｖｉｎ Ｊ． Ｈｉｎｉｃｈ ｅｔ． ａｌ， “Ｏｎ ｔｈｅ Ｐｒｉｎｃｉｐａｌ Ｄｏｍａｉｎ ｏｆ ｔｈｅ Ｄｉｓｃｒｅｔｅ Ｂｉｓｐｅｃｔｒｕｍ ｏｆ ａ Ｓｔａｔｉｏｎａｒｙ Ｓｉｇｎａｌ”， ＩＥＥＥ Ｔｒａｎｓ． Ｏｎ Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓｉｎｇ， Ｖｏｌ．４３， Ｎｏ．９， Ｓｅｐｔｅｍｂｅｒ １９９５

【発明の概要】

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

【0022】

上記した従来の方法は、観測データが非ガウス性又は非線形性を有すれば、観測データがある程度の歪度を有しているということに着目した処理である。しかしながら、観測データにガウス性又は線形性があっても、有限のデータを抽出すると、抽出の仕方によっては非ガウス性又は非線形性を有するように見え、そのため判定に誤りが生じる可能性がある。

【0023】

本発明は、上記の課題の解決のためのものであり、その目的は、観測データの大きさ、歪度の影響が抑制された処理結果を得ることを可能にすることにある。

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

【0024】

本発明の第１の態様のデータ解析装置は、
観測値を表す時系列データ中の、各々相連続するサンプル値で構成される複数個のデータブロックをフーリエ変換して、各データブロックについて複数の周波数成分を計算するフーリエ変換部と、
前記フーリエ変換部で計算された各データブロックについての周波数成分のうちの２つの周波数からなる周波数ペアの各々のバイコヒーレンスを計算するバイコヒーレンス計算部と、
前記周波数ペアの各々についての前記バイコヒーレンスに基づいて、前記時系列データについての検定量を計算する検定量計算部と、
前記検定量計算部で計算された前記検定量が予め定められた閾値よりも大きいか否かを判断する比較部とを有し、
前記比較部による判断結果が、解析の結果として出力され、
前記検定量計算部は、前記検定量の平均値及び分散が、各データブロックを構成するサンプル値の数、前記時系列データを構成するデータブロックの数、及び
前記データブロックのオーバーラップ率によらずに一定に保たれるようにするための調整量を用いて、前記検定量の計算を行うものである。
本発明の第２の態様のデータ解析装置は、
観測値を表す時系列データ中の、各々相連続するサンプル値で構成される複数個のデータブロックをフーリエ変換して、各データブロックについて複数の周波数成分を計算するフーリエ変換部と、
前記フーリエ変換部で計算された各データブロックについての周波数成分のうちの２つの周波数からなる周波数ペアの各々のバイコヒーレンスを計算するバイコヒーレンス計算部と、
前記周波数ペアの各々についての前記バイコヒーレンスに基づいて、前記時系列データについての検定量を計算する検定量計算部と、
前記検定量計算部で計算された前記検定量が予め定められた閾値よりも大きいか否かを判断する比較部とを有し、
前記比較部による判断結果が、解析の結果として出力され、
前記検定量計算部は、前記検定量をｚ_γで表すとき、

【数27】

（但し、
Ｑは前記周波数ペアの数、
γ_ｑは、前記周波数ペアの各々について、前記バイコヒーレンス計算部で計算されたバイコヒーレンス、
ｂｉａｓは、前記検定量の平均を０とするための調整量、
σ^２は、前記検定量の分散を補正するための調整量）
により前記検定量ｚ_γを計算するものである。

【発明の効果】

【0025】

本発明によれば、観測データの周波数間の位相的な相関関係を表す指標であるバイコヒーレンスを用い、バイコヒーレンスに基づいた検定量を用いて、観測データの非ガウス性及び非線形性を検出するので、観測データの大きさ、歪度の影響を抑制しつつ、非ガウス性又は非線形性の判定を正確に行うことができる。

【図面の簡単な説明】

【0026】

【図1】従来のデータ解析方法を実施するための装置の概略を示すブロック図である。

【図2】時系列データからのデータブロックの生成方法を示す図である。

【図3】データ解析に用いられる周波数ペアの２次元周波数平面上での配置の一例を示す図である。

【図4】本発明のデータ解析装置の概略を示すブロック図である。

【図5】図４のデータ解析装置の検出処理部の構成例を示すブロック図である。

【図6】図５の検出処理部のコヒーレンス計算部の構成例を示すブロック図である。

【図7】図５の検出処理部で用いられるパラメータを生成する装置の構成例を示すブロック図である。

【図8】（ａ）及び（ｂ）は、模擬した過渡的信号、及びその模擬信号に白色雑音を付加した波形を示す図である。

【図9】（ａ）及び（ｂ）は、図８（ｂ）に示される信号に対して、従来の方法及び本発明の方法を適用した結果を示す図である。

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

【0027】

図４は、本発明の時系列データ解析装置を示す。以下では、一例として、可聴域の音響信号を受信して、解析を行う場合について説明する。図示の時系列データ解析装置は、入力信号処理部２と、時系列データメモリ４と、検出処理部６とを有する。

【0028】

入力信号処理部２は、Ａ／Ｄ変換器３を有し、音響信号を電気信号に変換するマイクロフォン１に接続され、該マイクロフォン１から与えられるアナログ信号ｘ（ｔ）をサンプリングして、ディジタル化して、時系列データ（観測値を表す時系列データ）ｘ［ｎ］を生成する。Ａ／Ｄ変換器３のサンプリング周波数は、例えば５０ｋＨｚであり、Ａ／Ｄ変換器３から出力される時系列データはｘ［ｎ］で表される。ｎは、時刻を示すインデックスである。

【0029】

時系列データメモリ４は、入力信号処理部２から出力される時系列データｘ［ｎ］を記憶する。

【0030】

検出処理部６は、図５に示すように、データブロック生成部１２と、フーリエ変換部１４と、周波数成分メモリ１６と、周波数成分選択部１８と、バイコヒーレンス計算部２０と、検定量計算部３２と、比較部３４と、パラメータメモリ３６と、調整量計算部３８とを有する。図１と同じ符号は同様のブロックを示す。データブロック生成部１２、周波数成分メモリ１６及び周波数成分選択部１８は図１には示されていないが、従来の方法でも実際には同様のものが用いられる。

【0031】

データブロック生成部１２は、時系列データメモリ４に記憶された時系列データｘ［ｎ］から、各々複数のサンプル値で構成される、複数個即ちＰ個のデータブロックＤＢ_ｉ（ｉ＝１〜Ｐ）を生成して、フーリエ変換部１４に供給する。
データブロックＤＢ_ｉは、例えば、図２に示されるように、各々相連続するＮ個のサンプル値から成り、オーバーラップ率Ｌが５０％である。即ち、各データブロック（第ｉ番目のデータブロック）の後半のＮ／２個のサンプル値と、次のデータブロック（第ｉ＋１番目のデータブロック）の前半のＮ／２個のサンプル値とは同じものである。

【0032】

フーリエ変換部１４は、Ｐ個のデータブロックＤＢ_１〜ＤＢ_Ｐの各々のＮ個のサンプル値（観測データ）ｘ［０］〜ｘ［Ｎ−１］に対して離散フーリエ変換を行い、変換の結果として得られる周波数成分Ｘ_１［ｋ］〜Ｘ_Ｐ［ｋ］を出力する。ｋは周波数を表すインデックスである。各データブロック（ｉ番目のデータブロック）ＤＢ_ｉ（ｉは、１乃至Ｐのいずれか）についてのフーリエ変換の結果（周波数ｋの成分）をＸ_ｉ［ｋ］で表す。

【0033】

フーリエ変換の結果得られる周波数成分Ｘ_ｉ［ｋ］は、下記の式（１）で表される。式（１）は従来の方法に関して先に記載した式（１）と同じであるが、再度記載する。以下他の数式についても同様である。

【0034】

【数8】

【0035】

ここで、ｋはＭ個の周波数の値ｆ_１、ｆ_２、…ｆ_Ｍのいずれかを表すインデックスである。ｆ_１〜ｆ_Ｍを一般化してｆ_ｍ（ｍは１〜Ｍ）で表す。
周波数ｆ_ｍは、以下の値を有する。
ｍ＝１のとき、ｆ_ｍ＝０ （８Ａ）
２≦ｍ＜Ｍのとき、ｆ_ｍ＝ｆ_{（ｍ−１）}＋ΔＦ （８Ｂ）
但し、ΔＦ＝１５６２．５Ｈｚであり、これは５０ｋＨｚ（サンプリング周波数）をＮ＝３２（各データブロック中のサンプル値の数）で割った値に等しい。
Ｍは例えば２１である。

【0036】

フーリエ変換部１４で算出された、複数のデータブロックＤＢ_ｉ（ｉ＝１〜Ｐ）の各々についての周波数成分Ｘ_ｉ［ｋ］（ｋ＝ｆ_１〜ｆ_Ｍ）は、周波数成分メモリ１６に記憶される。

【0037】

周波数成分選択部１８は、周波数成分メモリ１６に記憶された複数のデータブロックＤＢ_ｉ（ｉ＝１〜Ｐ）の各々についての周波数成分Ｘ_ｉ［ｋ］（ｉ＝１〜Ｐ、ｋ＝ｆ_１〜ｆ_Ｍ）のうち、バイコヒーレンスの計算に用いられるものを逐次読み出して、バイコヒーレンス計算部２０に供給する。

【0038】

バイコヒーレンス計算部２０におけるバイコヒーレンスの計算は、複数個即ちＱ個の周波数ペアについて順に行われる。Ｑは例えば４９であり、４９個の周波数ペアの２次元周波数平面上における配置は例えば図３に「○」印で示される如くである。

【0039】

周波数ペア（ｋ_１，ｋ_２）についてバイコヒーレンスを計算する際、複数個即ちＰ個のデータブロックの各々ＤＢ_ｉについての、周波数ｋ_１、ｋ_２の周波数成分Ｘ_ｉ［ｋ_１］、Ｘ_ｉ［ｋ_２］のみならず、ｋ_１とｋ_２の和に等しい周波数ｋ_３の周波数成分Ｘ_ｉ［ｋ_３］も用いられる。そこで周波数成分選択部１８は、各データブロックＤＢ_ｉについての周波数ｋ_１、ｋ_２、ｋ_３の周波数成分Ｘ_ｉ［ｋ_１］、Ｘ_ｉ［ｋ_２］、Ｘ_ｉ［ｋ_３］をバイコヒーレンス計算部２０に供給する。
上記のように、
ｋ_３＝ｋ_１＋ｋ_２
の関係があるので、以下では、「ｋ_３」を、「ｋ_１＋ｋ_２」と表記する。
一方、任意の周波数についての説明、或いはｋ_１、ｋ_２、ｋ_３のすべてに当てはまる説明の際は、符号「ｋ」を用いる。

【0040】

バイコヒーレンス計算部２０は、各周波数ペア（ｋ_１，ｋ_２）について周波数成分選択部１８から供給された周波数成分Ｘ_ｉ［ｋ_１］、Ｘ_ｉ［ｋ_２］、Ｘ_ｉ［ｋ_１＋ｋ_２］を用いて、当該周波数ペア（ｋ_１，ｋ_２）のバイコヒーレンスγ［ｋ_１，ｋ_２］を計算する。以下に詳しく述べるように、Ｐ個のデータブロックから成る時系列データＤＳについて、Ｑ個のバイコヒーレンスの値が求められる。検定量計算部３２では、Ｑ個のバイコヒーレンスの値に基づいて時系列データＤＳについての検定量ｚ_γを求める。比較部３４では、求められた検定量ｚ_γを予め定められた閾値ｚｔｈと比較する。

【0041】

バイコヒーレンス計算部２０は、例えば図６に示すように、バイスペクトル計算部２２、パワースペクトル計算部２４、積自乗平均計算部２６、及び規格化部２８を有する。

【0042】

バイスペクトル計算部２２、パワースペクトル計算部２４、及び積自乗平均計算部２６は、周波数成分選択部１８から供給された周波数成分を受けて、それぞれバイスペクトル、パワースペクトル、及び積自乗平均を計算する。

【0043】

周波数ペア（ｋ_１，ｋ_２）のバイコヒーレンスを計算する際は、バイスペクトル計算部２２には、周波数成分Ｘ_ｉ［ｋ_１］、Ｘ_ｉ［ｋ_２］、Ｘ_ｉ［ｋ_１＋ｋ_２］（ｉ＝１〜Ｐ）が供給され、パワースペクトル計算部２４には、周波数成分Ｘ_ｉ［ｋ_１＋ｋ_２］（ｉ＝１〜Ｐ）が供給され、積自乗平均計算部２６には、周波数成分Ｘ_ｉ［ｋ_１］、Ｘ_ｉ［ｋ_２］（ｉ＝１〜Ｐ）が供給される。

【0044】

バイスペクトル計算部２２は、Ｐ個のデータブロックの各々についての、周波数ペア（ｋ_１，ｋ_２）のバイスペクトルを計算し、さらに該バイスペクトルの、Ｐ個のデータブロックについての平均を求め、該平均を、Ｐ個のデータブロックから成る時系列データＤＳについてのバイスペクトルとして出力する。即ち、まず各データブロック（ｉ番目のデータブロック）ＤＢ_ｉについてのバイスペクトルＢ_ｉ［ｋ_１，ｋ_２］を、次式（２）によって計算する。

【数9】

ここで、
Ｘ_ｉ［ｋ_１］はｉ番目のデータブロックの周波数ｋ_１の周波数成分であり、
Ｘ_ｉ［ｋ_２］はｉ番目のデータブロックの周波数ｋ_２の周波数成分であり、
Ｘ_ｉ^＊［ｋ_１＋ｋ_２］はＸ_ｉ［ｋ_１＋ｋ_２］の複素共役を表し、
Ｘ_ｉ［ｋ_１＋ｋ_２］はｉ番目のデータブロックの周波数ｋ_１＋ｋ_２の周波数成分である。

【0045】

式（２）で得られるバイスペクトルＢ_ｉ［ｋ_１，ｋ_２］は、観測データｘ［ｎ］の周波数間の相関の程度を示す指標である。

【0046】

バイスペクトル計算部２２は、式（２）の計算をＰ個のデータブロックＤＢ_ｉ（ｉ＝１〜Ｐ）の各々に対して行い、計算の結果得られたＰ個のバイスペクトルＢ_ｉ［ｋ_１，ｋ_２］の値を平均することによって、Ｐ個のデータブロックから成る時系列データＤＳについてのバイスペクトルＢａ［ｋ_１，ｋ_２］を得る。バイスペクトルＢａ［ｋ_１，ｋ_２］を求めるための計算式は、次式（３）で表される。

【数10】

【0047】

式（２）と式（３）を組合せて書き直すと、次の式（９）となる。

【数11】

【0048】

パワースペクトル計算部２４は、Ｐ個のデータブロックの各々についての、各周波数ペア（ｋ_１，ｋ_２）を構成する２つの周波数ｋ_１、ｋ_２の和に等しい周波数（ｋ_１＋ｋ_２）の周波数成分Ｘ_ｉ［ｋ_１＋ｋ_２］のパワースペクトルを計算し、さらに該パワースペクトルの、Ｐ個のデータブロックについての平均を求めることによって、Ｐ個のデータブロックから成る時系列データＤＳについてのパワースペクトルＳａ［ｋ_１＋ｋ_２］を得る。即ち、時系列データＤＳについてのパワースペクトルＳａ［ｋ_１＋ｋ_２］を求めるための計算は、次の式（１０）で表される。

【数12】

ここで、
Ｘ_ｉ［ｋ_１＋ｋ_２］はｉ番目のデータブロックの、周波数ｋ_１＋ｋ_２の周波数成分であり、
Ｘ_ｉ^＊［ｋ_１＋ｋ_２］は、Ｘ_ｉ［ｋ_１＋ｋ_２］の複素共役である。

【0049】

積自乗平均計算部２６は、Ｐ個のデータブロックの各々についての、各周波数ペア（ｋ_１，ｋ_２）を構成する２つの周波数ｋ_１、ｋ_２の周波数成分Ｘ_ｉ［ｋ_１］、Ｘ_ｉ［ｋ_２］の積の絶対値の自乗を求め、該自乗の、Ｐ個のデータブロックについての平均を求め、該平均をＰ個のデータブロックから成る時系列データＤＳについての積自乗平均Ｃａ［ｋ_１，ｋ_２］として出力する。積自乗平均Ｃａ［ｋ_１，ｋ_２］を求めるための計算は、次の式（１１）で表される。

【数13】

ここで、
Ｘ_ｉ［ｋ_１］はｉ番目のデータブロックの周波数ｋ_１の周波数成分であり、
Ｘ_ｉ［ｋ_２］はｉ番目のデータブロックの周波数ｋ_２の周波数成分である。

【0050】

規格化部２８は、バイスペクトル計算部２２で計算されたバイスペクトルＢａ［ｋ_１，ｋ_２］と、パワースペクトル計算部２４で計算されたパワースペクトルＳａ［ｋ_１＋ｋ_２］と、積自乗平均計算部２６で計算された積自乗平均Ｃａ［ｋ_１，ｋ_２］に基づいて、バイコヒーレンスγ［ｋ_１，ｋ_２］を計算する。その計算は次の式（１２）で表される。

【数14】

この計算は、バイスペクトルＢａ［ｋ_１，ｋ_２］を、積自乗平均Ｃａ［ｋ_１，ｋ_２］とパワースペクトルＳａ［ｋ_１＋ｋ_２］の積の平方根で割ることで規格化することを意味する。
規格化部２８で算出されたバイコヒーレンスγ［ｋ_１，ｋ_２］は、バイコヒーレンス計算部２０の出力となる。
式（１２）で求められるバイコヒーレンスγ［ｋ_１，ｋ_２］は、周波数ペアごとに計算され、周波数ペアを構成する２つの周波数ｋ_１及びｋ_２の成分間の位相的な相関関係を表す複素数であり、その振幅が０〜１の間になるように規格化されており、周波数ｋ_１とｋ_２の関数である。式（５）または（６）で求められる規格化バイスペクトルは、振幅の大きさが観測データの歪度に依存している。歪度は、観測データの抽出の仕方によって大きさが変化する場合があるので、結果的に検定量の大きさに影響を及ぼす。
一方、式（１２）で求められるバイコヒーレンスは、周波数ｋ_１及びｋ_２の成分間の位相的な相関関係のみを抽出し、観測データの抽出の仕方によって大きさが変化する可能性のある歪度の影響を抑制している。

【0051】

検定量計算部３２は、Ｑ個の周波数ペア（ｋ_１，ｋ_２）の各々についてバイコヒーレンス計算部２０で計算される、Ｐ個のデータブロックから成る時系列データＤＳについてのバイコヒーレンスγ［ｋ_１，ｋ_２］を、Ｑ個の周波数ペアのすべてについて積算することで、検定量ｚ_γを計算する。この計算は、算出される検定量ｚ_γの平均値及び分散が、各データブロックを構成するサンプル値の数Ｎ、時系列データＤＳを構成するデータブロックの数Ｐ、及びデータブロックのオーバーラップ率Ｌによらずに一定に保たれるようにするための調整量（下記のｂｉａｓ、σ^２）を用いて行われるものであり、次の式（１３）で表される。

【数15】

【0052】

式（１３）で、γ_ｑ（但しｑは１乃至Ｑのいずれか）は、Ｑ個の周波数ペア（ｋ_１，ｋ_２）の各々についてのγ［ｋ_１，ｋ_２］と同じものを表す。
式（１３）の分子のｂｉａｓ項は分子の平均値を０とする。従って、ｚ_γの平均値を０とするための調整項、分母のσ^２は検定量ｚ_γの分散を補正するための調整項であり、調整量計算部３８で、それぞれ以下の式（１４）、式（１５）によって求められる。

【0053】

【数16】

【0054】

式（１４）、式（１５）で、α及びβは定数であり、パラメータメモリ３６に記憶されている。
例えば、オーバーラップ率Ｌを５０％として、Ｐ＝３２個のデータブロックのＮ＝３２サンプル離散フーリエ変換を計算する場合には、α＝２、β＝６と設定すれば、観測データがガウス性を有し、且つ線形性を有する場合に、検定量ｚ_γは、平均０、分散１の正規分布に従う。
これによって検定量ｚ_γに対してある閾値を設定し、検定量ｚ_γが閾値よりも大きくなったときに、観測データが非ガウス性を有する、又は観測データが非線形性を有すると判断することが可能となる。

【0055】

比較部３４は、検定量計算部３２で計算された検定量ｚ_γが予め定められた閾値ｚｔｈよりも大きいか否かを判定し、判定結果ＤＪを出力する。即ち閾値ｚｔｈよりも大きければ、非ガウス性又は非線形性があるとの判断結果を出力し、閾値ｚｔｈ以下であれば、非ガウス性も非線形性もないとの判断結果を出力する。判定結果ＤＪは、データ解析装置による解析の結果として出力され、例えば図４に示される表示部８に解析の結果が表示される。

【0056】

上記のように、式（１３）の分子のｂｉａｓ項は分子の平均値を０とするための調整項、分母のσ^２は検定量ｚ_γの分散を補正するための調整項であり、この２つの調整項は、Ｎ（時間方向のサンプル数）及びＰ（データブロック数）を変えた場合でも、ガウス性及び線形性の信号に対する検定量ｚ_γの大きさを一定に保つためのものであり、上記の式（１４）、式（１５）によって求められる。
即ち、パラメータα及びβは、Ｐ、Ｎ、Ｌに応じて定める必要がある。

【0057】

式（１４）、式（１５）中のα、βは、以下のようにして定められる。即ち、パラメータα、βは、疑似的な白色雑音を用いることによって、実験的に求め、パラメータメモリ３６に記憶させておく。

【0058】

具体的には、疑似白色雑音に対して式（１２）でバイコヒーレンスγｗ［ｋ_１，ｋ_２］を計算し、これを用いて式（１３）と同様の計算式で計算した検定量（試験的検定量）ｚｗの平均値ｚｗａと分散ｚｗｖからパラメータα、βを求める。

【0059】

パラメータα、βの決定には、図７に示すパラメータ生成装置を用いる。このパラメータ生成装置のうち、図４、図５と同じ符号は、同じ機能を有するものであり、図４、図５の装置と共用することができる。

【0060】

（ａ） まず、白色雑音生成部４２により、疑似的な白色雑音を表す時系列データ（観測データと同じ符号ｘ［ｎ］で表す）を生成する。白色雑音の生成は、例えばソフトウェアにより行うことができる。生成された白色雑音を表す時系列データｘ［ｎ］を、図４に関連して説明した観測データ（観測値を表す時系列データ）の代わりに、時系列データメモリ４に記憶する。

【0061】

（ｂ） データブロック生成部１２により、時系列データメモリ４に記憶されている時系列データｘ［ｎ］を、各々Ｎサンプル値から成るＰ個のデータブロックを、オーバーラップ率Ｌで読み出し、フーリエ変換部１４に供給する。
フーリエ変換部１４、及びバイコヒーレンス計算部２０は、時系列データメモリ４から読み出された白色雑音を表すデータｘ［ｎ］に対して、図４を参照して説明した観測データに対する処理と同様の処理を行い、バイコヒーレンス計算部２０からは、計算されたバイコヒーレンスが出力される。出力されるバイコヒーレンスを、観測データに対する処理の場合と異なる記号γｗ［ｋ_１，ｋ_２］で表す。

【0062】

（ｃ） 検定量計算部３２は、バイコヒーレンスγｗ［ｋ_１，ｋ_２］に基づいて検定量（試験的検定量）ｚｗを計算する。この計算は次の式（１６）で表される。

【数17】

式（１６）でγｗ_ｑの添え字ｑは、それぞれ異なる周波数ペア（ｋ_１，ｋ_２）に対応するものであり、Ｑ個の周波数ペアを用いる場合、ｑは１からＱまでの値を取る。例えば、図３に示されるのと同じ配置の周波数ペアを用いる場合、Ｑは例えば４９である。

【0063】

（ｄ） 上記（ａ）〜（ｃ）の処理を複数回繰り返し、複数回、即ちＲ回の試行で得られた検定量ｚｗ_ｒ（ｒは１〜Ｒのいずれか）の平均及び分散を平均分散計算部４４に供給する。平均分散計算部４４は、複数回の試行で得られた検定量ｚｗ_ｒの平均ｚｗａ、分散ｚｗｖを計算する。この計算は、次の式（１７）、（１８）で表される。

【0064】

【数18】

【0065】

平均分散計算部４４は、このようにして求めた平均値ｚｗａ、分散ｚｗｖをパラメータ計算部４６に供給する。パラメータ計算部４６では、平均値ｚｗａ、分散ｚｗｖに基づいてパラメータα、βを計算する。この計算は下記の式（１９）、（２０）で表される。

【数19】

【0066】

ここで、Ｎ、Ｐはそれぞれ白色雑音を表す時系列データを読み出して検出処理部６に供給する際のデータブロック数及びサンプル数である。

【0067】

（ｅ） 手順（ａ）〜（ｄ）までの計算を、データブロック数Ｐ、サンプル数Ｎ及びオーバーラップ率Ｌの異なる値について行って、それぞれα、βの値を求め、パラメータメモリ３６に記憶する。例えば、Ｐ、Ｎ、Ｌに対応付けて、例えばルックアップテーブルの形式で記憶する。

【0068】

観測データに対する解析を行う際には、その時のＰ、Ｎ、Ｌに対応するα、βを読み出して、式（１４）、式（１５）の計算に用いる。

【0069】

なお、上記のように、それぞれの設定条件（Ｐ、Ｎ、Ｌ）に対して、α及びβをテーブル形式で記憶しておく代わりに、設定条件Ｐ、Ｎ、Ｌごとにｂｉａｓ、σ^２を計算してメモリに記憶させておくこととしても良い。その場合には、図５の調整量計算部３８は不要となる。

【0070】

以上のように、本発明では、従来の方法と同様に、観測データからＮサンプル×Ｐデータブロックのデータを抽出し、Ｐ個のＮサンプル離散フーリエ変換を計算する。その後、Ｐ個のデータブロックの各々のフーリエ変換の出力を用いて、パワースペクトル、バイスペクトルを計算するところまでは、従来の方法と同じである。

【0071】

バイコヒーレンスγ［ｋ_１，ｋ_２］（式（１２））と、従来の規格化バイスペクトルＢｎ［ｋ_１，ｋ_２］（式（６））との違いは、式中の分母にある。従来の規格化バイスペクトルＢｎ［ｋ_１，ｋ_２］の分母が

【数20】

であるのに対し、バイコヒーレンスγ［ｋ_１，ｋ_２］の分母は

【数21】

であり、両者の違いは、Ｓａ［ｋ_１］Ｓａ［ｋ_２］と、Ｃａ［ｋ_１，ｋ_２］の部分にある。この違いによって、バイコヒーレンスγ［ｋ_１，ｋ_２］は、式の上では信号Ｘ［ｋ_１］Ｘ［ｋ_２］と信号Ｘ［ｋ_１＋ｋ_２］間のコヒーレンスを求めているが、結果的には観測データの周波数ｋ_１及びｋ_２の成分間の位相的な相関関係を評価していることになる。

【0072】

なお、バイスペクトル自体も、周波数間の位相関係を表すが、バイスペクトルは、これとともに、観測データの振幅にも依存する。これに対して、本発明で用いるバイコヒーレンスは、式（１２）に示すように規格化を行うことで、振幅の大きさを０〜１の間になるようにして、２つの周波数成分間の位相に関する情報のみを抽出するようにしたものである。

【0073】

観測データが非ガウス性及び非線形性を有していれば、周波数ｋ_１とｋ_２の信号間に位相的なつながりがあるので、バイコヒーレンスは高い値を有することとなる。以上から、バイコヒーレンスによって観測データが非ガウス性を有し、かつ非線形性を有するか否かを確認することが可能となる。

【0074】

本発明の方法の効果を確認するため、模擬的に生成した過渡的信号を用いて実験を行った。過渡的な信号は、非ガウス性及び非線形性を有している。図８（ａ）及び（ｂ）に、模擬した過渡的信号Ｓｔｒ、及び該過渡的信号Ｓｔｒに白色雑音Ｎｗｔを付加した波形を示す。図８（ａ）は生成した過渡的信号であり、過渡的信号はある起点からの経過時間５、１０、１５、２０秒の時刻に現れている。図８（ｂ）は、生成した過渡的信号Ｓｔｒに白色雑音Ｎｗｔを付加した波形であり、この波形に対して従来の方法及び本発明の方法を適用した。

【0075】

図９（ａ）及び（ｂ）は、生成した信号（図８（ｂ）に示す信号）に対して、従来の方法及び本発明の方法を適用した結果を示す図である。図９（ａ）は、従来の方法を適用した結果を示し、図９（ｂ）は、本発明の方法を適用した結果を示す。なお、計算パラメータは以下のように設定した。
・１データブロックあたりのサンプル数Ｎ＝３２
・データブロック数Ｐ＝３２
・オーバーラップ率Ｌ＝５０％
・周波数ペア数Ｑ＝４９
・検定量ｚ_γの計算に用いられる調整量ｂｉａｓ＝０．０６２５
・検定量ｚ_γの計算に用いられる調整量σ^２＝０．１８７５
なお、サンプリング周波数は５０ｋＨｚであり、
周波数ペアの２次元平面上の配置は、図３に示す如くであり、図３に示される周波数ｆ_ｍ（ｍ＝１〜１４）はそれぞれ以下の値を有する。
ｍ＝１のとき、ｆ_ｍ＝０
１＜ｍ＜１４のとき、ｆ_ｍ＝ｆ_{（ｍ−１）}＋ΔＦ
但し、ΔＦ＝１５６２．５Ｈｚである。

【0076】

図９（ａ）及び（ｂ）の２つの処理の結果を比べると、白色雑音のみを観測しているときのレベルは同等である。一方、過渡的信号を観測しているときの応答レベルを見ると、本発明の方法は、従来の方法の概ね２倍以上の大きさである。したがって、本発明の方法は、従来の方法よりも出力ＳＮＲが高く、観測データの非ガウス性及び非線形性の検出能力が高いことが確認できる。

【0077】

以上本発明を、データ解析装置として説明したが、データ解析装置で実施されるデータ解析方法もまた本発明の一部を成す。さらに、データ解析装置における各部の処理、又はデータ解析方法における各処理をコンピュータに実行させるためのプログラム、及びそのようなプログラムを記録したコンピュータで読み取り可能な記録媒体もまた本発明の一部を成す。

【産業上の利用可能性】

【0078】

以上観測データが可聴域の音響信号である場合について説明したが、本発明はこれに限定されず、例えば、超音波信号の解析、機械振動の解析、脳波の解析にも適用できる。

【符号の説明】

【0079】

２ 入力信号処理部、 ４ 時系列データメモリ、 ６ 検出処理部、 １２ データブロック生成部、 １４ フーリエ変換部、 １６ 周波数成分メモリ、 １８ 周波数成分選択部、 ２０ バイコヒーレンス計算部、 ２２ バイスペクトル計算部、 ２４ パワースペクトル計算部、 ２６ 積自乗平均計算部、 ２８ 規格化部、 ３２ 検定量計算部、 ３４ 比較部、 ３６ パラメータメモリ、 ３８ 調整量計算部、 ４２ 白色雑音生成部、 ４４ 平均分散計算部、 ４６ パラメータ計算部。

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】