CN109472239B - Individual identification method of frequency hopping radio station - Google Patents

Individual identification method of frequency hopping radio station Download PDF

Info

Publication number
CN109472239B
CN109472239B CN201811328578.4A CN201811328578A CN109472239B CN 109472239 B CN109472239 B CN 109472239B CN 201811328578 A CN201811328578 A CN 201811328578A CN 109472239 B CN109472239 B CN 109472239B
Authority
CN
China
Prior art keywords
frequency
time
equal
energy spectrum
frequency hopping
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201811328578.4A
Other languages
Chinese (zh)
Other versions
CN109472239A (en
Inventor
李红光
齐子森
郭英
眭萍
许华
王少波
杨鑫
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Air Force Engineering University of PLA
Original Assignee
Air Force Engineering University of PLA
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Air Force Engineering University of PLA filed Critical Air Force Engineering University of PLA
Priority to CN201811328578.4A priority Critical patent/CN109472239B/en
Publication of CN109472239A publication Critical patent/CN109472239A/en
Application granted granted Critical
Publication of CN109472239B publication Critical patent/CN109472239B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/08—Feature extraction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
    • G06F18/2136—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on sparsity criteria, e.g. with an overcomplete basis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/214—Generating training patterns; Bootstrap methods, e.g. bagging or boosting
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G06F18/241—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches
    • G06F18/2411—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches based on the proximity to a decision surface, e.g. support vector machines
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/12—Classification; Matching

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Data Mining & Analysis (AREA)
  • Artificial Intelligence (AREA)
  • Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Signal Processing (AREA)
  • Mobile Radio Communication Systems (AREA)

Abstract

The method comprises the steps of firstly sparsely reconstructing a frequency hopping signal to obtain a time frequency spectrum, then partitioning the time frequency energy spectrum under different scales, and respectively extracting three characteristics of a Rayleigh entropy, a difference box dimension and a multi-fractal dimension of the time frequency energy spectrum. The classification and identification experiments show that the identification performance of the method is less influenced by the number of training samples, and the method has higher identification accuracy under the condition of less training samples.

Description

Individual identification method of frequency hopping radio station
Technical Field
The invention relates to wireless communication and signal processing technology, in particular to a frequency hopping radio station individual identification method.
Background
The sorting of the frequency hopping communication network station is the first premise of intercepting enemy communication and generating the best interference signal. The existing frequency hopping signal network station sorting mainly utilizes the duration, the azimuth information, the power, the signal time correlation and the like of a frequency hopping signal to realize the network station sorting identification of the frequency hopping signal. However, with the increase of the frequency hopping patterns, it is difficult to accurately sort the frequency hopping signals only by the above features. Due to the random discreteness of the component performance, the production process, the debugging and the like of each frequency hopping radio station, the frequency hopping signals radiated by the frequency hopping radio stations have individual characteristics which are different from those of other frequency hopping radio station signals, and therefore, the network station sorting and identification of the frequency hopping signals can be realized by utilizing the unique fine characteristics of different frequency hopping radio station signals, namely the fingerprint characteristics. At present, the commonly used fine feature extraction methods mainly include: the classification and identification are realized by extracting the characteristics of the frequency conversion transient signal of the frequency hopping communication equipment; sorting frequency hopping radio stations by using different characteristics of the transient response process of the transmitting power amplifier; the individual identification of the radiation source is realized by extracting the multi-dimensional characteristics of the time-frequency energy spectrum, but the identification accuracy is low due to more information redundancy of the characteristic set extracted by the method; the first moment, the second moment and the energy spectrum entropy characteristics of the time-frequency energy spectrum with fixed scale are extracted, but the time-frequency energy spectrum has different distribution information characteristics under different blocking conditions, so that the classification and identification accuracy is low. In addition, most of the existing algorithms for researching the fine feature recognition of the frequency hopping radio station are carried out under the conditions of large samples and high signal-to-noise ratio, and the adverse effects of the number of samples and the signal-to-noise ratio on the classification recognition rate cannot be effectively overcome.
Disclosure of Invention
Aiming at the problems in the prior art, the invention provides a frequency hopping radio station individual identification method, which comprises the following steps:
the first step is as follows: extracting frequency hopping signal time frequency energy spectrum
Assuming a finite set of frequencies for the frequency hopping signal as w, the received signal frequency is set as
Figure GSB0000194869010000021
M is more than or equal to 0 and less than or equal to P-1, a frequency hopping section is divided into P frequency points at equal intervals according to the requirement on time precision, a received frequency hopping signal y is divided into G sections at equal intervals, and each section of data ykHas a length of P, then
yk＝y(k·L：k·L+P-1) (1)
Wherein k is more than or equal to 1 and less than or equal to G, and L is a segment interval;
observation matrix Y ═ Y1，y2，...，yG]Can be expressed as
Y＝WX+V (2)
Wherein W is [ omega ]0，...，ωP-1]，ωi＝[ejω1，...，ejωP]X is a time-frequency matrix, V is a noise matrix, and X, V are belonged to CP×GC represents a complex matrix;
due to the sparsity of the frequency hopping signals in the time-frequency domain, time-frequency points in the X are sparse, non-zero points are all on the corresponding rows of each frequency hopping point, and the X is also sparse; therefore, an unconstrained optimization function with a penalty function is assumed:
Figure GSB0000194869010000022
wherein x iskWhich represents the k-th row of X,
Figure GSB0000194869010000023
denotes the value of X, μ, which minimizes L (X)1Is an X-point sparsity penalty factor, μ2Is an X row sparse penalty factor;
the AL0 algorithm introduces a gaussian function to approximate l0Norm of
Figure GSB0000194869010000024
Where s represents an arbitrary D-dimensional vector and σ represents a normal number close to zero; when σ is approximately 0, there are
Figure GSB0000194869010000025
Wherein s isiRepresents the ith row of vector s, D represents the dimension of vector s; approximating formula (4) to l0The norm represents the sparsity of the time-frequency domain of the frequency hopping signal, and then the sparse reconstruction problem can be converted into the minimization problem for solving the Gaussian sum function
Figure GSB0000194869010000031
Wherein Fσ(x) Indicating a frequency hopping signal0Norm, introducing a penalty factor lambda, and dividing l0The norm minimization problem translates into an optimization problem, i.e.
Figure GSB0000194869010000032
Wherein
Figure GSB0000194869010000033
Denotes the value of X which minimizes L (X);
the steepest descent direction is the direction of the conjugate gradient of the unconstrained optimization function L (x)
Figure GSB0000194869010000034
Figure GSB0000194869010000035
Wherein x*Is the conjugate of x; the definition of the conjugate gradient of the complex variable can be obtained
Figure GSB0000194869010000036
Wherein x isRRepresenting the real part, x, of a signal xIRepresenting the imaginary part of the signal x, i being the unit of an imaginary number, diagonal matrix Ai，i＝exp(-|xk|2/2σ2)/σ2；
According to
Figure GSB0000194869010000037
One can obtain, where the upper superscript H represents the conjugate matrix;
Figure GSB0000194869010000038
by substituting the formulae (7) and (8) into the formula (6)
Figure GSB0000194869010000039
Then the direction of the conjugate gradient of L (x) in equation (3) is calculated as
Figure GSB0000194869010000041
Wherein A is1＝fσ(X)/σ2；A2Is a diagonal matrix whose k-th diagonal element is
Figure GSB0000194869010000042
K is not less than 1 and not more than G, and X (k, k) represents the kth line of X;
the second step is that: fractal feature extraction
(1) Differential box dimension feature extraction
And (3) regarding the time-frequency axis of the frequency-hopping signal time-frequency energy spectrum subjected to sparse reconstruction as a plane, regarding the energy spectrum value as an image gray value, and calculating the difference box dimension of the time-frequency energy spectrum by the following steps:
step 1: assuming that the size of the time-frequency energy spectrum F is NxN, S is a grid of MxM, wherein M is more than or equal to 1 and less than or equal to N/2, and M belongs to Z+，Z+Representing a set of positive integers; f is divided into a plurality of sub-blocks by S, each grid S is provided with a box column with the size of M multiplied by M ', M' is more than or equal to 1 and less than or equal to N/2,
Figure GSB0000194869010000043
e represents the maximum energy spectrum value and,
Figure GSB0000194869010000044
represents a round-down operator;
step 2: setting a (mu, v) th grid, wherein mu is more than or equal to 1 and less than or equal to M, and v is more than or equal to 1 and less than or equal to M; the maximum and minimum energy values of the energy spectrum are denoted g, respectivelymaxAnd gminAt the scale M, the box count of the (μ, ν) th grid is:
Figure GSB0000194869010000045
step 3: by analogy, the number of all grid boxes under the M scale is calculated
Figure GSB0000194869010000046
Step 4: changing the scale M of the grid S, repeating the Step1-3, and calculating N under the condition of different block scales MM；
Step 5: estimating differential box dimensions
Figure GSB0000194869010000047
Wherein the content of the first and second substances,
Figure GSB0000194869010000051
(2) multi-fractal feature extraction
In order to more effectively represent the local characteristics of the time-frequency energy spectrum, extracting a multi-fractal dimension of the time-frequency energy spectrum as a second-dimensional characteristic;
the fractal MD can be expressed as:
Figure GSB0000194869010000052
in the formula, q represents a multi-fractal parameter, and in practical application, different choices of q values influence the stability of MD; according to difference boxStep1-3 of the dimension extraction algorithm, calculating N under the condition of different block scales MMThen, combining with the formula (14), obtaining the multi-fractal characteristic of the frequency-hopping signal time-frequency energy spectrum;
(3) time-frequency energy spectrum Rayleigh entropy feature extraction
The time-frequency energy distribution rules of the two frequency hopping radio station signals have certain difference, and the difference of the distribution can be measured by counting the Rayleigh entropy of the time-frequency energy spectrum, namely the difference is converted from the intuitive difference to the difference in numerical value;
let the probability density distribution of the random variable J be P ═ { P ═ P1，p2，...，pTWhere T represents the dimension of the random variable J, satisfying the condition
Figure GSB0000194869010000053
T is more than or equal to 1 and less than or equal to T, the Rayleigh entropy of J is defined as follows:
Figure GSB0000194869010000054
wherein alpha is the Rayleigh entropy order, and the Rayleigh entropy is defined by applying one-dimensional probability density distribution in the formula (15);
the continuous form of the rayleigh entropy of order α is as follows:
Figure GSB0000194869010000061
wherein f (χ, γ) is a continuous two-dimensional probability density distribution, fα(χ, γ) represents f (χ, γ) to the power α;
after the frequency hopping signal y is sparsely reconstructed, the obtained time-frequency energy spectrum P (t, f) has time-frequency edge characteristics and energy retention characteristics, as shown in formulas (17) and (18):
∫P(t，f)df＝|y(t)|2，∫P(t，f)dt＝|y(f)|2 (17)
Figure GSB0000194869010000062
wherein y (f) is the Fourier transform of y (t), and f represents the signal frequency; therefore, the two-dimensional probability density distribution f (χ, γ) of the time-frequency energy spectrum P (t, f) and the equation (16) has similar properties, so that the rayleigh entropy of the time-frequency energy spectrum of the frequency hopping signal can be defined by P (t, f) as shown in the formula:
Figure GSB0000194869010000063
in the formula Pα(t, f) represents P (t, f) to the power of alpha;
the stabilizing condition of formula (19) is ^ jeppeα(t, f) dtdf > 0; for convenience of calculation, the discrete expression of the rayleigh entropy of the time-frequency energy spectrum is as follows:
Figure GSB0000194869010000064
wherein X is a time-frequency matrix, k and k 'represent the number of sampling points of observation data, k is more than or equal to 1 and less than or equal to G, and k' is more than or equal to 1 and less than or equal to G;
psi and psi' denote frequency set dictionary, omega0≤ψ′≤ωP-1，ω0≤ψ≤ωP-1；
And finally, forming a characteristic vector V [ FD, MD (q), R from the three characteristics of the Rayleigh entropy, the difference box dimension and the multi-fractal dimension of the calculated time-frequency energy spectrumα]And carrying out classification and identification by using a support vector machine classifier.
In the individual identification method of the frequency hopping radio station, the AL0 algorithm comprises the following steps:
inputting: matrix W, measured value vector y;
step (1): let initial value x(0)＝WT(WWT)-1y；
Step (2): selection of descending sequence [ sigma ]1σ2...σJ]J ═ 6; let the convergence criterion be ε, ε representing a mean of 0 and a variance of σ2White gaussian noise of (1);
step (3): iterating an algorithm;
(1) j takes the values of 1, 2,.. J in sequence;
(2) for each value of j, let σ be σ ═ σj，
(3) When in use
Figure GSB0000194869010000071
Where norm represents the solution
Figure GSB0000194869010000072
Figure GSB0000194869010000072
2 norm of (d);
(4) step size u is determined such that
Figure GSB0000194869010000073
(5) Is updated in the gradient direction such that
Figure GSB0000194869010000074
(6) At this time, let x(j)Continuing to take the next value of j as x, and carrying out iterative operation until all values are taken in sequence;
step 4: outputting the result
Figure GSB0000194869010000075
In a specific embodiment of the present invention, σ takes on values of [2, 1, 0.5, 0.2, 0.1, 0.05] in sequence.
The method provided by the invention has the advantages that the recognition accuracy is not influenced by the change of the number of training samples and is always kept at a higher level. The main reason is that the algorithm respectively extracts the difference box dimension, the multi-fractal dimension and the Rayleigh entropy of the time-frequency energy spectrum as the feature vectors, quantitatively describes the energy change of the time-frequency energy spectrum, the complexity and the regularity of the time-frequency energy spectrum, and avoids the problem of misjudgment caused by the similarity of single features. The classification features of the method are only slightly overlapped, and the method has obvious clustering effect.
Drawings
FIG. 1 shows a flow chart of the method of the present invention;
FIG. 2 shows a time-frequency energy diagram for two frequency hopping stations; wherein 2(a) shows a first station time-frequency energy graph, and 2(b) shows a second station time-frequency energy graph;
FIG. 3 illustrates a fractal variation law;
FIG. 4 shows the law of change of time-frequency energy spectrum entropy with order;
FIG. 5 shows the recognition rate as a function of the number of training samples;
fig. 6 shows a classification effect diagram.
Detailed Description
The technical scheme and the implementation process of the invention are described in detail by combining specific examples.
The process flow of the method of the invention is shown in figure 1.
First by approximating l0And obtaining a time-frequency energy spectrum curved surface of the frequency hopping signal by a norm (AL0) algorithm, and then extracting a difference box dimension, a multi-fractal dimension and a time-frequency Rayleigh entropy under different block scales to form a feature vector. And finally, training, identifying and classifying through a support vector machine classifier.
The first step is as follows: extracting frequency hopping signal time frequency energy spectrum
Assuming a finite set of frequencies for the frequency hopping signal as w, the received signal frequency is set as
Figure GSB0000194869010000081
M is more than or equal to 0 and less than or equal to P-1, according to the requirement of the invention on time precision, the frequency hopping section is divided into P frequency points at equal intervals, the received frequency hopping signal y is divided into G sections at equal intervals, and each section of data ykHas a length of P, then
yk＝y(k·L：k·L+P-1) (1)
Wherein k is more than or equal to 1 and less than or equal to G, and L is a segmentation interval.
Observation matrix Y ═ Y1，y2，...，yG]Can be expressed as
Y＝WX+V (2)
Wherein W is [ omega ]0，...，ωP-1]，ωi＝[ejω1，...，ejωP]X is a time-frequency matrix, V is a noise matrix, and X, V are belonged to CP×GC represents a complex momentAnd (5) arraying.
Due to the sparsity of the frequency hopping signals in the time-frequency domain, time-frequency points in X are sparse, non-zero points are all on the corresponding rows of the frequency hopping points, and X is also sparse. An unconstrained optimization function with a penalty function can therefore be assumed:
Figure GSB0000194869010000091
wherein x iskWhich represents the k-th row of X,
Figure GSB0000194869010000092
denotes the value of X, μ, which minimizes L (X)1Is an X-point sparsity penalty factor, μ2Is an X-row sparse penalty factor.
The AL0 algorithm introduces a gaussian function to approximate l0Norm of
Figure GSB0000194869010000093
Where s represents an arbitrary D-dimensional vector and σ represents a normal number close to zero. When σ is approximately 0, there are
Figure GSB0000194869010000094
Wherein s isiThe ith row of the vector s is indicated and D indicates the dimension of the vector s. Can approximate the formula (4) to l0The norm represents the sparsity of the time-frequency domain of the frequency hopping signal, and then the sparse reconstruction problem can be converted into the minimization problem for solving the Gaussian sum function
Figure GSB0000194869010000095
Wherein Fσ(x) Indicating a frequency hopping signal0Norm, introducing a penalty factor lambda, and dividing l0The norm minimization problem translates into an optimization problem, i.e.
Figure GSB0000194869010000096
Wherein
Figure GSB0000194869010000097
Denotes the value of X which minimizes L (X).
The steepest descent direction is the direction of the conjugate gradient of the unconstrained optimization function L (x)
Figure GSB0000194869010000098
Figure GSB0000194869010000099
Wherein x*Is the conjugate of x. The definition of the conjugate gradient of the complex variable can be obtained
Figure GSB00001948690100000910
Wherein x isRRepresenting the real part, x, of a signal xIRepresenting the imaginary part of the signal x, i being the unit of an imaginary number, diagonal matrix Λi，i＝exp(-|xk|2/2σ2)/σ2。
According to
Figure GSB0000194869010000107
It can be seen that the superscript H represents the conjugate matrix.
Figure GSB0000194869010000101
By substituting the formulae (7) and (8) into the formula (6)
Figure GSB0000194869010000102
Then the direction of the conjugate gradient of L (x) in equation (3) is calculated as
Figure GSB0000194869010000103
Wherein, Λ1＝fσ(X)/σ2；Λ2Is a diagonal matrix whose k-th diagonal element is
Figure GSB0000194869010000104
K is not less than 1 and not more than G, and X (k,: represents the kth line of X. The detailed algorithm steps of the AL0 algorithm are as follows, where σ takes the values of [2, 1, 0.5, 0.2, 0.1, 0.05]：
Inputting: matrix W, measurement vector y.
Step 1: let initial value x(0)＝WT(WWT)-1y。
Step 2: selection of descending sequence [ sigma ]1σ2...σJ]In one embodiment of the invention is the sequence [2, 1, 0.5, 0.2, 0.1, 0.05]]And J is 6. Let the convergence criterion be ε, ε representing a mean of 0 and a variance of σ2White gaussian noise.
Step 3: and (6) iterating the algorithm.
(1) J takes the values of 1, 2,.. J in sequence;
(2) for each value of j, let σ be σ ═ σj，
(3) When in use
Figure GSB0000194869010000105
Where norm represents the solution
Figure GSB0000194869010000106
Figure GSB0000194869010000106
2 norm of (d);
(4) step size u is determined such that
Figure GSB0000194869010000111
(5) Is updated in the gradient direction such that
Figure GSB0000194869010000112
(6) At this time, let x(j)And continuing to take the next value of j and carrying out iterative operation until the values are sequentially subjected to the iterative operationAll values are taken out.
Step 4: outputting the result
Figure GSB0000194869010000113
The time-frequency energy spectrum of two synchronous networking frequency hopping radio stations obtained through sparse reconstruction of the AL0 algorithm is shown in figure 2.
As can be seen from fig. 2, two frequency hopping stations of the same type and operation carry respective fine feature information due to their different radiation sources. The individual identification of the frequency hopping radio station can be realized by extracting the characteristic rules of energy change, frequency point distribution and the like of the time-frequency energy spectrum of each radio station signal.
The second step is that: extracting fractal characteristics;
(1) differential box dimension feature extraction
And (3) regarding the time-frequency axis of the frequency-hopping signal time-frequency energy spectrum subjected to sparse reconstruction as a plane, regarding the energy spectrum value as an image gray value, and calculating the difference box dimension of the time-frequency energy spectrum by the following steps:
step 1: assuming that the size of the time-frequency energy spectrum F is NxN, S is a grid of MxM, wherein M is more than or equal to 1 and less than or equal to N/2, and M belongs to Z+，Z+Representing a set of positive integers. F is divided into a plurality of sub-blocks by S, each grid S is provided with a box column with the size of M multiplied by M ', M' is more than or equal to 1 and less than or equal to N/2,
Figure GSB0000194869010000114
e represents the maximum energy spectrum value and,
Figure GSB0000194869010000115
indicating the rounding-down operator.
Step 2: setting the (mu, v) th grid, wherein mu is more than or equal to 1 and less than or equal to M, and v is more than or equal to 1 and less than or equal to M. The maximum and minimum energy values of the energy spectrum are denoted g, respectivelymaxAnd gminAt the scale M, the box count of the (μ, ν) th grid is:
Figure GSB0000194869010000116
step 3: by analogy, the number of all grid boxes under the M scale is calculated
Figure GSB0000194869010000121
Step 4: changing the scale M of the grid S, repeating the Step1-3, and calculating N under the condition of different block scales MM。
Step 5: estimating differential box dimensions
Figure GSB0000194869010000122
Wherein the content of the first and second substances,
Figure GSB0000194869010000123
(2) multi-fractal feature extraction
Many scholars find that when single fractal describes most of objectively existing fractal objects, the complexity and nonlinear characteristics of the fractal objects cannot be completely measured, and if the overall characteristics of a system are described by using the single fractal dimension, all the characteristics of the fractal objects cannot be fully reflected to a certain extent. The multi-fractal can study the overall characteristics of the multi-fractal from local parts, and improves the fine degree of describing the geometric characteristics and local scale behaviors of the object. In order to more effectively represent the local characteristics of the time-frequency energy spectrum, the multi-fractal dimension of the time-frequency energy spectrum is extracted as a second-dimensional characteristic.
The fractal MD can be expressed as:
Figure GSB0000194869010000124
in the formula, q represents a multi-fractal parameter, and in practical application, different choices of q values influence the stability of MD. According to Step1-3 of the difference box dimension extraction algorithm, calculating N under the condition of different block scales MMThen, combining with the formula (14), obtaining the multi-fractal characteristic of the frequency-hopping signal time-frequency energy spectrum。
(3) Time-frequency energy spectrum Rayleigh entropy feature extraction
Entropy can be defined as a quantity that measures chaotic irregularities, uncertainty, etc. chaotic states. As can be seen from fig. 2, the time-frequency energy distribution rules of the two frequency hopping radio station signals have a certain difference, and the difference of the distribution can be measured by counting the rayleigh entropy of the time-frequency energy spectrum, that is, the difference is converted from the intuitive difference to the numerical difference.
Let the probability density distribution of the random variable J be P ═ { P ═ P1，p2，...，pTWhere T represents the dimension of the random variable J, satisfying the condition
Figure GSB0000194869010000131
T is more than or equal to 1 and less than or equal to T, the Rayleigh entropy of J is defined as follows:
Figure GSB0000194869010000132
where α is the rayleigh entropy order, and equation (15) defines rayleigh entropy using a one-dimensional probability density distribution.
The continuous form of the rayleigh entropy of order α is as follows:
Figure GSB0000194869010000133
wherein f (χ, γ) is a continuous two-dimensional probability density distribution, fα(χ, γ) represents f (χ, γ) to the α -power.
After the frequency hopping signal y is sparsely reconstructed, the obtained time-frequency energy spectrum P (t, f) has time-frequency edge characteristics and energy retention characteristics, as shown in formulas (17) and (18):
∫P(t，f)df＝|y(t)|2，∫P(t，f)dt＝|y(f)|2 (17)
Figure GSB0000194869010000134
where y (f) is the Fourier transform of y (t), and f represents the signal frequency. Therefore, the two-dimensional probability density distribution f (χ, γ) of the time-frequency energy spectrum P (t, f) and the equation (16) has similar properties, so that the rayleigh entropy of the time-frequency energy spectrum of the frequency hopping signal can be defined by P (t, f) as shown in the formula:
Figure GSB0000194869010000135
in the formula Pα(t, f) represents P (t, f) to the power of alpha.
The stabilizing condition of formula (19) is ^ jeppeα(t, f) dtdf > 0. For convenience of calculation, the discrete expression of the rayleigh entropy of the time-frequency energy spectrum is as follows:
Figure GSB0000194869010000141
wherein X is a time-frequency matrix, k and k 'represent the number of sampling points of observation data, k is more than or equal to 1 and less than or equal to G, and k' is more than or equal to 1 and less than or equal to G;
psi and psi' denote frequency set dictionary, omega0≤ψ′≤ωP-1，ω0≤ψ≤ωP-1。
And finally, forming a characteristic vector V [ FD, MD (q), R from the three characteristics of the Rayleigh entropy, the difference box dimension and the multi-fractal dimension of the calculated time-frequency energy spectrumα]And carrying out classification and identification by using a support vector machine classifier.
Example testing
The experimental data of the invention is collected from four frequency hopping radio stations with the same model, the working frequency of the radio stations is set to 150MHz, the frequency hopping bandwidth is 6.4MHz, the sampling rate is 600MHz, and the sampling duration is 5 seconds. Considering that when the time-frequency energy spectrum after sparse reconstruction is partitioned, if the value of the dimension L is too large, the boundary information of the energy spectrum cannot be completely utilized, so that the extracted features cannot completely reflect the essential information of the signal; if the value of the scale L is too small, the amount of calculation increases. Therefore, in the simulation experiment of the invention, the block scale L takes five values of 4, 8, 12, 16 and 20.
As can be seen from equation (14), the parameter q affects the stability of the md (q) values, and md (q) under different q values is calculated by using the data samples of the frequency hopping signal of one of the stations, and the change rule is shown in fig. 3.
As can be seen from FIG. 3, MD (q) is less affected by the value of q only when 10 < q. Through repeated experimental observation, the md (q) values of the hopping signals of the 4 stations all become stable when q is 13, and q is 13 in order to ensure the stability of the extracted features.
The rayleigh entropy of the time-frequency energy spectrum can well reflect the time-frequency energy spectrum rule and complexity of each frequency hopping signal, but as can be known from the formula (20), the extraction of the rayleigh entropy characteristics of the time-frequency energy spectrum is greatly influenced by the value of the order alpha. The invention selects the frequency hopping signal sample data of 4 radio stations, and respectively calculates the Rayleigh entropy of the energy spectrum of the integer order of 0-30 of the order alpha, as shown in figure 4.
As can be seen from fig. 4, as the order α increases, the rayleigh entropy of the time-frequency energy spectrum of the hopping signal of the 4 radio stations tends to become substantially smaller, and the rayleigh entropy of the time-frequency energy spectrum of the hopping signal of the 4 radio stations gradually converges to the same value. When alpha is 14, the rayleigh entropy of the time-frequency energy spectrum of the hopping signals of different stations has a certain discrimination, and when alpha is more than 16, the discrimination becomes more and more fuzzy. Fig. 4 shows that when α is 3, the discrimination is the best, and as the order α increases, the discrimination is continuously reduced, so the invention selects the rayleigh entropy of the 3-order time-frequency energy spectrum with higher discrimination as the characteristic of individual identification of the frequency hopping radio station.
Respectively intercepting 4 radio station frequency hopping signal data, starting from the head of the sampled data, intercepting 1024 sampling points as training samples at intervals of 200 sampling points, and intercepting 200 sections. Then, in the same manner, 200 pieces of data are cut out as test samples from the end of the sampling data. Under the condition of different training sample numbers, the algorithm is adopted to extract three characteristics of the time-frequency energy spectrum Rayleigh entropy, the difference box dimension and the multi-fractal dimension, then a support vector machine classifier is used for classification and identification, and the average identification accuracy of 50 experiments is counted as shown in figure 5.
As can be seen from FIG. 5, the recognition accuracy of the algorithm of the present invention is not affected by the variation of the number of training samples, and is always kept at 85.6%. The main reason is that the algorithm respectively extracts the difference box dimension, the multi-fractal dimension and the Rayleigh entropy of the time-frequency energy spectrum as the feature vectors, quantitatively describes the energy change of the time-frequency energy spectrum, the complexity and the regularity of the time-frequency energy spectrum, and avoids the problem of misjudgment caused by the similarity of single features.
In order to facilitate visualization of the classification effect, the number of the training samples is 200, and the classification effect of the three-dimensional features is drawn by using the algorithm of the invention and is shown in fig. 6.
As can be seen from FIG. 6, the classification features of the algorithm of the present invention are only rarely partially overlapped, and have an obvious clustering effect, which is in accordance with the recognition result of FIG. 5.

Claims (3)

1. A frequency hopping radio station individual identification method comprises the following steps:
the first step is as follows: extracting frequency hopping signal time frequency energy spectrum
Assuming a finite set of frequencies for the frequency hopping signal as w, the received signal frequency is set as
Figure FSB0000194829000000011
M is more than or equal to 0 and less than or equal to P-1, a frequency hopping section is divided into P frequency points at equal intervals according to the requirement on time precision, a received frequency hopping signal y is divided into G sections at equal intervals, and each section of data ykHas a length of P, then
yk＝y(k·L：k·L+P-1) (1)
Wherein k is more than or equal to 1 and less than or equal to G, and L is a segment interval;
observation matrix Y ═ Y1，y2，...，yG]Can be expressed as
Y＝WX+V (2)
Wherein W is [ omega ]0，...，ωP-1]，ωi＝[ejω1，...，ejωP]X is a time-frequency matrix, V is a noise matrix, and X, V are belonged to CP ×GC represents a complex matrix;
due to the sparsity of the frequency hopping signals in the time-frequency domain, time-frequency points in the X are sparse, non-zero points are all on the corresponding rows of each frequency hopping point, and the X is also sparse; therefore, an unconstrained optimization function with a penalty function is assumed:
Figure FSB0000194829000000012
Figure FSB0000194829000000013
wherein x iskWhich represents the k-th row of X,
Figure FSB0000194829000000014
denotes the value of X, μ, which minimizes L (X)1Is an X-point sparsity penalty factor, μ2Is an X row sparse penalty factor;
the AL0 algorithm introduces a gaussian function to approximate l0Norm of
Figure FSB0000194829000000015
Where s represents an arbitrary D-dimensional vector and σ represents a normal number close to zero; when σ is approximately 0, there are
Figure FSB0000194829000000016
Wherein s isiRepresents the ith row of vector s, D represents the dimension of vector s; approximating formula (4) to l0The norm represents the sparsity of the time-frequency domain of the frequency hopping signal, and then the sparse reconstruction problem can be converted into the minimization problem for solving the Gaussian sum function
Figure FSB0000194829000000021
Wherein Fσ(x) Indicating a frequency hopping signal0Norm, introducing a penalty factor lambda, and dividing l0The norm minimization problem translates into an optimization problem, i.e.
Figure FSB0000194829000000022
Wherein
Figure FSB0000194829000000023
Denotes the value of X which minimizes L (X);
the steepest descent direction is the direction of the conjugate gradient of the unconstrained optimization function L (x)
Figure FSB0000194829000000024
Figure FSB0000194829000000025
Wherein x*Is the conjugate of x; the definition of the conjugate gradient of the complex variable can be obtained
Figure FSB0000194829000000026
Wherein x isRRepresenting the real part, x, of a signal xIRepresenting the imaginary part of the signal x, i being the unit of an imaginary number, diagonal matrix Λi，i＝exp(-|xk|2/2σ2)/σ2；
According to
Figure FSB0000194829000000027
One can obtain, where the upper superscript H represents the conjugate matrix;
Figure FSB0000194829000000028
by substituting the formulae (7) and (8) into the formula (6)
Figure FSB0000194829000000029
Then the direction of the conjugate gradient of L (x) in equation (3) is calculated as
Figure FSB0000194829000000031
Wherein, Λ1＝fσ(X)/σ2；Λ2Is a diagonal matrix whose k-th diagonal element is
Figure FSB0000194829000000032
K is not less than 1 and not more than G, and X (k, k) represents the kth line of X;
the second step is that: fractal feature extraction
(1) Differential box dimension feature extraction
And (3) regarding the time-frequency axis of the frequency-hopping signal time-frequency energy spectrum subjected to sparse reconstruction as a plane, regarding the energy spectrum value as an image gray value, and calculating the difference box dimension of the time-frequency energy spectrum by the following steps:
step 1: assuming that the size of the time-frequency energy spectrum F is NxN, S is a grid of MxM, wherein M is more than or equal to 1 and less than or equal to N/2, and M belongs to Z+，Z+Representing a set of positive integers; f is divided into a plurality of sub-blocks by S, each grid S is provided with a box column with the size of M multiplied by M ', M' is more than or equal to 1 and less than or equal to N/2,
Figure FSB0000194829000000033
e represents the maximum energy spectrum value and,
Figure FSB0000194829000000034
represents a round-down operator;
step 2: setting the (mu, v) grid, wherein mu is more than or equal to 1 and less than or equal to M, and v is more than or equal to 1 and less than or equal to M; the maximum and minimum energy values of the energy spectrum are denoted g, respectivelymaxAnd gminAt scale M, the number of boxes in the (μ, v) th grid is:
Figure FSB0000194829000000035
step 3: by analogy, the number of all grid boxes under the M scale is calculated
Figure FSB0000194829000000036
Step 4: changing the scale M of the grid S, repeating the Step1-3, and calculating N under the condition of different block scales MM；
Step 5: estimating differential box dimensions
Figure FSB0000194829000000037
Wherein the content of the first and second substances,
Figure FSB0000194829000000041
(2) multi-fractal feature extraction
In order to more effectively represent the local characteristics of the time-frequency energy spectrum, extracting a multi-fractal dimension of the time-frequency energy spectrum as a second-dimensional characteristic;
the fractal MD can be expressed as:
Figure FSB0000194829000000042
in the formula, q represents a multi-fractal parameter, and in practical application, different choices of q values influence the stability of MD; according to Step1-3 of the difference box dimension extraction algorithm, calculating N under the condition of different block scales MMThen, combining with the formula (14), obtaining the multi-fractal characteristic of the frequency-hopping signal time-frequency energy spectrum;
(3) time-frequency energy spectrum Rayleigh entropy feature extraction
The time-frequency energy distribution rules of the two frequency hopping radio station signals have certain difference, and the difference of the distribution can be measured by counting the Rayleigh entropy of the time-frequency energy spectrum, namely the difference is converted from the intuitive difference to the difference in numerical value;
let the probability density distribution of the random variable J be P ═ { P ═ P1，p2，...，pTWhere T represents the dimension of the random variable J, satisfying the condition
Figure FSB0000194829000000043
T is more than or equal to 1 and less than or equal to T, the Rayleigh entropy of J is defined as follows:
Figure FSB0000194829000000044
wherein alpha is the Rayleigh entropy order, and the Rayleigh entropy is defined by applying one-dimensional probability density distribution in the formula (15);
the continuous form of the rayleigh entropy of order α is as follows:
Figure FSB0000194829000000051
wherein f (χ, γ) is a continuous two-dimensional probability density distribution, fα(χ, γ) represents f (χ, γ) to the power α;
after the frequency hopping signal y is sparsely reconstructed, the obtained time-frequency energy spectrum P (t, f) has time-frequency edge characteristics and energy retention characteristics, as shown in formulas (17) and (18):
∫P(t，f)df＝|y(t)|2，∫P(t，f)dt＝|y(f)|2 (17)
Figure FSB0000194829000000052
wherein y (f) is the Fourier transform of y (t), and f represents the signal frequency; therefore, the two-dimensional probability density distribution f (χ, γ) of the time-frequency energy spectrum P (t, f) and the equation (16) has similar properties, so that the rayleigh entropy of the time-frequency energy spectrum of the frequency hopping signal can be defined by P (t, f) as shown in the formula:
Figure FSB0000194829000000053
in the formula Pα(t, f) represents P (t, f) to the power of alpha;
the stabilizing condition of formula (19) is ^ jeppeα(t, f) dtdf > 0; for convenience of calculation, the discrete expression of the rayleigh entropy of the time-frequency energy spectrum is as follows:
Figure FSB0000194829000000054
wherein X is a time-frequency matrix, k and k 'represent the number of sampling points of observation data, k is more than or equal to 1 and less than or equal to G, and k' is more than or equal to 1 and less than or equal to G; psi and psi' denote frequency set dictionary, omega0≤ψ′≤ωP-1，ω0≤ψ≤ωP-1；
And finally, forming a characteristic vector V [ FD, MD (q), R from the three characteristics of the Rayleigh entropy, the difference box dimension and the multi-fractal dimension of the calculated time-frequency energy spectrumα]And carrying out classification and identification by using a support vector machine classifier.
2. The individual identification method of frequency hopping radio of claim 1, wherein AL0 algorithm steps are as follows:
inputting: matrix W, measured value vector y;
step (1): let initial value x(0)＝WT(WWT)-1y；
Step (2): selection of descending sequence [ sigma ]1σ2...σJ]J ═ 6; let the convergence criterion be ε, ε representing a mean of 0 and a variance of σ2White gaussian noise of (1);
step (3): iterating an algorithm;
(1) j takes the values of 1, 2,.. J in sequence;
(2) for each value of j, let σ be σ ═ σj，
(3) When in use
Figure FSB0000194829000000061
Where norm represents the solution
Figure FSB0000194829000000062
2 norm of (d);
(4) step size u is determined such that
Figure FSB0000194829000000063
(5) Is updated in the gradient direction such that
Figure FSB0000194829000000064
(6) At this time, let x(j)Continuing to take the next value of j as x, and carrying out iterative operation until all values are taken in sequence;
step 4: outputting the result
Figure FSB0000194829000000065
3. The individual identification method of a frequency hopping radio station as claimed in claim 2, wherein σ takes values of [2, 1, 0.5, 0.2, 0.1, 0.05] in sequence.
CN201811328578.4A 2018-10-28 2018-10-28 Individual identification method of frequency hopping radio station Active CN109472239B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811328578.4A CN109472239B (en) 2018-10-28 2018-10-28 Individual identification method of frequency hopping radio station

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811328578.4A CN109472239B (en) 2018-10-28 2018-10-28 Individual identification method of frequency hopping radio station

