JP7040668B2

JP7040668B2 - 異常検出装置、方法、プログラム

Info

Publication number: JP7040668B2
Application number: JP2021511011A
Authority: JP
Inventors: ムルトゥザペトラードワラー; 充仙洞田; 裕子太田; 亮太鈴木
Original assignee: NEC Corp
Current assignee: NEC Corp
Priority date: 2018-08-30
Filing date: 2018-08-30
Publication date: 2022-03-23
Anticipated expiration: 2038-08-30
Also published as: WO2020044495A1; JP2021535393A; US20210247427A1

Description

本発明は、異常検出を行う装置、方法、プログラムに関する。

教師付き学習では、ラベル付けされたデータセット、すなわち、入力の値と期待される出力の値を準備し、入力と出力の間の関係に関する学習モデルを作成する。新しい入力が与えられると、学習したモデルを用いることで出力を予測することができる。教師なし学習では、入力は提供されるが出力は提供されず、入力データセットはグループ又はクラスタに分類（cluster）される。新しい入力が与えられると、クラスタモデルを用いることで出力を予測することが可能である。

クラスタリング手法では、入力データに最も適合するクラスタの数を探索して特徴量をそれぞれのクラスタに分類することが好ましい。例えば、ドメインの知識や測定データに関する事前情報を用いて、最適なクラスタ数を見つけるための試行が行われる。

測定装置で測定された電気機器の消費電流などの電気信号の波形データにより、電気機器の正常動作と異常動作を学習するためには、２組のクラスタ（すなわち、正常クラスタと異常クラスタ）の中心が互いに十分に離れていることが好ましい。特徴抽出は、ドメイン知識又は次元削減技術に基づいて行われる。主成分分析（Principal Component Analysis, ＰＣＡ）、すなわち線形変換法が一般的な特徴抽出法として利用される。また、オートエンコーダ法、すなわち非線形変換法も、特徴抽出に利用されている。オートエンコーダ法は、入力よりも少ない数の隠れユニットをネットワークに配置するなどの制約をネットワークに与えることで、ネットワークへの入力が出力と同じとなるニューラルネットワークである。オートエンコーダ法は、ラベルを必要とせず、入力をコンパクトに符号化した表現を学習する。

信号処理に基づく異常検出について、特許文献１には、可動要素を有する装置の異常診断方法として、装置の状態が正常か否（異常）かを容易に把握可能とした方法が開示されている。装置の振動に関する生波形を正規化して正規化波形を求める。正規化された波形のデータを高速フーリエ変換（FFT）して周波数成分のデータ、すなわちパワースペクトルを求める。正常な条件下での装置の振動に関するデータも同様の方法で事前に求めておく。診断におけるパワースペクトルから正常な条件におけるパワースペクトルを差し引いた残りのパワースペクトルに基づいて装置の異常診断が行われる。

特開２００１－３１８０３１号公報

教師なし学習手法では、特徴変換を行っても正常電流波形と異常電流波形から判別可能な特徴を生成できない場合がある。教師なし学習手法では、教師あり分類(supervised classification)法や半教師付き１クラス（semi-supervised one-class）法に比べて、決定境界（decision boundary）が不明であったり、十分に分離されていなかったりする場合がある。特徴変換に基づく教師なしクラスタリング方式では、決定境界を推定するためのラベルや情報が入手できないかあるいは知られていないため、正常波形信号と異常波形信号の良好な分離の生成に失敗する可能性がある。

したがって、本発明の目的は、正常波形信号と異常波形信号を良好に分離したクラスタを作成できる装置、方法、プログラムを提供することにある。

本発明の１つの側面によれば、電気機器の周期的な電気波形信号の複数セットの波形データであって、個々の波形データの長さが、前記周期的な電気波形信号の１サイクル時間に設定された、複数セットの波形データを取得する入力部と、前記複数の個々の波形データセットのそれぞれについて周波数スペクトルを計算する周波数スペクトル処理部と、前記複数の個々の波形データセットの周波数スペクトルのそれぞれから、前記１サイクル時間の逆数である基本周波数の奇数次及び偶数次高調波を抽出して分離し、奇数次及び偶数次の高調波行列を作成する分離部と、前記奇数次及び偶数次の高調波行列に対して正準相関分析（Canonical correlation Analysis、ＣＣＡ）を行い、ＣＣＡ特徴量を取得する正準相関分析（ＣＣＡ）部と、前記ＣＣＡ特徴量に対して線形変換を行い、線形変換された特徴量を得る線形変換部と、前記線形変換された特徴量に基づいてモデルを生成するモデル生成部と、前記複数セットの波形データの前記複数の周波数スペクトルをマグニチュード量子化して、正常な波形信号と異常な波形信号とを識別する識別部と、少なくとも識別結果を出力する出力部と、を備えた異常検出装置が提供される。

本発明の１つの側面によれば、電気機器の周期的な電気波形信号の複数セットの波形データであって、個々の波形データの長さが、前記周期的な電気波形信号の１サイクル時間に設定された、複数セットの波形データを取得し、
前記複数の個々の波形データセットのそれぞれについて周波数スペクトルを計算し、
前記複数の個々の波形データセットの周波数スペクトルのそれぞれから、前記１サイクル時間の逆数である基本周波数の奇数次及び偶数次高調波を抽出して分離し、奇数次及び偶数次の高調波行列を作成し、
前記奇数次及び偶数次の高調波行列に対して正準相関分析（Canonical correlation Analysis、ＣＣＡ）を行い、ＣＣＡ特徴量を取得し、
前記ＣＣＡ特徴量に対して線形変換を行い、線形変換された特徴量を取得し、
前記線形変換された特徴量に基づいてモデルを生成し、
前記複数セットの波形データの前記複数の周波数スペクトルをマグニチュード量子化して、正常な波形信号と異常な波形信号とを識別し、
少なくとも識別結果を出力する、異常検出方法が提供される。

