JP6567242B2 - 無線機識別装置及び無線機識別方法 - Google Patents

Description

この発明は、無線機を識別する無線機識別装置及び無線機識別方法に関するものである。

無線機から送信される無線信号の立ち上がりなどの過渡応答は、個体差があることが知られており、近年では、無線機から送信される無線信号を解析することで、無線機を識別する無線機識別装置が提案されている。
例えば、無線機識別装置は、無線機から送信される無線信号の立ち上がりの特徴量を算出し、算出した特徴量に基づいて、無線機を識別する。
無線信号の立ち上がりの特徴量としては、立ち上がり時の無線信号における振幅等の分散値、歪度及び尖度などがある。

ただし、信号対雑音比 （Ｓｉｇｎａｌ−ｔｏ−Ｎｏｉｓｅ Ｒａｔｉｏ：ＳＮＲ）が低い環境下では、立ち上がりの特徴量の算出精度が低下して、無線機の識別精度が低下する。このため、無線機から送信される無線信号を短時間フーリエ変換することで、ＳＮＲを改善してから、立ち上がりの特徴量を算出する無線機識別装置が以下の特許文献１に開示されている。
この無線機識別装置では、立ち上がりの特徴量として、立ち上がりの時定数を算出している。

特開２００２−２６８２６号公報

従来の無線機識別装置は、立ち上がりの時定数に基づく無線機の識別方法を採用している。立ち上がりの時定数に基づく無線機の識別方法では、無線機の種類については識別することが可能である。しかし、同一種類の無線機は、異なる無線機であっても、立ち上がりの時定数の差異が小さいため、無線機の個体を識別することができないという課題があった。

この発明は上記のような課題を解決するためになされたもので、無線機の個体識別を行うことができる無線機識別装置及び無線機識別方法を得ることを目的とする。

この発明に係る無線機識別装置は、識別対象の無線機から送信された無線信号を受信する信号受信部と、信号受信部により受信された無線信号を短時間フーリエ変換し、無線信号の短時間フーリエ変換結果を示すフーリエ変換信号を出力するフーリエ変換部と、フーリエ変換部から出力されたフーリエ変換信号のスペクトログラムにおけるエネルギーの時間波形を算出し、エネルギーの時間波形から無線信号の立ち上がり時刻を検出する立ち上がり検出部と、フーリエ変換部から出力されたフーリエ変換信号の中から、立ち上がり検出部により検出された立ち上がり時刻を含む時間帯のフーリエ変換信号を抽出し、抽出したフーリエ変換信号に含まれている設定周波数成分の複素時間信号から、抽出したフーリエ変換信号の変動を特徴量として算出する特徴量算出部とを設け、無線機識別部が、特徴量算出部により算出された特徴量に基づいて、無線機を識別するようにしたものである。

この発明によれば、フーリエ変換部から出力されたフーリエ変換信号のスペクトログラムにおけるエネルギーの時間波形を算出し、エネルギーの時間波形から無線信号の立ち上がり時刻を検出する立ち上がり検出部と、フーリエ変換部から出力されたフーリエ変換信号の中から、立ち上がり検出部により検出された立ち上がり時刻を含む時間帯のフーリエ変換信号を抽出し、抽出したフーリエ変換信号に含まれている設定周波数成分の複素時間信号から、抽出したフーリエ変換信号の変動を特徴量として算出する特徴量算出部とを設け、無線機識別部が、特徴量算出部により算出された特徴量に基づいて、無線機を識別するように構成したので、無線機の個体を識別することができる効果がある。

この発明の実施の形態１による無線機識別装置を示す構成図である。 この発明の実施の形態１による無線機識別装置を示すハードウェア構成図である。 無線機識別装置における信号受信部１を除く構成要素がソフトウェア又はファームウェアなどで実現される場合のコンピュータのハードウェア構成図である。 この発明の実施の形態１による無線機識別方法における学習処理を示すフローチャートである。 フーリエ変換部２による学習用無線信号ｓ_ｇ（ｎ）及び識別用無線信号ｓ_ｄ（ｎ）のＳＴＦＴを示す説明図である。 立ち上がり検出部３による立ち上がり時刻Ｔ_０の検出処理を示す説明図である。 特徴量算出部７による特徴量の算出処理を示す説明図である。 図８Ａは、特徴量算出用のパラメータが適正でないために、種類（Ａ）の無線機に係る特徴量ベクトル〇と、種類（Ｂ）の無線機に係る特徴量ベクトル◇と、種類（Ｃ）の無線機に係る特徴量ベクトル△とが一部重なっている例を示す説明図、図８Ｂは、特徴量算出用のパラメータが適正であるために、種類（Ａ）の無線機に係る特徴量ベクトル〇と、種類（Ｂ）の無線機に係る特徴量ベクトル◇と、種類（Ｃ）の無線機に係る特徴量ベクトル△とが重なっていない例を示す説明図である。 第１のデータベース部９に格納されているシリアルナンバーが付加されている学習用無線信号を示す説明図である。 第２のデータベース部１０に格納される複数の既知の無線機に係る特徴量の一例を示す説明図である。 この発明の実施の形態１による無線機識別方法における識別処理を示すフローチャートである。 無線機識別部１１による無線機の識別処理を示す説明図である。

以下、この発明をより詳細に説明するために、この発明を実施するための形態について、添付の図面に従って説明する。

実施の形態１．
図１は、この発明の実施の形態１による無線機識別装置を示す構成図である。
図２は、この発明の実施の形態１による無線機識別装置を示すハードウェア構成図である。
図１及び図２において、信号受信部１は、例えば図２に示す信号受信回路２１で実現される。
信号受信部１は、既知の無線機から送信された無線信号を学習用無線信号として受信し、受信した学習用無線信号をアナログ信号からデジタル信号に変換することで、デジタルの学習用無線信号をフーリエ変換部２及び第１のデータベース部９に出力する。
また、信号受信部１は、識別対象の無線機から送信された無線信号を識別用無線信号として受信し、受信した識別用無線信号をアナログ信号からデジタル信号に変換することで、デジタルの識別用無線信号をフーリエ変換部２に出力する。

フーリエ変換部２は、例えば図２に示すフーリエ変換回路２２で実現される。
フーリエ変換部２は、信号受信部１から出力されたデジタルの学習用無線信号を短時間フーリエ変換（ＳＴＦＴ：Ｓｈｏｒｔ Ｔｉｍｅ Ｆｏｕｒｉｅｒ Ｔｒａｎｓｆｏｒｍ）し、学習用無線信号の短時間フーリエ変換結果であるＳＴＦＴ結果を示すフーリエ変換信号（以下、学習用フーリエ変換信号と称する）を立ち上がり検出部３及び特徴量算出部７に出力する処理を実施する。
また、フーリエ変換部２は、信号受信部１から出力されたデジタルの識別用無線信号をＳＴＦＴし、識別用無線信号のＳＴＦＴ結果を示すフーリエ変換信号（以下、識別用フーリエ変換信号と称する）を立ち上がり検出部３及び特徴量算出部７に出力する処理を実施する。
ＳＴＦＴは、一定時間毎に異なる時刻の無線信号をそれぞれ切り出し、切り出した各々の無線信号をそれぞれ高速フーリエ変換（ＦＦＴ：Ｆａｓｔ Ｆｏｕｒｉｅｒ Ｔｒａｎｓｆｏｒｍ）する手法である。