Publications (2)

Publication Number Publication Date
CN109472239A CN109472239A (en) 2019-03-15
CN109472239B true CN109472239B (en) 2021-10-01

Family

ID=65672201

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811328578.4A Active CN109472239B (en) 2018-10-28 2018-10-28 Individual identification method of frequency hopping radio station

Country Status (1)

Country Link
CN (1) CN109472239B (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP4227691A1 (en) * 2022-02-10 2023-08-16 Rohde & Schwarz GmbH & Co. KG Method of classifying a radio frequency signal

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110472584A (en) * 2019-08-16 2019-11-19 四川九洲电器集团有限责任公司 A kind of communication equipment personal identification method, electronic equipment and computer program product
CN112994740B (en) * 2021-04-23 2021-07-23 成都天锐星通科技有限公司 Frequency hopping signal parameter estimation method and device, electronic equipment and readable storage medium

Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102437984A (en) * 2011-11-07 2012-05-02 哈尔滨工程大学 Modulation signal identification method based on complexity characteristic under low signal-to-noise ratio condition
CN203151487U (en) * 2013-04-11 2013-08-21 中国人民解放军装甲兵工程学院 Performance test system for frequency hopping radio station
CN103340625A (en) * 2013-06-18 2013-10-09 中国人民解放军第四军医大学 Regularization method of fast optimization in electrical impedance tomography
CN103959839A (en) * 2011-08-12 2014-07-30 黑莓有限公司 Methods of channel state information feedback and transmission in coordinated multi-point wireless communications system
CN104485979A (en) * 2014-12-09 2015-04-01 西安电子科技大学 Blind estimation method for underdetermined hybrid frequency hopping parameters based on time frequency diagram correction
CN105891825A (en) * 2016-03-29 2016-08-24 西安电子科技大学 Multiple-input multiple-output array radar staring imaging method based on tensor compression perception
CN106908774A (en) * 2017-01-06 2017-06-30 南京航空航天大学 Based on the sparse one-dimensional range profile recognition methods for keeping projecting of multiple dimensioned core
CN107526074A (en) * 2017-07-19 2017-12-29 上海无线电设备研究所 A kind of distance of sparse Frequency Hopping Signal and Speed Two Dimensions high resolution processing method
CN107894591A (en) * 2017-09-30 2018-04-10 沈阳航空航天大学 Through-wall radar diffraction tomography method based on compressed sensing
CN108700652A (en) * 2015-12-09 2018-10-23 欧利景无线有限公司 The methods, devices and systems for detecting and monitoring for wireless event

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20140181166A1 (en) * 2012-12-26 2014-06-26 Industrial Technology Research Institute Apparatus for low complexity sub-nyquist sampling of sparse wideband signals

Patent Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103959839A (en) * 2011-08-12 2014-07-30 黑莓有限公司 Methods of channel state information feedback and transmission in coordinated multi-point wireless communications system
CN102437984A (en) * 2011-11-07 2012-05-02 哈尔滨工程大学 Modulation signal identification method based on complexity characteristic under low signal-to-noise ratio condition
CN203151487U (en) * 2013-04-11 2013-08-21 中国人民解放军装甲兵工程学院 Performance test system for frequency hopping radio station
CN103340625A (en) * 2013-06-18 2013-10-09 中国人民解放军第四军医大学 Regularization method of fast optimization in electrical impedance tomography
CN104485979A (en) * 2014-12-09 2015-04-01 西安电子科技大学 Blind estimation method for underdetermined hybrid frequency hopping parameters based on time frequency diagram correction
CN108700652A (en) * 2015-12-09 2018-10-23 欧利景无线有限公司 The methods, devices and systems for detecting and monitoring for wireless event
CN105891825A (en) * 2016-03-29 2016-08-24 西安电子科技大学 Multiple-input multiple-output array radar staring imaging method based on tensor compression perception
CN106908774A (en) * 2017-01-06 2017-06-30 南京航空航天大学 Based on the sparse one-dimensional range profile recognition methods for keeping projecting of multiple dimensioned core
CN107526074A (en) * 2017-07-19 2017-12-29 上海无线电设备研究所 A kind of distance of sparse Frequency Hopping Signal and Speed Two Dimensions high resolution processing method
CN107894591A (en) * 2017-09-30 2018-04-10 沈阳航空航天大学 Through-wall radar diffraction tomography method based on compressed sensing

Non-Patent Citations (7)

* Cited by examiner, † Cited by third party
Title
Frequency hopping radio individual identification based on energy spectrum blended subtle characteristics;Yang Xin 等;《Journal of Physics: Conference Series》;20190707;第1325卷;第1-7页 *
Specific emitter identification based on normalized frequency spectrum;Degang Sun 等;《2016 2nd IEEE International Conference on Computer and Communications》;20161017;第1192-1205页 *
Specific Emitter Identification via Hilbert–Huang Transform in Single-Hop and Relaying Scenarios;Jingwen Zhang 等;《IEEE Transactions on Information Forensics and Security》;20160630;第11卷(第6期);第1192-1205页 *
Structure-Aware Bayesian Compressive Sensing for Frequency-Hopping Spectrum Estimation With Missing Observations;Shengheng Liu 等;《IEEE Transactions on Signal Processing》;20180415;第66卷(第8期);第2153-2166页 *
块稀疏贝叶斯模型下的跳频信号时频分析;李雷;《信号处理》;20180131;第34卷(第01期);第107-113页 *
基于ST-RFT算法的信号调制方式识别与参数估计方法研究;刘丹;《中国优秀硕士学位论文全文数据库 信息科技辑》;20171015(第10期);第I136-60页 *
基于时频能量谱特征的跳频电台个体识别;杨鑫 等;《信号处理》;20191031;第35卷(第10期);第1671-1679页 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP4227691A1 (en) * 2022-02-10 2023-08-16 Rohde & Schwarz GmbH & Co. KG Method of classifying a radio frequency signal

Also Published As

Publication number Publication date
CN109472239A (en) 2019-03-15

Similar Documents

Publication Publication Date Title
CN109472239B (en) Individual identification method of frequency hopping radio station
Jamali et al. Identification of optimal features for fast and accurate classification of power quality disturbances
CN109307862A (en) A kind of target radiation source individual discrimination method
Wang et al. Fractal complexity-based feature extraction algorithm of communication signals
Iorkyase et al. Improving RF-based partial discharge localization via machine learning ensemble method
Liu et al. Deep learning and recognition of radar jamming based on CNN
CN105913081B (en) SAR image classification method based on improved PCAnet
CN112213688B (en) Feature extraction method for identifying low-altitude airspace low-small slow aircraft target individuals
CN108090462B (en) Radiation source fingerprint feature extraction method based on box dimensions
Ghadimi et al. Deep learning-based approach for low probability of intercept radar signal detection and classification
CN108734228A (en) The polarimetric SAR image random forest classification method of comprehensive multiple features
CN111680737B (en) Radar radiation source individual identification method under differential signal-to-noise ratio condition
CN110244275A (en) The reconstruct of sea clutter bispectrum and emulation mode
CN116047427B (en) Small sample radar active interference identification method
CN115034261A (en) Method and equipment for extracting features between pulses of radar radiation source signal and storage medium
CN109446910A (en) A kind of communication emitter Signals classifying identification method
Charalampidis et al. Removal of nonprecipitation echoes in weather radar using multifractals and intensity
CN114841195B (en) Avionics space signal modeling method and system
Shreyamsha Kumar et al. Target identification using harmonic wavelet based ISAR imaging
Melinda et al. Implementation of segmentation scheme based on wavelet transform in multi-spectral fluctuation patterns
CN106125073B (en) Scattering mechanism identification and extracting method based on adaptive Gauss expression
Shan et al. Wavelet based recognition for pulsar signals
Saha et al. Curvelet entropy for facial expression recognition
Yang et al. Frequency hopping radio individual identification based on energy spectrum blended subtle characteristics
CN111083632A (en) Ultra-wideband indoor positioning method based on support vector machine

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant