JP6567242B2

JP6567242B2 - Radio identification device and radio identification method

Info

Publication number: JP6567242B2
Application number: JP2019530259A
Authority: JP
Inventors: 貴文松田; 正資大島; 網嶋　武; 武網嶋; 信弘鈴木
Original assignee: Mitsubishi Electric Corp
Current assignee: Mitsubishi Electric Corp
Priority date: 2017-07-18
Filing date: 2017-07-18
Publication date: 2019-08-28
Anticipated expiration: 2037-07-18
Also published as: JPWO2019016853A1; WO2019016853A1

Description

この発明は、無線機を識別する無線機識別装置及び無線機識別方法に関するものである。 The present invention relates to a radio identification device and a radio identification method for identifying a radio.

無線機から送信される無線信号の立ち上がりなどの過渡応答は、個体差があることが知られており、近年では、無線機から送信される無線信号を解析することで、無線機を識別する無線機識別装置が提案されている。
例えば、無線機識別装置は、無線機から送信される無線信号の立ち上がりの特徴量を算出し、算出した特徴量に基づいて、無線機を識別する。
無線信号の立ち上がりの特徴量としては、立ち上がり時の無線信号における振幅等の分散値、歪度及び尖度などがある。It is known that transient responses such as the rise of a radio signal transmitted from a radio device have individual differences. In recent years, a radio signal that identifies a radio device by analyzing a radio signal transmitted from the radio device. A machine identification device has been proposed.
For example, the wireless device identification device calculates the feature amount of the rising edge of the wireless signal transmitted from the wireless device, and identifies the wireless device based on the calculated feature amount.
The feature amount of the rising edge of the wireless signal includes a dispersion value such as an amplitude, skewness, and kurtosis in the wireless signal at the time of rising.

ただし、信号対雑音比（Ｓｉｇｎａｌ−ｔｏ−ＮｏｉｓｅＲａｔｉｏ：ＳＮＲ）が低い環境下では、立ち上がりの特徴量の算出精度が低下して、無線機の識別精度が低下する。このため、無線機から送信される無線信号を短時間フーリエ変換することで、ＳＮＲを改善してから、立ち上がりの特徴量を算出する無線機識別装置が以下の特許文献１に開示されている。
この無線機識別装置では、立ち上がりの特徴量として、立ち上がりの時定数を算出している。However, in an environment where the signal-to-noise ratio (Signal-to-Noise Ratio: SNR) is low, the calculation accuracy of the rising feature amount decreases, and the identification accuracy of the radio device decreases. For this reason, Patent Document 1 below discloses a wireless device identification apparatus that calculates a feature quantity of rising after improving SNR by short-time Fourier transform of a wireless signal transmitted from a wireless device.
In this wireless device identification device, a rising time constant is calculated as a rising feature amount.

特開２００２−２６８２６号公報JP 2002-26826 A

従来の無線機識別装置は、立ち上がりの時定数に基づく無線機の識別方法を採用している。立ち上がりの時定数に基づく無線機の識別方法では、無線機の種類については識別することが可能である。しかし、同一種類の無線機は、異なる無線機であっても、立ち上がりの時定数の差異が小さいため、無線機の個体を識別することができないという課題があった。 A conventional radio identification device employs a radio identification method based on a rising time constant. In the wireless device identification method based on the rise time constant, the type of the wireless device can be identified. However, even if the same type of wireless device is a different wireless device, there is a problem that the individual wireless devices cannot be identified because the difference in the time constant of the rise is small.

この発明は上記のような課題を解決するためになされたもので、無線機の個体識別を行うことができる無線機識別装置及び無線機識別方法を得ることを目的とする。 The present invention has been made to solve the above-described problems, and an object of the present invention is to provide a wireless device identification apparatus and a wireless device identification method capable of performing individual identification of a wireless device.

この発明に係る無線機識別装置は、識別対象の無線機から送信された無線信号を受信する信号受信部と、信号受信部により受信された無線信号を短時間フーリエ変換し、無線信号の短時間フーリエ変換結果を示すフーリエ変換信号を出力するフーリエ変換部と、フーリエ変換部から出力されたフーリエ変換信号のスペクトログラムにおけるエネルギーの時間波形を算出し、エネルギーの時間波形から無線信号の立ち上がり時刻を検出する立ち上がり検出部と、フーリエ変換部から出力されたフーリエ変換信号の中から、立ち上がり検出部により検出された立ち上がり時刻を含む時間帯のフーリエ変換信号を抽出し、抽出したフーリエ変換信号に含まれている設定周波数成分の複素時間信号から、抽出したフーリエ変換信号の変動を特徴量として算出する特徴量算出部とを設け、無線機識別部が、特徴量算出部により算出された特徴量に基づいて、無線機を識別するようにしたものである。 A radio identification device according to the present invention includes a signal receiving unit that receives a radio signal transmitted from a radio to be identified, and a short-time Fourier transform of the radio signal received by the signal receiving unit, thereby A Fourier transform unit that outputs a Fourier transform signal indicating a Fourier transform result, and calculates a time waveform of energy in a spectrogram of the Fourier transform signal output from the Fourier transform unit, and detects a rise time of the radio signal from the time waveform of the energy The Fourier transform signal in the time zone including the rise time detected by the rise detection unit is extracted from the rise detection unit and the Fourier transform signal output from the Fourier transform unit, and is included in the extracted Fourier transform signal. Fluctuations in the Fourier transform signal extracted from the complex time signal of the set frequency component are used as feature quantities. A feature quantity calculating unit that calculates provided, the wireless device identification unit, based on the feature amount calculated by the feature calculation unit, in which so as to identify the radio.

この発明によれば、フーリエ変換部から出力されたフーリエ変換信号のスペクトログラムにおけるエネルギーの時間波形を算出し、エネルギーの時間波形から無線信号の立ち上がり時刻を検出する立ち上がり検出部と、フーリエ変換部から出力されたフーリエ変換信号の中から、立ち上がり検出部により検出された立ち上がり時刻を含む時間帯のフーリエ変換信号を抽出し、抽出したフーリエ変換信号に含まれている設定周波数成分の複素時間信号から、抽出したフーリエ変換信号の変動を特徴量として算出する特徴量算出部とを設け、無線機識別部が、特徴量算出部により算出された特徴量に基づいて、無線機を識別するように構成したので、無線機の個体を識別することができる効果がある。 According to this invention, the time waveform of energy in the spectrogram of the Fourier transform signal output from the Fourier transform unit is calculated, the rise detection unit that detects the rise time of the radio signal from the time waveform of energy, and the output from the Fourier transform unit The Fourier transform signal in the time zone that includes the rise time detected by the rise detection unit is extracted from the Fourier transform signal that has been detected, and extracted from the complex time signal of the set frequency component contained in the extracted Fourier transform signal And a feature amount calculation unit that calculates a variation of the Fourier transform signal as a feature amount, and the wireless device identification unit is configured to identify the wireless device based on the feature amount calculated by the feature amount calculation unit. There is an effect that an individual wireless device can be identified.

この発明の実施の形態１による無線機識別装置を示す構成図である。It is a block diagram which shows the radio | wireless machine identification device by Embodiment 1 of this invention. この発明の実施の形態１による無線機識別装置を示すハードウェア構成図である。It is a hardware block diagram which shows the radio | wireless machine identification device by Embodiment 1 of this invention. 無線機識別装置における信号受信部１を除く構成要素がソフトウェア又はファームウェアなどで実現される場合のコンピュータのハードウェア構成図である。It is a hardware block diagram of a computer in case the components except the signal receiving part 1 in a radio | wireless machine identification device are implement | achieved by software or firmware. この発明の実施の形態１による無線機識別方法における学習処理を示すフローチャートである。It is a flowchart which shows the learning process in the radio | wireless machine identification method by Embodiment 1 of this invention. フーリエ変換部２による学習用無線信号ｓ_ｇ（ｎ）及び識別用無線信号ｓ_ｄ（ｎ）のＳＴＦＴを示す説明図である。It is explanatory drawing which shows STFT of the radio signal for learning _sg (n) by the Fourier-transform part 2, and the radio signal for identification _sd (n). 立ち上がり検出部３による立ち上がり時刻Ｔ_０の検出処理を示す説明図である。It is an explanatory view showing the process of detecting the rising time T ₀ by the rising edge detection unit 3. 特徴量算出部７による特徴量の算出処理を示す説明図である。It is explanatory drawing which shows the calculation process of the feature-value by the feature-value calculation part. 図８Ａは、特徴量算出用のパラメータが適正でないために、種類（Ａ）の無線機に係る特徴量ベクトル〇と、種類（Ｂ）の無線機に係る特徴量ベクトル◇と、種類（Ｃ）の無線機に係る特徴量ベクトル△とが一部重なっている例を示す説明図、図８Ｂは、特徴量算出用のパラメータが適正であるために、種類（Ａ）の無線機に係る特徴量ベクトル〇と、種類（Ｂ）の無線機に係る特徴量ベクトル◇と、種類（Ｃ）の無線機に係る特徴量ベクトル△とが重なっていない例を示す説明図である。FIG. 8A shows that since the parameter for calculating the feature quantity is not appropriate, the feature quantity vector 0 related to the type (A) radio, the feature quantity vector ◇ related to the type (B) radio, and the type (C) FIG. 8B is an explanatory diagram showing an example in which the feature quantity vector Δ related to the other radio device partially overlaps, and FIG. 8B shows the feature quantity related to the type (A) radio device because the parameter for calculating the feature value is appropriate. It is explanatory drawing which shows the example in which the vector (circle), the feature-value vector ◇ which concerns on the type (B) radio | wireless machine, and the feature-value vector (DELTA) which concerns on the type (C) radio | wireless machine do not overlap. 第１のデータベース部９に格納されているシリアルナンバーが付加されている学習用無線信号を示す説明図である。It is explanatory drawing which shows the radio signal for learning to which the serial number stored in the 1st database part 9 is added. 第２のデータベース部１０に格納される複数の既知の無線機に係る特徴量の一例を示す説明図である。FIG. 6 is an explanatory diagram illustrating an example of feature amounts related to a plurality of known wireless devices stored in a second database unit 10; この発明の実施の形態１による無線機識別方法における識別処理を示すフローチャートである。It is a flowchart which shows the identification process in the radio | wireless machine identification method by Embodiment 1 of this invention. 無線機識別部１１による無線機の識別処理を示す説明図である。It is explanatory drawing which shows the identification process of the radio | wireless machine by the radio | wireless machine identification part.

以下、この発明をより詳細に説明するために、この発明を実施するための形態について、添付の図面に従って説明する。 Hereinafter, in order to explain the present invention in more detail, modes for carrying out the present invention will be described with reference to the accompanying drawings.

実施の形態１．
図１は、この発明の実施の形態１による無線機識別装置を示す構成図である。
図２は、この発明の実施の形態１による無線機識別装置を示すハードウェア構成図である。
図１及び図２において、信号受信部１は、例えば図２に示す信号受信回路２１で実現される。
信号受信部１は、既知の無線機から送信された無線信号を学習用無線信号として受信し、受信した学習用無線信号をアナログ信号からデジタル信号に変換することで、デジタルの学習用無線信号をフーリエ変換部２及び第１のデータベース部９に出力する。
また、信号受信部１は、識別対象の無線機から送信された無線信号を識別用無線信号として受信し、受信した識別用無線信号をアナログ信号からデジタル信号に変換することで、デジタルの識別用無線信号をフーリエ変換部２に出力する。Embodiment 1 FIG.
FIG. 1 is a block diagram showing a radio identification device according to Embodiment 1 of the present invention.
FIG. 2 is a hardware configuration diagram showing the radio identification device according to Embodiment 1 of the present invention.
1 and 2, the signal receiving unit 1 is realized by, for example, a signal receiving circuit 21 shown in FIG.
The signal receiving unit 1 receives a radio signal transmitted from a known radio as a learning radio signal, and converts the received learning radio signal from an analog signal to a digital signal, thereby converting the digital learning radio signal into a digital signal. The data is output to the Fourier transform unit 2 and the first database unit 9.
The signal receiving unit 1 receives a radio signal transmitted from a radio device to be identified as an identification radio signal, and converts the received identification radio signal from an analog signal to a digital signal, thereby enabling digital identification. The radio signal is output to the Fourier transform unit 2.

フーリエ変換部２は、例えば図２に示すフーリエ変換回路２２で実現される。
フーリエ変換部２は、信号受信部１から出力されたデジタルの学習用無線信号を短時間フーリエ変換（ＳＴＦＴ：ＳｈｏｒｔＴｉｍｅＦｏｕｒｉｅｒＴｒａｎｓｆｏｒｍ）し、学習用無線信号の短時間フーリエ変換結果であるＳＴＦＴ結果を示すフーリエ変換信号（以下、学習用フーリエ変換信号と称する）を立ち上がり検出部３及び特徴量算出部７に出力する処理を実施する。
また、フーリエ変換部２は、信号受信部１から出力されたデジタルの識別用無線信号をＳＴＦＴし、識別用無線信号のＳＴＦＴ結果を示すフーリエ変換信号（以下、識別用フーリエ変換信号と称する）を立ち上がり検出部３及び特徴量算出部７に出力する処理を実施する。
ＳＴＦＴは、一定時間毎に異なる時刻の無線信号をそれぞれ切り出し、切り出した各々の無線信号をそれぞれ高速フーリエ変換（ＦＦＴ：ＦａｓｔＦｏｕｒｉｅｒＴｒａｎｓｆｏｒｍ）する手法である。The Fourier transform unit 2 is realized by, for example, a Fourier transform circuit 22 shown in FIG.
The Fourier transform unit 2 performs a short time Fourier transform (STFT) on the digital learning radio signal output from the signal receiving unit 1 and outputs an STFT result that is a short time Fourier transform result of the learning radio signal. A process of outputting the illustrated Fourier transform signal (hereinafter referred to as a learning Fourier transform signal) to the rise detection unit 3 and the feature amount calculation unit 7 is performed.
Further, the Fourier transform unit 2 performs STFT on the digital identification radio signal output from the signal reception unit 1 and generates a Fourier transform signal (hereinafter referred to as identification Fourier transform signal) indicating the STFT result of the identification radio signal. A process of outputting to the rise detection unit 3 and the feature amount calculation unit 7 is performed.
The STFT is a technique in which radio signals at different times are cut out at regular intervals, and each cut out radio signal is subjected to fast Fourier transform (FFT).

立ち上がり検出部３は、例えば図２に示す立ち上がり検出回路２３で実現される。
立ち上がり検出部３は、フーリエ変換部２から出力された学習用フーリエ変換信号のスペクトログラムにおけるエネルギーの時間波形（以下、学習用のエネルギー時間波形と称する）を算出し、学習用のエネルギー時間波形から学習用無線信号における立ち上がり時刻を検出する処理を実施する。
また、立ち上がり検出部３は、フーリエ変換部２から出力された識別用フーリエ変換信号のスペクトログラムにおけるエネルギーの時間波形（以下、識別用のエネルギー時間波形と称する）を算出し、識別用のエネルギー時間波形から識別用無線信号における立ち上がり時刻を検出する処理を実施する。The rise detection unit 3 is realized by, for example, a rise detection circuit 23 shown in FIG.
The rise detection unit 3 calculates a time waveform of energy (hereinafter referred to as a learning energy time waveform) in a spectrogram of the learning Fourier transform signal output from the Fourier transform unit 2, and learns from the learning energy time waveform. The process which detects the rise time in the radio signal for operation is performed.
The rise detection unit 3 calculates a time waveform of energy (hereinafter referred to as an identification energy time waveform) in a spectrogram of the identification Fourier transform signal output from the Fourier transform unit 2, and performs an identification energy time waveform. To detect the rise time in the identification radio signal.

パラメータ設定部４は、第１のパラメータ設定部５及び第２のパラメータ設定部６を備えており、例えば図２に示すパラメータ設定回路２４で実現される。
第１のパラメータ設定部５は、特徴量算出用のパラメータとして、立ち上がり検出部３により検出された立ち上がり時刻を含む時間帯の長さである時間範囲Δｔを示すパラメータと、設定周波数成分である設定周波数ビンを示すパラメータとを設定する処理を実施する。
第２のパラメータ設定部６は、特徴量算出部７により算出される第１から第９の特徴量の重みを示す重み付けパラメータを設定する処理を実施する。
また、第２のパラメータ設定部６は、フーリエ変換部２によりＳＴＦＴが行われる際に用いられる窓関数のサンプル数及び窓関数をシフトさせるサンプル数を設定する処理を実施する。The parameter setting unit 4 includes a first parameter setting unit 5 and a second parameter setting unit 6, and is realized by, for example, a parameter setting circuit 24 shown in FIG.
The first parameter setting unit 5 includes a parameter indicating a time range Δt that is a length of a time zone including a rising time detected by the rising detection unit 3 and a setting frequency component as parameters for calculating a feature amount. A process for setting a parameter indicating a frequency bin is performed.
The second parameter setting unit 6 performs a process of setting weighting parameters indicating the weights of the first to ninth feature amounts calculated by the feature amount calculation unit 7.
In addition, the second parameter setting unit 6 performs a process of setting the number of samples of the window function and the number of samples for shifting the window function used when the STFT is performed by the Fourier transform unit 2.

特徴量算出部７は、例えば図２に示す特徴量算出回路２５で実現される。
特徴量算出部７は、パラメータ設定部４により設定された特徴量算出用のパラメータから、立ち上がり時刻を含む時間帯及び設定周波数ビンを認識する。
特徴量算出部７は、フーリエ変換部２から出力された学習用フーリエ変換信号の中から、立ち上がり検出部３により検出された学習用無線信号における立ち上がり時刻を含む時間帯の学習用フーリエ変換信号を抽出する処理を実施する。
そして、特徴量算出部７は、抽出した学習用フーリエ変換信号に含まれている設定周波数ビンの複素時間信号から、抽出した学習用フーリエ変換信号の変動を特徴量として算出する処理を実施する。The feature amount calculation unit 7 is realized by, for example, a feature amount calculation circuit 25 shown in FIG.
The feature amount calculation unit 7 recognizes the time zone including the rise time and the set frequency bin from the feature amount calculation parameters set by the parameter setting unit 4.
The feature amount calculation unit 7 calculates the learning Fourier transform signal in the time zone including the rise time in the learning radio signal detected by the rise detection unit 3 from the learning Fourier transform signal output from the Fourier transform unit 2. Perform the extraction process.
Then, the feature amount calculation unit 7 performs a process of calculating the variation of the extracted learning Fourier transform signal as a feature amount from the complex time signal of the set frequency bin included in the extracted learning Fourier transform signal.

また、特徴量算出部７は、フーリエ変換部２から出力された識別用フーリエ変換信号の中から、立ち上がり検出部３により算出された識別用無線信号における立ち上がり時刻を含む時間帯の識別用フーリエ変換信号を抽出する処理を実施する。
そして、特徴量算出部７は、抽出した識別用フーリエ変換信号に含まれている設定周波数成分の複素時間信号から、抽出した識別用フーリエ変換信号の変動を特徴量として算出する処理を実施する。Further, the feature amount calculation unit 7 performs an identification Fourier transform for a time zone including the rising time in the identification wireless signal calculated by the rising detection unit 3 from among the identification Fourier transform signals output from the Fourier transform unit 2. A process of extracting a signal is performed.
Then, the feature amount calculation unit 7 performs a process of calculating, as a feature amount, the variation of the extracted identification Fourier transform signal from the complex time signal of the set frequency component included in the extracted identification Fourier transform signal.

データベース部８は、第１のデータベース部９及び第２のデータベース部１０を備えており、例えば図２に示すデータベース回路２６で実現される。
第１のデータベース部９は、信号受信部１から出力されたデジタルの学習用無線信号を格納する。
第２のデータベース部１０は、既知の無線機に係る特徴量として、特徴量算出部７により算出された学習用フーリエ変換信号の変動を格納する。The database unit 8 includes a first database unit 9 and a second database unit 10, and is realized by, for example, the database circuit 26 shown in FIG.
The first database unit 9 stores the digital learning radio signal output from the signal receiving unit 1.
The second database unit 10 stores the variation of the learning Fourier transform signal calculated by the feature amount calculation unit 7 as a feature amount related to a known wireless device.

無線機識別部１１は、例えば図２に示す無線機識別回路２７で実現される。
無線機識別部１１は、特徴量算出部７により算出された識別対象の無線機に係る特徴量と、第２のデータベース部１０により格納されている既知の無線機に係る特徴量とを比較する処理を実施する。
無線機識別部１１は、識別対象の無線機に係る特徴量と既知の無線機に係る特徴量との比較結果に基づいて、識別対象の無線機を識別する処理を実施する。The radio identification unit 11 is realized by, for example, a radio identification circuit 27 illustrated in FIG.
The wireless device identification unit 11 compares the feature value related to the identification target wireless device calculated by the feature value calculation unit 7 with the feature value related to the known wireless device stored by the second database unit 10. Perform the process.
The wireless device identification unit 11 performs a process of identifying the wireless device to be identified based on the comparison result between the feature value related to the wireless device to be identified and the feature value related to the known wireless device.

図１では、無線機識別装置の構成要素である信号受信部１、フーリエ変換部２、立ち上がり検出部３、パラメータ設定部４、特徴量算出部７、データベース部８及び無線機識別部１１のそれぞれが、図２に示すような専用のハードウェアで実現されるものを想定している。即ち、信号受信回路２１、フーリエ変換回路２２、立ち上がり検出回路２３、パラメータ設定回路２４、特徴量算出回路２５、データベース回路２６及び無線機識別回路２７で実現されるものを想定している。 In FIG. 1, each of a signal receiving unit 1, a Fourier transform unit 2, a rise detection unit 3, a parameter setting unit 4, a feature amount calculation unit 7, a database unit 8, and a radio unit identification unit 11, which are components of the radio unit identification device. However, it is assumed to be realized by dedicated hardware as shown in FIG. That is, it is assumed that the signal receiving circuit 21, the Fourier transform circuit 22, the rising edge detection circuit 23, the parameter setting circuit 24, the feature amount calculation circuit 25, the database circuit 26, and the wireless device identification circuit 27 are realized.