立ち上がり検出部３は、例えば図２に示す立ち上がり検出回路２３で実現される。
立ち上がり検出部３は、フーリエ変換部２から出力された学習用フーリエ変換信号のスペクトログラムにおけるエネルギーの時間波形（以下、学習用のエネルギー時間波形と称する）を算出し、学習用のエネルギー時間波形から学習用無線信号における立ち上がり時刻を検出する処理を実施する。
また、立ち上がり検出部３は、フーリエ変換部２から出力された識別用フーリエ変換信号のスペクトログラムにおけるエネルギーの時間波形（以下、識別用のエネルギー時間波形と称する）を算出し、識別用のエネルギー時間波形から識別用無線信号における立ち上がり時刻を検出する処理を実施する。

パラメータ設定部４は、第１のパラメータ設定部５及び第２のパラメータ設定部６を備えており、例えば図２に示すパラメータ設定回路２４で実現される。
第１のパラメータ設定部５は、特徴量算出用のパラメータとして、立ち上がり検出部３により検出された立ち上がり時刻を含む時間帯の長さである時間範囲Δｔを示すパラメータと、設定周波数成分である設定周波数ビンを示すパラメータとを設定する処理を実施する。
第２のパラメータ設定部６は、特徴量算出部７により算出される第１から第９の特徴量の重みを示す重み付けパラメータを設定する処理を実施する。
また、第２のパラメータ設定部６は、フーリエ変換部２によりＳＴＦＴが行われる際に用いられる窓関数のサンプル数及び窓関数をシフトさせるサンプル数を設定する処理を実施する。

特徴量算出部７は、例えば図２に示す特徴量算出回路２５で実現される。
特徴量算出部７は、パラメータ設定部４により設定された特徴量算出用のパラメータから、立ち上がり時刻を含む時間帯及び設定周波数ビンを認識する。
特徴量算出部７は、フーリエ変換部２から出力された学習用フーリエ変換信号の中から、立ち上がり検出部３により検出された学習用無線信号における立ち上がり時刻を含む時間帯の学習用フーリエ変換信号を抽出する処理を実施する。
そして、特徴量算出部７は、抽出した学習用フーリエ変換信号に含まれている設定周波数ビンの複素時間信号から、抽出した学習用フーリエ変換信号の変動を特徴量として算出する処理を実施する。

また、特徴量算出部７は、フーリエ変換部２から出力された識別用フーリエ変換信号の中から、立ち上がり検出部３により算出された識別用無線信号における立ち上がり時刻を含む時間帯の識別用フーリエ変換信号を抽出する処理を実施する。
そして、特徴量算出部７は、抽出した識別用フーリエ変換信号に含まれている設定周波数成分の複素時間信号から、抽出した識別用フーリエ変換信号の変動を特徴量として算出する処理を実施する。

データベース部８は、第１のデータベース部９及び第２のデータベース部１０を備えており、例えば図２に示すデータベース回路２６で実現される。
第１のデータベース部９は、信号受信部１から出力されたデジタルの学習用無線信号を格納する。
第２のデータベース部１０は、既知の無線機に係る特徴量として、特徴量算出部７により算出された学習用フーリエ変換信号の変動を格納する。

無線機識別部１１は、例えば図２に示す無線機識別回路２７で実現される。
無線機識別部１１は、特徴量算出部７により算出された識別対象の無線機に係る特徴量と、第２のデータベース部１０により格納されている既知の無線機に係る特徴量とを比較する処理を実施する。
無線機識別部１１は、識別対象の無線機に係る特徴量と既知の無線機に係る特徴量との比較結果に基づいて、識別対象の無線機を識別する処理を実施する。

図１では、無線機識別装置の構成要素である信号受信部１、フーリエ変換部２、立ち上がり検出部３、パラメータ設定部４、特徴量算出部７、データベース部８及び無線機識別部１１のそれぞれが、図２に示すような専用のハードウェアで実現されるものを想定している。即ち、信号受信回路２１、フーリエ変換回路２２、立ち上がり検出回路２３、パラメータ設定回路２４、特徴量算出回路２５、データベース回路２６及び無線機識別回路２７で実現されるものを想定している。

ここで、データベース回路２６は、例えば、ＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）、ＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）、フラッシュメモリ、ＥＰＲＯＭ（Ｅｒａｓａｂｌｅ Ｐｒｏｇｒａｍｍａｂｌｅ Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）、ＥＥＰＲＯＭ（Ｅｌｅｃｔｒｉｃａｌｌｙ Ｅｒａｓａｂｌｅ Ｐｒｏｇｒａｍｍａｂｌｅ Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒy）などの不揮発性又は揮発性の半導体メモリ、磁気ディスク、フレキシブルディスク、光ディスク、コンパクトディスク、ミニディスク、あるいは、ＤＶＤ（Ｄｉｇｉｔａｌ Ｖｅｒｓａｔｉｌｅ Ｄｉｓｃ）が該当する。
また、信号受信回路２１、フーリエ変換回路２２、立ち上がり検出回路２３、パラメータ設定回路２４、特徴量算出回路２５及び無線機識別回路２７は、例えば、単一回路、複合回路、プログラム化したプロセッサ、並列プログラム化したプロセッサ、ＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ Ｓｐｅｃｉｆｉｃ Ｉｎｔｅｇｒａｔｅｄ Ｃｉｒｃｕｉｔ）、ＦＰＧＡ（Ｆｉｅｌｄ−Ｐｒｏｇｒａｍｍａｂｌｅ Ｇａｔｅ Ａｒｒａｙ）、または、これらを組み合わせたものが該当する。

無線機識別装置における信号受信部１を除く構成要素は、専用のハードウェアで実現されるものに限るものではなく、無線機識別装置における信号受信部１を除く構成要素が、ソフトウェア、ファームウェア、または、ソフトウェアとファームウェアとの組み合わせで実現されるものであってもよい。
ソフトウェア又はファームウェアはプログラムとして、コンピュータのメモリに格納される。コンピュータは、プログラムを実行するハードウェアを意味し、例えば、ＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）、中央処理装置、処理装置、演算装置、マイクロプロセッサ、マイクロコンピュータ、プロセッサ、あるいは、ＤＳＰ（Ｄｉｇｉｔａｌ Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓｏｒ）が該当する。

図３は、無線機識別装置における信号受信部１を除く構成要素がソフトウェア又はファームウェアなどで実現される場合のコンピュータのハードウェア構成図である。
無線機識別装置における信号受信部１を除く構成要素がソフトウェア又はファームウェアなどで実現される場合、データベース部８をコンピュータのメモリ３１上に構成するとともに、フーリエ変換部２、立ち上がり検出部３、パラメータ設定部４、特徴量算出部７及び無線機識別部１１の処理手順をコンピュータに実行させるためのプログラムをメモリ３１に格納し、コンピュータのプロセッサ３２がメモリ３１に格納されているプログラムを実行するようにすればよい。

また、図２では、無線機識別装置の構成要素のそれぞれが専用のハードウェアで実現される例を示し、図３では、無線機識別装置がソフトウェアやファームウェアなどで実現される例を示しているが、無線機識別装置における一部の構成要素が専用のハードウェアで実現され、残りの構成要素がソフトウェアやファームウェアなどで実現されるものであってもよい。

