CN103064063A - Poly-phase code radar signal waveform automatic identification method based on continuous wave Doppler (CWD) feature - Google Patents

Abstract

The invention relates to a poly-phase code radar signal waveform automatic identification method based on continuous wave Doppler (CWD) feature and belongs to the technical field of information countermeasures. According to the poly-phase code radar signal waveform automatic identification method based on CWD feature, a discrete sampling type Choi-Williams conversion is used as a basic tool, a CWD image of poly-phase code pulse compression radar signal is used as a feature extraction object, a Pseudo-Zernike moment of the CWD image, a target number of an image, a time position of a peak power in a CWD and a poly-phase code waveform symmetrical property are taken as features to identify a poly-phase code pulse compression radar waveform, and a neural network which is composed of 10 perceptrons and stops in advance of overall average is set up for automatic waveform identification. The poly-phase code radar signal waveform automatic identification method based on CWD feature has the advantages that the identification precision of the poly-phase code radar signal waveform is improved, the requirement of a signal to noise ratio is further decreased, and a new way can be developed for the design of radar signal identification and sorting by popularized and used in a poly-phase code continuous wave radar signal.

Description

Heterogeneous code radar signal waveform automatic identifying method based on the CWD feature

Technical field

The present invention relates to a kind of heterogeneous code radar signal waveform automatic identifying method based on the CWD feature, belong to the information countermeasure technical field.

Background technology

In recent years, along with deepening continuously that heterogeneous code radar signal is studied, Frank code and P1, P2, P3 and P4 coded signal have appearred.In a sense, Frank code and P1～P4 coded signal all belongs to heterogeneous coded signal (Polyphase Coded).Traditional melodeon can't be realized the effective identification to it, and how research effectively identified this class signal has important theoretical significance and engineering application is worth.

The concept of Fractional Fourier Transform namely was suggested as far back as nineteen twenty-nine, applied to optical field in the eighties in 20th century, became one of study hotspot of signal process field from the nineties.Fractional Fourier Transform can be understood as the Chirp base decomposes, and therefore is particularly suitable for processing Chirp class signal.Utilize linear frequency modulation (LFM) signal (namely Chirp signal) to present the characteristic of different energy accumulatings at the Fractional Fourier Domain of different orders, just can realize detection and parameter estimation to the LFM signal by do the peak value two-dimensional search at Fractional Fourier Domain.And the triangle linear frequency modulation continuous wave can be thought to be comprised of many LFM Signal, and Fractional Fourier Transform is linear transformation, does not have cross term to disturb, and has more advantage in the situation that have additive noise.

Present method for parameter estimation to the triangle linear frequency modulation continuous wave, mostly can only estimate the partial parameters of signal, what have requires signal to noise ratio (S/N ratio) higher, and the requirement that also has just in time is a modulation period to signal sampling, and these methods all do not have well to solve the characteristic parameter extraction problem of this type of signal.

Summary of the invention

The present invention relates to a kind of heterogeneous code radar signal waveform automatic identifying method based on the CWD feature, belong to the information countermeasure technical field.The present invention is transformed to basic tool with discrete sampling type Choi-Williams, with the CWD image of heterogeneous coded pulse compression radar signal as the feature extraction object, proposition is with the time location of peak power in the Pseudo-Zernike square of CWD image, the target number, CWD in the image and the polyphase code waveform symmetry character feature as the heterogeneous coded pulse compression radar waveform of identification, utilize the difference of polyphase code signal on feature, set up the neural network that a population mean that is comprised of 10 perceptrons stops in advance and carry out automatic waveform recognition.The method that the present invention proposes has improved the recognition accuracy of heterogeneous code radar signal waveform, further reduced the requirement of signal to noise ratio (S/N ratio), and can be through promoting the use in heterogeneous coded continous wave radar signal, for the design of Radar Signal Recognition and sorting provides a new approach.

Based on the heterogeneous code radar signal waveform automatic identifying method of CWD feature, be divided into normalized, signal characteristic abstraction and 3 parts of signal classifier design of heterogeneous code radar signal discrete CWD image, amount to 4 treatment steps (Fig. 1).

Beneficial effect

1. the present invention propose based on the heterogeneous code radar signal waveform automatic identifying method of CWD feature can pin-point accuracy the identification of the heterogeneous code radar signal of realization.

2. the present invention does not have specific (special) requirements substantially to the sampling time of signal, can begin from any time of signal to analyze, and still has higher recognition accuracy under low signal-to-noise ratio.

3. the parallel multilayer neural network waveform recognition device of the present invention's proposition can apply to the identification of heterogeneous code radar signal, can also apply to the identification of other radar signal, and key is choosing of signal characteristic.

Description of drawings

Fig. 1 is based on the automatic identifying schemes schematic diagram of heterogeneous code radar signal waveform of CWD feature

Five kinds of waveform recognition accuracy of Fig. 2 and total recognition correct rate

Embodiment

Suppose that radar signal that melodeon receives mixes and additive white Gaussian noise arranged (AWGN), and signal is as follows through processing the complex envelope that becomes the radar signal y (t) that baseband signal then receives:

y(t)＝x(t)+ω(t) (1)