ここで、データベース回路２６は、例えば、ＲＡＭ（ＲａｎｄｏｍＡｃｃｅｓｓＭｅｍｏｒｙ）、ＲＯＭ（ＲｅａｄＯｎｌｙＭｅｍｏｒｙ）、フラッシュメモリ、ＥＰＲＯＭ（ＥｒａｓａｂｌｅＰｒｏｇｒａｍｍａｂｌｅＲｅａｄＯｎｌｙＭｅｍｏｒｙ）、ＥＥＰＲＯＭ（ＥｌｅｃｔｒｉｃａｌｌｙＥｒａｓａｂｌｅＰｒｏｇｒａｍｍａｂｌｅＲｅａｄＯｎｌｙＭｅｍｏｒy）などの不揮発性又は揮発性の半導体メモリ、磁気ディスク、フレキシブルディスク、光ディスク、コンパクトディスク、ミニディスク、あるいは、ＤＶＤ（ＤｉｇｉｔａｌＶｅｒｓａｔｉｌｅＤｉｓｃ）が該当する。
また、信号受信回路２１、フーリエ変換回路２２、立ち上がり検出回路２３、パラメータ設定回路２４、特徴量算出回路２５及び無線機識別回路２７は、例えば、単一回路、複合回路、プログラム化したプロセッサ、並列プログラム化したプロセッサ、ＡＳＩＣ（ＡｐｐｌｉｃａｔｉｏｎＳｐｅｃｉｆｉｃＩｎｔｅｇｒａｔｅｄＣｉｒｃｕｉｔ）、ＦＰＧＡ（Ｆｉｅｌｄ−ＰｒｏｇｒａｍｍａｂｌｅＧａｔｅＡｒｒａｙ）、または、これらを組み合わせたものが該当する。Here, the database circuit 26 includes, for example, a RAM (Random Access Memory), a ROM (Read Only Memory), a flash memory, an EPROM (Erasable Programmable Read Only Memory), or an EEPROM (Electrically Erasable Memory). A volatile semiconductor memory, a magnetic disk, a flexible disk, an optical disk, a compact disk, a mini disk, or a DVD (Digital Versatile Disc) is applicable.
The signal receiving circuit 21, the Fourier transform circuit 22, the rising edge detection circuit 23, the parameter setting circuit 24, the feature amount calculation circuit 25, and the wireless device identification circuit 27 are, for example, a single circuit, a composite circuit, a programmed processor, a parallel processor A programmed processor, an application specific integrated circuit (ASIC), a field-programmable gate array (FPGA), or a combination thereof is applicable.

無線機識別装置における信号受信部１を除く構成要素は、専用のハードウェアで実現されるものに限るものではなく、無線機識別装置における信号受信部１を除く構成要素が、ソフトウェア、ファームウェア、または、ソフトウェアとファームウェアとの組み合わせで実現されるものであってもよい。
ソフトウェア又はファームウェアはプログラムとして、コンピュータのメモリに格納される。コンピュータは、プログラムを実行するハードウェアを意味し、例えば、ＣＰＵ（ＣｅｎｔｒａｌＰｒｏｃｅｓｓｉｎｇＵｎｉｔ）、中央処理装置、処理装置、演算装置、マイクロプロセッサ、マイクロコンピュータ、プロセッサ、あるいは、ＤＳＰ（ＤｉｇｉｔａｌＳｉｇｎａｌＰｒｏｃｅｓｓｏｒ）が該当する。The components other than the signal receiving unit 1 in the radio identification device are not limited to those realized by dedicated hardware, and the components other than the signal receiving unit 1 in the radio identification device are software, firmware, or It may be realized by a combination of software and firmware.
Software or firmware is stored as a program in the memory of a computer. The computer means hardware that executes a program, for example, a CPU (Central Processing Unit), a central processing unit, a processing unit, a processing unit, a microprocessor, a microcomputer, a processor, or a DSP (Digital Signal Processor) To do.

図３は、無線機識別装置における信号受信部１を除く構成要素がソフトウェア又はファームウェアなどで実現される場合のコンピュータのハードウェア構成図である。
無線機識別装置における信号受信部１を除く構成要素がソフトウェア又はファームウェアなどで実現される場合、データベース部８をコンピュータのメモリ３１上に構成するとともに、フーリエ変換部２、立ち上がり検出部３、パラメータ設定部４、特徴量算出部７及び無線機識別部１１の処理手順をコンピュータに実行させるためのプログラムをメモリ３１に格納し、コンピュータのプロセッサ３２がメモリ３１に格納されているプログラムを実行するようにすればよい。FIG. 3 is a hardware configuration diagram of a computer when components other than the signal receiving unit 1 in the wireless device identification apparatus are realized by software or firmware.
When components other than the signal receiving unit 1 in the radio identification device are realized by software or firmware, the database unit 8 is configured on the memory 31 of the computer, and the Fourier transform unit 2, the rise detection unit 3, and parameter setting The program for causing the computer to execute the processing procedure of the unit 4, the feature amount calculation unit 7, and the wireless device identification unit 11 is stored in the memory 31, and the processor 32 of the computer executes the program stored in the memory 31. do it.

また、図２では、無線機識別装置の構成要素のそれぞれが専用のハードウェアで実現される例を示し、図３では、無線機識別装置がソフトウェアやファームウェアなどで実現される例を示しているが、無線機識別装置における一部の構成要素が専用のハードウェアで実現され、残りの構成要素がソフトウェアやファームウェアなどで実現されるものであってもよい。 2 shows an example in which each component of the radio identification device is realized by dedicated hardware, and FIG. 3 shows an example in which the radio identification device is realized by software, firmware, or the like. However, some components in the wireless device identification device may be realized by dedicated hardware, and the remaining components may be realized by software, firmware, or the like.

次に動作について説明する。
＜学習時の処理＞
最初に、既知の無線機から送信された無線信号を学習する際の処理内容を説明する。
図４は、この発明の実施の形態１による無線機識別方法における学習処理を示すフローチャートである。
信号受信部１は、既知の無線機から送信された無線信号を学習用無線信号として受信し、受信した学習用無線信号をアナログ信号からデジタル信号に変換する（図４のステップＳＴ１）。
学習用無線信号には、既知の無線機を識別する個体情報であるシリアルナンバーが付加されているものとする。
信号受信部１は、デジタルの学習用無線信号ｓ_ｇ（ｎ）をフーリエ変換部２及び第１のデータベース部９に出力する。ｓ_ｇ（ｎ）におけるｎは、学習用無線信号のサンプル番号である。例えば、図９には、サンプル番号ｎが１〜４の学習用無線信号が示されている。
第１のデータベース部９には、シリアルナンバーが付加されている学習用無線信号ｓ_ｇ（ｎ）が格納される。
図９は、第１のデータベース部９に格納されているシリアルナンバーが付加されている学習用無線信号を示す説明図である。
図９の例では、シリアルナンバーが同じ無線機であっても、受信時のＳＮＲ又はサンプリング周波数などの受信条件が異なっている複数の学習用無線信号が格納されている。Next, the operation will be described.
<Process during learning>
First, processing contents when learning a wireless signal transmitted from a known wireless device will be described.
FIG. 4 is a flowchart showing a learning process in the wireless device identification method according to Embodiment 1 of the present invention.
The signal receiving unit 1 receives a radio signal transmitted from a known radio as a learning radio signal, and converts the received learning radio signal from an analog signal to a digital signal (step ST1 in FIG. 4).
It is assumed that a serial number, which is individual information for identifying a known radio device, is added to the learning radio signal.
The signal receiving unit 1 outputs the digital learning radio signal s _g (n) to the Fourier transform unit 2 and the first database unit 9. n in s _g (n) is a sample number of the learning radio signal. For example, FIG. 9 shows learning radio signals whose sample numbers n are 1 to 4.
The first database unit 9 stores a learning radio signal s _g (n) to which a serial number is added.
FIG. 9 is an explanatory diagram showing the learning radio signal to which the serial number stored in the first database unit 9 is added.
In the example of FIG. 9, a plurality of learning radio signals having different reception conditions such as SNR or sampling frequency at the time of reception are stored even if the radios have the same serial number.

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

As shown in the following formula (1), the Fourier transform unit 2 performs STFT on the digital learning radio signal s _g (n) output from the signal receiving unit 1 to obtain the learning radio signal s _g (n). The learning Fourier transform signal S _g (m, k), which is a Fourier transform signal indicating the STFT result, is output to the rise detection unit 3 and the feature amount calculation unit 7 (step ST2 in FIG. 4).

式（１）において、ｗ（ｎ−ｍ）は、窓関数であり、例えば、矩形窓又はハミング窓が考えられる。
ｍは、窓関数ｗ（ｎ−ｍ）をシフトさせるサンプル数、ｍ_０は、窓関数ｗ（ｎ−ｍ）の適用区間の初めのサンプル番号、Ｍは、窓関数ｗ（ｎ−ｍ）のサンプル数、ｋは、周波数ビン番号である。
なお、窓関数ｗ（ｎ−ｍ）のサンプル数Ｍ及び窓関数ｗ（ｎ−ｍ）をシフトさせるサンプル数ｍのそれぞれは、第２のパラメータ設定部６によって事前に設定される。In Expression (1), w (n−m) is a window function, and for example, a rectangular window or a Hamming window can be considered.
m is the number of samples to shift the window function w (n−m), m ₀ is the first sample number of the application section of the window function w (n−m), and M is the window function w (n−m). The number of samples, k is a frequency bin number.
Note that the number of samples M of the window function w (nm) and the number of samples m of shifting the window function w (nm) are set in advance by the second parameter setting unit 6.

ここで、図５は、フーリエ変換部２による学習用無線信号ｓ_ｇ（ｎ）及び識別用無線信号ｓ_ｄ（ｎ）のＳＴＦＴを示す説明図である。
学習用無線信号ｓ_ｇ（ｎ）及び識別用無線信号ｓ_ｄ（ｎ）は、時間と振幅との関係を示す２次元信号である。識別用無線信号ｓ_ｄ（ｎ）については後述する。
ＳＴＦＴ結果を示す学習用フーリエ変換信号Ｓ_ｇ（ｍ，ｋ）及び識別用フーリエ変換信号Ｓ_ｄ（ｍ，ｋ）は、時間と周波数及び電力との関係を示す３次元信号である。識別用フーリエ変換信号Ｓ_ｄ（ｍ，ｋ）については後述する。
フーリエ変換部２により学習用無線信号ｓ_ｇ（ｎ）及び識別用無線信号ｓ_ｄ（ｎ）のそれぞれがＳＴＦＴされることで、学習用無線信号ｓ_ｇ（ｎ）及び識別用無線信号ｓ_ｄ（ｎ）のそれぞれが、コヒーレント積分されるので、ＳＮＲが改善される効果が得られる。Here, FIG. 5 is an explanatory diagram showing STFTs of the learning radio signal s _g (n) and the identification radio signal s _d (n) by the Fourier transform unit 2.
The learning radio signal s _g (n) and the identification radio signal s _d (n) are two-dimensional signals indicating the relationship between time and amplitude. The identification radio signal s _d (n) will be described later.
The learning Fourier transform signal S _g (m, k) indicating the STFT result and the identification Fourier transform signal S _d (m, k) are three-dimensional signals indicating the relationship between time, frequency, and power. The Fourier transform signal S _d (m, k) for identification will be described later.
Each of the learning radio signal s _g (n) and the identification radio signal s _d (n) is subjected to STFT by the Fourier transform unit 2, so that the learning radio signal s _g (n) and the identification radio signal s _d ( Since each of n) is coherently integrated, the effect of improving the SNR is obtained.