本発明の１つの側面によれば、電気機器の周期的な電気波形信号の複数セットの波形データであって、個々の波形データの長さが、前記周期的な電気波形信号の１サイクル時間に設定された、複数セットの波形データを取得し、
前記複数の個々の波形データセットのそれぞれについて周波数スペクトルを計算し、
前記複数の個々の波形データセットの周波数スペクトルのそれぞれから、前記１サイクル時間の逆数である基本周波数の奇数次及び偶数次高調波を抽出して分離し、奇数次及び偶数次の高調波行列を作成し、
前記奇数次及び偶数次の高調波行列に対して正準相関分析（Canonical correlation Analysis、ＣＣＡ）を行い、ＣＣＡ特徴量を取得し、
前記ＣＣＡ特徴量に対して線形変換を行い、線形変換された特徴量を取得し、
前記線形変換された特徴量に基づいてモデルを生成し、
前記複数セットの波形データの前記複数の周波数スペクトルをマグニチュード量子化して、正常な波形信号と異常な波形信号とを識別し、
少なくとも前記識別結果を出力する処理をコンピュータに実行させるプログラムが提供される。本発明によれば、該プログラムはこれを記録した記録媒体で提供してもよい。

記録媒体は、半導体メモリ（ＲＡＭ（Random Access Memory）、ＲＯＭ（Read Only Memory）、ＥＥＰＲＯＭ（Electrically Erasable and Programmable Read Only Memory）、フラッシュメモリ等）、ハードディスクドライブ（ＨＤＤ）、ソリッドステートドライブ（ＳＳＤ）、コンパクトディスク、デジタルバーサタイルディスク等の非一時的（non-transitory）コンピュータ可読記録媒体であってもよい。

本発明によれば、正常波形信号と異常波形信号との良好に分離されたクラスタを作成することが可能となる。本発明の他の利点及び特徴は、以下の説明、図面及び特許請求の範囲により明らかであろう。

図１は、第１の実施形態の動作を説明するフローチャートである。図２は、第１の実施形態の構成を示す図である。図３は、第１の実施形態の動作を説明するフローチャートである。図４は、第１の実施形態の動作を説明するフローチャートである。図５は、図４のＲＭＳ量子化の構成を示す図である。図６は、第１の実施形態を説明するための図である。図７は、第１の実施形態を説明するための図である。図８は、第１の実施形態を説明するための図である。図９は、第２の実施形態の構成を説明するための図である。図１０（Ａ）及び図１０（Ｂ）は、実験結果（精度及びＮＭＩ）を示す。

以下に本発明の実施形態について述べる。図１は、本発明の一実施形態による異常検出方法を説明するフローチャートである。この方法は、プロセッサによって実行可能なコンピュータプログラム命令及びプロセッサによって処理するためのデータを記憶するためのメモリに結合されたプロセッサによって好ましくは実行されてもよい。ここで、プロセッサによって実行される、メモリに記憶されたプログラム命令は、図面を参照して後述される異常検出装置を構成する。

＜ステップＳ１０１＞
異常検出装置は、電気機器の交流（Alternating-Current, ＡＣ）電流信号（消費電流信号など）の複数セットの波形データを取得する。個々の波形データセットは、交流電流信号（交流電源）の１サイクル時間に設定されている。異常検出装置は、測定装置から交流電流信号の複数セットの波形データを受信するように構成されていてもよい。あるいは、異常検出装置は、交流電流信号の複数サイクルの波形データを受信し、波形データを複数セットの波形データにスライスし、各波形データは交流電流信号の１サイクル時間に設定された長さを有するように構成されていてもよい。

＜ステップＳ１０２＞
異常検出装置は、個々の波形データセットのそれぞれの周波数スペクトルを算出して処理する。異常検出装置は、個々の波形データセットの周波数スペクトルの集合を、行数＝個々の波形データセットの数、列数＝周波数スペクトルの周波数成分（ビン）の数の行列（matrix）（周波数スペクトル行列）の形式で格納してもよい。

＜ステップＳ１０３＞
異常検出装置は、交流電流信号の周期（１サイクル時間）の逆数である基本周波数（基本周波数は交流電源周波数に等しい）の奇数次高調波と偶数次高調波を、周波数スペクトル行列の複数の周波数スペクトルのそれぞれから分離し、奇数次高調波行列と偶数次高調波行列とを作成する。奇数次高調波行列及び偶数次高調波行列の行数は波形データセットの数であり、奇数次高調波行列及び偶数次高調波行列の列数はそれぞれ、各波形データセットの周波数スペクトルにおける奇数次高調波の数及び偶数次高調波の数に設定される。

＜ステップＳ１０４＞
異常検出装置は、奇数次及び偶数次高調波行列の特徴量(FFT特徴量)に対して正準相関分析(Canonical Correlation Analysis, ＣＣＡ)を行う。ＣＣＡは奇数次高調波行列と偶数次高調波行列を特徴空間に変換し、変換された各徴次元は、２つの行列間の対応する次元（dimensions）と相関性が高い。

＜ステップＳ１０５＞
異常検出装置は、ＣＣＡ特徴量に対して線形変換を行い、線形変換された特徴量を得る。

＜ステップＳ１０６＞
異常検出装置は、線形変換された特徴量に基づいて、クラスタモデルなどの教師なし学習モデルを生成する。

＜ステップＳ１０７＞
異常検出装置は、個々の波形データセットの周波数スペクトルの集合（sets）（周波数スペクトル行列）の高調波に対して、反復的にマグニチュード量子化（算出された振幅（magnitude）を２ビンヒストグラムとして符号化）を行い、正常波形信号と異常波形信号（異常）を識別する。

＜ステップＳ１０８＞
異常検出装置は、少なくとも識別された結果を出力する。

本実施形態によれば、正常波形信号と異常波形信号のよく分離されたクラスタ（複数）を作成することができる。

関連技術では、どのクラスタが正常で、どのクラスタが異常であるかを自動的に識別することは困難である。本実施形態によれば、信号解析に教師なし手法を適用することで、クラスタが正常か異常かを自動的に識別することが可能となる。

図２は、第１の実施形態による異常検出装置１００の構成を例示する図である。図２を参照すると、異常検出装置１００は、入力部１０１と、記憶部１０２と、出力部１０３と、周波数スペクトル処理部１０４と、分離部１０５と、正準相関分析（ＣＣＡ）部１０６と、線形変換部１０７と、モデル生成部１０８と、識別部１０９とを含む。

入力部１０１は、測定装置（図示せず）から、電気機器の交流電流信号について、個々の波形データの各セットの長さを交流電流信号の周期（１サイクル時間）（すなわち、交流電源の１サイクル時間）に設定された複数セットの波形データ（サンプリングされたデジタルデータ）を受信する。測定装置としては、スマートタップ、ＣＴ（Current Transformer）等の電流センサ等を利用してもよい（ただし、これらに制限されない）。測定装置（図示せず）は、個々の波形データセットの各々は、位相が、交流電圧信号のサイクル開始（ゼロクロス）に整合された複数セットの波形データを入力部１０１に供給してもよい。入力部１０１は、測定装置（図示せず）とネットワークを介して通信するためのネットワーク・インタフェース・カードを含んでもよい。異常検出装置１００は、測定装置（図示せず）に、交流電流信号の複数セットの波形データを異常検出装置１００に送信するように指示するように構成されていてもよい。入力部１０１は、複数セットの波形データを記憶部１０２に記憶する。