次に動作について説明する。
＜学習時の処理＞
最初に、既知の無線機から送信された無線信号を学習する際の処理内容を説明する。
図４は、この発明の実施の形態１による無線機識別方法における学習処理を示すフローチャートである。
信号受信部１は、既知の無線機から送信された無線信号を学習用無線信号として受信し、受信した学習用無線信号をアナログ信号からデジタル信号に変換する（図４のステップＳＴ１）。
学習用無線信号には、既知の無線機を識別する個体情報であるシリアルナンバーが付加されているものとする。
信号受信部１は、デジタルの学習用無線信号ｓ_ｇ（ｎ）をフーリエ変換部２及び第１のデータベース部９に出力する。ｓ_ｇ（ｎ）におけるｎは、学習用無線信号のサンプル番号である。例えば、図９には、サンプル番号ｎが１〜４の学習用無線信号が示されている。
第１のデータベース部９には、シリアルナンバーが付加されている学習用無線信号ｓ_ｇ（ｎ）が格納される。
図９は、第１のデータベース部９に格納されているシリアルナンバーが付加されている学習用無線信号を示す説明図である。
図９の例では、シリアルナンバーが同じ無線機であっても、受信時のＳＮＲ又はサンプリング周波数などの受信条件が異なっている複数の学習用無線信号が格納されている。

フーリエ変換部２は、以下の式（１）に示すように、信号受信部１から出力されたデジタルの学習用無線信号ｓ_ｇ（ｎ）をＳＴＦＴし、学習用無線信号ｓ_ｇ（ｎ）のＳＴＦＴ結果を示すフーリエ変換信号である学習用フーリエ変換信号Ｓ_ｇ（ｍ，ｋ）を立ち上がり検出部３及び特徴量算出部７に出力する（図４のステップＳＴ２）。

式（１）において、ｗ（ｎ−ｍ）は、窓関数であり、例えば、矩形窓又はハミング窓が考えられる。
ｍは、窓関数ｗ（ｎ−ｍ）をシフトさせるサンプル数、ｍ_０は、窓関数ｗ（ｎ−ｍ）の適用区間の初めのサンプル番号、Ｍは、窓関数ｗ（ｎ−ｍ）のサンプル数、ｋは、周波数ビン番号である。
なお、窓関数ｗ（ｎ−ｍ）のサンプル数Ｍ及び窓関数ｗ（ｎ−ｍ）をシフトさせるサンプル数ｍのそれぞれは、第２のパラメータ設定部６によって事前に設定される。

ここで、図５は、フーリエ変換部２による学習用無線信号ｓ_ｇ（ｎ）及び識別用無線信号ｓ_ｄ（ｎ）のＳＴＦＴを示す説明図である。
学習用無線信号ｓ_ｇ（ｎ）及び識別用無線信号ｓ_ｄ（ｎ）は、時間と振幅との関係を示す２次元信号である。識別用無線信号ｓ_ｄ（ｎ）については後述する。
ＳＴＦＴ結果を示す学習用フーリエ変換信号Ｓ_ｇ（ｍ，ｋ）及び識別用フーリエ変換信号Ｓ_ｄ（ｍ，ｋ）は、時間と周波数及び電力との関係を示す３次元信号である。識別用フーリエ変換信号Ｓ_ｄ（ｍ，ｋ）については後述する。
フーリエ変換部２により学習用無線信号ｓ_ｇ（ｎ）及び識別用無線信号ｓ_ｄ（ｎ）のそれぞれがＳＴＦＴされることで、学習用無線信号ｓ_ｇ（ｎ）及び識別用無線信号ｓ_ｄ（ｎ）のそれぞれが、コヒーレント積分されるので、ＳＮＲが改善される効果が得られる。

立ち上がり検出部３は、学習用無線信号ｓ_ｇ（ｎ）の立ち上がり時刻Ｔ_０を検出する処理を実施する（図４のステップＳＴ３）。
図６は、立ち上がり検出部３による立ち上がり時刻Ｔ_０の検出処理を示す説明図である。
以下、立ち上がり検出部３による立ち上がり時刻Ｔ_０の検出処理を具体的に説明する。

まず、立ち上がり検出部３は、以下の式（２）に示すように、フーリエ変換部２から出力された学習用フーリエ変換信号Ｓ_ｇ（ｍ，ｋ）の２乗値｜Ｓ_ｇ（ｍ，ｋ）｜^２をスペクトログラムＳＰ_ｇ（ｍ，ｋ）として算出する。

次に、立ち上がり検出部３は、以下の式（３）に示すように、スペクトログラムＳＰ_ｇ（ｍ，ｋ）におけるエネルギーの時間波形である学習用のエネルギー時間波形Ｅ_ｇ（ｍ）を算出する。
学習用のエネルギー時間波形Ｅ_ｇ（ｍ）は、スペクトログラムＳＰ_ｇ（ｍ，ｋ）に含まれている複数の周波数における電力のうち、同一時刻の電力同士がそれぞれ合計されることで算出される各々の時刻における電力の合計値である。

式（３）において、ｋ_１は、スペクトログラムＳＰ_ｇ（ｍ，ｋ）に含まれている複数の周波数のうち、電力を合計する周波数の下限周波数を示す周波数ビン番号、ｋ_ｎは、スペクトログラムＳＰ_ｇ（ｍ，ｋ）に含まれている複数の周波数のうち、電力を合計する周波数の上限周波数を示す周波数ビン番号である。
ｋ＝ｋ_１，・・・，ｋ_ｎである。

次に、立ち上がり検出部３は、学習用のエネルギー時間波形Ｅ_ｇ（ｍ）と閾値Ｅ_ｔｈとを比較する。
閾値Ｅ_ｔｈとしては、例えば、学習用のエネルギー時間波形Ｅ_ｇ（ｍ）の最大値から設定値Ｅ_０を減算した値が考えられる。
立ち上がり検出部３は、学習用のエネルギー時間波形Ｅ_ｇ（ｍ）と閾値Ｅ_ｔｈとの比較結果を参照することで、最初は、ノイズレベルのエネルギーを示していた学習用のエネルギー時間波形Ｅ_ｇ（ｍ）が、閾値Ｅ_ｔｈに到達したか否かを判定する。
立ち上がり検出部３は、学習用のエネルギー時間波形Ｅ_ｇ（ｍ）が閾値Ｅ_ｔｈに到達した時刻を、学習用無線信号ｓ_ｇ（ｎ）の立ち上がり時刻Ｔ_０として検出する。
この実施の形態１では、立ち上がり検出部３が、学習用無線信号ｓ_ｇ（ｎ）の立ち上がり時刻Ｔ_０を検出する際、スペクトログラムＳＰ_ｇ（ｍ，ｋ）に含まれている複数の周波数における電力のうち、同一時刻の電力同士をそれぞれ合計するようにしている。これにより、立ち上がり時に電力が大きく変動する周波数が含まれている場合でも、当該電力の変動の影響が低減されるため、立ち上がり時刻Ｔ_０の検出精度が向上する。

特徴量算出部７は、第１のパラメータ設定部５により設定された特徴量算出用のパラメータから、立ち上がり時刻Ｔ_０を含む時間帯（ｔ_１〜ｔ_２）及び設定周波数ビンｋ_ｂを認識する。
第１のパラメータ設定部５による特徴量算出用のパラメータの設定処理については後述する。
次に、特徴量算出部７は、図６に示すように、フーリエ変換部２から出力された学習用フーリエ変換信号Ｓ_ｇ（ｍ，ｋ）の中から、立ち上がり検出部３により検出された学習用無線信号ｓ_ｇ（ｎ）の立ち上がり時刻Ｔ_０を含む時間帯（ｔ_１〜ｔ_２）の学習用フーリエ変換信号ｓ_ｇ（ｎ）_{ｔ１〜ｔ２}を抽出する。