Wherein x (t) is the complex envelope (including only a code-element period) of radar emission signal, and ω (t) is circulation additivity white complex gaussian noise.

The phase encoding complex signal of radar emission is expressed as follows:

$x\left(t\right)={\mathrm{Ae}}^{j(2\mathrm{\π}{f}_{c}t+{\mathrm{\φ}}_i)}---\left(2\right)$

Wherein, A is signal amplitude, f _cSignal(-) carrier frequency, φ _iBe the signal discrete phase sequence, each phase place has the identical duration.The below provides the mathematical model of 5 kinds of heterogeneous code radar signals will studying.

The Frank coded signal is that a kind of stepping to the LFM signal approaches, and it has adopted N step frequency, and carries out N sampling at each Frequency point.Therefore, total hits of a Frank code is N ²The phase place of i sampling of j frequency of a Frank code is as follows:

${\mathrm{\φ}}_{i,j}=\frac{2\mathrm{\π}}{N}(i-1)(j-1)---\left(3\right)$

Wherein, i=1,2 ..., N, j=1,2 ..., N.Then the pulse compression ratio of this Frank code is N ²

The P1 coded signal also is that the stepping to the LFM signal approaches, and it has adopted N step frequency, and carries out N sampling at each Frequency point.The phase place of i sampling of j frequency of a P1 code is as follows:

${\mathrm{\φ}}_{i,j}=-\frac{\mathrm{\π}}{N}[N-(2j-1)][(j-1)N+(i-1)]---\left(4\right)$

Wherein, i=1,2 ..., N, j=1,2 ..., N.Then the pulse compression ratio of this P1 code is N ²

The P2 code has the characteristics of the palindrome, and namely the positive and negative of P2 code is identical.The phase place of i sampling of j frequency of a P2 code is as follows:

${\mathrm{\φ}}_{i,j}=[\frac{\mathrm{\π}}{2}\·\frac{N-1}{N}\·\frac{\mathrm{\π}}{N}(i-1)](N+1-2j)---\left(5\right)$

Wherein, i=1,2 ..., N, j=1,2 ..., N.Then the pulse compression ratio of this P2 code is N ²It is pointed out that N is necessary for even number in the P2 code, if N is odd number, then the numerical value of the autocorrelation sidelobe of this P2 code can be too high.

The P3 code is from evolution that a LFM signal is sampled.The phase place of i sampling of a P3 code is as follows:

${\mathrm{\φ}}_i=\frac{\mathrm{\π}}{\mathrm{\ρ}}{(i-1)}^2---\left(6\right)$

I=1 wherein, 2 ..., ρ, ρ are pulse compression ratio.

The P4 code is from evolution that the identical signal of P3 code is sampled.The phase place of i sampling of a P4 code is as follows:

${\mathrm{\φ}}_i=\frac{\mathrm{\π}}{\mathrm{\ρ}}{(i-1)}^2-\mathrm{\π}(i-1)---\left(7\right)$

I=1 wherein, 2 ..., ρ, ρ are pulse compression ratio.

The step that realizes heterogeneous code radar signal waveform automatic identifying method is as follows:

Step 1, the observation signal y (t) of any one section heterogeneous code radar signal=x (t)+ω (t) is sampled, obtain discrete form y (n)=x (n)+ω (n), sample frequency f _s, sampling time T _sThen calculate the discrete CWD conversion of signal y (n).

$W(t,\mathrm{\ω})=\∫\∫\frac{1}{\sqrt{4\mathrm{\π}{\mathrm{\τ}}^2/\mathrm{\σ}}}\exp (-\frac{{(\mathrm{\μ}-t)}^2}{4{\mathrm{\τ}}^2/\mathrm{\σ}})\·y(\mathrm{\μ}+\frac{\mathrm{\τ}}{2})\·y^{*}(\mathrm{\μ}-\frac{\mathrm{\τ}}{2})\exp (-\mathrm{j\ω\τ})\mathrm{d\μd\τ}---\left(8\right)$

Wherein, σ (σ＞0) is scale factor.

Step 2, for the impact that minimum signal bandwidth and sample frequency are brought signal CWD image, need to carry out normalized to the CWD image of signal.The step of normalized is as follows:

(1) the CWD image of signal carried out the threshold test processing;

(2) image after the threshold test processing is carried out time gated and frequency domain filtering, namely from the image border, remove the zone that does not contain signal;

(3) final bianry image depth-width ratio is normalized to 1.

The threshold value gating is processed vital effect for whole normalization algorithm.Should only comprise signal content and not include any independently noise spot through the view data after the processing of threshold value gating, because the result of second step is very responsive to these noise spots.Concerning the threshold value gating was processed, the output of choosing for the first step of threshold value played a crucial role.Adopt iterative algorithm that global threshold T is found the solution herein, specific algorithm is as follows:

(1) select the initial estimate of a T, this numerical value how obtains by the maximum grey level of CWD image and minimal gray level are averaging;

(2) according to selected T value the CWD image is divided into G ₁And G ₂Two parts, wherein G ₁Comprise the point that all grey levels are higher than T, G ₂Comprise the point that all grey levels are equal to or less than T;

(3) calculate respectively G ₁And G ₂Average intensity level μ in two parts ₁And μ ₂