測定装置（図示せず）は、交流電流信号の波形データをスライス又はセグメント化して、個々の波形データセットの各々が交流電流信号の１周期（１サイクル時間）（すなわち、交流電源サイクル）に設定された長さを有し、電気機器に供給される交流電圧信号のサイクルの開始点（ゼロクロス）に整合された位相を有する複数セットの波形データとし、複数セットの波形データを入力部１０１に送信するように構成されていてもよい。あるいは、入力部１０１は、電気機器を流れる交流電流信号の複数サイクルの波形データと、電気機器に供給される交流電圧信号の複数サイクルの波形データとを取得してもよく、交流電圧信号のサンプリング（時間と周波数）は、交流電流信号のサンプリングと同じ（同期している）ように構成されている。入力部１０１は、交流電流信号の複数サイクルの波形データを、個々の波形データセットの各々が交流電流信号の１周期（１サイクル時間）（すなわち、交流電圧信号の１サイクル時間）に設定された長さを有し、交流電流信号のサイクルの開始点（ゼロクロス）に整合された位相を有する複数セットの波形データにスライスしてもよい。

周波数スペクトル処理部１０４は、記憶部１０２から複数セットの波形データを読み出す。周波数スペクトル処理部１０４は、ＦＦＴ（Fast Fourier Transform）を用いて、個々の波形データセットのそれぞれの周波数スペクトル（フレームともいう）を算出し、個々の波形データセットのそれぞれの周波数スペクトルの正規化やフィルタリングなどの処理を、個別に、すなわちフレーム単位で行う。周波数スペクトル処理部１０４は、複数のフレームの周波数処理結果を、行数＝個別波形データセットの数、列数＝周波数スペクトルに含まれる周波数成分の数の周波数スペクトル行列の形式で記憶部１０２に格納する。

分離部１０５は、複数の個別波形データセットの複数の周波数スペクトル（フレーム）の各々から、交流電流信号の基本周波数（＝交流電源周波数＝５０Ｈｚ又は６０Ｈｚ）の奇数次及び偶数次の高調波を分離して、行数＝個別波形データセット（フレーム）の数、列数＝周波数スペクトルの奇数次及び偶数次の高調波の数を有する奇数次及び偶数次周波数の高調波行列を作成する。交流電流信号の１周期（１サイクル時間）に設定された長さの個別波形データセットの各々は、周期的な交流電流信号を１サイクルごとにスライスして得られ、個別波形データセットに対してＦＦＴを行うことにより、各個別波形データセットの周波数スペクトルにおける周波数成分は、基本周波数（周期の逆数）の高調波であることから、関連情報を失うことなく、奇数次高調波と偶数次高調波の分離を行うことができる。測定装置（図示せず）がモータ等の回転機械を有する電気機器の交流電流信号を測定する場合には、高調波は、周波数スペクトルにピークを有する。電気機器が異常状態にある場合、個々の波形データの周波数スペクトルに周波数成分の余分な（extra）ピークが観測されることがある。

正準相関分析部１０６は、２つのデータセット、すなわち奇数次高調波行列と偶数次高調波行列との間で正準相関分析(canonical correlation analysis)を行い、奇数次高調波行列と偶数次高調波行列の潜在特徴ベクトル(latent feature vectors)間の相関を最大化する射影方向(projective directions)を求める。

線形変換部１０７は、正準相関分析を受けた特徴量（ＣＣＡ変換された特徴量という）に対して線形変換を行い、変換された特徴量を得る。変換された特徴量の次元は削減されている。

モデル生成部１０８は、線形変換部１０７により次元削減された変換特徴量に基づいてモデルを生成する。また、ＣＣＡ変換された特徴量は、奇数次高調波のみ、偶数次高調波のみ、奇数次高調波と偶数次高調波との組み合わせでモデル化することも可能である。

識別部１０９は、正常波形信号と異常波形信号を識別するために、周波数スペクトルマトリクスの高調波に対して反復的にマグニチュード量子化（ここでは、マグニチュードが２ビンヒストグラムとして符号化されている）を実行する。

出力部１０３は、少なくとも識別された結果又は異常検出処理を出力してもよい。出力部１０３は、表示装置、プリンタ、データを記憶装置に記憶するデータ記憶装置、通信部、又はそれらの任意の組み合わせを含んでいてもよい。

図３は、周波数スペクトル処理部１０４の動作を説明するフローチャートである。入力部１０１が、電気機器の交流電流信号の波形データをＭセット（所定の整数）受信したとする（ステップＳ２０１）。

周波数スペクトル処理部１０４は、Ｍセットの波形データ（時間領域データ）から、Ｍセットの波形データにそれぞれＦＦＴ（Fast Fourie Transform：高速フーリエ変換）を適用して、Ｍセットの周波数スペクトル（スペクトラ）を算出する（ステップＳ２０２）。なお、以下では、簡略化と説明の容易さのために、ＤＦＴ（Discrete Fourier Transform：離散フーリエ変換）のみを用いている。

個々の波形データ（１フレーム）が x = [ x(0), x(T), x(2*T), ..., x((N-1)*T) ] で構成されているものとする。Ｔはサンプリング間隔である。x(n*T)をx(n)と表記すると、波形データのＤＦＴ（Discrete Fourier Transform)は、以下で与えられる。

… (1)
ここで、X(k)はk番目の周波数ビン(周波数: k/(T*N))における複素DFT係数である。

… (2)

ここで、R(k)はDFT係数X(k)の実部、I(k)は虚部をそれぞれ表し、j² =-1である。

ｋ番目の周波数ビンにおける振幅（amplitude）（大きさ： magnitude）と位相は次のように与えられる。

…(3)

… (4)

Ｍセットの波形データを時間順にＸ_１、Ｘ_２、Ｘ_Ｍと表記すると、ｉ番目の波形データ（ｉ番目のフレーム）に対応するｉ番目の周波数スペクトルは以下で与えられる。

…(5)

Ｍ掛けるＮの行列（M*N行列、又は（M, N）の行列と表記する）F（周波数スペクトル行列とも称する）は、以下で与えられる。