次に、特徴量算出部７は、図７に示すように、抽出した学習用フーリエ変換信号ｓ_ｇ（ｎ）_{ｔ１〜ｔ２}に含まれている複素時間信号の中から、設定周波数ビンｋ_ｂの複素時間信号ｓ_ｇハット（ｎ）_{ｔ１〜ｔ２}の切り出しを行う。
電子出願の都合上、明細書の文章中では、ｓの文字の上に“＾”の記号を付することができないので、ｓ_ｇハット（ｎ）_{ｔ１〜ｔ２}のように表記している。
図７は、特徴量算出部７による特徴量の算出処理を示す説明図である。
特徴量算出部７は、切り出した複素時間信号ｓ_ｇハット（ｎ）_{ｔ１〜ｔ２}から、抽出した学習用フーリエ変換信号ｓ_ｇ（ｎ）_{ｔ１〜ｔ２}の変動を特徴量として算出する（図４のステップＳＴ４）。

以下、特徴量算出部７による特徴量の算出処理を具体的に説明する。
特徴量算出部７は、切り出した複素時間信号ｓ_ｇハット（ｎ）_{ｔ１〜ｔ２}から、以下の式（４）に示すように、複素時間信号ｓ_ｇハット（ｎ）_{ｔ１〜ｔ２}の瞬時振幅ａ_ｇ（ｎ）を算出する。

式（４）において、Ｉ_ｇ（ｎ）は、複素時間信号ｓ_ｇハット（ｎ）_{ｔ１〜ｔ２}の実部、Ｑ_ｇ（ｎ）は、複素時間信号ｓ_ｇハット（ｎ）_{ｔ１〜ｔ２}の虚部である。

また、特徴量算出部７は、切り出した複素時間信号ｓ_ｇハット（ｎ）_{ｔ１〜ｔ２}から、以下の式（５）に示すように、複素時間信号ｓ_ｇハット（ｎ）_{ｔ１〜ｔ２}の瞬時位相φ_ｇ（ｎ）を算出する。

また、特徴量算出部７は、切り出した複素時間信号ｓ_ｇハット（ｎ）_{ｔ１〜ｔ２}から、以下の式（６）に示すように、複素時間信号ｓ_ｇハット（ｎ）_{ｔ１〜ｔ２}の瞬時周波数ｆ_ｇ（ｎ）を算出する。

次に、特徴量算出部７は、以下の式（７）に示すように、第１の特徴量として、複素時間信号ｓ_ｇハット（ｎ）_{ｔ１〜ｔ２}における瞬時振幅ａ_ｇ（ｎ）の分散値ｐ_ｇ１を算出する。

特徴量算出部７は、以下の式（９）に示すように、第２の特徴量として、複素時間信号ｓ_ｇハット（ｎ）_{ｔ１〜ｔ２}における瞬時振幅ａ_ｇ（ｎ）の歪度ｐ_ｇ２を算出する。

特徴量算出部７は、以下の式（１１）に示すように、第３の特徴量として、複素時間信号ｓ_ｇハット（ｎ）_{ｔ１〜ｔ２}における瞬時振幅ａ_ｇ（ｎ）の尖度ｐ_ｇ３を算出する。

次に、特徴量算出部７は、以下の式（１２）に示すように、第４の特徴量として、複素時間信号ｓ_ｇハット（ｎ）_{ｔ１〜ｔ２}における瞬時位相φ_ｇ（ｎ）の分散値ｐ_ｇ４を算出する。

特徴量算出部７は、以下の式（１４）に示すように、第５の特徴量として、複素時間信号ｓ_ｇハット（ｎ）_{ｔ１〜ｔ２}における瞬時位相φ_ｇ（ｎ）の歪度ｐ_ｇ５を算出する。

特徴量算出部７は、以下の式（１６）に示すように、第６の特徴量として、複素時間信号ｓ_ｇハット（ｎ）_{ｔ１〜ｔ２}における瞬時位相φ_ｇ（ｎ）の尖度ｐ_ｇ６を算出する。

次に、特徴量算出部７は、以下の式（１７）に示すように、第７の特徴量として、複素時間信号ｓ_ｇハット（ｎ）_{ｔ１〜ｔ２}における瞬時周波数ｆ_ｇ（ｎ）の分散値ｐ_ｇ７を算出する。

特徴量算出部７は、以下の式（１９）に示すように、第８の特徴量として、複素時間信号ｓ_ｇハット（ｎ）_{ｔ１〜ｔ２}における瞬時周波数ｆ_ｇ（ｎ）の歪度ｐ_ｇ８を算出する。

特徴量算出部７は、以下の式（２１）に示すように、第９の特徴量として、複素時間信号ｓ_ｇハット（ｎ）_{ｔ１〜ｔ２}における瞬時周波数ｆ_ｇ（ｎ）の尖度ｐ_ｇ９を算出する。

特徴量算出部７は、第２のパラメータ設定部６により設定された第１から第９の特徴量ｐ_ｇｉ（ｉ＝１，２，・・・，９）の重みｗ_ｉを示す重み付けパラメータを取得する。
第１から第９の特徴量ｐ_ｇｉの中には、無線機の識別に有効な特徴量と、識別に有効でない特徴量とが存在している。このため、識別に有効な特徴量には大きな重みが設定され、識別に有効でない特徴量には小さな重みが設定される。
例えば、ｗ_ｉ（ｉ＝１，２，・・・，９）は、０以上１以下の値であり、ｗ_１＋ｗ_２＋・・・＋ｗ_９＝１である。
第２のパラメータ設定部６による重み付けパラメータの設定方法は、特に問わないが、例えば、第１から第９の特徴量ｐ_ｇｉの分散値で、第１から第９の特徴量ｐ_ｇｉを規格化し、第１から第９の特徴量ｐ_ｇｉを規格化した値を重みとして用いる方法が考えられる。

特徴量算出部７は、以下の式（２２）に示すように、第１から第９の特徴量ｐ_ｇｉと重みｗ_ｉとをそれぞれ乗算することで、第１から第９の特徴量ｐ_ｇｉの重み付けを行う。
ｐ’_ｇｉ＝ｐ_ｇｉ×ｗ_ｉ （２２）
ｉ＝１，２，・・・，９
特徴量算出部７は、重み付け後の第１から第９の特徴量ｐ’_ｇｉを並べたベクトルを、既知の無線機の特徴量ベクトルｐ_ｇとして算出する。
特徴量算出部７は、算出した既知の無線機の特徴量ベクトルｐ_ｇを、学習用無線信号に付加されているシリアルナンバーと一緒に、第２のデータベース部１０に格納する（図４のステップＳＴ５）。
また、特徴量算出部７は、既知の無線機の特徴量ベクトルｐ_ｇを第２のデータベース部１０に格納した旨を示す情報を第１のパラメータ設定部５に出力する。