(4) according to T=0.5 (μ ₁+ μ ₂) the new threshold value T of calculating;

(5) repeat the step that (2) arrive (4), until the numerical value of T reaches the requirement of convergence.

Therefore the numerical value of only setting global threshold can not guarantee removing fully of independent noise point, before and the frequency domain filtering time gated at second step, at first will process the CWD image, to reject harmful noise spot.Hazardous noise point is rejected and is finished by following steps:

(1) at first to bianry image is corroded and expansion process, bianry image is corroded and expansion process level and smooth CWD image not only, and can under Low SNR, remove because the flat or vertical spectral line of Choi-Williams nuclear formula (8);

(2) then the image of processing is carried out mark, with have in the bianry image obvious differentiation object carry out mark;

(3) reject at last in the target-marking target less than certain threshold value (as, reject all targets less than maximum target in the image 10%).If reject enough height of the Threshold of little target, also can be used as little target for the non-principal ingredient in P1, P2 and the P4 signal CWD image so and weed out.

In the second step of normalized, the image-region that does not contain signal content is removed from whole CWD image, the 3rd step was contained residue the image-region depth-width ratio normalization of signal content.Bianry image size after treatment is M * M, and wherein M is the minimum dimension of image after normalization process second step is processed.

Extract the signal characteristic with obvious differentiation in step 3, the signal CWD bianry image after normalized, be used for realizing the waveform recognition of dissimilar heterogeneous coded signal.CWD signal characteristic and the account form thereof extracted are as follows:

(1) Pseudo-Zernike square: the Pseudo-Zernike square has translation invariance, convergent-divergent unchangeability, rotational invariance and mirror invariant performance.

The p+q rank geometric moment of a digital picture f (x, y) is defined as follows:

$m_{\mathrm{pq}}=\underset{x}{\mathrm{\Σ}}\underset{y}{\mathrm{\Σ}}f(x,y)x^p{y}^q---\left(9\right)$

Translation invariant and convergent-divergent permanent center geometric moment are defined as follows:

$G_{\mathrm{pq}}=\frac{1}{m_{00}^{(p+q+2)/2}}\underset{x}{\mathrm{\Σ}}\underset{y}{\mathrm{\Σ}}f(x,y){(x-\stackrel{\‾}{x})}^p{(y-\stackrel{\‾}{y})}^q---\left(10\right)$

Wherein,

$\stackrel{\‾}{x}=\frac{m_{10}}{m_{00}},\stackrel{\‾}{y}=\frac{m_{01}}{m_{00}}$

The constant radially geometric moment of translation invariant and convergent-divergent is defined as follows:

$G_{\mathrm{pq}}=\frac{1}{m_{00}^{(p+q+3)/2}}\underset{x}{\mathrm{\Σ}}\underset{y}{\mathrm{\Σ}}f(x,y){({\tilde{x}}^2-{\tilde{y}}^2)}^{1/2}{\tilde{x}}^p{\tilde{y}}^q---\left(11\right)$

Wherein, $\tilde{x}=x-\stackrel{\‾}{x},\tilde{y}=y-\stackrel{\‾}{y}.$

The n rank Pseudo-Zernike square of m circulation can calculate by translation invariant and convergent-divergent permanent center geometric moment and translation invariant and the constant radially geometric moment of convergent-divergent, and specific algorithm is as follows:

$Z_{\mathrm{nm}}=\frac{n+1}{\mathrm{\π}}\underset{\underset{n-s-m=\mathrm{even}}{s=0}}{\overset{n-\left|m\right|}{\mathrm{\Σ}}}{D}_{\mathrm{nms}}$

$\underset{a=0}{\overset{k}{\mathrm{\Σ}}}\underset{b=0}{\overset{m}{\mathrm{\Σ}}}{(-j)}^b\left(\begin{array}{c}k\\ a\end{array}\right)\left(\begin{array}{c}m\\ b\end{array}\right)B_{\mathrm{nms}}{G}_{2k-2a+m-b,2a+b}---\left(12\right)$

$+\frac{n+1}{\mathrm{\π}}\underset{\underset{n-s-m=\mathrm{odd}}{s=0}}{\overset{n-\left|m\right|}{\mathrm{\Σ}}}{D}_{\mathrm{nms}}\underset{a=0}{\overset{d}{\mathrm{\Σ}}}\underset{b=0}{\overset{m}{\mathrm{\Σ}}}{(-j)}^b\left(\begin{array}{c}d\\ a\end{array}\right.)]\left(\begin{array}{c}m\\ b\end{array}\right)B_{\mathrm{nms}}{G}_{2d-2a+m-b,2a+b}$

Wherein,

$k=\frac{n-s-m}{2},$

$d=\frac{n-s-m-1}{2}$

$B_{\mathrm{nms}}={(-1)}^s\frac{(n-s)!}{s!(\frac{n+\left|m\right|}{2}-s)!(\frac{n-\left|m\right|}{2}-s)!}$

$D_{\mathrm{nms}}={(-1)}^s\frac{(2n+1-s)!}{s!(n-|m|-s)!(n+|m|+1-s)!}$