立ち上がり検出部３は、学習用無線信号ｓ_ｇ（ｎ）の立ち上がり時刻Ｔ_０を検出する処理を実施する（図４のステップＳＴ３）。
図６は、立ち上がり検出部３による立ち上がり時刻Ｔ_０の検出処理を示す説明図である。
以下、立ち上がり検出部３による立ち上がり時刻Ｔ_０の検出処理を具体的に説明する。The rise detection unit 3 performs a process of detecting the rise time T ₀ of the learning radio signal s _g (n) (step ST3 in FIG. 4).
FIG. 6 is an explanatory diagram showing detection processing of the rising time T ₀ by the rising detection unit 3.
Hereinafter, the detection process of the rising time T ₀ by the rising detection unit 3 will be specifically described.

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

First, as shown in the following formula (2), the rising edge detection unit 3 squares the learning Fourier transform signal S _g (m, k) output from the Fourier transform unit 2 | S _g (m, k ) | ² is calculated as spectrogram SP _g (m, k).

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

式（３）において、ｋ_１は、スペクトログラムＳＰ_ｇ（ｍ，ｋ）に含まれている複数の周波数のうち、電力を合計する周波数の下限周波数を示す周波数ビン番号、ｋ_ｎは、スペクトログラムＳＰ_ｇ（ｍ，ｋ）に含まれている複数の周波数のうち、電力を合計する周波数の上限周波数を示す周波数ビン番号である。
ｋ＝ｋ_１，・・・，ｋ_ｎである。Next, the rising edge detection unit 3 calculates a learning energy time waveform E _g (m), which is a time waveform of energy in the spectrogram SP _g (m, k), as shown in the following equation (3).
The energy time waveform E _g (m) for learning is calculated by summing powers at the same time among powers at a plurality of frequencies included in the spectrogram SP _g (m, k). This is the total value of power at the time.

In the formula (3), _{k 1,} among a plurality of frequencies contained in the spectrogram _SP g (m, k), the frequency bin number that indicates the lower limit frequency of the frequency summing power, _{k n} is spectrogram SP _g It is a frequency bin number indicating the upper limit frequency of the frequency for summing power among a plurality of frequencies included in (m, k).
k ₌ k 1, ···, a _{k n.}

次に、立ち上がり検出部３は、学習用のエネルギー時間波形Ｅ_ｇ（ｍ）と閾値Ｅ_ｔｈとを比較する。
閾値Ｅ_ｔｈとしては、例えば、学習用のエネルギー時間波形Ｅ_ｇ（ｍ）の最大値から設定値Ｅ_０を減算した値が考えられる。
立ち上がり検出部３は、学習用のエネルギー時間波形Ｅ_ｇ（ｍ）と閾値Ｅ_ｔｈとの比較結果を参照することで、最初は、ノイズレベルのエネルギーを示していた学習用のエネルギー時間波形Ｅ_ｇ（ｍ）が、閾値Ｅ_ｔｈに到達したか否かを判定する。
立ち上がり検出部３は、学習用のエネルギー時間波形Ｅ_ｇ（ｍ）が閾値Ｅ_ｔｈに到達した時刻を、学習用無線信号ｓ_ｇ（ｎ）の立ち上がり時刻Ｔ_０として検出する。
この実施の形態１では、立ち上がり検出部３が、学習用無線信号ｓ_ｇ（ｎ）の立ち上がり時刻Ｔ_０を検出する際、スペクトログラムＳＰ_ｇ（ｍ，ｋ）に含まれている複数の周波数における電力のうち、同一時刻の電力同士をそれぞれ合計するようにしている。これにより、立ち上がり時に電力が大きく変動する周波数が含まれている場合でも、当該電力の変動の影響が低減されるため、立ち上がり時刻Ｔ_０の検出精度が向上する。Next, the rising detection unit 3 compares the energy time waveform E _g (m) for learning with the threshold value E _th .
As the threshold value E _th , for example, a value obtained by subtracting the set value E ₀ from the maximum value of the energy time waveform E _g (m) for learning can be considered.
The rising edge detection unit 3 refers to the comparison result between the learning energy time waveform E _g (m) and the threshold value E _th , so that the learning energy time waveform E _g that initially indicates the noise level energy is used. It is determined whether (m) has reached the threshold value _Eth .
The rise detection unit 3 detects the time when the learning energy time waveform E _g (m) reaches the threshold E _th as the rise time T ₀ of the learning radio signal s _g (n).
In the first embodiment, when the rise detection unit 3 detects the rise time T ₀ of the learning radio signal s _g (n), the power at a plurality of frequencies included in the spectrogram SP _g (m, k). Among them, the powers at the same time are summed up. Thus, even if it contains frequencies vary greatly power at the rise time, the influence of the variation of the power is reduced, thereby improving the detection accuracy of the rise time T _0.

特徴量算出部７は、第１のパラメータ設定部５により設定された特徴量算出用のパラメータから、立ち上がり時刻Ｔ_０を含む時間帯（ｔ_１〜ｔ_２）及び設定周波数ビンｋ_ｂを認識する。
第１のパラメータ設定部５による特徴量算出用のパラメータの設定処理については後述する。
次に、特徴量算出部７は、図６に示すように、フーリエ変換部２から出力された学習用フーリエ変換信号Ｓ_ｇ（ｍ，ｋ）の中から、立ち上がり検出部３により検出された学習用無線信号ｓ_ｇ（ｎ）の立ち上がり時刻Ｔ_０を含む時間帯（ｔ_１〜ｔ_２）の学習用フーリエ変換信号ｓ_ｇ（ｎ）_{ｔ１〜ｔ２}を抽出する。Feature amount calculation unit 7 recognizes from the parameter for the set characteristic quantity calculation, a time zone including the rise time T ₀ the (t 1 _~t ₂₎ and the set frequency bin k _b of the first parameter setting unit 5 .
The parameter setting process for calculating the feature amount by the first parameter setting unit 5 will be described later.
Next, as shown in FIG. 6, the feature amount calculation unit 7 learns the learning detected by the rising detection unit 3 from the learning Fourier transform signal S _g (m, k) output from the Fourier transform unit 2. The learning Fourier transform signals s _g (n) _t1 to _t2 in the time zone (t _{1 to} t ₂ ) including the rising time T ₀ of the radio signal s _g (n) for use are extracted.

次に、特徴量算出部７は、図７に示すように、抽出した学習用フーリエ変換信号ｓ_ｇ（ｎ）_{ｔ１〜ｔ２}に含まれている複素時間信号の中から、設定周波数ビンｋ_ｂの複素時間信号ｓ_ｇハット（ｎ）_{ｔ１〜ｔ２}の切り出しを行う。
電子出願の都合上、明細書の文章中では、ｓの文字の上に“＾”の記号を付することができないので、ｓ_ｇハット（ｎ）_{ｔ１〜ｔ２}のように表記している。
図７は、特徴量算出部７による特徴量の算出処理を示す説明図である。
特徴量算出部７は、切り出した複素時間信号ｓ_ｇハット（ｎ）_{ｔ１〜ｔ２}から、抽出した学習用フーリエ変換信号ｓ_ｇ（ｎ）_{ｔ１〜ｔ２}の変動を特徴量として算出する（図４のステップＳＴ４）。Next, the feature amount calculation unit 7, as shown in FIG. 7, from the complex time signal included in the extracted learning Fourier transform signal _s _{g (n)} t1~t2, the set frequency bin k _b The complex time signal s _g hat (n) _{t1 to t2} is cut out.
For the convenience of electronic filing, in the text of the specification, since the symbol “^” cannot be put on the letter s, it is represented as s _g hat (n) _{t1 to t2} .
FIG. 7 is an explanatory diagram illustrating a feature amount calculation process by the feature amount calculation unit 7.
The feature amount calculation unit 7 calculates, as a feature amount, a variation in the extracted learning Fourier transform signal s _g (n) _t1 to _t2 from the extracted complex time signal s _g hat (n) _t1 to _t2 (FIG. 4). Step ST4).

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

式（４）において、Ｉ_ｇ（ｎ）は、複素時間信号ｓ_ｇハット（ｎ）_{ｔ１〜ｔ２}の実部、Ｑ_ｇ（ｎ）は、複素時間信号ｓ_ｇハット（ｎ）_{ｔ１〜ｔ２}の虚部である。Hereinafter, the feature amount calculation processing by the feature amount calculation unit 7 will be described in detail.
Feature amount calculation unit 7, the complex time signal _{s g} hat _{(n) t1 to t2} cut, as shown in the following equation (4), the complex time signal _{s g} hat _{(n) t1 to t2} of the instantaneous amplitude a _g (n) is calculated.

In Equation (4), I _g (n) is the real part of the complex time signal s _g hat (n) _{t1 to t2} , and Q _g (n) is the imaginary of the complex time signal s _g hat (n) _{t1 to t2} . Part.

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

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

The feature quantity calculation unit 7, the complex time signal _{s g} hat _{(n) t1 to t2} cut, as shown in the following equation (5), the instantaneous complex time signal _{s g} hat _{(n) t1 to t2} The phase φ _g (n) is calculated.

The feature quantity calculation unit 7, the complex time signal _{s g} hat _{(n) t1 to t2} cut, as shown in the following equation (6), the instantaneous complex time signal _{s g} hat _{(n) t1 to t2} The frequency f _g (n) is calculated.

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

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

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

Next, as shown in the following formula (7), the feature amount calculation unit 7 uses the variance of the instantaneous amplitude a _g (n) in the complex time signal s _g hat (n) _{t1 to t2} as the first feature amount. The value _pg1 is calculated.

Feature amount calculation unit 7, as shown in the following equation (9), as a second feature quantity, skewness _{p g2} of the complex time signal _{s g} hat _(n) instantaneous amplitude _a g in _{t1 to t2} (n) Is calculated.

As shown in the following formula (11), the feature amount calculation unit 7 uses, as the third feature amount, the kurtosis _{pg 3} of the instantaneous amplitude a _g (n) in the complex time signal s _g hat (n) _{t1 to t2} . Is calculated.

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

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

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

Next, as shown in the following formula (12), the feature amount calculation unit 7 uses the dispersion of the instantaneous phase φ _g (n) in the complex time signal s _g hat (n) _{t1 to t2} as the fourth feature amount. The value _pg4 is calculated.

As shown in the following formula (14), the feature amount calculation unit 7 uses, as the fifth feature amount, the skewness p _g5 of the instantaneous phase φ _g (n) in the complex time signal s _g hat (n) _{t1 to t2} . Is calculated.

As shown in the following formula (16), the feature quantity calculation unit 7 uses, as the sixth feature quantity, the kurtosis p _g6 of the instantaneous phase φ _g (n) in the complex time signal s _g hat (n) _{t1 to t2} . Is calculated.

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

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

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

Next, as shown in the following equation (17), the feature amount calculation unit 7 uses the variance of the instantaneous frequency f _g (n) in the complex time signal s _g hat (n) _{t1 to t2} as the seventh feature amount. The value _pg7 is calculated.

Feature amount calculation unit 7, as shown in the following equation (19), as the feature amount of the 8, skewness of the complex time signal _{s g} hat _(n) instantaneous frequency _f g of _{t1 to t2} (n) _{p g8} Is calculated.

Feature amount calculation unit 7, as shown in the following equation (21), as the feature quantity of the ninth, kurtosis of the complex time signal _{s g} hat _(n) instantaneous frequency _f g of _{t1 to t2} (n) _{p g9} Is calculated.

特徴量算出部７は、第２のパラメータ設定部６により設定された第１から第９の特徴量ｐ_ｇｉ（ｉ＝１，２，・・・，９）の重みｗ_ｉを示す重み付けパラメータを取得する。
第１から第９の特徴量ｐ_ｇｉの中には、無線機の識別に有効な特徴量と、識別に有効でない特徴量とが存在している。このため、識別に有効な特徴量には大きな重みが設定され、識別に有効でない特徴量には小さな重みが設定される。
例えば、ｗ_ｉ（ｉ＝１，２，・・・，９）は、０以上１以下の値であり、ｗ_１＋ｗ_２＋・・・＋ｗ_９＝１である。
第２のパラメータ設定部６による重み付けパラメータの設定方法は、特に問わないが、例えば、第１から第９の特徴量ｐ_ｇｉの分散値で、第１から第９の特徴量ｐ_ｇｉを規格化し、第１から第９の特徴量ｐ_ｇｉを規格化した値を重みとして用いる方法が考えられる。The feature amount calculation unit 7 sets weighting parameters indicating the weights w _i of the first to ninth feature amounts p _gi (i = 1, 2,..., 9) set by the second parameter setting unit 6. get.
Among the first to ninth feature amounts p _gi , there are feature amounts that are effective for identifying a radio device and feature amounts that are not effective for identification. For this reason, a large weight is set for a feature quantity effective for identification, and a small weight is set for a feature quantity not effective for identification.
For example, w _i (i = 1, 2,..., 9) is a value of 0 or more and 1 or less, and w ₁ + w ₂ +... + W ₉ = 1.
The method for setting the weighting parameter by the second parameter setting unit 6 is not particularly limited. For example, the first to ninth feature values p _gi are normalized by the variance values of the first to ninth feature values p _gi. A method of using the normalized values of the first to ninth feature values p _gi as weights is conceivable.

特徴量算出部７は、以下の式（２２）に示すように、第１から第９の特徴量ｐ_ｇｉと重みｗ_ｉとをそれぞれ乗算することで、第１から第９の特徴量ｐ_ｇｉの重み付けを行う。
ｐ’_ｇｉ＝ｐ_ｇｉ×ｗ_ｉ（２２）
ｉ＝１，２，・・・，９
特徴量算出部７は、重み付け後の第１から第９の特徴量ｐ’_ｇｉを並べたベクトルを、既知の無線機の特徴量ベクトルｐ_ｇとして算出する。
特徴量算出部７は、算出した既知の無線機の特徴量ベクトルｐ_ｇを、学習用無線信号に付加されているシリアルナンバーと一緒に、第２のデータベース部１０に格納する（図４のステップＳＴ５）。
また、特徴量算出部７は、既知の無線機の特徴量ベクトルｐ_ｇを第２のデータベース部１０に格納した旨を示す情報を第１のパラメータ設定部５に出力する。The feature quantity calculation unit 7 multiplies the first to ninth feature quantities p _gi and the weights w _i , respectively, as shown in the following formula (22), thereby _obtaining the first to ninth feature quantities p _gi. Weighting is performed.
p ′ _gi = p _gi × w _i (22)
i = 1, 2,..., 9
Feature amount calculation unit 7, a vector from the first the weighted arranged feature quantities p _'gi ninth, calculated as a feature vector p _g of the known radio.
Feature amount calculation unit 7, a feature quantity vector p _g of the calculated known radio, together with the serial number is added to the learning radio signal is stored in the second database section 10 (step of FIG. 4 ST5).
The feature quantity calculation unit 7 outputs information indicating that the feature vector p _g of the known radio stored in a second database section 10 to the first parameter setting unit 5.

ここで、第１のパラメータ設定部５による特徴量算出用のパラメータの設定処理について説明する。
図８は、特徴量算出用のパラメータと複数の既知の無線機に係る特徴量ベクトルとの関係を示す説明図である。
図８の例では、〇は、種類（Ａ）の既知の無線機に係る特徴量ベクトルｐ_ｇ、◇は、種類（Ｂ）の既知の無線機に係る特徴量ベクトルｐ_ｇ、△は、種類（Ｃ）の既知の無線機に係る特徴量ベクトルｐ_ｇである。
図８の例では、種類（Ａ）に属する無線機の個体の数が１０であるため、１０個の特徴量ベクトル〇が得られている。
また、種類（Ｂ）に属する無線機の個体の数が１０であるため、１０個の特徴量ベクトル◇が得られており、種類（Ｃ）に属する無線機の個体の数が１０であるため、１０個の特徴量ベクトル△が得られている。Here, the parameter setting process for calculating the feature amount by the first parameter setting unit 5 will be described.
FIG. 8 is an explanatory diagram showing the relationship between the feature quantity calculation parameters and the feature quantity vectors related to a plurality of known wireless devices.
In the example of FIG. 8, ◯ is a feature vector p _g related to a known type radio (A), ◇ is a feature vector p _g related to a known type (B) radio, and Δ is a type a feature vector p _g according to known radios (C).
In the example of FIG. 8, since the number of wireless devices belonging to the type (A) is 10, ten feature vectors ◯ are obtained.
Further, since the number of wireless devices belonging to type (B) is 10, ten feature quantity vectors ◇ are obtained, and the number of wireless devices belonging to type (C) is 10. Ten feature quantity vectors Δ are obtained.

図８Ａは、特徴量算出用のパラメータが適正でないために、種類（Ａ）の無線機に係る特徴量ベクトル〇と、種類（Ｂ）の無線機に係る特徴量ベクトル◇と、種類（Ｃ）の無線機に係る特徴量ベクトル△とが一部重なっている例を示している。
既知の無線機に係る特徴量ベクトル〇と特徴量ベクトル◇と特徴量ベクトル△とが重なっている場合、これらの特徴量ベクトル〇，◇，△と、識別対象の無線機に係る特徴量ベクトルとを比較しても、特徴量ベクトル〇，◇，△の中で、識別対象の無線機に係る特徴量ベクトルと最も類似している特徴量ベクトルを高精度に特定することができない。
したがって、特徴量ベクトル〇，◇，△の中で、識別対象の無線機に係る特徴量ベクトルと最も類似している特徴量ベクトルを高精度に特定する上で、既知の無線機に係る特徴量ベクトル〇と特徴量ベクトル◇と特徴量ベクトル△とが重なっていないことが望ましい。
図８Ｂは、特徴量算出用のパラメータが適正であるために、種類（Ａ）の無線機に係る特徴量ベクトル〇と、種類（Ｂ）の無線機に係る特徴量ベクトル◇と、種類（Ｃ）の無線機に係る特徴量ベクトル△とが重なっていない例を示している。FIG. 8A shows that since the parameter for calculating the feature quantity is not appropriate, the feature quantity vector 0 related to the type (A) radio, the feature quantity vector ◇ related to the type (B) radio, and the type (C) In this example, the feature vector Δ related to the wireless device partially overlaps.
When the feature quantity vector 〇, feature quantity vector ◇, and feature quantity vector △ related to a known wireless device overlap, these feature quantity vectors 〇, ◇, △, and the feature quantity vector related to the wireless device to be identified are In the feature quantity vectors ◯, ◇, △, the feature quantity vector that is most similar to the feature quantity vector related to the wireless device to be identified cannot be identified with high accuracy.
Therefore, among the feature quantity vectors ◯, ◇, △, the feature quantity vector that is most similar to the feature quantity vector associated with the radio device to be identified is identified with high accuracy. It is desirable that the vector 0, the feature vector ◇, and the feature vector Δ do not overlap.
FIG. 8B shows that since the parameter for calculating the feature quantity is appropriate, the feature quantity vector 0 related to the type (A) radio, the feature quantity vector ◇ related to the type (B) radio, and the type (C ) Shows an example in which the feature quantity vector Δ related to the wireless device does not overlap.

第１のパラメータ設定部５は、特徴量算出用のパラメータの初期設定として、時間帯（ｔ_１〜ｔ_２）の長さである時間範囲Δｔを示すパラメータと、設定周波数ビンｋ_ｂを示すパラメータとを任意の値に設定する。
時間帯（ｔ_１〜ｔ_２）については、学習用のエネルギー時間波形Ｅ_ｇ（ｍ）に基づいて設定する方法が考えられる。
例えば、学習用無線信号ｓ_ｇ（ｎ）の立ち上がり時刻Ｔ_０よりも前の時刻において、学習用のエネルギー時間波形Ｅ_ｇ（ｍ）の最大値よりもエネルギーが３０ｄＢ下がっている時刻をｔ_１とする。また、学習用無線信号ｓ_ｇ（ｎ）の立ち上がり時刻Ｔ_０よりも後の時刻において、学習用のエネルギー時間波形Ｅ_ｇ（ｍ）の最大値よりもエネルギーが３０ｄＢ下がっている時刻をｔ_２とする。
また、ｔ_２については、学習用フーリエ変換信号ｓ_ｇ（ｎ）の周波数方向の変動に基づいて設定する方法が考えられる。
例えば、学習用フーリエ変換信号ｓ_ｇ（ｎ）の過渡応答時の周波数変動が、特定の周波数範囲内に収まった時刻をｔ_２とする。The first parameter setting unit 5, as the initial setting of the parameters for the feature quantity calculation, a parameter indicating a length of time ranging Δt time period (t 1 _~t _2), parameters indicating the set frequency bin k _b And are set to arbitrary values.
A method of setting the time zone (t _{1 to} t ₂ ) based on the learning energy time waveform E _g (m) is conceivable.
For example, the time when the energy is lower by 30 dB than the maximum value of the learning energy time waveform E _g (m) at the time before the rising time T _{0 of} the learning radio signal s _g (n) is denoted by t ₁ . To do. Further, the time when the energy is 30 dB lower than the maximum value of the learning energy time waveform E _g (m) at the time after the rising time T _{0 of} the learning radio signal s _g (n) is t ₂ . To do.
Also, a method of setting t ₂ based on the fluctuation in the frequency direction of the learning Fourier transform signal s _g (n) is conceivable.
For example, the time when the frequency variation at the time of the transient response of the learning Fourier transform signal s _g (n) falls within a specific frequency range is defined as t ₂ .

設定周波数ビンｋ_ｂについては、学習用のエネルギー時間波形Ｅ_ｇ（ｍ）において、平均電力が最大値をとる周波数ビンを採用する方法が考えられる。
また、設定周波数ビンｋ_ｂについては、学習用のエネルギー時間波形Ｅ_ｇ（ｍ）において、平均電力が閾値以上の周波数ビンを検出して、検出した１以上の周波数ビンに係る周波数の平均を求め、平均の周波数の周波数ビンを採用する方法が考えられる。For the set frequency bin k _b , a method of adopting a frequency bin in which the average power takes a maximum value in the energy time waveform E _g (m) for learning can be considered.
For the set frequency bin k _b , in the energy time waveform E _g (m) for learning, a frequency bin having an average power equal to or higher than a threshold is detected, and an average of the frequencies related to the detected one or more frequency bins is obtained. A method that employs frequency bins having an average frequency can be considered.

第１のパラメータ設定部５は、特徴量算出部７から既知の無線機の特徴量ベクトルｐ_ｇを第２のデータベース部１０に格納した旨を示す情報を受けると、第２のデータベース部１０から、複数の既知の無線機の特徴量ベクトルｐ_ｇを取得する。The first parameter setting unit 5 receives the information indicating that stored from the feature amount calculation unit 7 feature vector p _g of the known radio in the second database section 10, from the second database 10 , it acquires the feature quantity vector p _g plurality of known radio.