… (6)

周波数スペクトル行列Ｆの各行は、個別に正規化されるのが好ましい（ステップＳ２０３）。正規化の一例として、周波数スペクトル行列Ｆの行の各要素を正規化係数ＮＦ_ｉ（ｉ＝１，...，Ｍ）で除算して行列ＦＮを生成する。

… (7)

特に制限されないが、ｉ番目の正規化係数ＮＦ_ｉ（ｉ＝１，...,Ｍ）は、例えば、周波数スペクトル行列Ｆのｉ番目の行の要素の振幅の和として定義される。

… (8)

周波数スペクトル行列Ｆの各行の要素の値を０と１の間の範囲に設定する行毎の正規化は、行列ＦＮの各列の波形データの時間領域の特徴（characteristics）又は特性（properties）を変更するために、実行されるのが好ましい。

周波数スペクトル処理部１０４は、行列ＦＮの各行にハイパスフィルタリングを施す（ステップＳ２０４）。ハイパスフィルタリングは、基本周波数（交流電源周波数ｆ₀：５０／６０Ｈｚ）以下の直流成分及びその第２高調波未満の周波数成分をカットして、基本周波数の高調波（第２、第３,・・・）を抽出するために用いられる。ハイパスフィルタは、周波数領域において、ハイパスフィルタの周波数スペクトルと行列ＦＮの各行（正規化された周波数スペクトル）とを乗算して行ってもよい。

周波数スペクトル処理部１０４は、ハイパスフィルタリングを受けた行列ＦＮの各行について対数変換（例えば、１０を底とする対数）を行って行列ＦＬを作成してもよい（ステップＳ２０５）。行列ＦＬの各要素は、振幅（magnitude）スペクトルのデシベル表現（decibel representation）に対応する。

… (9)
ここで、(1/NF₁)X₁'、...、(1/NF_M)X_M'は、第１～Ｍ番目のフレーム（波形データセット）のハイパスフィルタリングされた正規化周波数スペクトルである。

時間領域の波形データ（フレーム）は、以下のようにフーリエ展開で表現（近似）することができる。

… (10)
Ａ_kは振幅、
ｆ_k =k/(T*N) (k=0,...,N-1)は周波数、ただしTはサンプリング間隔、
φ_kは位相である。

電気機器において、電流波形データの基本周波数は交流電源の周波数（50Hz/60Hz）である。基本周波数は、電気機器の負荷に応じて振幅（amplitude）（magnitude：振幅スペクトル）が変化することがある。一つの波形は、上記式に示すように、サンプリング周波数ｆ_s (=1/T)のＮ個の周波数成分の正弦波信号で表すことができる。サンプリング周波数は、信号に含まれる最高周波数ｆ_ｃの少なくとも２倍のｆ_ｓ > 2*ｆ_ｃ(ここでは、ｆ_ｃ＝ｆ_Ｎ-１)であることが好ましい。振幅、周波数、又は位相は、異常の影響を受けやすい変数である。この変数の変化を検出することで、現在の波形データの異常を識別することができる。

周波数スペクトルは、振幅変数と周波数変数の情報を含むものとする。また、Ｍフレームの波形データについて基本周波数の位相が一定であるとしてもよい。

位相整合は、周波数スペクトルを解析する前に実行される。各フレームの電流信号の波形データは、交流電圧に関して一定の位相関係、すなわち交流電圧サイクルの開始点又はゼロクロスに対する一定の位相遅れを有するようにサンプリングされてもよい。

奇数・偶数分離部１０５は、周波数スペクトル行列ＦＬの基本周波数（交流電源周波数：５０Ｈｚ／６０Ｈｚ）の奇数次高調波と偶数次高調波を分離して、奇数次高調波行列Ｏと偶数次高調波行列Ｅを作成する。

本実施形態では、特に限定されないが、サンプリング周波数ｆ_ｓは、交流電源周波数ｆ_０の倍数（ｆ_ｓ＝ｈ*ｆ_０：ｈは正の整数）とする。例えば、サンプリング周波数ｆ_ｓが１ＫＨｚ、交流電源周波数ｆ_０が５０Ｈｚの場合、サンプリング周波数ｆ_ｓは交流電源周波数ｆ_０の２０倍となる。ｈ＝ｆ_ｓ／ｆ_０とし、第１ビン（ビン０）が直流成分となるように周波数成分が配置されているとすると、Ｍ*Ｎの行列ＦＬの各行の周波数スペクトルの奇数次高調波は、周波数ビン（列）：｛２＊（ｊ-１）＋１｝＊ｈ＋１（ｊ＝１，２，...）に存在し、行列ＦＬの各行の周波数スペクトルの偶数次高調波は、周波数ビン（列）：2 * j * h + 1（ｊ＝１，２，...）に存在する。ここで、最大高調波周波数はナイキスト周波数ｆ_ｓ／２以下である。