ここで、第１のパラメータ設定部５による特徴量算出用のパラメータの設定処理について説明する。
図８は、特徴量算出用のパラメータと複数の既知の無線機に係る特徴量ベクトルとの関係を示す説明図である。
図８の例では、〇は、種類（Ａ）の既知の無線機に係る特徴量ベクトルｐ_ｇ、◇は、種類（Ｂ）の既知の無線機に係る特徴量ベクトルｐ_ｇ、△は、種類（Ｃ）の既知の無線機に係る特徴量ベクトルｐ_ｇである。
図８の例では、種類（Ａ）に属する無線機の個体の数が１０であるため、１０個の特徴量ベクトル〇が得られている。
また、種類（Ｂ）に属する無線機の個体の数が１０であるため、１０個の特徴量ベクトル◇が得られており、種類（Ｃ）に属する無線機の個体の数が１０であるため、１０個の特徴量ベクトル△が得られている。

図８Ａは、特徴量算出用のパラメータが適正でないために、種類（Ａ）の無線機に係る特徴量ベクトル〇と、種類（Ｂ）の無線機に係る特徴量ベクトル◇と、種類（Ｃ）の無線機に係る特徴量ベクトル△とが一部重なっている例を示している。
既知の無線機に係る特徴量ベクトル〇と特徴量ベクトル◇と特徴量ベクトル△とが重なっている場合、これらの特徴量ベクトル〇，◇，△と、識別対象の無線機に係る特徴量ベクトルとを比較しても、特徴量ベクトル〇，◇，△の中で、識別対象の無線機に係る特徴量ベクトルと最も類似している特徴量ベクトルを高精度に特定することができない。
したがって、特徴量ベクトル〇，◇，△の中で、識別対象の無線機に係る特徴量ベクトルと最も類似している特徴量ベクトルを高精度に特定する上で、既知の無線機に係る特徴量ベクトル〇と特徴量ベクトル◇と特徴量ベクトル△とが重なっていないことが望ましい。
図８Ｂは、特徴量算出用のパラメータが適正であるために、種類（Ａ）の無線機に係る特徴量ベクトル〇と、種類（Ｂ）の無線機に係る特徴量ベクトル◇と、種類（Ｃ）の無線機に係る特徴量ベクトル△とが重なっていない例を示している。

第１のパラメータ設定部５は、特徴量算出用のパラメータの初期設定として、時間帯（ｔ_１〜ｔ_２）の長さである時間範囲Δｔを示すパラメータと、設定周波数ビンｋ_ｂを示すパラメータとを任意の値に設定する。
時間帯（ｔ_１〜ｔ_２）については、学習用のエネルギー時間波形Ｅ_ｇ（ｍ）に基づいて設定する方法が考えられる。
例えば、学習用無線信号ｓ_ｇ（ｎ）の立ち上がり時刻Ｔ_０よりも前の時刻において、学習用のエネルギー時間波形Ｅ_ｇ（ｍ）の最大値よりもエネルギーが３０ｄＢ下がっている時刻をｔ_１とする。また、学習用無線信号ｓ_ｇ（ｎ）の立ち上がり時刻Ｔ_０よりも後の時刻において、学習用のエネルギー時間波形Ｅ_ｇ（ｍ）の最大値よりもエネルギーが３０ｄＢ下がっている時刻をｔ_２とする。
また、ｔ_２については、学習用フーリエ変換信号ｓ_ｇ（ｎ）の周波数方向の変動に基づいて設定する方法が考えられる。
例えば、学習用フーリエ変換信号ｓ_ｇ（ｎ）の過渡応答時の周波数変動が、特定の周波数範囲内に収まった時刻をｔ_２とする。

設定周波数ビンｋ_ｂについては、学習用のエネルギー時間波形Ｅ_ｇ（ｍ）において、平均電力が最大値をとる周波数ビンを採用する方法が考えられる。
また、設定周波数ビンｋ_ｂについては、学習用のエネルギー時間波形Ｅ_ｇ（ｍ）において、平均電力が閾値以上の周波数ビンを検出して、検出した１以上の周波数ビンに係る周波数の平均を求め、平均の周波数の周波数ビンを採用する方法が考えられる。

第１のパラメータ設定部５は、特徴量算出部７から既知の無線機の特徴量ベクトルｐ_ｇを第２のデータベース部１０に格納した旨を示す情報を受けると、第２のデータベース部１０から、複数の既知の無線機の特徴量ベクトルｐ_ｇを取得する。

第１のパラメータ設定部５は、複数の既知の無線機に係る特徴量ベクトルｐ_ｇの間に重なりがあるか否かを判定し、例えば、図８Ｂに示すように、複数の既知の無線機に係る特徴量ベクトルｐ_ｇの間に重なりがなければ、先に設定している特徴量算出用のパラメータを有効として、パラメータの設定処理を終了する。
第１のパラメータ設定部５は、例えば、図８Ａに示すように、複数の既知の無線機に係る特徴量ベクトルｐ_ｇの間に重なりがあれば、特徴量算出用のパラメータとして、先に設定している時間範囲Δｔ又は設定周波数ビンｋ_ｂのうち、少なくとも１つを更新する。
特徴量算出部７は、第１のパラメータ設定部５により更新された特徴量算出用のパラメータを用いて、再度、複数の既知の無線機の特徴量ベクトルｐ_ｇをそれぞれ算出する。
第１のパラメータ設定部５は、特徴量算出部７により再度算出された複数の既知の無線機に係る特徴量ベクトルｐ_ｇの間に重なりがあるか否かを判定する。
複数の既知の無線機に係る特徴量ベクトルｐ_ｇの間に重なりがなくなるまで、第１のパラメータ設定部５による特徴量算出用のパラメータの更新処理と、特徴量算出部７による特徴量ベクトルｐ_ｇの算出処理とが繰り返される。
なお、パラメータの更新処理と、特徴量ベクトルｐ_ｇの算出処理とを繰り返しても、複数の既知の無線機に係る重なりがなくならない場合には、重なりが最も小さくなったパラメータを適正なパラメータとして、パラメータの設定処理を終了する。

第１のパラメータ設定部５により特徴量算出用のパラメータが更新されることで、特徴量算出用のパラメータが適正な値になると、図８Ｂに示すように、第２のデータベース部１０に格納される複数の既知の無線機に係る特徴量ベクトルｐ_ｇの間の重なりがなくなる、もしくは、重なりが最も小さくなる。
特徴量算出用のパラメータが更新されたのち、特徴量算出部７により再度算出された複数の既知の無線機に係る特徴量ベクトルｐ_ｇが第２のデータベース部１０に格納される際、不適正な特徴量算出用のパラメータに従って算出されている複数の既知の無線機に係る特徴量ベクトルｐ_ｇは、第２のデータベース部１０から削除される。
図１０は、第２のデータベース部１０に格納される複数の既知の無線機に係る特徴量の一例を示す説明図である。
図１０では、既知の無線機の種類が、種類（Ａ）、種類（Ｂ）の例を示しており、また、種類（Ａ）に属する無線機が、シリアルナンバー（１）の個体、シリアルナンバー（２）の個体の例を示している。
また、図１０には、特徴量算出用のパラメータ及び重み付けパラメータの一例を示している。

この実施の形態１では、第１のパラメータ設定部５によって、特徴量算出用のパラメータが適正な値になるまで、特徴量算出用のパラメータが更新される例を示しているが、特徴量算出用のパラメータの更新タイミングは、これに限るものではない。
例えば、信号受信部１によって、既知の無線機から新たに送信された学習用無線信号が受信されることで、特徴量算出部７によって、新たに特徴量が算出されたタイミングで、第１のパラメータ設定部５が、特徴量算出用のパラメータを更新するようにしてもよい。