第１のパラメータ設定部５は、複数の既知の無線機に係る特徴量ベクトルｐ_ｇの間に重なりがあるか否かを判定し、例えば、図８Ｂに示すように、複数の既知の無線機に係る特徴量ベクトルｐ_ｇの間に重なりがなければ、先に設定している特徴量算出用のパラメータを有効として、パラメータの設定処理を終了する。
第１のパラメータ設定部５は、例えば、図８Ａに示すように、複数の既知の無線機に係る特徴量ベクトルｐ_ｇの間に重なりがあれば、特徴量算出用のパラメータとして、先に設定している時間範囲Δｔ又は設定周波数ビンｋ_ｂのうち、少なくとも１つを更新する。
特徴量算出部７は、第１のパラメータ設定部５により更新された特徴量算出用のパラメータを用いて、再度、複数の既知の無線機の特徴量ベクトルｐ_ｇをそれぞれ算出する。
第１のパラメータ設定部５は、特徴量算出部７により再度算出された複数の既知の無線機に係る特徴量ベクトルｐ_ｇの間に重なりがあるか否かを判定する。
複数の既知の無線機に係る特徴量ベクトルｐ_ｇの間に重なりがなくなるまで、第１のパラメータ設定部５による特徴量算出用のパラメータの更新処理と、特徴量算出部７による特徴量ベクトルｐ_ｇの算出処理とが繰り返される。
なお、パラメータの更新処理と、特徴量ベクトルｐ_ｇの算出処理とを繰り返しても、複数の既知の無線機に係る重なりがなくならない場合には、重なりが最も小さくなったパラメータを適正なパラメータとして、パラメータの設定処理を終了する。The first parameter setting unit 5 determines whether there is overlap between the feature vector p _g according to several known radio, for example, as shown in FIG. 8B, a plurality of known radio without overlap between the feature vector p _g according to, as an effective parameter for feature calculation is set to above, and ends the parameter setting processing.
The first parameter setting unit 5 sets, for example, as shown in Figure 8A, if there is overlap between the feature vector p _g according to several known radio, as a parameter for the feature amount calculation, first of to have a time range Δt or sets frequency bin k _b, updating at least one.
Feature amount calculation unit 7, using the parameters for the updated feature quantity calculation by the first parameter setting unit 5, again, it is calculated several known feature quantity vectors p _g radios respectively.
The first parameter setting unit 5 determines whether there is overlap between the feature vector p _g according to several known radio calculated again by the feature amount calculation unit 7.
Until the overlap is eliminated between the feature vector p _g according to several known radio, and updating the parameters for feature calculation by the first parameter setting unit 5, feature vector p by the feature amount calculating section 7 The calculation process of _g is repeated.
Note that, even if the parameter update process and the feature vector _pg calculation process are repeated, if the overlap related to a plurality of known wireless devices is not eliminated, the parameter with the smallest overlap is set as an appropriate parameter. Then, the parameter setting process is terminated.

第１のパラメータ設定部５により特徴量算出用のパラメータが更新されることで、特徴量算出用のパラメータが適正な値になると、図８Ｂに示すように、第２のデータベース部１０に格納される複数の既知の無線機に係る特徴量ベクトルｐ_ｇの間の重なりがなくなる、もしくは、重なりが最も小さくなる。
特徴量算出用のパラメータが更新されたのち、特徴量算出部７により再度算出された複数の既知の無線機に係る特徴量ベクトルｐ_ｇが第２のデータベース部１０に格納される際、不適正な特徴量算出用のパラメータに従って算出されている複数の既知の無線機に係る特徴量ベクトルｐ_ｇは、第２のデータベース部１０から削除される。
図１０は、第２のデータベース部１０に格納される複数の既知の無線機に係る特徴量の一例を示す説明図である。
図１０では、既知の無線機の種類が、種類（Ａ）、種類（Ｂ）の例を示しており、また、種類（Ａ）に属する無線機が、シリアルナンバー（１）の個体、シリアルナンバー（２）の個体の例を示している。
また、図１０には、特徴量算出用のパラメータ及び重み付けパラメータの一例を示している。When the parameter for feature quantity calculation is updated by the first parameter setting section 5 and the parameter for feature quantity calculation has an appropriate value, it is stored in the second database section 10 as shown in FIG. 8B. that overlap between the plurality of known feature vector p _g of the radio is eliminated, or, overlap smallest.
After the parameters for feature calculation is updated, when a feature vector p _g according to several known radio calculated again by the feature amount calculation unit 7 is stored in the second database section 10, improper feature vector p _g according to several known radio that is calculated in accordance with the parameters for the feature quantity calculation is deleted from the second database 10.
FIG. 10 is an explanatory diagram illustrating an example of feature amounts related to a plurality of known wireless devices stored in the second database unit 10.
FIG. 10 shows an example in which the types of known radios are type (A) and type (B), and the radios belonging to type (A) are serial number (1) individuals, serial numbers. The example of the individual of (2) is shown.
FIG. 10 shows an example of the feature amount calculation parameter and the weighting parameter.

この実施の形態１では、第１のパラメータ設定部５によって、特徴量算出用のパラメータが適正な値になるまで、特徴量算出用のパラメータが更新される例を示しているが、特徴量算出用のパラメータの更新タイミングは、これに限るものではない。
例えば、信号受信部１によって、既知の無線機から新たに送信された学習用無線信号が受信されることで、特徴量算出部７によって、新たに特徴量が算出されたタイミングで、第１のパラメータ設定部５が、特徴量算出用のパラメータを更新するようにしてもよい。In the first embodiment, an example is shown in which the first parameter setting unit 5 updates the feature amount calculation parameter until the feature amount calculation parameter reaches an appropriate value. However, the update timing of the parameters for use is not limited to this.
For example, when a learning radio signal newly transmitted from a known wireless device is received by the signal receiving unit 1, the first calculation is performed at the timing when the feature amount is newly calculated by the feature amount calculating unit 7. The parameter setting unit 5 may update the feature amount calculation parameter.

＜識別時の処理＞
次に、識別対象の無線機から送信された無線信号を用いて、無線機を識別する際の処理内容を説明する。
図１１は、この発明の実施の形態１による無線機識別方法における識別処理を示すフローチャートである。
信号受信部１は、識別対象の無線機から送信された無線信号を識別用無線信号として受信し、受信した識別用無線信号をアナログ信号からデジタル信号に変換する（図１１のステップＳＴ１１）。
信号受信部１は、デジタルの識別用無線信号ｓ_ｄ（ｎ）をフーリエ変換部２に出力する。ｓ_ｄ（ｎ）におけるｎは、識別用無線信号のサンプル番号である。<Process during identification>
Next, processing contents when identifying a wireless device using a wireless signal transmitted from the wireless device to be identified will be described.
FIG. 11 is a flowchart showing identification processing in the radio identification method according to Embodiment 1 of the present invention.
The signal receiving unit 1 receives the radio signal transmitted from the identification target radio as an identification radio signal, and converts the received identification radio signal from an analog signal to a digital signal (step ST11 in FIG. 11).
The signal receiving unit 1 outputs a digital identification radio signal s _d (n) to the Fourier transform unit 2. n in s _d (n) is the sample number of the radio signal for identification.

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

立ち上がり検出部３は、識別用無線信号ｓ_ｄ（ｎ）の立ち上がり時刻Ｔ_０を検出する処理を実施する（図１１のステップＳＴ１３）。
以下、立ち上がり検出部３による立ち上がり時刻Ｔ_０の検出処理を具体的に説明する。The Fourier transform unit 2 performs STFT on the digital identification radio signal s _d (n) output from the signal receiving unit 1 as shown in the following equation (23) (step ST12 in FIG. 11).
The Fourier transform unit 2 outputs the identification Fourier transform signal S _d (m, k), which is a Fourier transform signal indicating the STFT result of the identification radio signal s _d (n), to the rise detection unit 3 and the feature amount calculation unit 7. (Step ST12 in FIG. 11). FIG. 5 also shows the STFT of the identification radio signal s _d (n).

The rise detection unit 3 performs a process of detecting the rise time T ₀ of the identification radio signal s _d (n) (step ST13 in FIG. 11).
Hereinafter, the detection process of the rising time T ₀ by the rising detection unit 3 will be specifically described.

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

First, as shown in the following formula (24), the rising edge detection unit 3 calculates the square value | S _d (m, k) of the Fourier transform signal S _d (m, k) for identification output from the Fourier transform unit 2. ) | ² is calculated as spectrogram SP _d (m, k).

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

式（２５）において、ｋ_１は、スペクトログラムＳＰ_ｄ（ｍ，ｋ）に含まれている複数の周波数のうち、電力を合計する周波数の下限周波数を示す周波数ビン番号、ｋ_ｎは、スペクトログラムＳＰ_ｄ（ｍ，ｋ）に含まれている複数の周波数のうち、電力を合計する周波数の上限周波数を示す周波数ビン番号である。
ｋ＝ｋ_１，・・・，ｋ_ｎである。Next, the rising edge detection unit 3 calculates an energy time waveform E _d (m) for identification, which is a time waveform of energy in the spectrogram SP _d (m, k), as shown in the following formula (25).
The energy time waveform E _d (m) for identification is calculated by summing powers at the same time among powers at a plurality of frequencies included in the spectrogram SP _d (m, k). This is the total value of power at the time.

In the formula (25), _{k 1,} among a plurality of frequencies contained in the spectrogram _SP d (m, k), the frequency bin number that indicates the lower limit frequency of the frequency summing power, _{k n} is spectrogram SP _d It is a frequency bin number indicating the upper limit frequency of the frequency for summing power among a plurality of frequencies included in (m, k).
k ₌ k 1, ···, a _{k n.}

次に、立ち上がり検出部３は、識別用のエネルギー時間波形Ｅ_ｄ（ｍ）と閾値Ｅ_ｔｈとを比較する。
立ち上がり検出部３は、識別用のエネルギー時間波形Ｅ_ｄ（ｍ）と閾値Ｅ_ｔｈとの比較結果を参照することで、最初は、ノイズレベルのエネルギーを示していた識別用のエネルギー時間波形Ｅ_ｄ（ｍ）が、閾値Ｅ_ｔｈに到達したか否かを判定する。
立ち上がり検出部３は、識別用のエネルギー時間波形Ｅ_ｄ（ｍ）が閾値Ｅ_ｔｈに到達した時刻を、識別用無線信号ｓ_ｄ（ｎ）の立ち上がり時刻Ｔ_０として検出する。
この実施の形態１では、立ち上がり検出部３が、識別用無線信号ｓ_ｄ（ｎ）の立ち上がり時刻Ｔ_０を検出する際、スペクトログラムＳＰ_ｄ（ｍ，ｋ）に含まれている複数の周波数における電力のうち、同一時刻の電力同士をそれぞれ合計するようにしている。これにより、立ち上がり時に電力が大きく変動する周波数が含まれている場合でも、当該電力の変動の影響が低減されるため、立ち上がり時刻Ｔ_０の検出精度が向上する。Next, the rise detection unit 3 compares the energy time waveform E _d (m) for identification with the threshold value E _th .
The rising edge detection unit 3 refers to the comparison result between the identification energy time waveform E _d (m) and the threshold value E _th , so that the identification energy time waveform E _d that initially indicates the noise level energy is used. It is determined whether (m) has reached the threshold value _Eth .
The rise detection unit 3 detects the time when the identification energy time waveform E _d (m) reaches the threshold E _th as the rise time T ₀ of the identification radio signal s _d (n).
In the first embodiment, when the rise detection unit 3 detects the rise time T ₀ of the identification radio signal s _d (n), the power at a plurality of frequencies included in the spectrogram SP _d (m, k). Among them, the powers at the same time are summed up. Thus, even if it contains frequencies vary greatly power at the rise time, the influence of the variation of the power is reduced, thereby improving the detection accuracy of the rise time T _0.

特徴量算出部７は、第１のパラメータ設定部５により設定された特徴量算出用のパラメータから、立ち上がり時刻Ｔ_０を含む時間帯（ｔ_１〜ｔ_２）及び設定周波数ビンｋ_ｂを認識する。
この特徴量算出用のパラメータは、第１のパラメータ設定部５によって、複数の既知の無線機に係る特徴量ベクトルｐ_ｇの間に重なりがなくなるように、もしくは、重なりが最も小さくなるように設定されたパラメータである。
次に、特徴量算出部７は、図６に示すように、フーリエ変換部２から出力された識別用フーリエ変換信号Ｓ_ｄ（ｍ，ｋ）の中から、立ち上がり検出部３により検出された識別用無線信号ｓ_ｄ（ｎ）の立ち上がり時刻Ｔ_０を含む時間帯（ｔ_１〜ｔ_２）の識別用フーリエ変換信号ｓ_ｄ（ｎ）_{ｔ１〜ｔ２}を抽出する。Feature amount calculation unit 7 recognizes from the parameter for the set characteristic quantity calculation, a time zone including the rise time T ₀ the (t 1 _~t ₂₎ and the set frequency bin k _b of the first parameter setting unit 5 .
Parameters for the feature amount calculated by the first parameter setting unit 5, as the overlap between the feature vector p _g according to several known radio is eliminated, or set so that the overlap is minimized Parameter.
Next, as shown in FIG. 6, the feature quantity calculation unit 7 performs the identification detected by the rising detection unit 3 from the identification Fourier transform signal S _d (m, k) output from the Fourier transform unit 2. The Fourier transform signals for identification s _d (n) _t1 to _t2 in the time zone (t _{1 to} t ₂ ) including the rising time T ₀ of the radio signal s _d (n) for use are extracted.

次に、特徴量算出部７は、図７に示すように、抽出した識別用フーリエ変換信号ｓ_ｄ（ｎ）_{ｔ１〜ｔ２}に含まれている複素時間信号の中から、設定周波数ビンｋ_ｂの複素時間信号ｓ_ｄハット（ｎ）_{ｔ１〜ｔ２}の切り出しを行う。
電子出願の都合上、明細書の文章中では、ｓの文字の上に“＾”の記号を付することができないので、ｓ_ｄハット（ｎ）_{ｔ１〜ｔ２}のように表記している。
特徴量算出部７は、切り出した複素時間信号ｓ_ｄハット（ｎ）_{ｔ１〜ｔ２}から、抽出した識別用フーリエ変換信号ｓ_ｄ（ｎ）_{ｔ１〜ｔ２}の変動を特徴量として算出する（図１１のステップＳＴ１４）。Next, the feature amount calculation unit 7, as shown in FIG. 7, from the complex time signal included in the extracted identification for Fourier transform signal _s _{d (n)} t1~t2, the set frequency bin k _b The complex time signal s _d hat (n) _{t1 to t2} is cut out.
For the convenience of electronic filing, in the text of the specification, since the symbol “^” cannot be added on the letter s, it is represented as s _d hat (n) _{t1 to t2} .
The feature amount calculation unit 7 calculates, from the extracted complex time signal s _d hat (n) _t1 to _t2 , the fluctuation of the extracted Fourier transform signal s _d (n) _{t1 to t2} as a feature amount (see FIG. 11). Step ST14).

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

式（２６）において、Ｉ_ｄ（ｎ）は、複素時間信号ｓ_ｄハット（ｎ）_{ｔ１〜ｔ２}の実部、Ｑ_ｄ（ｎ）は、複素時間信号ｓ_ｄハット（ｎ）_{ｔ１〜ｔ２}の虚部である。Hereinafter, the feature amount calculation processing by the feature amount calculation unit 7 will be described in detail.
Feature amount calculation unit 7, the complex time signal _{s d} hat _{(n) t1 to t2} cut, as shown in the following equation (26), the complex time signal _{s d} hat _{(n) t1 to t2} of the instantaneous amplitude a _d (n) is calculated.

In Equation (26), I _d (n) is the real part of the complex time signal s _d hat (n) _{t1 to t2} , and Q _d (n) is the imaginary of the complex time signal s _d hat (n) _{t1 to t2} . Part.

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

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

The feature quantity calculation unit 7, the complex time signal _{s d} hat _{(n) t1 to t2} cut, as shown in the following equation (27), the instantaneous complex time signal _{s d} hat _{(n) t1 to t2} The phase φ _d (n) is calculated.

The feature quantity calculation unit 7, the cut complex time signal _{s d} hat _{(n) t1 to t2,} as shown in the following equation (28), the instantaneous complex time signal _{s d} hat _{(n) t1 to t2} The frequency f _d (n) is calculated.

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

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

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

Next, as shown in the following formula (29), the feature amount calculation unit 7 uses the variance of the instantaneous amplitude a _d (n) in the complex time signal s _d hat (n) _{t1 to t2} as the first feature amount. The value p _d1 is calculated.

As shown in the following equation (31), the feature amount calculation unit 7 uses the skewness p _d2 of the instantaneous amplitude a _d (n) in the complex time signal s _d hat (n) _{t1 to t2} as the second feature amount. Is calculated.

As shown in the following formula (33), the feature amount calculation unit 7 uses the kurtosis p _d3 of the instantaneous amplitude a _d (n) in the complex time signal s _d hat (n) _{t1 to t2} as the third feature amount. Is calculated.

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

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

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

Next, as shown in the following formula (34), the feature amount calculation unit 7 uses the dispersion of the instantaneous phase φ _d (n) in the complex time signal s _d hat (n) _{t1 to t2} as the fourth feature amount. The value _pd4 is calculated.

As shown in the following formula (36), the feature quantity calculation unit 7 uses, as the fifth feature quantity, the skewness p _d5 of the instantaneous phase φ _d (n) in the complex time signal s _d hat (n) _{t1 to t2} . Is calculated.

As shown in the following formula (38), the feature quantity calculation unit 7 uses the kurtosis p _d6 of the instantaneous phase φ _d (n) in the complex time signal s _d hat (n) _{t1 to t2} as the sixth feature quantity. Is calculated.

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

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

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

Next, as shown in the following formula (39), the feature amount calculation unit 7 uses the variance of the instantaneous frequency f _d (n) in the complex time signal s _d hat (n) _{t1 to t2} as the seventh feature amount. The value _pd7 is calculated.

As shown in the following equation (41), the feature amount calculation unit 7 uses the skewness p _d8 of the instantaneous frequency f _d (n) in the complex time signal s _d hat (n) _{t1 to t2} as the eighth feature amount. Is calculated.

As shown in the following equation (43), the feature amount calculation unit 7 uses, as the ninth feature amount, the kurtosis p _d9 of the instantaneous frequency f _d (n) in the complex time signal s _d hat (n) _{t1 to t2} . Is calculated.

特徴量算出部７は、第２のパラメータ設定部６により設定された第１から第９の特徴量ｐ_ｄｉ（ｉ＝１，２，・・・，９）の重みｗ_ｉを示す重み付けパラメータを取得する。
特徴量算出部７は、以下の式（４４）に示すように、第１から第９の特徴量ｐ_ｄｉと重みｗ_ｉとをそれぞれ乗算することで、第１から第９の特徴量ｐ_ｄｉの重み付けを行う。
ｐ’_ｄｉ＝ｐ_ｄｉ×ｗ_ｉ（４４）
ｉ＝１，２，・・・，９
特徴量算出部７は、重み付け後の第１から第９の特徴量ｐ’_ｄｉを並べたベクトルを、識別対象の無線機の特徴量ベクトルｐ_ｄとして算出する。
特徴量算出部７は、算出した識別対象の無線機の特徴量ベクトルｐ_ｄを無線機識別部１１に出力する。The feature amount calculation unit 7 sets weighting parameters indicating the weights w _i of the first to ninth feature amounts p _di (i = 1, 2,..., 9) set by the second parameter setting unit 6. get.
As shown in the following formula (44), the feature amount calculation unit 7 multiplies the first to ninth feature amounts p _di and the weights w _i , respectively, so that the first to ninth feature amounts p _{di are obtained.} Weighting is performed.
p ′ _di = p _di × w _i (44)
i = 1, 2,..., 9
Feature amount calculation unit 7, a vector from the first the weighted arranged feature quantities p _'di ninth, calculated as a feature vector p _d of the identification target radio.
The feature amount calculation unit 7 outputs the calculated feature amount vector _pd of the wireless device to be identified to the wireless device identification unit 11.

無線機識別部１１は、特徴量算出部７により算出された識別対象の無線機に係る特徴量ベクトルｐ_ｄと、第２のデータベース部１０により格納されている複数の既知の無線機に係る特徴量ベクトルｐ_ｇとをそれぞれ比較する。
そして、無線機識別部１１は、識別対象の無線機に係る特徴量ベクトルｐ_ｄと、複数の既知の無線機に係る特徴量ベクトルｐ_ｇとの比較結果に基づいて、識別対象の無線機を識別する（図１１のステップＳＴ１５）。Radio identification unit 11, a feature quantity vector p _d of the characteristic amount calculation unit 7 identifies target wireless device calculated by the feature according to the plurality of known radio stored by the second database section 10 comparing the amount vector p _g respectively.
The radio identification unit 11 includes a feature vector p _d according to the identification target radio, based on a result of comparison between the feature vector p _g according to several known radio, to be identified radios Identify (step ST15 in FIG. 11).

以下、無線機識別部１１による無線機の識別処理を具体的に説明する。
図１２は、無線機識別部１１による無線機の識別処理を示す説明図である。
図１２では、説明の簡単化のため、特徴量算出部７により算出される第１から第９の特徴量のうち、第１から第３の特徴量だけを例示している。Hereinafter, the wireless device identification processing by the wireless device identification unit 11 will be described in detail.
FIG. 12 is an explanatory diagram showing a radio identification process performed by the radio identification unit 11.
In FIG. 12, only the first to third feature amounts of the first to ninth feature amounts calculated by the feature amount calculation unit 7 are illustrated for simplicity of explanation.

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

式（４５）において、ｘ_ｉ，ｈは、ｎ次元の特徴量空間における既知の無線機ｈ（ｈ＝１，・・・，Ｈ）に係る特徴量ベクトルｐ_ｇの位置、ｘ_ｊは、ｎ次元の特徴量空間における識別対象の無線機に係る特徴量ベクトルｐ_ｄの位置である。
（ｘ_ｉ，ｈ−ｘ_ｊ）は、ｎ次元の特徴量空間における２点間のユークリッド距離、Ｓは、既知の無線機ｈから送信される学習用無線信号ｓ_ｇ（ｎ）の特徴量空間における分布の分散共分散行列である。
マハラノビス距離ｄ_ｉｊ，ｈは、学習用無線信号ｓ_ｇ（ｎ）の特徴量空間における分布の分散による影響が正規化された空間での距離となる。Radio identification unit 11, as shown in the following equation (45), a feature amount vector p _d of the identification target radio, a plurality of known radio stored by the second database section 10 Mahalanobis distance _{d ij} between the feature vector _{p g} _{according, h} a is calculated.
Here, for convenience of explanation, the known radio is located H pieces, it is assumed that the feature vector p _g of the H-number of radio is respectively stored in the second database section 10.
Among the H radio devices, for example, a radio device belonging to type (A), a radio device belonging to type (B), or a radio device belonging to type (C) is included.