行列ＦＬ、奇数次高調波行列Ｏ、偶数次高調波行列Ｅのｉ番目の行、ｊ番目の列の要素をFL(i,j)、O(i,j)、E(i,j)で表すと
O(i,j)＝FL(i, {2*(j-1)+1}*h+1)，
E(i,j)= FL(i, 2*j*h+1)
(i=1,...,M, j=1,....,j_max, ここで (2*j_max*h+1)*f₀<=f_s/2
… (11)

正準相関分析部１０６は、奇数次高調波行列Ｏと偶数次高調波行列Ｅとの間で正準相関分析（ＣＣＡ）を行う。ＣＣＡは、２つの特徴量セット（集合）の特徴ベクトル間の相関を最大とする射影方向（projective directions）を求める。

n対の特徴ベクトルを持つデータセットについて、
{(a_i, b_i) i=1,...,n}, ここで a_i, b_i は p 次元ベクトル、q 次元ベクトルである。

ただし

… (12)

ＣＣＡは次のように表現できる。

… (13)
ここで argmax_ga,gbは g_a ^TAB^Tg_bの関数値が最大となる射影方向g_a、g_bを与える。上付き文字 T は転置演算子を表す。

上記の条件を満たす一対の射影方向（g_ａ１，g_ｂ１）を求めた後、次のような最適化問題を解くことで、第２の対の射影方向を求めることができる。

… (14)

以上の処理を繰り返すことで、対になっている固有値問題を解くことができる。

… (15)

固有ベクトル（g_a1，g_b1）(i=1,…,l)、すなわち、最大のl個の固有値に対応するものは、ＣＣＡの射影方向の対である。

… (16)
行列Ａ、ＢからＣＣＡで抽出した特徴量セットを合成する。

正準相関分析部１０６は、行列Ｏと行列Ｅを行列Ａと行列Ｂとし、行列Ｏと行列Ｅの２つの特徴量セットの特徴ベクトル間の相関が最大となる射影方向を求めることにより、上記処理を算出する。

… (17)
この結果、行列ＯとＥの次元が削減される。

本実施形態では、好ましくは、行列Ｏと行列Ｅの２つの特徴量セットの特徴ベクトル間の相関が最大となるように射影方向を求めるＣＣＡが用いられる。しかし、本実施形態は、ＣＣＡに限定されるものではない。例えば、ＰＣＡ（主成分分析）を用いてもよい。

Ｙをm*pの行列とした場合、共分散行列は次のように計算できる。

… (18)

ＰＣＡは、共分散行列Cの固有ベクトルと共分散行列Cの固有値を求めるために、共分散行列Cの固有値分解を行う。

… (19)

Uが直交行列ならば

Σはm個の対角要素の共分散行列の固有値λを降順に持ち、他の要素は値0を持つ対角行列である。

… (20)

ｍ個の固有値のうち、最大の固有値がｌ個の固有値に対応する最初のｌ番目のＰＣＡ成分を選択することで、次元が削減される。ＰＣＡを用いる場合、式（１８）の行列Ｙを行列Ｏとし、行列Ｅとしてもよい。またＳＶＤ（Singular Value Decomposition, 特異値分解）を用いてもよい。

線形変換部１０７は、正準相関分析により抽出された特徴ベクトルに対して、平均、ＰＣＡ等の線形変換を行い、変換された特徴量を得る。より具体的には、ＣＣＡ特徴量、すなわち行列Ｏ及び行列Ｅに対して線形変換を行う。例えば、行列Ｏの第１列と行列Ｅの第１列を連結して、サイズ（Ｍ，２）（Ｍ行、２列）の新しい行列を作成してもよい（Ｍはフレーム数）。サイズ（Ｍ，２）の新しい行列にＰＣＡのような線形変換が適用され、サイズ（Ｍ，１）の変換された行列（すなわち、Ｍ行の列ベクトル）が得られる。ＰＣＡの代わりに、平均のような任意の線形変換を使用して、サイズ（Ｍ，２）の行列を（Ｍ，１）に変換するために適用することができる。同様に、行列Ｏと行列Ｅの各列を結合し、それらに線形変換を適用してもよい。線形変換部１０７は、ＣＣＡ変換された特徴量の空間をそのまま保持する。ＣＣＡ特徴量に対して線形変換を適用することで、変換された特徴量を変換してＣＣＡの元の空間に戻すことができる。

奇数行列や偶数行列の特徴ベクトルの次元削減には，高次元空間から低次元空間への写像を行うＰＣＡを用いてもよい。

モデル生成部１０８は、線形変換部１０７が出力する変換された特徴量に基づいて、クラスタモデルなどのモデルを生成する。モデル生成部１０８は、ラベルが利用可能であれば、当然のことながら、教師付き学習モデルを生成してもよい。
また、モデル生成部１０８は、ラベルが利用可能であれば、半教師付き学習モデルを生成してもよいことは勿論である。

線形変換にＰＣＡを用いた場合、最初の３つの主成分を選択すると、モデルの特徴空間が３次元空間になる。

識別部１０９は、Ｍセットの波形データ（Ｍフレーム）にＦＦＴを適用して得られた周波数スペクトル行列Ｆ（式（６）で定義される）、又は行列ＦＮ（式（７）で定義される）を受け取り、正規化処理を行う。識別部１０９は、反復的なマグニチュード量子化処理を行い、正常波形信号と異常波形信号を識別する。図４は、識別部１０９の動作を説明するフローチャートである。図４において、ステップＳ３０１が図３のステップＳ２０１に対応し、ステップＳ３０２がステップＳ２０２、すなわちＳ２０２～Ｓ２０４に対応する。なお、周波数スペクトル処理部１０４は、算出された周波数スペクトル行列Ｆ又はＦＮを記憶部１０２に記憶し、識別部１０９は、記憶部１０２から周波数スペクトル行列Ｆ又はＦＮを受信（読み出し）するように構成してもよい。ここで、周波数スペクトル行列Ｆ又はＦＮを行列Ｆと表記する。

forループのループ変数mを1に初期化し、リストビン＿カウント（bin_count）をヌル（null）に初期化する（ステップＳ３０３）。

forループの本体（Ｓ３０４～Ｓ３０７）が実行される。すなわち、m=1から周波数スペクトル行列Ｆの周波数成分の数（交流電源周波数の高調波の数）まで、ステップＳ３０５～Ｓ３０７（forループの本体）が繰り返される。

＜ステップＳ３０５＞
引数にmを設定して、サブルーチン：マグニチュード量子化（magnitude_quantization）（ＲＭＳ量子化（RMS_quantization）ともいう）を呼び出す。

＜ステップＳ３０６＞
マグニチュード量子化(m)のサブルーチンで返された2つのビンのうち、最大のカウント値（又は大きい方のカウント値）をビン＿カウント（bin_count）に付加する。

＜ステップＳ３０７＞
mを1増加させる。

識別部１０９は、繰り返し（m=1から周波数成分の最大数（周波数スペクトルの高調波の数））ごとに、最大のカウントをビン＿カウント（bin_count）に付加する。ここで、ビン＿カウント（bin_count）は、それぞれの周波数成分で得られた最大のカウントから構成される線形リストであり、リストの長さは、周波数成分の数（すなわち、ナイキスト周波数未満の高調波周波数成分の数）である。また、本実施形態では、ビン＿カウント（bin_count）は、周波数成分の昇順に配列された最大カウント値のリストで構成されている（ただし、これに制限されない）。

図５は、図４のマグニチュード（ＲＭＳ）量子化（サブルーチン）の動作を説明するフローチャートである。

＜ステップＳ４０１＞
周波数スペクトル行列Ｆ（Ｍ行の各行は１フレームの振幅スペクトルを持つ）を入力する。入力された周波数スペクトル行列Ｆのコピーは、元の周波数スペクトル行列Ｆの内容が振幅（ＲＭＳ）量子化（サブルーチン）によって変更されないように、マグニチュード（ＲＭＳ）量子化（サブルーチン）によって局所的に書き込み／読み出しアクセスが可能なメモリに行列Ｆ'として格納されていてもよい。

＜ステップＳ４０２＞
周波数スペクトル行列Ｆ'のＭ行の各行にハイパスフィルタリングを適用してm個の周波数成分（ｍ番目以降の周波数成分）を抽出する。値mはサブルーチン：マグニチュード量子化（magnitude quantization）へ入力引数（周波数成分の数）：mとして引き渡される。なお、ハイパスフィルタリングの代わりに、周波数スペクトル行列F'の各行にバンドパスフィルタリングを適用してもよい。

＜ステップＳ４０３＞
ハイパスフィルタリングされた周波数スペクトル行列Ｆ'のＭ行のそれぞれに逆ＦＦＴを適用して、ハイパスフィルタリングされた周波数スペクトル行列Ｆ'の各行に対応する各時間領域波形データのＭセットの時間領域波形データを得る。

＜ステップＳ４０４＞
ハイパスフィルタリングされた周波数スペクトル行列Ｆ'に逆ＦＦＴを適用して得られた時間領域波形データのＭセットごとのＲＭＳ２乗（RMS square）を算出する。

＜ステップＳ４０５＞２ビンヒストグラムを作成する。
Ｍセットの時間領域の波形データについて、Ｍ個のＲＭＳ値を計算したＭ次元のＲＭＳベクトルを作成する。次いで、このＲＭＳベクトルに対して２つのビン数のヒストグラムを計算する。すなわち、２ビン数のヒストグラムでは、横軸は、下位レンジ（lower_range）と上位レンジ（upper range）からなる２つのビンに分割したＲＭＳ値であり、縦軸はヒストグラムカウント（各ビンにおけるカウント）である。ＲＭＳベクトルの値は２つのビンのうちのいずれかに分類(categorized)され、対応するビンのカウントがインクリメントされる。マグニチュード（ＲＭＳ）量子化サブルーチンは、各ビンのカウント値（ヒストグラムカウント）に加えて、ビンレンジ（bin_range）を戻り値として返しても良い。

図４を参照して、ステップＳ３０４では、ｍが周波数成分の数より大である場合には、forループを終了してステップＳ３０８にジャンプする処理が実行される。識別部１０９は、以下の処理を実行する。

＜ステップＳ３０８＞
ビンカウント（bin_count）が最小値をとる周波数成分m_minを取得する。

＜ステップＳ３０９＞
m_minを引数に設定してサブルーチン：マグニチュード量子化（magnitude_quantization）を呼び出す。

＜ステップＳ３１０＞
サブルーチン：マグニチュード量子化（magnitude_quantization(m_min)）で得られた値をラベル：ｑに変換する。ＲＭＳベクトル内のＲＭＳ値がどのビンレンジ（bin_range）に該当するかを、ビンレンジ（bin_range）から算出する。次に、低いビンレンジを０、高いビンレンジを１としてタグ付けを行う。これは、低いビンレンジの値が正常、高いビンレンジの値が異常インスタンス（anomaly instance）であるとの仮定に基づいている。ステップＳ３１０では、各インスタンスにラベルを０、１としてタグ付けする。０は正常、１は異常である。これらのタグ又はラベルは、ＲＭＳ量子化ラベルと命名される。

＜ステップＳ３１１＞
予測ラベルを取得し、ｐに設定する。例えば、クラスタリングが適用される場合、各インスタンスはいずれかのクラスタに属する。クラスタが２個の場合、各インスタンスはクラスタ番号０又はクラスタ番号１のいずれかに属すると予測される。これらのクラスタ番号を予測ラベルと呼ぶ。

＜ステップＳ３１２＞
量子化されたラベルベクトルと予測されたラベルベクトルのワンホット（one hot）エンコーディングラベルをＰとＱに設定する。ワンホットエンコーディングとは、カテゴリ変数を、予測のために機械学習アルゴリズムに提供できる形に変換する処理である。２つのクラスタの場合、ＱとＰはＭ個の行と２つの列を持ち、Ｍはサンプル数である。

量子化されたラベルベクトルと予測されたラベルベクトルのワンホット（one_hot）ラベルをq = [1 1 0 0 1] と p = [1 1 0 0 1]とする（ここで、 M=5 ）。

＜ステップＳ３１３＞
行列Ｄ＝ドット積（Ｑ，Ｐ）
ＱとＰは、サイズ (Ｎ，２) の行列であり、ここでＮは標本（samples）数，2 は列の数（0 と 1 のラベル）である。行列 Q と P のドット積（内積）をとることは、行列 Q と P の共分散を計算することになる。ドット積の結果である行列 D は、行が量子化ラベル、列が予測ラベルである 2×2 の行列である。

＜ステップＳ３１４＞
正常＝行列Ｄの１行目が最大値、異常＝行列Ｄの２行目が最大値。
本実施形態では、予測されたクラスタが正常であるか異常であるかを識別することは、量子化されたラベルが０の場合は正常であり、量子化されたラベルが１の場合は異常であるという仮定に基づいて行われる。

行列Ｄの各行で、最大値を持つ値は、予測されたラベル（列内）が正常か異常かを識別する。
行列Ｄを

とする。

量子化されたラベルが０のものを正常とし、予測ラベルの名前（正常又は異常）を求めるために、行列Ｄの１行目及び２行目の最大値を名前として選択する。例えば、行列Ｄがグレースケールで示されている図６に示すように、行列Ｄの１行目では、最大値が９０であり、これは予測ラベル１が正常クラスタであることを意味する。行列Ｄの２行目では、最大値は８０であり、これは予測されたラベル０が異常クラスタであることを意味する。このようにして、量子化されたラベルに基づいて、クラスタが正常クラスタであるか異常クラスタであるか（どのクラスタが正常で、どのクラスタが異常であるか）を識別することが可能となる。

図７及び図８は、本実施形態のシーケンスの一例を模式的に示す図であり、図２の出力部１０３が出力する実験結果（画像データ）の一部を示す図である。

時系列波形信号（２０１）は、電気機器の消費電流のサンプル波形データである。画像３０１は、波形の一例を示す図である。時系列信号をＭセット（それぞれＮ個のサンプルを有する）の波形データにスライスし、Ｍセットの波形データのそれぞれを時間順に一列に並べることで、２次元行列Ｘ（２０２）を構成する。電流信号をスライスした複数セットの波形データ（３０２）の各々は、電気機器に供給される交流電圧信号の１周期の開始点（ゼロクロス）に位相合わせされている。

複数セットのスライス波形データは、ＦＦＴ又はＤＦＴ（Ｎが２のべき乗の場合、ＦＦＴを行うようにしてもよい）により周波数領域に変換され（２０３）、Ｍ×Ｎの周波数スペクトル行列Ｆが得られる（３０３）。ＦＦＴにより得られた周波数スペクトル（２０４）は、データからＭセットｆ波の周波数スペクトルの各周波数ビンにおける強度（magnitude）（信号強度：intensity）が、輝度、すなわち明るさ、に比例して出力（表示）された画像データ３０４として表示される。

画像データ３０５は、周波数スペクトルの行毎の正規化（２０５）の結果（表示装置に表示される画像データ）である。画像データ３０６は、ハイパスフィルタリング（２０６）の結果である。

画像データ３０７は、ハイパスフィルタリング（２０７）後の周波数スペクトル行列の対数をとった結果である

画像データ３０８Ｅと３０８Ｏは偶数次高調波と奇数次高調波とを分離した結果である（３０８）。

画像データ３０９Ｅと３０９ＯはＣＣＡの結果である（２０９）。

画像データ３１０は、第１の２次元の線形変換（２１０）の結果である。

画像データ３１１は、モデル（クラスタモデル）（２１１）を用いたクラスタリング結果である。画像データ３１１では、３次元軸上に２クラスのクラスタリングが表示される。すなわち、線形変換部１０７により特徴空間の次元が３に縮小された３次元特徴空間上で、２クラス（正常・異常）クラスタリングが行われる。

画像データ３１２は、正常／異常の識別（２１２）を行う識別部１０９が出力する正常波形信号と異常波形信号の各々について２ビンヒストグラムを例示したものである。画像データ３１２において、横軸はＲＭＳ（root mean square: 二乗平均平方根）値、縦軸は各ビンのカウント（ヒストグラムカウント）である。異常セクション（anomaly section）は、ＲＭＳ値の小さい範囲よりもＲＭＳ値の大きい範囲でより大きなカウントを有し、正常セクション(normal section)は、ＲＭＳ値の大きい範囲よりもＲＭＳ値の小さい範囲でより大きなカウントを有する。

画像データ３１３は、識別部１０９による波形データの予測結果の出力である。画像データ３１３に示すように、画像データ３１１のクラスタリング結果におけるクラスタリングモデルから、予測ラベルに基づく異常クラスタと正常クラスタが入れ替わっている。

上記の実施形態で説明した異常検出装置１００（又はシステム）は、例えば図９に示すように、サーバシステム（又はクラウドシステム）のようなコンピュータシステム上で実装されてもよい。図９を参照すると、サーバシステムなどのコンピュータシステム４００は、プロセッサ（中央処理装置）４０１と、例えば半導体メモリ（例えば、ランダムアクセスメモリ（ＲＡＭ）、リードオンリーメモリ（ＲＯＭ）、電気的に消去可能でプログラム可能なＲＯＭ（ＥＥＰＲＯＭ）などを含むことができるメモリ４０２、及び／又は、ハードディスクドライブ（ＨＤＤ）、コンパクトディスク（ＣＤ）、デジタルバーサタイルディスク（ＤＶＤ）等の少なくとも１つを含むメモリ４０２、入出力部（表示端末）４０３、記憶部４０４、通信部４０５とを含む。

コンピュータシステム４００は、１つ又は複数のネットワークインターフェースコントローラ（カード）（ＮＩＣ）を含んでもよい通信部４０５を介して、ＬＡＮ（ＬｏｃａｌＡｒｅａＮｅｔｗｏｒｋ）及び／又はＷＡＮ（ＷｉｄｅＡｒｅａＮｅｔｗｏｒｋ）などのネットワーク４０６に接続することができる。図２の異常検出装置１００の処理を実行するためのプログラム（命令及びデータ）が記憶部４０３に格納されており、プロセッサ４０１は、記憶部４０３から、メモリ４０２に設けられたメインメモリにプログラムを読み出して、実施形態に従った異常検知機能を実現するためのプログラムを実行する。

業務用コーヒーマシンの正常・異常電流波形（電流シグネチャ（Electric Current Signature, ＥＣＳ）とも称する）を測定するために、制御実験を行った。

検証のため、コーヒーマシンが正常な状態であった時間を記載している。製造されたコーヒーの味を変えるために、コーヒー豆の量を 60%に減らした。特に制限されないが、異常とは、モータによる粉砕工程では消費電力が変化してＥＣＳが変化するため、コーヒーマシンのモータの動作が正常状態から変化した場合と定義される。正常状態と異常状態のそれぞれ１６回のコーヒーをいれるシーケンスを測定し、モータによる粉砕サブシーケンスを抽出する。

アプローチ：
主成分分析（ＰＣＡ）、オートエンコーダ（ＡＥ）、及びＣＣＡの３つの方法が比較され、図１０のＣＣＡは、上述した本実施形態の方法を表している。

主成分分析(ＰＣＡ)やオートエンコーダ（ＡＥ）方法の入力は、結果比較のためにＣＣＡの入力との類似が維持された前処理された周波数特徴量である。

k-meansにおけるランダムな開始点の問題を避けるために、同じ選択された特徴量を用いたクラスタリングを１０回行い、平均結果を２個のクラスタ数だけ再コード化する。

トレーニング段階では正常状態と異常状態の各１０個のカップを選択し、テスト段階ではトレーニング段階では図示しない残りの６個のカップを選択する。

結果：
性能を測定するために、クラスタリング精度と正規化相互情報（Normalized Mutual Information, ＮＭＩ）を含む2つの指標（metrics）を採用した。どちらの指標も値が大きいほど性能が良いことを示している。正規化相互情報（ＮＭＩ）は、相互情報（Mutual Information, ＭＩ）スコアを正規化したもので、0（相互情報なし）から1（完全相関）の間で結果をスケーリングする。すなわち、ＮＭＩ値は、0.0から1.0までの範囲であり、値が高いほどクラスタリング性能が良いことを意味する。精度は、クラスタリング結果のうち正しいと思われる部分の割合を測定したもので、0(=0%)から1(=100%)までの範囲であり、値が高いほどクラスタリング性能が良いことを示している。図１０は、実験結果を示しており、図１０（Ａ）、（Ｂ）はクラスタリング精度（Accuracy）と正規化相互情報（ＮＭＩ）を示しており、横軸は次元性（dimensionality）を示している。

図１０に示した性能は、テストデータに基づいて計算されたものである。ＮＭＩはＣＣＡの値が高いことを示している。これは、他の方法と比べて、両方のクラスタが均質（homogeneous）であり、完全(complete)であり、特徴空間のオーバーラップが少ないことを意味している。

なお、上記の特許文献１の開示を、本書に引用をもって繰り込むものとする。本発明の全開示（請求の範囲を含む）の枠内において、さらにその基本的技術思想に基づいて、実施の形態ないし実施例の変更・調整が可能である。また、本発明の請求の範囲の枠内において種々の開示要素（各請求項の各要素、各実施例の各要素、各図面の各要素等を含む）の多様な組み合わせ乃至選択が可能である。すなわち、本発明は、請求の範囲を含む全開示、技術的思想にしたがって当業者であればなし得るであろう各種変形、修正を含むことは勿論である。

１００異常検出装置
１０１入力部
１０２記憶部
１０３出力部
１０４周波数スペクトル処理部
１０５分離部
１０６ＣＣＡ（正準相関分析）部
１０７線形変換部
１０８モデル生成部
１０９識別部
２０１－２１３処理
３０１－３１３表示イメージ
４００コンピュータシステム
４０１プロセッサ
４０２メモリ
４０３記憶部
４０４入出力部（表示端末）
４０５通信部
４０６ネットワーク

Claims

電気機器の周期的な電気波形信号の複数セットの波形データであって、個々の波形データの長さが、前記周期的な電気波形信号の１サイクル時間に設定された、複数セットの波形データを取得する入力部と、
前記複数の個々の波形データセットのそれぞれについて周波数スペクトルを計算する周波数スペクトル処理部と、
前記複数の個々の波形データセットの周波数スペクトルのそれぞれから、前記１サイクル時間の逆数である基本周波数の奇数次及び偶数次高調波を抽出して分離し、奇数次及び偶数次の高調波行列を作成する分離部と、
前記奇数次及び偶数次の高調波行列に対して正準相関分析（Canonical correlation Analysis、ＣＣＡ）を行い、ＣＣＡ特徴量を取得する正準相関分析（ＣＣＡ）部と、
前記ＣＣＡ特徴量に対して線形変換を行い、線形変換された特徴量を得る線形変換部と、
前記線形変換された特徴量に基づいてモデルを生成するモデル生成部と、
前記複数セットの波形データの前記複数の周波数スペクトルをマグニチュード量子化して、正常な波形信号と異常な波形信号とを識別する識別部と、
少なくとも識別結果を出力する出力部と、
を備える、異常検出装置。
前記周期的な電気波形信号は、交流（Alternating Current：ＡＣ）電流信号であり、前記交流電流信号の前記複数セットの波形データはそれぞれ交流電圧信号に位相整合されている、ことを特徴とする、請求項１記載の異常検出装置。
前記モデル生成部は、
奇数次高調波のみ、
偶数次高調波のみ、または、
奇数次高調波と偶数次高調波の任意の組み合わせを用いて、前記ＣＣＡ特徴量のモデル化を行う、ことを特徴とする、請求項１又は請求項２記載の異常検出装置。
前記周波数スペクトル処理部は、前記複数セットの波形データのそれぞれの周波数ビンの各々について振幅スペクトルを計算する、請求項１乃至３のいずれか一に記載の異常検出装置。
前記周波数スペクトル処理部は、前記複数セットの波形データのそれぞれについて、前記周波数スペクトルの正規化処理を個別に行う、請求項４記載の異常検出装置。
前記周波数スペクトル処理部は、前記複数セットの波形データのそれぞれの前記正規化された周波数スペクトルにハイパスフィルタリングを適用して前記基本周波数以下の周波数成分を遮断する、請求項５記載の異常検出装置。
前記線形変換部は、前記奇数次の高調波行列と前記偶数次の高調波行列の線形結合を生成する、請求項１乃至６のいずれか一に記載の異常検出装置。
前記識別部は、前記マグニチュード量子化として、前記複数セットの波形データの前記周波数スペクトルにおける周波数成分の数だけ反復的に、ＲＭＳ（Root Mean Square: 二乗平均平方根）量子化を行いＲＭＳ値に基づいて２つのビンのヒストグラムを作成し、
前記２つのビンのうち、より大きなヒストグラムカウント値で構成されたリストまたはベクトルを作成する、請求項１乃至７のいずれか一に記載の異常検出装置。
電気機器の周期的な電気波形信号の複数セットの波形データであって、個々の波形データの長さが、前記周期的な電気波形信号の１サイクル時間に設定された、複数セットの波形データを取得し、
前記複数の個々の波形データセットのそれぞれについて周波数スペクトルを計算し、
前記複数の個々の波形データセットの周波数スペクトルのそれぞれから、前記１サイクル時間の逆数である基本周波数の奇数次及び偶数次高調波を抽出して分離し、奇数次及び偶数次の高調波行列を作成し、
前記奇数次及び偶数次の高調波行列に対して正準相関分析（Canonical correlation Analysis、ＣＣＡ）を行い、ＣＣＡ特徴量を取得し、
前記ＣＣＡ特徴量に対して線形変換を行い、線形変換された特徴量を取得し、
前記線形変換された特徴量に基づいてモデルを生成し、
前記複数セットの波形データの前記複数の周波数スペクトルをマグニチュード量子化して、正常な波形信号と異常な波形信号とを識別し、
少なくとも識別結果を出力する、異常検出方法。
電気機器の周期的な電気波形信号の複数セットの波形データであって、個々の波形データの長さが、前記周期的な電気波形信号の１サイクル時間に設定された、複数セットの波形データを取得し、
前記複数の個々の波形データセットのそれぞれについて周波数スペクトルを計算し、
前記複数の個々の波形データセットの周波数スペクトルのそれぞれから、前記１サイクル時間の逆数である基本周波数の奇数次及び偶数次高調波を抽出して分離し、奇数次及び偶数次の高調波行列を作成し、
前記奇数次及び偶数次の高調波行列に対して正準相関分析（Canonical correlation Analysis、ＣＣＡ）を行い、ＣＣＡ特徴量を取得し、
前記ＣＣＡ特徴量に対して線形変換を行い、線形変換された特徴量を取得し、
前記線形変換された特徴量に基づいてモデルを生成し、
前記複数セットの波形データの前記複数の周波数スペクトルをマグニチュード量子化して、正常な波形信号と異常な波形信号とを識別し、
少なくとも識別結果を出力する、
処理をコンピュータに実行させるプログラム。