＜識別時の処理＞
次に、識別対象の無線機から送信された無線信号を用いて、無線機を識別する際の処理内容を説明する。
図１１は、この発明の実施の形態１による無線機識別方法における識別処理を示すフローチャートである。
信号受信部１は、識別対象の無線機から送信された無線信号を識別用無線信号として受信し、受信した識別用無線信号をアナログ信号からデジタル信号に変換する（図１１のステップＳＴ１１）。
信号受信部１は、デジタルの識別用無線信号ｓ_ｄ（ｎ）をフーリエ変換部２に出力する。ｓ_ｄ（ｎ）におけるｎは、識別用無線信号のサンプル番号である。

フーリエ変換部２は、以下の式（２３）に示すように、信号受信部１から出力されたデジタルの識別用無線信号ｓ_ｄ（ｎ）をＳＴＦＴする（図１１のステップＳＴ１２）。
フーリエ変換部２は、識別用無線信号ｓ_ｄ（ｎ）のＳＴＦＴ結果を示すフーリエ変換信号である識別用フーリエ変換信号Ｓ_ｄ（ｍ，ｋ）を立ち上がり検出部３及び特徴量算出部７に出力する（図１１のステップＳＴ１２）。図５は、識別用無線信号ｓ_ｄ（ｎ）のＳＴＦＴについても表している。

立ち上がり検出部３は、識別用無線信号ｓ_ｄ（ｎ）の立ち上がり時刻Ｔ_０を検出する処理を実施する（図１１のステップＳＴ１３）。
以下、立ち上がり検出部３による立ち上がり時刻Ｔ_０の検出処理を具体的に説明する。

まず、立ち上がり検出部３は、以下の式（２４）に示すように、フーリエ変換部２から出力された識別用フーリエ変換信号Ｓ_ｄ（ｍ，ｋ）の２乗値｜Ｓ_ｄ（ｍ，ｋ）｜^２をスペクトログラムＳＰ_ｄ（ｍ，ｋ）として算出する。

次に、立ち上がり検出部３は、以下の式（２５）に示すように、スペクトログラムＳＰ_ｄ（ｍ，ｋ）におけるエネルギーの時間波形である識別用のエネルギー時間波形Ｅ_ｄ（ｍ）を算出する。
識別用のエネルギー時間波形Ｅ_ｄ（ｍ）は、スペクトログラムＳＰ_ｄ（ｍ，ｋ）に含まれている複数の周波数における電力のうち、同一時刻の電力同士がそれぞれ合計されることで算出される各々の時刻における電力の合計値である。

式（２５）において、ｋ_１は、スペクトログラムＳＰ_ｄ（ｍ，ｋ）に含まれている複数の周波数のうち、電力を合計する周波数の下限周波数を示す周波数ビン番号、ｋ_ｎは、スペクトログラムＳＰ_ｄ（ｍ，ｋ）に含まれている複数の周波数のうち、電力を合計する周波数の上限周波数を示す周波数ビン番号である。
ｋ＝ｋ_１，・・・，ｋ_ｎである。

次に、立ち上がり検出部３は、識別用のエネルギー時間波形Ｅ_ｄ（ｍ）と閾値Ｅ_ｔｈとを比較する。
立ち上がり検出部３は、識別用のエネルギー時間波形Ｅ_ｄ（ｍ）と閾値Ｅ_ｔｈとの比較結果を参照することで、最初は、ノイズレベルのエネルギーを示していた識別用のエネルギー時間波形Ｅ_ｄ（ｍ）が、閾値Ｅ_ｔｈに到達したか否かを判定する。
立ち上がり検出部３は、識別用のエネルギー時間波形Ｅ_ｄ（ｍ）が閾値Ｅ_ｔｈに到達した時刻を、識別用無線信号ｓ_ｄ（ｎ）の立ち上がり時刻Ｔ_０として検出する。
この実施の形態１では、立ち上がり検出部３が、識別用無線信号ｓ_ｄ（ｎ）の立ち上がり時刻Ｔ_０を検出する際、スペクトログラムＳＰ_ｄ（ｍ，ｋ）に含まれている複数の周波数における電力のうち、同一時刻の電力同士をそれぞれ合計するようにしている。これにより、立ち上がり時に電力が大きく変動する周波数が含まれている場合でも、当該電力の変動の影響が低減されるため、立ち上がり時刻Ｔ_０の検出精度が向上する。

特徴量算出部７は、第１のパラメータ設定部５により設定された特徴量算出用のパラメータから、立ち上がり時刻Ｔ_０を含む時間帯（ｔ_１〜ｔ_２）及び設定周波数ビンｋ_ｂを認識する。
この特徴量算出用のパラメータは、第１のパラメータ設定部５によって、複数の既知の無線機に係る特徴量ベクトルｐ_ｇの間に重なりがなくなるように、もしくは、重なりが最も小さくなるように設定されたパラメータである。
次に、特徴量算出部７は、図６に示すように、フーリエ変換部２から出力された識別用フーリエ変換信号Ｓ_ｄ（ｍ，ｋ）の中から、立ち上がり検出部３により検出された識別用無線信号ｓ_ｄ（ｎ）の立ち上がり時刻Ｔ_０を含む時間帯（ｔ_１〜ｔ_２）の識別用フーリエ変換信号ｓ_ｄ（ｎ）_{ｔ１〜ｔ２}を抽出する。

次に、特徴量算出部７は、図７に示すように、抽出した識別用フーリエ変換信号ｓ_ｄ（ｎ）_{ｔ１〜ｔ２}に含まれている複素時間信号の中から、設定周波数ビンｋ_ｂの複素時間信号ｓ_ｄハット（ｎ）_{ｔ１〜ｔ２}の切り出しを行う。
電子出願の都合上、明細書の文章中では、ｓの文字の上に“＾”の記号を付することができないので、ｓ_ｄハット（ｎ）_{ｔ１〜ｔ２}のように表記している。
特徴量算出部７は、切り出した複素時間信号ｓ_ｄハット（ｎ）_{ｔ１〜ｔ２}から、抽出した識別用フーリエ変換信号ｓ_ｄ（ｎ）_{ｔ１〜ｔ２}の変動を特徴量として算出する（図１１のステップＳＴ１４）。

以下、特徴量算出部７による特徴量の算出処理を具体的に説明する。
特徴量算出部７は、切り出した複素時間信号ｓ_ｄハット（ｎ）_{ｔ１〜ｔ２}から、以下の式（２６）に示すように、複素時間信号ｓ_ｄハット（ｎ）_{ｔ１〜ｔ２}の瞬時振幅ａ_ｄ（ｎ）を算出する。

式（２６）において、Ｉ_ｄ（ｎ）は、複素時間信号ｓ_ｄハット（ｎ）_{ｔ１〜ｔ２}の実部、Ｑ_ｄ（ｎ）は、複素時間信号ｓ_ｄハット（ｎ）_{ｔ１〜ｔ２}の虚部である。

また、特徴量算出部７は、切り出した複素時間信号ｓ_ｄハット（ｎ）_{ｔ１〜ｔ２}から、以下の式（２７）に示すように、複素時間信号ｓ_ｄハット（ｎ）_{ｔ１〜ｔ２}の瞬時位相φ_ｄ（ｎ）を算出する。

また、特徴量算出部７は、切り出した複素時間信号ｓ_ｄハット（ｎ）_{ｔ１〜ｔ２}から、以下の式（２８）に示すように、複素時間信号ｓ_ｄハット（ｎ）_{ｔ１〜ｔ２}の瞬時周波数ｆ_ｄ（ｎ）を算出する。

次に、特徴量算出部７は、以下の式（２９）に示すように、第１の特徴量として、複素時間信号ｓ_ｄハット（ｎ）_{ｔ１〜ｔ２}における瞬時振幅ａ_ｄ（ｎ）の分散値ｐ_ｄ１を算出する。

特徴量算出部７は、以下の式（３１）に示すように、第２の特徴量として、複素時間信号ｓ_ｄハット（ｎ）_{ｔ１〜ｔ２}における瞬時振幅ａ_ｄ（ｎ）の歪度ｐ_ｄ２を算出する。

特徴量算出部７は、以下の式（３３）に示すように、第３の特徴量として、複素時間信号ｓ_ｄハット（ｎ）_{ｔ１〜ｔ２}における瞬時振幅ａ_ｄ（ｎ）の尖度ｐ_ｄ３を算出する。

次に、特徴量算出部７は、以下の式（３４）に示すように、第４の特徴量として、複素時間信号ｓ_ｄハット（ｎ）_{ｔ１〜ｔ２}における瞬時位相φ_ｄ（ｎ）の分散値ｐ_ｄ４を算出する。

特徴量算出部７は、以下の式（３６）に示すように、第５の特徴量として、複素時間信号ｓ_ｄハット（ｎ）_{ｔ１〜ｔ２}における瞬時位相φ_ｄ（ｎ）の歪度ｐ_ｄ５を算出する。

特徴量算出部７は、以下の式（３８）に示すように、第６の特徴量として、複素時間信号ｓ_ｄハット（ｎ）_{ｔ１〜ｔ２}における瞬時位相φ_ｄ（ｎ）の尖度ｐ_ｄ６を算出する。

次に、特徴量算出部７は、以下の式（３９）に示すように、第７の特徴量として、複素時間信号ｓ_ｄハット（ｎ）_{ｔ１〜ｔ２}における瞬時周波数ｆ_ｄ（ｎ）の分散値ｐ_ｄ７を算出する。

特徴量算出部７は、以下の式（４１）に示すように、第８の特徴量として、複素時間信号ｓ_ｄハット（ｎ）_{ｔ１〜ｔ２}における瞬時周波数ｆ_ｄ（ｎ）の歪度ｐ_ｄ８を算出する。

特徴量算出部７は、以下の式（４３）に示すように、第９の特徴量として、複素時間信号ｓ_ｄハット（ｎ）_{ｔ１〜ｔ２}における瞬時周波数ｆ_ｄ（ｎ）の尖度ｐ_ｄ９を算出する。

特徴量算出部７は、第２のパラメータ設定部６により設定された第１から第９の特徴量ｐ_ｄｉ（ｉ＝１，２，・・・，９）の重みｗ_ｉを示す重み付けパラメータを取得する。
特徴量算出部７は、以下の式（４４）に示すように、第１から第９の特徴量ｐ_ｄｉと重みｗ_ｉとをそれぞれ乗算することで、第１から第９の特徴量ｐ_ｄｉの重み付けを行う。
ｐ’_ｄｉ＝ｐ_ｄｉ×ｗ_ｉ （４４）
ｉ＝１，２，・・・，９
特徴量算出部７は、重み付け後の第１から第９の特徴量ｐ’_ｄｉを並べたベクトルを、識別対象の無線機の特徴量ベクトルｐ_ｄとして算出する。
特徴量算出部７は、算出した識別対象の無線機の特徴量ベクトルｐ_ｄを無線機識別部１１に出力する。

無線機識別部１１は、特徴量算出部７により算出された識別対象の無線機に係る特徴量ベクトルｐ_ｄと、第２のデータベース部１０により格納されている複数の既知の無線機に係る特徴量ベクトルｐ_ｇとをそれぞれ比較する。
そして、無線機識別部１１は、識別対象の無線機に係る特徴量ベクトルｐ_ｄと、複数の既知の無線機に係る特徴量ベクトルｐ_ｇとの比較結果に基づいて、識別対象の無線機を識別する（図１１のステップＳＴ１５）。

以下、無線機識別部１１による無線機の識別処理を具体的に説明する。
図１２は、無線機識別部１１による無線機の識別処理を示す説明図である。
図１２では、説明の簡単化のため、特徴量算出部７により算出される第１から第９の特徴量のうち、第１から第３の特徴量だけを例示している。

無線機識別部１１は、以下の式（４５）に示すように、識別対象の無線機に係る特徴量ベクトルｐ_ｄと、第２のデータベース部１０により格納されている複数の既知の無線機に係る特徴量ベクトルｐ_ｇとの間のマハラノビス距離ｄ_ｉｊ，ｈをそれぞれ算出する。
ここでは、説明の便宜上、既知の無線機がＨ個あり、Ｈ個の無線機に係る特徴量ベクトルｐ_ｇが第２のデータベース部１０にそれぞれ格納されているものとする。
Ｈ個の無線機の中には、例えば、種類（Ａ）に属する無線機、種類（Ｂ）に属する無線機、あるいは、種類（Ｃ）に属する無線機が含まれている。

式（４５）において、ｘ_ｉ，ｈは、ｎ次元の特徴量空間における既知の無線機ｈ（ｈ＝１，・・・，Ｈ）に係る特徴量ベクトルｐ_ｇの位置、ｘ_ｊは、ｎ次元の特徴量空間における識別対象の無線機に係る特徴量ベクトルｐ_ｄの位置である。
（ｘ_ｉ，ｈ−ｘ_ｊ）は、ｎ次元の特徴量空間における２点間のユークリッド距離、Ｓは、既知の無線機ｈから送信される学習用無線信号ｓ_ｇ（ｎ）の特徴量空間における分布の分散共分散行列である。
マハラノビス距離ｄ_ｉｊ，ｈは、学習用無線信号ｓ_ｇ（ｎ）の特徴量空間における分布の分散による影響が正規化された空間での距離となる。

無線機識別部１１は、算出したＨ個のマハラノビス距離ｄ_ｉｊ，ｈを相互に比較し、Ｈ個のマハラノビス距離ｄ_ｉｊ，ｈの中で、最小のマハラノビス距離ｄ_ｉｊ，ｈを特定する。
無線機識別部１１は、Ｈ個の無線機の中で、最小のマハラノビス距離ｄ_ｉｊ，ｈに係る既知の無線機が、識別対象の無線機と最も類似していると判断する。
無線機識別部１１は、識別対象の無線機の識別結果として、最小のマハラノビス距離ｄ_ｉｊ，ｈに係る既知の無線機のシリアルナンバーを出力する。これにより、無線機の種類だけでなく、無線機の個体識別が可能になる。
図１２の例では、３つの既知の無線機に係る特徴量ベクトルｐ_ｇの位置が示されている。３つの既知の無線機の中で、特徴量ベクトルｐ_ｇの位置が●である既知の無線機とのマハラノビス距離が最小であるため、識別対象の無線機は、特徴量ベクトルｐ_ｇの位置が●の無線機であると識別される。

ここでは、無線機識別部１１が、マハラノビス距離ｄ_ｉｊ，ｈを用いて、無線機を識別する例を示しているが、これに限るものではなく、例えば、最小二乗法による線形識別、サポートベクターマシン（ＳＶＭ：Ｓｕｐｐｏｒｔ Ｖｅｃｔｏｒ Ｍａｃｈｉｎｅ）、あるいは、ニューラルネットワークを用いて、無線機を識別するようにしてもよい。

以上で明らかなように、この実施の形態１によれば、フーリエ変換部２から出力されたフーリエ変換信号のスペクトログラムにおけるエネルギーの時間波形を算出し、エネルギーの時間波形から無線信号の立ち上がり時刻を検出する立ち上がり検出部３と、フーリエ変換部２から出力されたフーリエ変換信号の中から、立ち上がり検出部３により検出された立ち上がり時刻を含む時間帯のフーリエ変換信号を抽出し、抽出したフーリエ変換信号に含まれている設定周波数成分の複素時間信号から、抽出したフーリエ変換信号の変動を特徴量として算出する特徴量算出部７とを設け、無線機識別部１１が、特徴量算出部７により算出された特徴量に基づいて、無線機を識別するように構成したので、無線機の個体を識別することができる効果を奏する。

なお、本願発明はその発明の範囲内において、実施の形態の任意の構成要素の変形、もしくは実施の形態の任意の構成要素の省略が可能である。

この発明は、無線機を識別する無線機識別装置及び無線機識別方法に適している。

１ 信号受信部、２ フーリエ変換部、３ 立ち上がり検出部、４ パラメータ設定部、５ 第１のパラメータ設定部、６ 第２のパラメータ設定部、７ 特徴量算出部、８ データベース部、９ 第１のデータベース部、１０ 第２のデータベース部、１１ 無線機識別部、２１ 信号受信回路、２２ フーリエ変換回路、２３ 立ち上がり検出回路、２４ パラメータ設定回路、２５ 特徴量算出回路、２６ データベース回路、２７ 無線機識別回路、３１ メモリ、３２ プロセッサ。

Claims (8)

1. 識別対象の無線機から送信された無線信号を受信する信号受信部と、
   前記信号受信部により受信された無線信号を短時間フーリエ変換し、前記無線信号の短時間フーリエ変換結果を示すフーリエ変換信号を出力するフーリエ変換部と、
   前記フーリエ変換部から出力されたフーリエ変換信号のスペクトログラムにおけるエネルギーの時間波形を算出し、前記エネルギーの時間波形から前記無線信号の立ち上がり時刻を検出する立ち上がり検出部と、
   前記フーリエ変換部から出力されたフーリエ変換信号の中から、前記立ち上がり検出部により検出された立ち上がり時刻を含む時間帯のフーリエ変換信号を抽出し、抽出したフーリエ変換信号に含まれている設定周波数成分の複素時間信号から、抽出したフーリエ変換信号の変動を特徴量として算出する特徴量算出部と、
   前記特徴量算出部により算出された特徴量に基づいて、前記無線機を識別する無線機識別部と
   を備えた無線機識別装置。
2. 前記信号受信部は、既知の無線機から送信された無線信号を受信し、
   前記フーリエ変換部は、前記信号受信部により受信された既知の無線機の無線信号を短時間フーリエ変換し、既知の無線機の無線信号の短時間フーリエ変換結果を示すフーリエ変換信号を出力し、
   前記立ち上がり検出部は、前記フーリエ変換部から出力された既知の無線機に係るフーリエ変換信号のスペクトログラムにおけるエネルギーの時間波形を算出し、当該エネルギーの時間波形から既知の無線機の無線信号の立ち上がり時刻を検出し、
   前記特徴量算出部は、前記フーリエ変換部から出力された既知の無線機に係るフーリエ変換信号の中から、前記立ち上がり検出部により検出された立ち上がり時刻を含む時間帯のフーリエ変換信号を抽出し、抽出した既知の無線機に係るフーリエ変換信号に含まれている設定周波数成分の複素時間信号から、抽出した既知の無線機に係るフーリエ変換信号の変動を特徴量として算出し、算出した既知の無線機に係る特徴量をデータベース部に格納し、
   前記無線機識別部は、前記特徴量算出部により算出された識別対象の無線機に係る特徴量と、前記データベース部により格納されている既知の無線機に係る特徴量とを比較し、特徴量の比較結果に基づいて、識別対象の無線機を識別することを特徴とする請求項１記載の無線機識別装置。
3. 特徴量算出用のパラメータとして、前記時間帯の長さを示すパラメータと、前記設定周波数成分を示すパラメータとを設定する第１のパラメータ設定部を備え、
   前記特徴量算出部は、前記第１のパラメータ設定部により設定された特徴量算出用のパラメータを用いて、識別対象の無線機に係る特徴量及び既知の無線機に係る特徴量のそれぞれを算出することを特徴とする請求項２記載の無線機識別装置。
4. 前記第１のパラメータ設定部は、互いに異なる既知の無線機から送信された無線信号のそれぞれが前記信号受信部に受信された際に、前記特徴量算出部によりそれぞれ算出される特徴量が重ならないように、あるいは、前記特徴量算出部によりそれぞれ算出される特徴量の重なりが最も小さくなるように、前記特徴量算出用のパラメータを更新することを特徴とする請求項３記載の無線機識別装置。
5. 前記立ち上がり検出部は、前記スペクトログラムに含まれている複数の周波数における電力のうち、同一時刻の電力同士をそれぞれ合計することで、前記エネルギーの時間波形として、各々の時刻における電力の合計値を算出し、無線信号の立ち上がり時刻として、前記電力の合計値が閾値に到達した時刻を検出することを特徴とする請求項１記載の無線機識別装置。
6. 前記特徴量算出部は、
   第１から第３の特徴量として、前記複素時間信号における振幅の分散値、歪度及び尖度をそれぞれ算出し、第４から第６の特徴量として、前記複素時間信号における位相の分散値、歪度及び尖度をそれぞれ算出し、第７から第９の特徴量として、前記複素時間信号における周波数の分散値、歪度及び尖度をそれぞれ算出することを特徴とする請求項１記載の無線機識別装置。
7. 前記第１から第９の特徴量の重みを示す重み付けパラメータを設定する第２のパラメータ設定部を備え、
   前記特徴量算出部は、前記第２のパラメータ設定部により設定された重み付けパラメータに従って前記第１から第９の特徴量の重み付けを実施することを特徴とする請求項６記載の無線機識別装置。
8. 信号受信部が、識別対象の無線機から送信された無線信号を受信し、
   フーリエ変換部が、前記信号受信部により受信された無線信号を短時間フーリエ変換し、前記無線信号の短時間フーリエ変換結果を示すフーリエ変換信号を出力し、
   立ち上がり検出部が、前記フーリエ変換部から出力されたフーリエ変換信号のスペクトログラムにおけるエネルギーの時間波形を算出し、前記エネルギーの時間波形から前記無線信号の立ち上がり時刻を検出し、
   特徴量算出部が、前記フーリエ変換部から出力されたフーリエ変換信号の中から、前記立ち上がり検出部により検出された立ち上がり時刻を含む時間帯のフーリエ変換信号を抽出し、抽出したフーリエ変換信号に含まれている設定周波数成分の複素時間信号から、抽出したフーリエ変換信号の変動を特徴量として算出し、
   無線機識別部が、前記特徴量算出部により算出された特徴量に基づいて、前記無線機を識別する
   無線機識別方法。

Images