In the formula _{(45), x i, h} is the known radio h in the feature quantity space of n dimensions (h = 1, ···, H ) position of the feature vector _{p g} according to, _{x j} is n is the position of the feature vector p _d of the identification target radio in feature space dimensions.
(X _{i, h} −x _j ) is the Euclidean distance between two points in the n-dimensional feature amount space, and S is the feature amount space of the learning radio signal s _g (n) transmitted from the known radio h. Is the variance-covariance matrix of the distribution at.
The Mahalanobis distance d _{ij, h} is a distance in a space where the influence of the distribution of the distribution in the feature amount space of the learning radio signal s _g (n) is normalized.

無線機識別部１１は、算出したＨ個のマハラノビス距離ｄ_ｉｊ，ｈを相互に比較し、Ｈ個のマハラノビス距離ｄ_ｉｊ，ｈの中で、最小のマハラノビス距離ｄ_ｉｊ，ｈを特定する。
無線機識別部１１は、Ｈ個の無線機の中で、最小のマハラノビス距離ｄ_ｉｊ，ｈに係る既知の無線機が、識別対象の無線機と最も類似していると判断する。
無線機識別部１１は、識別対象の無線機の識別結果として、最小のマハラノビス距離ｄ_ｉｊ，ｈに係る既知の無線機のシリアルナンバーを出力する。これにより、無線機の種類だけでなく、無線機の個体識別が可能になる。
図１２の例では、３つの既知の無線機に係る特徴量ベクトルｐ_ｇの位置が示されている。３つの既知の無線機の中で、特徴量ベクトルｐ_ｇの位置が●である既知の無線機とのマハラノビス距離が最小であるため、識別対象の無線機は、特徴量ベクトルｐ_ｇの位置が●の無線機であると識別される。The radio equipment identification unit 11 compares the calculated H Mahalanobis distances d _{ij, h} with each other, and specifies the minimum Mahalanobis distance d _{ij, h} among the H Mahalanobis distances _{dij, h} .
The wireless device identification unit 11 determines that the known wireless device having the smallest Mahalanobis distance _{dij, h} among the H wireless devices is most similar to the wireless device to be identified.
The wireless device identification unit 11 outputs the serial number of a known wireless device related to the minimum Mahalanobis distance _{dij, h} as the identification result of the wireless device to be identified. As a result, not only the type of the wireless device but also the individual identification of the wireless device becomes possible.
In the example of FIG. 12, the position of the feature vector p _g of the three known radio is shown. Among the three known radio, since the Mahalanobis distance between the feature vector p position of _g is known radio is ● is the minimum, radio identification target, the position of the feature vector p _g ● It is identified as a radio.

ここでは、無線機識別部１１が、マハラノビス距離ｄ_ｉｊ，ｈを用いて、無線機を識別する例を示しているが、これに限るものではなく、例えば、最小二乗法による線形識別、サポートベクターマシン（ＳＶＭ：ＳｕｐｐｏｒｔＶｅｃｔｏｒＭａｃｈｉｎｅ）、あるいは、ニューラルネットワークを用いて、無線機を識別するようにしてもよい。Here, an example is shown in which the wireless device identification unit 11 identifies a wireless device using the Mahalanobis distance _{dij, h} . However, the present invention is not limited to this. For example, linear identification by the least square method, support vector A radio may be identified using a machine (SVM: Support Vector Machine) or a neural network.

以上で明らかなように、この実施の形態１によれば、フーリエ変換部２から出力されたフーリエ変換信号のスペクトログラムにおけるエネルギーの時間波形を算出し、エネルギーの時間波形から無線信号の立ち上がり時刻を検出する立ち上がり検出部３と、フーリエ変換部２から出力されたフーリエ変換信号の中から、立ち上がり検出部３により検出された立ち上がり時刻を含む時間帯のフーリエ変換信号を抽出し、抽出したフーリエ変換信号に含まれている設定周波数成分の複素時間信号から、抽出したフーリエ変換信号の変動を特徴量として算出する特徴量算出部７とを設け、無線機識別部１１が、特徴量算出部７により算出された特徴量に基づいて、無線機を識別するように構成したので、無線機の個体を識別することができる効果を奏する。 As apparent from the above, according to the first embodiment, the time waveform of energy in the spectrogram of the Fourier transform signal output from the Fourier transform unit 2 is calculated, and the rising time of the radio signal is detected from the time waveform of energy. The rise time detection unit 3 and the Fourier transform signal output from the Fourier transform unit 2 extract a Fourier transform signal in a time zone including the rise time detected by the rise time detection unit 3, and the extracted Fourier transform signal A feature amount calculation unit 7 is provided that calculates, as a feature amount, fluctuations in the extracted Fourier transform signal from the complex time signal of the set frequency component included, and the wireless device identification unit 11 is calculated by the feature amount calculation unit 7. Since the wireless device is identified based on the feature amount, it is possible to identify the individual wireless device. To.

なお、本願発明はその発明の範囲内において、実施の形態の任意の構成要素の変形、もしくは実施の形態の任意の構成要素の省略が可能である。 In the present invention, any constituent element of the embodiment can be modified or any constituent element of the embodiment can be omitted within the scope of the invention.

この発明は、無線機を識別する無線機識別装置及び無線機識別方法に適している。 The present invention is suitable for a wireless device identification device and a wireless device identification method for identifying a wireless device.

１信号受信部、２フーリエ変換部、３立ち上がり検出部、４パラメータ設定部、５第１のパラメータ設定部、６第２のパラメータ設定部、７特徴量算出部、８データベース部、９第１のデータベース部、１０第２のデータベース部、１１無線機識別部、２１信号受信回路、２２フーリエ変換回路、２３立ち上がり検出回路、２４パラメータ設定回路、２５特徴量算出回路、２６データベース回路、２７無線機識別回路、３１メモリ、３２プロセッサ。 DESCRIPTION OF SYMBOLS 1 Signal receiving part, 2 Fourier transform part, 3 Rise detection part, 4 Parameter setting part, 5 1st parameter setting part, 6 2nd parameter setting part, 7 Feature-value calculation part, 8 Database part, 9 1st Database unit, 10 Second database unit, 11 Radio identification unit, 21 Signal reception circuit, 22 Fourier transform circuit, 23 Rise detection circuit, 24 Parameter setting circuit, 25 Feature quantity calculation circuit, 26 Database circuit, 27 Radio identification Circuit, 31 memory, 32 processor.

Claims

A signal receiving unit that receives a radio signal transmitted from a radio to be identified; and
A Fourier transform unit that performs a short-time Fourier transform on the radio signal received by the signal receiving unit and outputs a Fourier transform signal indicating a short-time Fourier transform result of the radio signal;
Calculating a time waveform of energy in a spectrogram of a Fourier transform signal output from the Fourier transform unit, and detecting a rise time of the radio signal from the time waveform of the energy;
From the Fourier transform signal output from the Fourier transform unit, the Fourier transform signal in the time zone including the rise time detected by the rise detection unit is extracted, and the set frequency component included in the extracted Fourier transform signal A feature amount calculation unit that calculates, as a feature amount, a variation of the Fourier transform signal extracted from the complex time signal of
A radio identification apparatus comprising: a radio identification unit that identifies the radio based on the feature quantity calculated by the feature quantity calculation unit.

The signal receiving unit receives a radio signal transmitted from a known radio,
The Fourier transform unit performs a short-time Fourier transform on the radio signal of the known radio received by the signal receiver, and outputs a Fourier transform signal indicating a short-time Fourier transform result of the radio signal of the known radio.
The rise detection unit calculates a time waveform of energy in a spectrogram of a Fourier transform signal related to a known radio output from the Fourier transform unit, and a rise time of the radio signal of the known radio from the time waveform of the energy Detect
The feature amount calculation unit extracts a Fourier transform signal in a time zone including a rise time detected by the rise detection unit from Fourier transform signals related to a known radio output from the Fourier transform unit, From the complex time signal of the set frequency component included in the Fourier transform signal related to the extracted known radio device, the variation of the Fourier transform signal related to the extracted known radio device is calculated as a feature amount, and the calculated known radio The feature quantity related to the machine is stored in the database part,
The wireless device identification unit compares the feature amount related to the identification target wireless device calculated by the feature amount calculation unit with the feature amount related to the known wireless device stored in the database unit, and the feature amount 2. The radio identification device according to claim 1, wherein a radio to be identified is identified based on the comparison result.

A first parameter setting unit that sets a parameter indicating the length of the time zone and a parameter indicating the set frequency component as a parameter for calculating the feature amount;
The feature amount calculation unit calculates each of the feature amount related to the identification target wireless device and the feature amount related to the known wireless device using the feature amount calculation parameters set by the first parameter setting unit. The wireless device identification device according to claim 2, wherein:

The first parameter setting unit does not overlap the feature amounts calculated by the feature amount calculation unit when radio signals transmitted from different known radios are received by the signal reception unit. 4. The wireless device identification apparatus according to claim 3, wherein the parameter for calculating the feature amount is updated so that the overlap of the feature amounts respectively calculated by the feature amount calculation unit is minimized. .

The rise detection unit calculates a total value of power at each time as a time waveform of the energy by summing powers at the same time among powers at a plurality of frequencies included in the spectrogram. The wireless device identification device according to claim 1, wherein a time when the total value of the power reaches a threshold is detected as a rising time of the wireless signal.

The feature amount calculation unit includes:
As the first to third feature quantities, the variance value, skewness and kurtosis of the amplitude in the complex time signal are calculated, respectively, and as the fourth to sixth feature quantities, the phase variance value in the complex time signal, 2. The radio according to claim 1, wherein skewness and kurtosis are calculated, respectively, and a frequency dispersion value, skewness and kurtosis in the complex time signal are respectively calculated as seventh to ninth feature amounts. Machine identification device.

A second parameter setting unit configured to set a weighting parameter indicating a weight of the first to ninth feature amounts;
7. The radio identification apparatus according to claim 6, wherein the feature amount calculation unit performs weighting of the first to ninth feature amounts according to a weighting parameter set by the second parameter setting unit.

The signal receiving unit receives a wireless signal transmitted from the wireless device to be identified,
The Fourier transform unit performs a short-time Fourier transform on the radio signal received by the signal reception unit, and outputs a Fourier transform signal indicating a short-time Fourier transform result of the radio signal,
The rise detection unit calculates the time waveform of energy in the spectrogram of the Fourier transform signal output from the Fourier transform unit, detects the rise time of the radio signal from the time waveform of the energy,
The feature amount calculation unit extracts a Fourier transform signal in a time zone including the rise time detected by the rise detection unit from the Fourier transform signal output from the Fourier transform unit, and is included in the extracted Fourier transform signal From the complex time signal of the set frequency component, the fluctuation of the extracted Fourier transform signal is calculated as a feature value.
A wireless device identification method in which a wireless device identification unit identifies the wireless device based on the feature amount calculated by the feature amount calculation unit.