By to Pseudo-Zernike square Z _NmTaking absolute value obtains rotational invariance, and by it being taken the logarithm to carry out dynamic range compression.More than comprehensive, can draw final extraction and be characterized as:

${\hat{Z}}_{\mathrm{nm}}={\log }_e\left|Z_{\mathrm{nm}}\right|---\left(13\right)$

Wherein, choose the following Pseudo-Zernike square of bianry image after the normalized as the feature of waveform recognition:

With

(2) the target number in the bianry image after the normalized: on the basis that normalization is successfully finished, 2 signal target components are arranged in the two-value CWD image of Frank code and P3 coded signal, and only have 1 signal target component in the two-value CWD image of other 3 heterogeneous coded signals.In order to improve the robustness of feature, will Remove All less than the component of signal of maximum target component of signal 20%.

(3) peak power of the time location of peak power: P1, P2 and P4 coded signal is relatively near with the distance at coding center among the CWD, and Frank code and P3 code have the highest peak power at the coding end.This feature does not calculate from bianry image, so do not need the complete normalization to original CWD image, only just can realize by time gated processing.The method of calculating this feature in the CWD image is as follows:

$t_{\max }=\frac{1}{N-1}\underset{x}{\arg \max }\left\{W_{\mathrm{CW}}(x,y)\right\}---\left(14\right)$

Wherein, x represents time shaft, and y represents frequency axis, and N represents W _CW(x, y) length on time shaft.

Purpose be that the numerical value of this feature is normalized between 0～1.

(4) identification for the block structure of the CWD image of Frank code, P1 and P2 coded signal proposes P3 and P4 coded signal to be distinguished with other 3 kinds of coded signals by calculating the standard deviation of target component width in the bianry image, and circular is as follows:

After the signal object component in mark CWD image, process separately for each component of signal object, namely all component of signals outside the required component to be processed are removed at every turn.Major component among the bianry image B (x, y) is the proper vector of its covariance matrix, and its computing method are as follows:

$C=\underset{x=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{y=0}{\overset{N-1}{\mathrm{\Σ}}}(z-c){(z-c)}^{T}B(x,y)---\left(15\right)$

Wherein, the bianry image size is N * N,

With

, ordinate horizontal for the center of image can calculate by through type (11).

Bianry image is rotated, so that primary axis is parallel with the vertical or abscissa axis of image.Because the discrete coordinates characteristics of image, this rotary course need to carry out interpolation calculation, adopts the arest neighbors differential technique.The standard deviation of image object width can be calculated by postrotational binary image data.Suppose in the image that first primary axis corresponding with the highest major component of energy rotates to parallel with the longitudinal axis of image, then need the computed image row and value

$r\left(x\right)=\underset{y=0}{\overset{N-1}{\mathrm{\Σ}}}\hat{B}(x,y),x=\mathrm{0,1},...,N-1---\left(16\right)$

Wherein,

Represent postrotational bianry image.

Normalized r (x) is expressed as follows, and it is limited in 0～1 interval.

$\hat{r}\left(x\right)=\frac{r\left(x\right)}{\max r\left(x\right)}---\left(17\right)$

More than comprehensive, the standard deviation computing formula of target component width is as follows in the bianry image:

${\mathrm{\σ}}_{\mathrm{obj}}=\sqrt{\frac{1}{M}\underset{x}{\mathrm{\Σ}}{\hat{r}}^2\left(x\right)-{\left(\frac{1}{M}\underset{x}{\mathrm{\Σ}}\hat{r}\left(x\right)\right)}^2}---\left(18\right)$

Wherein, M represents

The quantity of non-small and weak sampling, the summation in the following formula is for non-small and weak sampling summation.Because small and weak sampling has a strong impact on the quality that standard deviation is estimated, especially when the row or column (depending on sense of rotation) of image when not containing any signal content, so before calculating, it will be got rid of.Set T _ObjValue be 0.3, be about to all

Small and weak sampling get rid of outside read group total.

The numerical value of final feature is exactly the standard deviation sigma of all component of signals in the image _ObjMean value.

(5) the 5th signal characteristic characteristic use the different coding symmetry characteristic of coded signal, the cross correlation function of the time energizing signal by compute sign rate sampling pulse and this pulse obtains.The peaked time delay size of this cross correlation function is a key character of the above-mentioned different coding signal of identification, and its computing method are as follows:

${\hat{r}}_y\left(\mathrm{\&tau;}\right)=\left\{\begin{array}{cc}\underset{n=0}{\overset{N-\mathrm{\&tau;}-1}{\mathrm{\&Sigma;}}}y(n+\mathrm{\&tau;})y^{*}(N-1-n),& \mathrm{\&tau;}\&GreaterEqual;0\\ \underset{n=-\mathrm{\&tau;}}{\overset{N-1}{\mathrm{\&Sigma;}}}y(n+\mathrm{\&tau;})y^{*}(N-1-n),& \mathrm{\&tau;}<0\end{array}\right.---\left(19\right)$

Wherein, y (n), n=0,1 ..., N-1 is the discrete time complex envelope signal of single symbol rate sampling, | τ |≤N-1.Then the peaked time delay size of final cross correlation function is expressed as:

${\mathrm{\&tau;}}_{\max }=\underset{\mathrm{\&tau;}}{\arg \max }\left({\hat{r}}_y\left(\mathrm{\&tau;}\right)\right)---\left(20\right)$

This feature have constant rotation (be n=0 among the y (n), 1 ..., the rotation amplitude is identical during N-1) and unchangeability is in order to obtain the feature identical with time-reversal signal, and can also be with | τ _Max| identify as feature.

Step 4, design a kind of combined classifier of neural network, in order to improve the recognition performance of neural network from accuracy of identification and efficient.Designed neural network classifier is as follows:

Employing is to the posterior probability weighting, by voting to determine modulation type.

If the classification number of signal to be sorted is K, the sorter number is N, and for input feature vector vector X, then the k of n sorter is output as

O _nk(X)＝P(c _k|X)+e _nk(X) (21)

Wherein, P (c _k| X) expression is judged as the posterior probability of k class, e when being input as X _Nk(X) output error of k node of n sorter of expression.Weight vector ω _k={ ω _1k, ω _2k..., ω _NkBe the output weights of k node of n sorter.Then each sorter judges that other output weighted sum of same class can be expressed as follows:

$S_k\left(X\right)=\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_{\mathrm{nk}}{O}_{\mathrm{nk}}\left(X\right)---\left(22\right)$

Wherein, k=1,2 ..., K.

Increase constraint condition $\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_{\mathrm{nk}}=1$ With $\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_{\mathrm{nk}}{e}_{\mathrm{nk}}\left(X\right)=1,$ Then have

S _k(X)＝P(c _k|X) (23)

When

The time, judge that k class signal exists.

Below in conjunction with example the present invention is done the emulation explanation:

When SNR be-during 2dB～31dB, each modulation signal interval 3dB is produced altogether 1000 feature samples, carry out Monte Carlo emulation emulation, and computation of mean values is as test result.Each training set and checking collection are that verification msg is carried out different cutting apart to set, and wherein the checking collection comprises the data of original training set 10%.

Fig. 2 points out that this sorter carried out the waveform recognition function reliably.Overall correct classification rate is higher than 97%, and the correct recognition rata that belongs to independent modulation type when signal to noise ratio (S/N ratio) is 3dB is higher than 91%.Yet, trickle obscuring still appear in the situation that signal to noise ratio (S/N ratio) is higher, and there is about 1%～2% P2 coded signal to be identified as the P1 coded signal by mistake.This obscures the error that occurs when coming from the estimate symbol rate.In addition, about 1% P1 coded signal is identified as the P4 coded signal by mistake.Table 1 has recorded the classification number percent when signal to noise ratio (S/N ratio) is 3dB.

Table 1 unlike signal identification crossing-over rate

Above-mentioned simulation result shows, the present invention is higher to the recognition correct rate of heterogeneous coded signal; Under Low SNR, heterogeneous code radar signal still had higher correct identification probability.

Claims (1)

1. based on the heterogeneous code radar signal waveform automatic identifying method of CWD feature, be divided into normalized, signal characteristic abstraction and 3 parts of signal classifier design of heterogeneous code radar signal discrete CWD image, amount to 4 treatment steps.

Suppose that radar signal that melodeon receives mixes and additive white Gaussian noise arranged (AWGN), and signal is as follows through processing the complex envelope that becomes the radar signal y (t) that baseband signal then receives:

y(t)＝x(t)+ω(t) (1)

Wherein x (t) is the complex envelope (including only a code-element period) of radar emission signal, and ω (t) is circulation additivity white complex gaussian noise.

The phase encoding complex signal of radar emission is expressed as follows:

$x\left(t\right)={\mathrm{Ae}}^{j(2{\mathrm{\&pi;f}}_{c}t+{\mathrm{\&phi;}}_i)}---\left(2\right)$

Wherein, A is signal amplitude, f _cSignal(-) carrier frequency, φ _iBe the signal discrete phase sequence, each phase place has the identical duration.The below provides the mathematical model of 5 kinds of heterogeneous code radar signals will studying.

The Frank coded signal is that a kind of stepping to the LFM signal approaches, and it has adopted N step frequency, and carries out N sampling at each Frequency point.Therefore, total hits of a Frank code is N ²The phase place of i sampling of j frequency of a Frank code is as follows:

${\mathrm{\&phi;}}_{i,j}=\frac{2\mathrm{\&pi;}}{N}(i-1)(j-1)---\left(3\right)$

Wherein, i=1,2 ..., N, j=1,2 ..., N.Then the pulse compression ratio of this Frank code is N ²

The P1 coded signal also is that the stepping to the LFM signal approaches, and it has adopted N step frequency, and carries out N sampling at each Frequency point.The phase place of i sampling of j frequency of a P1 code is as follows:

${\mathrm{\&phi;}}_{i,j}=-\frac{\mathrm{\&pi;}}{N}[N-(2j-1)][(j-1)N+(i-1)]---\left(4\right)$

Wherein, i=1,2 ..., N, j=1,2 ..., N.Then the pulse compression ratio of this P1 code is N ²

The P2 code has the characteristics of the palindrome, and namely the positive and negative of P2 code is identical.The phase place of i sampling of j frequency of a P2 code is as follows:

${\mathrm{\&phi;}}_{i,j}=[\frac{\mathrm{\&pi;}}{2}\&CenterDot;\frac{N-1}{N}\&CenterDot;\frac{\mathrm{\&pi;}}{N}(i-1)](N+1-2j)---\left(5\right)$

Wherein, i=1,2 ..., N, j=1,2 ..., N.Then the pulse compression ratio of this P2 code is N ²It is pointed out that N is necessary for even number in the P2 code, if N is radix, then the numerical value of the autocorrelation sidelobe of this P2 code can be too high.

The P3 code is from evolution that a LFM signal is sampled.The phase place of i sampling of a P3 code is as follows:

${\mathrm{\&phi;}}_i=\frac{\mathrm{\&pi;}}{\mathrm{\&rho;}}{(i-1)}^2---\left(6\right)$

I=1 wherein, 2 ..., ρ, ρ are pulse compression ratio.

The P4 code is from evolution that the identical signal of P3 code is sampled.The phase place of i sampling of a P4 code is as follows:

${\mathrm{\&phi;}}_i=\frac{\mathrm{\&pi;}}{\mathrm{\&rho;}}{(i-1)}^2-\mathrm{\&pi;}(i-1)---\left(7\right)$

I=1 wherein, 2 ..., ρ, ρ are pulse compression ratio.

The step that realizes heterogeneous code radar signal waveform automatic identifying method is as follows:

Step 1, the observation signal y (t) of any one section heterogeneous code radar signal=x (t)+ω (t) is sampled, obtain discrete form y (n)=x (n)+ω (n), sample frequency f _s, sampling time T _sThen calculate the discrete CWD conversion of signal y (n).

$W(t,\mathrm{\&omega;})=\&Integral;\&Integral;\frac{1}{\sqrt{4\mathrm{\&pi;}{\mathrm{\&tau;}}^2/\mathrm{\&sigma;}}}\exp (-\frac{{(\mathrm{\&mu;}-t)}^2}{4{\mathrm{\&tau;}}^2/\mathrm{\&sigma;}})\&CenterDot;y(\mathrm{\&mu;}+\frac{\mathrm{\&tau;}}{2})\&CenterDot;y^{*}(\mathrm{\&mu;}-\frac{\mathrm{\&tau;}}{2})\exp (-\mathrm{j\&omega;\&tau;})\mathrm{d\&mu;d\&tau;}---\left(8\right)$

Wherein, σ (σ＞0) is scale factor.

Step 2, for the impact that minimum signal bandwidth and sample frequency are brought signal CWD image, need to carry out normalized to the CWD image of signal.The step of normalized is as follows:

(1) the CWD image of signal carried out the threshold test processing;

(2) image after the threshold test processing is carried out time gated and frequency domain filtering, namely from the image border, remove the zone that does not contain signal;

(3) final bianry image depth-width ratio is normalized to 1.

The threshold value gating is processed vital effect for whole normalization algorithm.Should only comprise signal content and not include any independently noise spot through the view data after the processing of threshold value gating, because the result of second step is very responsive to these noise spots.Concerning the threshold value gating was processed, the output of choosing for the first step of threshold value played a crucial role.Adopt iterative algorithm that global threshold T is found the solution herein, specific algorithm is as follows:

(1) select the initial estimate of a T, this numerical value how obtains by the maximum grey level of CWD image and minimal gray level are averaging;

(2) according to selected T value the CWD image is divided into G ₁And G ₂Two parts, wherein G ₁Comprise the point that all grey levels are higher than T, G ₂Comprise the point that all grey levels are equal to or less than T;

(3) calculate respectively G ₁And G ₂Average intensity level μ in two parts ₁And μ ₂

(4) according to T=0.5 (μ ₁+ μ ₂) the new threshold value T of calculating;

(5) repeat the step that (2) arrive (4), until the numerical value of T reaches the requirement of convergence.

Therefore the numerical value of only setting global threshold can not guarantee removing fully of independent noise point, before and the frequency domain filtering time gated at second step, at first will process the CWD image, to reject harmful noise spot.Hazardous noise point is rejected and is finished by following steps:

(1) at first to bianry image is corroded and expansion process, bianry image is corroded and expansion process level and smooth CWD image not only, and can under Low SNR, remove because the flat or vertical spectral line of Choi-Williams nuclear formula (8);

(2) then the image of processing is carried out mark, with have in the bianry image obvious differentiation object carry out mark;

(3) reject at last in the target-marking target less than certain threshold value (as, reject all targets less than maximum target in the image 10%).If reject enough height of the Threshold of little target, also can be used as little target for the non-principal ingredient in P1, P2 and the P4 signal CWD image so and weed out.

In the second step of normalized, the image-region that does not contain signal content is removed from whole CWD image, the 3rd step was contained residue the image-region depth-width ratio normalization of signal content.Bianry image size after treatment is M * M, and wherein M is the minimum dimension of image after normalization process second step is processed.

Extract the signal characteristic with obvious differentiation in step 3, the signal CWD bianry image after normalized, be used for realizing the waveform recognition of dissimilar heterogeneous coded signal.CWD signal characteristic and the account form thereof extracted are as follows:

(1) Pseudo-Zernike square: the Pseudo-Zernike square has translation invariance, convergent-divergent unchangeability, rotational invariance and mirror invariant performance.

The p+q rank geometric moment of a digital picture f (x, y) is defined as follows:

$m_{\mathrm{pq}}=\underset{x}{\mathrm{\&Sigma;}}\underset{y}{\mathrm{\&Sigma;}}f(x,y)x^p{y}^q---\left(9\right)$

Translation invariant and convergent-divergent permanent center geometric moment are defined as follows:

$G_{\mathrm{pq}}=\frac{1}{m_{00}^{(p+q+2)/2}}\underset{x}{\mathrm{\&Sigma;}}\underset{y}{\mathrm{\&Sigma;}}f(x,y){(x-\stackrel{\&OverBar;}{x})}^p{(y-\stackrel{\&OverBar;}{y})}^q---\left(10\right)$

Wherein,

$\stackrel{\&OverBar;}{x}=\frac{m_{10}}{m_{00}},$ $\stackrel{\&OverBar;}{y}=\frac{m_{01}}{m_{00}}$

The constant radially geometric moment of translation invariant and convergent-divergent is defined as follows:

$G_{\mathrm{pq}}=\frac{1}{m_{00}^{(p+q+3)/2}}\underset{x}{\mathrm{\&Sigma;}}\underset{y}{\mathrm{\&Sigma;}}f(x,y){({\tilde{x}}^2-{\tilde{y}}^2)}^{1/2}{\tilde{x}}^p{\tilde{y}}^q---\left(11\right)$

Wherein, $\tilde{x}=x-\stackrel{\&OverBar;}{x},$ $\tilde{y}=y-\stackrel{\&OverBar;}{y}.$

The n rank Pseudo-Zernike square of m circulation can calculate by translation invariant and convergent-divergent permanent center geometric moment and translation invariant and the constant radially geometric moment of convergent-divergent, and specific algorithm is as follows:

$Z_{\mathrm{nm}}=\frac{n+1}{\mathrm{\&pi;}}\underset{\underset{n-s-m=\mathrm{even}}{s=0}}{\overset{n-\left|m\right|}{\mathrm{\&Sigma;}}}{D}_{\mathrm{nms}}$

$\underset{a=0}{\overset{k}{\mathrm{\&Sigma;}}}\underset{b=0}{\overset{m}{\mathrm{\&Sigma;}}}{(-j)}^b\left(\begin{array}{c}k\\ a\end{array}\right)\left(\begin{array}{c}m\\ b\end{array}\right)B_{\mathrm{nms}}{G}_{2k-2a+m-b,2a+b}---\left(12\right)$

$+\frac{n+1}{\mathrm{\&pi;}}\underset{\underset{n-s-m=\mathrm{odd}}{s=0}}{\overset{n-\left|m\right|}{\mathrm{\&Sigma;}}}{D}_{\mathrm{nms}}\underset{a=0}{\overset{d}{\mathrm{\&Sigma;}}}\underset{b=0}{\overset{m}{\mathrm{\&Sigma;}}}{(-j)}^b\left(\begin{array}{c}d\\ a\end{array}\right)\left(\begin{array}{c}m\\ b\end{array}\right)B_{\mathrm{nms}}{G}_{2d-2a+m-b,2a+b}$

Wherein,

$k=\frac{n-s-m}{2},$

$d=\frac{n-s-m-1}{2}$

$B_{\mathrm{nms}}={(-1)}^s\frac{(n-s)!}{s!(\frac{n+\left|m\right|}{2}-s)!(\frac{n-\left|m\right|}{2}-s)!}$

$D_{\mathrm{nms}}={(-1)}^s\frac{(2n+1-s)!}{s!(n-\left|m\right|-s)!(n+\left|m\right|+1-s)!}$

By to Pseudo-Zernike square Z _NmTaking absolute value obtains rotational invariance, and by it being taken the logarithm to carry out dynamic range compression.More than comprehensive, can draw final extraction and be characterized as:

Wherein, choose the following Pseudo-Zernike square of bianry image after the normalized as the feature of waveform recognition:

With

(2) the target number in the bianry image after the normalized: on the basis that normalization is successfully finished, 2 signal target components are arranged in the two-value CWD image of Frank code and P3 coded signal, and only have 1 signal target component in the two-value CWD image of other 3 heterogeneous coded signals.In order to improve the robustness of feature, will Remove All less than the component of signal of maximum target component of signal 20%.

(3) peak power of the time location of peak power: P1, P2 and P4 coded signal is relatively near with the distance at coding center among the CWD, and Frank code and P3 code have the highest peak power at the coding end.This feature does not calculate from bianry image, so do not need the complete normalization to original CWD image, only just can realize by time gated processing.The method of calculating this feature in the CWD image is as follows:

$t_{\max }=\frac{1}{N-1}\underset{x}{\arg \max }\left\{W_{\mathrm{CW}}(x,y)\right\}---\left(14\right)$

Wherein, x represents time shaft, and y represents frequency axis, and N represents W _CW(x, y) length on time shaft. Purpose be that the numerical value of this feature is normalized between 0～1.

(4) identification for the block structure of the CWD image of Frank code, P1 and P2 coded signal proposes P3 and P4 coded signal to be distinguished with other 3 kinds of coded signals by calculating the standard deviation of target component width in the bianry image, and circular is as follows:

After the signal object component in mark CWD image, process separately for each component of signal object, namely all component of signals outside the required component to be processed are removed at every turn.Major component among the bianry image B (x, y) is the proper vector of its covariance matrix, and its computing method are as follows:

$C=\underset{x=0}{\overset{N-1}{\mathrm{\&Sigma;}}}\underset{y=0}{\overset{N-1}{\mathrm{\&Sigma;}}}(z-c){(z-c)}^{T}B(x,y)---\left(15\right)$

Wherein, the bianry image size is N * N, z=(x, y) ^T,

With

, ordinate horizontal for the center of image can calculate by through type (11).

Bianry image is rotated, so that primary axis is parallel with the vertical or abscissa axis of image.Because the discrete coordinates characteristics of image, this rotary course need to carry out interpolation calculation, adopts the arest neighbors differential technique.The standard deviation of image object width can be calculated by postrotational binary image data.Suppose in the image that first primary axis corresponding with the highest major component of energy rotates to parallel with the longitudinal axis of image, then need the computed image row and value

Wherein,

Represent postrotational bianry image.

Normalized r (x) is expressed as follows, and it is limited in 0～1 interval.

$\hat{r}\left(x\right)=\frac{r\left(x\right)}{\max r\left(x\right)}---\left(17\right)$

More than comprehensive, the standard deviation computing formula of target component width is as follows in the bianry image:

${\mathrm{\&sigma;}}_{\mathrm{obj}}=\sqrt{\frac{1}{M}\underset{x}{\mathrm{\&Sigma;}}{\hat{r}}^2\left(x\right)-{\left(\frac{1}{M}\underset{x}{\mathrm{\&Sigma;}}\hat{r}\left(x\right)\right)}^2}---\left(18\right)$

Wherein, M represents

The quantity of non-small and weak sampling, the summation in the following formula is for non-small and weak sampling summation.Because small and weak sampling has a strong impact on the quality that standard deviation is estimated, especially when the row or column (depending on sense of rotation) of image when not containing any signal content, so before calculating, it will be got rid of.Set T _ObjValue be 0.3, be about to all

Small and weak sampling get rid of outside read group total.

The numerical value of final feature is exactly the standard deviation sigma of all component of signals in the image _ObjMean value.

(5) the 5th signal characteristic characteristic use the different coding symmetry characteristic of coded signal, the cross correlation function of the time energizing signal by compute sign rate sampling pulse and this pulse obtains.The peaked time delay size of this cross correlation function is a key character of the above-mentioned different coding signal of identification, and its computing method are as follows:

${\hat{r}}_y\left(\mathrm{\&tau;}\right)=\left\{\begin{array}{cc}\underset{n=0}{\overset{N-\mathrm{\&tau;}-1}{\mathrm{\&Sigma;}}}y(n+\mathrm{\&tau;})y^{*}(N-1-n),& \mathrm{\&tau;}\&GreaterEqual;0\\ \underset{n=-\mathrm{\&tau;}}{\overset{N-1}{\mathrm{\&Sigma;}}}y(n+\mathrm{\&tau;})y^{*}(N-1-n),& \mathrm{\&tau;}<0\end{array}\right.---\left(19\right)$

Wherein, y (n), n=0,1 ..., N-1 is the discrete time complex envelope signal of single symbol rate sampling, | τ |≤N-1.Then the peaked time delay size of final cross correlation function is expressed as:

${\mathrm{\&tau;}}_{\max }=\underset{\mathrm{\&tau;}}{\arg \max }\left({\hat{r}}_y\left(\mathrm{\&tau;}\right)\right)---\left(20\right)$

This feature have constant rotation (be n=0 among the y (n), 1 ..., the rotation amplitude is identical during N-1) and unchangeability is in order to obtain the feature identical with time-reversal signal, and can also be with | τ _Max| identify as feature.

Step 4, design a kind of combined classifier of neural network, in order to improve the recognition performance of neural network from accuracy of identification and efficient.Designed neural network classifier is as follows:

Employing is to the posterior probability weighting, by voting to determine modulation type.

If the classification number of signal to be sorted is K, the sorter number is N, and for input feature vector vector X, then the k of n sorter is output as

O _nk(X)＝P(c _k|X)+e _nk(X) (21)

Wherein, P (c _k| X) expression is judged as the posterior probability of k class, e when being input as X _Nk(X) output error of k node of n sorter of expression.Weight vector ω _k={ ω _1k, ω _2k..., ω _NkBe the output weights of k node of n sorter.Then each sorter judges that other output weighted sum of same class can be expressed as follows:

$S_k\left(X\right)=\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_{\mathrm{nk}}{O}_{\mathrm{nk}}\left(X\right)---\left(22\right)$

Wherein, k=1,2 ..., K.

Increase constraint condition $\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_{\mathrm{nk}}=1$ With $\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_{\mathrm{nk}}{e}_{\mathrm{nk}}\left(X\right)=1,$ Then have

S _k(X)＝P(c _k|X) (23)

When

The time, judge that k class signal exists.

Images