CN100521670C - Detecting and analyzing method for multi system frequency shift key control signal - Google Patents

Abstract

This invention relates to automatic identification to communication signals and analysis technology aiming at testing and analyzing multi-bit frequency shift key controlled signals characterizing in comparing the mean value of the absolute value of the variance rate sequence of the normalized instant frequencies with the set threshold value to judge MFSK signal, when analyzing, Smooth filtration is used in reducing shake of the MFSK signal transient frequency then the entropy principle is used in analyzing the histogram for the distribution of transient frequencies to estimate the state number of code elements and frequency sphere with an equal probability distribution principle then to use the equal frequency distance to get more accurate code element state number, frequency sphere and central frequency.

Description

A kind of detection of multi system frequency shift key control signal and analytical method

Technical field

The invention belongs to the signal of communication analysis field, the particularly automatic identification and the analysis of multi system frequency shift key control (MFSK) signal in the digital modulation mode.

Background technology

Modulation system identification automatically is the important step that electronic warfare (signal reconnaissance), radio monitoring, adaptive communications etc. are used, and is an important topic in non-collaboration communication field, and the various theories of signal of communication analysis have obtained extensive use in this direction.Frequency shift keying modulation system (FSK) is being played the part of important role in the shortwave/ultra short wave communication of army, government and amateur radio communication, being a kind of important digital signal modulation mode, is the important content of radio signal monitoring and military communication scouting to the detection of such signal and analysis.The present invention relates generally to the detection and the analysis of frequency shift keyed signals, the i.e. analysis of the identification of frequency shift keyed signals and order of modulation.

In recent years, at multi system frequency shift key control (MFSK) input and analysis, following several method has been proposed:

1. the method for analyzing based on signal envelope analysis and instantaneous frequency

People such as E.E.Azzouz in Automatic modulation recognition (Journal Of The FranklinInstitute-Engineering And Applied Mathematics, Mar 1997), proposed to use signal envelope spectrum peak γ _MaxPermanent envelope (frequency modulated signal) signal and non-constant envelope signal are classified, use the standard deviation of normalization instantaneous frequency absolute value _AfTo BFSK signal and the classification of QFSK signal, step is as follows:

1). the signal that receives is converted into analytic signal, removes the weak signal section;

2). calculate zero center normalization instantaneous amplitude acn;

3). the spectrum peak γ of signal calculated envelope _Max, should be worth with predefine thresholding τ _γCompare, if γ _Max＜τ _γThen received signal is judged to be frequency modulated signal;

${\mathrm{\γ}}_{\max }=\frac{\max \left(\mathrm{DFT}\left(a_{\mathrm{cn}}\left(n\right)\right)\right)}{N_s}$

4). calculate its zero center normalization instantaneous frequency f by the instantaneous phase of analytic signal _Cn

5). calculate the standard deviation of normalization instantaneous frequency absolute value _Af, should be worth with the predefine thresholding Compare, if σ _Af＜τ _σThen be judged to be the BFSK signal, otherwise be judged to be the QPSK signal.

${\mathrm{\σ}}_{\mathrm{af}}=\sqrt{\frac{1}{N_s}\mathrm{\Σ}{f}_{\mathrm{cn}}^2\left(n\right)-{\left[\frac{1}{N_s}\mathrm{\Σ}{f}_{\mathrm{cn}}\right]}^2}$

The problem of this method is:

● because factor affecting such as frequency selective fadings, the MFSK sample of signal of many actual reception is not a constant envelope signal;

● this method can only be discerned BFSK and QFSK signal.

Cao Zhi has just waited the people " to need not the digital signal automatic Modulation Recognition method commonly used of priori " and also related in (patent No. 02123627.5) detection and the analysis of MFSK signal in relevant patent.Relevant step is as follows:

1). the R parameter of signal calculated sample, if this parameter less than the predefine thresholding, then is classified as constant envelope signal with this signal;

2). to constant envelope signal, the instantaneous phase of signal calculated s (n)=x (n)+iy (n):

3). go volume folded to instantaneous phase, add the phasing sequence for the phase sequence of mould 2 π:

4). the instantaneous frequency f of signal calculated ₁(n):

$f_1\left(n\right)=f_{11}\left(n\right)-\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{f}_{11}\left(n\right),f_{11}\left(n\right)=\mathrm{\θ}\left(n\right)-\mathrm{\θ}(n-1);$

5). to f ₁(n) carry out the statistics of absolute value, obtain statistics number C greater than pi/2 _p, satisfy C _p＜T _CpSignal determining be the signal that does not have phase hit, i.e. MFSK signal;

6). to f ₁(n) carry out interpolation processing, obtain f ₂(n):

$f_2\left(n\right)=\left\{\begin{array}{cc}{f}_1\left(n\right)& ,\left|f_1\left(n\right)\right|<\frac{\mathrm{\&pi;}}{2}\\ \frac{1}{2}[f_1(n+1)+f_1(n-1)]& ,\left|f_1\left(n\right)\right|\&GreaterEqual;\frac{\mathrm{\&pi;}}{2}\end{array}\right.;$

7). to f ₂(n) carry out normalized and obtain normalized frequency f (n):

$f\left(n\right)=\frac{f_2\left(n\right)}{\max \left\{f_2\left(n\right)\right\}};$

8). the distribution to f (n) is added up, and obtains the peak value number N of f (n) _f, if N _f=2, then be judged to be the BFSK signal, otherwise be judged to be the FM signal.

The problem of this method is:

● because factor affecting such as frequency selective fadings, the MFSK sample of signal of many actual reception is not a constant envelope signal;

● this patent does not provide the peak value number N of the Distribution Statistics of calculating f (n) _fMethod, and can only discern the BFSK signal.

2. based on the method for zero passage sampling

People such as S.-Z.Hsue are at paper Automatic Modulation Classification Using Zero Crossing (Radar andSignal Processing, IEE Proceedings For, vol.37, No.6, pp.459-464, Dec.1990.K.) proposed in to use the zero passage sampler to extract the inverse of the instantaneous frequency of signal, detect the MFSK signal by analyzing this variance reciprocal then, analyze the method for MFSK signal by analyzing the instantaneous frequency distribution histogram, step is as follows:

1). received signal, the time point sequence x (i) of the zero passage of usefulness zero passage sampler tracer signal;

2). by time x (i) calculated value sequences y (i)=x (the i+1)-x (i) and value sequence z (i)=y (the i+1)-y (i) of signal zero passage;

3). set thresholding τ _z, will satisfy z (i)〉and τ _zSample point be defined as " intersymbol saltus step (IST) sample ", remove corresponding y (i), obtain revised value sequence y _a(i);

4). calculate y _a(i) variance G, and set thresholding T satisfies G〉sample of signal of T is judged as the MFSK signal;

5). with the y between the IST sample point _a(i) get average, calculate its distribution histogram;

6). calculate the order of modulation of MFSK signal by analyzing this histogrammic peak value number.

The problem of this method is:

● the accuracy of zero passage Frequency Estimation is very responsive to SNR, and in relevant report, the author has only provided the simulation result of CNR 〉=15dB;

● zero passage method is measured instantaneous frequency needs very high sample rate;

● under the low signal-to-noise ratio condition, calculate accurately relatively difficulty of order of modulation by the peak value number of analyzing the resulting frequency distribution block diagram of zero passage Frequency Estimation.

3. based on the method for signal discrete Fourier transform

People such as Zaihe Yu are at paper A practical classification algorithm for M-ary frequency shift keyingsignals (Military Communications Conference, 2004.MILCOM 2004.IEEE Volume 2,31 Oct.-3Nov.2004 Page (s): proposed a kind of method that obtains MFSK signal order of modulation by direct analysis MFSK power spectrum signal 1123-1128 Vol.2).Do not have the detection method of MFSK signal in this report, but suppose that signal to be analyzed is the MFSK signal.This method step is as follows:

1). the discrete Fourier transform (DFT) of signal calculated sample (DFT);

2). according to thresholding rule, local extremum rule, equidistant rule, perfect information rule, rule of classification, adaptive threshold rule and confidence level rule analysis sample of signal DFT peak value number, obtain order of modulation M.

The problem of this method is:

● this method supposes that the distribution of each code element on power equates, still, in the test environment of reality, because factors such as frequency selective fadings, each symbol power distributes and often do not equate;

● when calculating DFT, the author states the FFT length that use is long as far as possible, and this is unfavorable for the realization of real system, and under the signal bandwidth condition of unknown, estimates that a rational FFT length also is very difficult.

Summary of the invention

Main purpose of the present invention is identification and analyzes multi system frequency shift key control (MFSK) signal, comprises a kind of detection method of MFSK signal and a kind of analytical method of MFSK signal.

The invention is characterized in that described method realizes successively according to the following steps in computer:

Step (1) receives pending data, obtains real signal x (n);

Step (2) converts real signal x (n) to analytic signal s (n) according to the following steps:

Step (2.1) is calculated FFT after real signal x (n) zero padding, obtains the discrete Fourier transform (DFT) F of signal _x(n), the length N of zero padding _zObtain by following formula:

N _z=N _FFT-N _x, N _xBe the length of real signal x (n),

Expression is greater than the smallest positive integral of " ", N _FFTBe the length of carrying out the FFT computing;

Step (2.2) is got described F _x(n) first half, n=1 ..., N _FFT/ 2, obtain the discrete Fourier transform (DFT) F of signal s (n) _s(n);

Step (2.3) is F _s(n) first half, n=1 ..., N _FFT/ 4 and latter half n=N _FFT/ 4+1 ..., N _FFT/ 2 transposings;

The F that step (2.4) calculation procedure (2.3) obtains _s(n) IFFT obtains analytic signal s (n), s (n)=I _FFT(F _s(n));

The analytic signal s's (n) that step (2.5) obtains step (2.4)

Part is clipped, and remainder is the analytic signal that subsequent analysis is used, and its effective length is

Step (3) is according to the predefine thresholding, and segmentation detects the amplitude of s (n), removes the weak signal section, and its steps in sequence is as follows:

The signal s (n) that step (3.1) obtains step (2.5)=x (n)+jy (n) is divided into isometric N _SegSection is with a sequence S _i, i=1,2 ..., N _SegExpression, N _Seg=5～20;

Each the segment signal amplitude that obtains in step (3.2) calculation procedure (3.1) and m _i:

$m_i=\underset{s\left(n\right)\&Element;S_i}{\mathrm{\&Sigma;}}\left|s\left(n\right)\right|;$

Step (3.3) is set thresholding τ _m:

${\mathrm{\&tau;}}_m=\frac{1}{2}\max \{m_i,i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;,N_{\mathrm{seg}}\};$

Step (3.4) is with m _i＜τ _mSignal segment be judged to be the weak signal section, remove all weak signal sections, remaining signal segment is connected into new analytic signal sequence s (n), length is N;

The normalization instantaneous frequency sequence f of step (4) calculation procedure (3.4) described analytic signal sequence s (n) _n(n), realize according to the following steps:

Step (4.1) is calculated as follows analytic signal sequence s (n)=x (n)+jy (n), n=1, and 2 ..., the instantaneous phase sequence of N

The instantaneous phase sequence that step (4.2) obtains step (4.1)

Do and get the surplus instantaneous frequency sequence f (n) that obtains signal behind the calculus of differences:

Step (4.3) is normalized to zero-mean, unit variance sequence f to this instantaneous frequency sequence f (n) _n(n):

$f_n\left(n\right)=\frac{f\left(n\right)-m_f}{{\mathrm{\&sigma;}}_f},$

Wherein, $m_f=\frac{1}{N-1}\underset{n=1}{\overset{N-1}{\mathrm{\&Sigma;}}}f\left(n\right),$ ${\mathrm{\&sigma;}}_f=\sqrt{\frac{1}{N-2}\underset{n=1}{\overset{N-1}{\mathrm{\&Sigma;}}}{(f\left(n\right)-m_f)}^2};$

The f that step (5) calculation procedure (4.3) obtains _nThe rate of change sequence Δ f of normalization instantaneous frequency (n) _n(n) and the average of absolute value

Wherein:

Δf _n(n)＝|f _n(n+2)-f _n(n)|，n＝1，2，……，N-3， $m_{\mathrm{\&Delta;}{f}_n}=\frac{1}{N-3}\underset{n=1}{\overset{N-3}{\mathrm{\&Sigma;}}}\mathrm{\&Delta;}{f}_n\left(n\right)$

Step (6) is described

With predefine thresholding τ _mRelatively: if $m_{\mathrm{\&Delta;}{f}_n}<{\mathrm{\&tau;}}_m,$ Then described signal x (n) is judged to be the MFSK signal, and execution in step (7), otherwise, judge that receiving sequence is not the MFSK signal, flow process finishes;

Step (7) is carried out following steps successively and is reduced the shake of instantaneous frequency to utilize mean filter:

The instantaneous frequency sequence f that step (7.1) utilizes mean filter that step (4.3) is obtained _n(n) carry out smothing filtering, the instantaneous frequency sequence f that obtains upgrading _n(n), filter length l _fBy the possible maximum symbol rate fd of analyzed signal _MaxSetting sample frequency f with signal _sObtain by following formula:

$l_f=m\frac{f_s}{f_{d\max }},m=\frac{1}{3}~\frac{1}{2};$

The instantaneous frequency sequence of the renewal that step (7.2) obtains according to step (7.1) recomputates the normalization instantaneous frequency rate of change sequence Δ f of signal _n(n): Δ f _n(n)=| f _n(n+2)-f _n(n) |;

The variance of the normalization instantaneous frequency rate of change sequence that step (7.3) obtains according to following formula calculation procedure (7.2) _{τ Δ fn}:

${\mathrm{\&tau;}}_{\mathrm{\&Delta;}{f}_n}=\sqrt{\frac{1}{N-4}\underset{n=1}{\overset{N-3}{\mathrm{\&Sigma;}}}{(\mathrm{\&Delta;}{f}_n\left(n\right)-m_{\mathrm{\&Delta;}{f}_n})}^2},$ $m_{\mathrm{\&Delta;}{f}_n}=\frac{1}{N-3}\underset{n=1}{\overset{N-3}{\mathrm{\&Sigma;}}}\mathrm{\&Delta;}{f}_n\left(n\right);$

Step (7.4) search Δ f _n(n) sequence, if $\left|\mathrm{\&Delta;}{f}_n\left(n\right)\right|<{\mathrm{\&tau;}}_{\mathrm{\&Delta;}{f}_n},$ F then _n(n+1) belong to symbol stable region Γ _Si, i=1,2 ..., N _s, N wherein _sBe the number of symbol stable region, if $\left|\mathrm{\&Delta;}{f}_n\left(n\right)\right|>{\mathrm{\&tau;}}_{\mathrm{\&Delta;}{f}_n}$ F then _n(n+1) belong to the interval Γ of symbol saltus step _Ti, i=1,2 ..., N _t, N wherein _tNumber for symbol saltus step interval;

Step (7.5) changes the instantaneous frequency value of each symbol stable region in the step (7.4) average of all instantaneous frequency values of this interval into, gives up the instantaneous frequency value in all symbol saltus step intervals, constitutes new instantaneous frequency sequence f _n(n);

Step (8) basis " entropy principle " is the histogram of the described instantaneous frequency sequential value of analytical procedure (7.5) according to the following steps:

Step (8.1) is according to normalization instantaneous frequency sequence f _n(n) set up distribution histogram

With this as f _nThe estimation of probability density function (n):

$\hat{p}\left(f_i\right)=\frac{\mathrm{count}\{f_n\left(n\right),f_i\&le;f_n\left(n\right)<f_i+\mathrm{\&Delta;f}\}}{N_{f_n}},$ $f_i=\min \left\{f_n\left(n\right)\right\}+\mathrm{i\&Delta;f},i=\mathrm{0,1},\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;,N_{\hat{p}}-1,$

Wherein,

Be the length of f (n), $\mathrm{\&Delta;f}=(\max \left\{f_n\left(n\right)\right\}-\min \left\{f_n\left(n\right)\right\})/N_{\hat{p}},$

Length for distribution histogram $N_{\hat{p}}=512~4096,$ Count{x, y} represent the to satisfy condition number of x value of y

Step (8.2) is determined thresholding τ _pThe hunting zone be $0~\max \left\{\hat{p}\left(f_i\right)\right\}/2,$ Step-size in search $\mathrm{\&Delta;}{\mathrm{\&tau;}}_p=\max \left\{\hat{p}\left(f_i\right)\right\}/\left(2K\right),$ K＝10～100；

Step (8.3) is to each τ in the hunting zone _pValue is carried out following steps:

Step (8.3.1) search is satisfied $\hat{p}\left(f_i\right)>{\mathrm{\&tau;}}_p$ F _i, continuous f _iBe classified as one group;

Step (8.3.2) is calculated the probability density function of each group and is estimated

Sum obtains the estimated value of the distribution probability of each code element of MFSK signal

Step (8.3.3) is calculated should τ _pEntropy e:

$e=-\underset{k=1}{\overset{\hat{M}}{\mathrm{\&Sigma;}}}{\hat{P}}_k\log \left({\hat{P}}_k\right),$

Be continuous f _iThe number of the group that value is formed;

Step (8.4) is got the τ of corresponding entropy maximum _p, with corresponding f _iThe value of the number of group

As the initial estimate of code element state number M, every group minimum instantaneous frequency f _MinkWith the maximum instantaneous frequency f _MaxkFrequency range (f as each code element state _Mink, f _Maxk), k=1,2 ...,

Step (9) basis " equiprobability distribution principle ", the exact value of accurate estimating code element state number M according to the following steps:

Step (9.1) is set thresholding $0<{\mathrm{\&tau;}}_{P_{\min }}<1;$

Step (9.2) is removed ${\hat{P}}_k<{\mathrm{\&tau;}}_{P_{\min }}\&times;\max \{{\hat{P}}_k,k=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;,\hat{M}\}$ State, herein, Be a span at 0～1 thresholding, $\max \{{\hat{P}}_k,k=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;,\hat{M}\}$ Be

In maximum, the two multiplies each other and obtains a value that is used to screen rational code element state;

Step (9.3) recomputates the estimated value of code element state number M Frequency range with each code element state

$(f_{\min k},f_{\max k}),k=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;,\hat{M};$

Step (10) is according to accurately estimating code element state number M and each state centre frequency such as " frequency difference principle ", and its step is as follows:

Step (10.1) is calculated the centre frequency of each code element state

${\hat{f}}_{\mathrm{mk}}=\frac{\underset{i=i_{\min k}}{\overset{i_{\max k}}{\mathrm{\&Sigma;}}}\hat{p}\left(f_i\right)f_i}{\underset{i=i_{\min k}}{\overset{i_{\max k}}{\mathrm{\&Sigma;}}}\hat{p}\left(f_i\right)}$

Wherein, i _MinkAnd i _MaxkBe the minimum of every set of symbols state, the label of maximum instantaneous frequency,

It is the estimation of the probability density function of the instantaneous frequency sequence that calculates of step (8.1);

Step (10.2) is calculated the frequency difference between each state ${\hat{f}}_{\mathrm{dk},l}:{\hat{f}}_{\mathrm{dk},l}=|{\hat{f}}_{\mathrm{mk}}-{\hat{f}}_{\mathrm{ml}}|;$

Step (10.3) is calculated the weight w of each frequency difference _Kl:

Weight equals to obtain the smaller value that the probability distribution of two code element states of this frequency difference is estimated, the probability distribution that deducts all the code element states between these two code element states is estimated sum, if difference less than 0, then weight is 0:

$w_{\mathrm{kl}}=\left\{\begin{array}{cc}\max \{0,\min \{{\hat{P}}_k,{\hat{P}}_l\}-\underset{n=k+1}{\overset{l-1}{\mathrm{\&Sigma;}}}{\hat{P}}_n\}& ,k<l-1\\ \max \{0,\min \{{\hat{P}}_k,{\hat{P}}_l\}\}& ,k=l-1\end{array}\right.;$

Step (10.4) obtains the estimation of frequency difference to each frequency difference weighted average:

${\hat{f}}_d=\frac{\underset{k,l}{\mathrm{\&Sigma;}}{w}_{\mathrm{kl}}{\hat{f}}_{\mathrm{dk},l}}{\underset{k,l}{\mathrm{\&Sigma;}}{w}_{\mathrm{kl}}}$

Step (10.5) is set thresholding ${\mathrm{\&tau;}}_{\min },{\mathrm{\&tau;}}_{\max }:{\mathrm{\&tau;}}_{\min }=0.4{\hat{f}}_d~0.8{\hat{f}}_d,$ ${\mathrm{\&tau;}}_{\max }=1.2{\hat{f}}_d~1.6{\hat{f}}_d;$

Step (10.6) for ${\hat{f}}_{\mathrm{dk},k+1}<{\mathrm{\&tau;}}_{\min }$ Adjacent states, merge code element state k and k+1, wherein,

Adjacent code element state k that obtains in the expression step (10.2) and the frequency difference between the k+1;

Step (10.7) for ${\hat{f}}_{\mathrm{dk},k+1}>{\mathrm{\&tau;}}_{\max }$ Adjacent states, between code element state k and k+1, insert

Individual code element state, wherein,

Expression is greater than the smallest positive integral of " ",

Adjacent code element state k that obtains in the expression step (10.2) and the frequency difference between the k+1;

Step (10.8) is upgraded the estimation of code element state number

Centre frequency with each code element state

Step (11) subsequent treatment:

Be calculated as follows the final estimation that obtains MFSK signal order of modulation

${\hat{M}}_0=2^{\left[{\log }_2^{\hat{M}}\right]},$

Wherein [] computing is represented to get from " " nearest integer.

The present invention has following advantage:

● do not rely on the hypothesized model that the MFSK signal is a constant envelope signal, allow the decline of signal frequency of occurrences selectivity, the scope of application is greater than other existing algorithms;

● caught the MFSK signal to have the substantive characteristics of instantaneous frequency saltus step, discerned and analyzed the MFSK signal as essential characteristic, had robustness preferably with instantaneous frequency;

● at the test of computer mould analog signal with at the test shows of actual signal, the algorithm that the present invention proposes has good recognition effect and stronger practicality.

The effect that the present invention reaches:

Multi system frequency shift key control (MFSK) is the important digital signal modulation mode of a class, and the present invention relates generally to the detection and the analysis of the signal of communication of such modulation system.Computer simulation experiment shows, to order of modulation is 2～32 MFSK signal, and under the condition of signal to noise ratio greater than 5dB, the detection accuracy of institute of the present invention extracting method has reached more than 99%, under the condition of signal to noise ratio greater than 7.5dB, analyze accuracy and reached more than 95%.To the test shows of the sample of signal of actual reception, the detection accuracy of institute of the present invention extracting method has reached more than 95%, analyzes accuracy and has reached more than 90%.

Description of drawings

MFSK input and analysis process figure that Fig. 1 proposes for the present invention.

Fig. 2 is the hardware elementary diagram that adopts when actual signal is tested among the embodiment.

Fig. 3 is the identification test result of emulation MFSK signal, and wherein dotted line is represented fsk signal

Value, solid line are represented non-fsk signal

Value.

Fig. 4 is the analytical test result of emulation MFSK signal, and every curve is represented the analysis accuracy of a kind of modulation system under different signal to noise ratio conditions.

Fig. 5 for the identification test result of actual reception MFSK signal wherein dotted line represent fsk signal

Value, solid line are represented non-fsk signal

Value.

Embodiment

MFSK input proposed by the invention and analytical method overall procedure may further comprise the steps as shown in Figure 1:

1) receives pending data, obtain real signal x (n);

2) real signal x (n) is become analytic signal s (n);

3) according to the predefine thresholding, segmentation detects s (n) amplitude, removes the weak signal section;

4) the normalization instantaneous frequency sequence f of calculating s (n) _n(n);

5) the rate of change sequence Δ f of calculating normalization instantaneous frequency _n(n) and the average of absolute value

6) will

With predefine thresholding τ _mCompare, if $m_{\mathrm{\&Delta;}{f}_n}<{\mathrm{\&tau;}}_m,$ Be the MFSK signal then, otherwise judge that this signal is not the MFSK signal this signal determining;

7) utilize mean filter to reduce the shake of instantaneous frequency;

8) analyze the instantaneous frequency distribution histogram according to " entropy principle ", obtain code element state number

Initial estimation;

9) according to " equiprobability distribution principle " accurate estimating code element state number

10) according to accurate estimating code element state number such as " frequency difference principle "

With each state centre frequency;

11) subsequent treatment.

Here provide at the test of Computer Simulation signal with at two embodiment of test of actual reception signal.

1. at the test of Computer Simulation signal

A) identification of MFSK signal

Under different signal to noise ratio conditions, generate two groups of signals and be used for test.One group is fsk signal, and one group is non-fsk signal:

● the fsk signal group comprises 2FSK, 4FSK, 8FSK, 16FSK, 32FSK signal, 256 code elements of signal length, and 100 sampled points of each code element, frequency difference 100Hz between the code element, sample rate is 4 times of maximum modulating frequency;

● non-fsk signal group comprises BPSK, QPSK, 8PSK, 16QAM, 64QAM, OQPSK, CW signal, and 256 code elements of signal length, sample frequency are 8 times of signal baud rate, adopts the raised cosine window shaping pulse of rolloff-factor 0.8.

The signal that receives in the actual application environment, signal to noise ratio are unknown, also are unsettled.For the situation of more realistic application,, added the white noise of 5dB to 30dB at random to two groups of test signals that generate.Obtain two groups of simulate signals at last, respectively comprise 500 segment signal samples, calculate every segment signal respectively

Value is estimated rational thresholding τ _m, obtain test result as accompanying drawing 3.Transverse axis is the sample of signal numbering in the accompanying drawing 3, and the longitudinal axis is Value.As can be seen, the fsk signal group is obviously separated with non-fsk signal group, and reasonably the threshold value scope is τ _m=0.8～1.2.Work as τ _m=0.84 o'clock, the accuracy of identification reached 100%.

B) analysis of MFSK signal

Under 0～20dB signal to noise ratio condition, generation 2FSK, 4FSK, 8FSK, 16FSK, 32FSK signal are used for test.Signal length is 1024 code elements, and 100 sampled points of each code element, frequency difference 200Hz between the code element, signal sampling rate are 4 times of maximum modulating frequency.With the algorithm that the present invention proposes the signal that generates is analyzed, every kind of order of modulation/signal to noise ratio combination is repeated 100 tests, obtained the test result as accompanying drawing 4.Accompanying drawing 4 is that this algorithm of utilization is analyzed the curve of the analysis accuracy that obtains to the test signal of different modulating exponent number under different signal to noise ratio conditions, and transverse axis is the signal to noise ratio of test signal, and the longitudinal axis is to analyze accuracy, and every curve is represented a kind of modulation system.As can be seen, when SNR 〉=7.5dB, the analysis accuracy of 2FSK, 4FSK, 8FSK, 16FSK, 32FSK modulation signal has all been reached more than 95%.

2. at the test of actual reception signal (test is seen accompanying drawing 2 with hardware elementary diagram)

A) identification of MFSK signal

With seemingly, test data is divided into two groups at the test class of simulate signal:

● the MFSK signal, 2FSK (58), 4FSK (16), 6FSK (3), 8FSK (17), 12FSK (3), 13FSK (2), 16FSK (1), 20FSK (2), 32FSK (5) signal (being the sample number of this order of modulation signal in the bracket) have been comprised, the sample rate scope is 6000～8000Hz, chip rate scope 13Bd～1000Bd.

● non-MFSK signal, comprised BPSK (8), QPSK (2), 8PSK (9) signal (being the sample number of this modulation system signal in the bracket), the sample rate scope is 6000～8000Hz, the chip rate scope is 30Bd～2400Bd.

From two groups of samples, the random extraction time span is the signal segment of 1s (6000～8000 sampled point), calculates it

Value is estimated rational thresholding τ _m, obtain test result as accompanying drawing 5.In the accompanying drawing 5, transverse axis is the sample signal segment number, and the longitudinal axis is

Value.As can be seen from the figure, two groups of samples are well separated basically, and reasonably the threshold value scope is τ _m=0.8～1.2, conform to simulation result.Get τ _mIn the time of=1.17, the accuracy of separation has reached 95.6%.

B) analysis of MFSK signal

We utilize the actual reception signal that the parser of the MFSK signal of the present invention's proposition is tested, because sufficiently long 16FSK signal and 32FSK signal are difficult to obtain, test has included only 2FSK, 4FSK, 8FSK signal with set of signals.In the test every segment signal sample is divided into the data segment that length is 1s (6000～8000 sample points), analyzes with the MFSK signal analysis algorithm that the present invention proposes, the test result that obtains is as follows:

Signal type	The data hop count	Correct number of results	Accuracy (%)
				2FSK	1078	988	91.7
4FSK	127	115	90.6
				8FSK	169	155	91.8
Amount to	1374	1258	91.56

Test result from last table as can be seen, this algorithm has reached than higher accuracy the analysis of actual signal.

Claims (1)

1. the detection of a multi system frequency shift key control signal and analytical method is characterized in that, described method be in computer successively by

Following steps realize:

Step (1) receives pending data, obtains real signal x (n);

Step (2) converts real signal x (n) to analytic signal s (n) according to the following steps:

Step (2.1) is calculated FFT after real signal x (n) zero padding, obtains the discrete Fourier transform (DFT) F of signal _x(n), the length N of zero padding _zObtain by following formula:

N _z=N _FFT-N _x, N _xBe the length of real signal x (n),

Expression is greater than the smallest positive integral of " ", N _FFTBe the length of carrying out the FFT computing;

Step (2.2) is got described F _x(n) first half, n=1 ..., N _FFT/ 2, obtain the discrete Fourier transform (DFT) F of signal s (n) _s(n);

Step (2.3) is F _s(n) first half, n=1 ..., N _FFT/ 4 and latter half n=N _FFT/ 4+1 ..., N _FFT/ 2 transposings;

The F that step (2.4) calculation procedure (2.3) obtains _s(n) IFFT obtains analytic signal s (n), s (n)=IFFT (F _s(n));

The analytic signal s's (n) that step (2.5) obtains step (2.4)

Part is clipped, and remainder is the analytic signal that subsequent analysis is used, and its effective length is

Step (3) is according to the predefine thresholding, and segmentation detects the amplitude of s (n), removes the weak signal section, and its steps in sequence is as follows:

The signal s (n) that step (3.1) obtains step (2.5)=x (n)+jy (n) is divided into isometric N _SegSection is with a sequence S _i, i=1,2 ..., N _SegExpression, N _Seg=5～20;

Each the segment signal amplitude that obtains in step (3.2) calculation procedure (3.1) and m _i:

$m_i=\underset{s\left(n\right)\&Element;S_i}{\mathrm{\&Sigma;}}\left|s\left(n\right)\right|;$

Step (3.3) is set thresholding τ _m:

${\mathrm{\&tau;}}_m=\frac{1}{2}\max \{m_i,i=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;,N_{\mathrm{seg}}\};$

Step (3.4) is with m _i＜τ _mSignal segment be judged to be the weak signal section, remove all weak signal sections, remaining signal segment is connected into new analytic signal sequence s (n), length is N;

The normalization instantaneous frequency sequence f of step (4) calculation procedure (3.4) described analytic signal sequence s (n) _n(n), realize according to the following steps:

Step (4.1) is calculated as follows analytic signal sequence s (n)=x (n)+jy (n), n=1, and 2 ..., the instantaneous phase sequence of N

n＝1，2，……，N；

The instantaneous phase sequence that step (4.2) obtains step (4.1) Do and get the surplus instantaneous frequency sequence f (n) that obtains signal behind the calculus of differences:

n＝1，2，……，N-1；

Step (4.3) is normalized to zero-mean, unit variance sequence f to this instantaneous frequency sequence f (n) _n(n):

$f_n\left(n\right)=\frac{f\left(n\right)-m_f}{{\mathrm{\&sigma;}}_f},$

Wherein, $m_f=\frac{1}{N-1}\underset{n=1}{\overset{N-1}{\mathrm{\&Sigma;}}}f\left(n\right),$ ${\mathrm{\&sigma;}}_f=\sqrt{\frac{1}{N-2}\underset{n=1}{\overset{N-1}{\mathrm{\&Sigma;}}}{(f\left(n\right)-m_f)}^2};$

The f that step (5) calculation procedure (4.3) obtains _nThe rate of change sequence Δ f of normalization instantaneous frequency (n) _n(n) and the average of absolute value

Wherein:

${\mathrm{\&Delta;f}}_n\left(n\right)=|f_n(n+2)-f_n\left(n\right)|,$ n＝1，2，……，N-3， $m_{{\mathrm{\&Delta;f}}_n}=\frac{1}{N-3}\underset{n=1}{\overset{N-3}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;f}}_n\left(n\right)$

Step (6) is described

With predefine thresholding τ _mRelatively: if $M_{\mathrm{\&Delta;}{f}_n}<{\mathrm{\&tau;}}_m,$ Then described signal x (n) is judged to be the MFSK signal, and execution in step (7), otherwise, judge that receiving sequence is not the MFSK signal, flow process finishes;

Step (7) is carried out following steps successively and is reduced the shake of instantaneous frequency to utilize mean filter:

The instantaneous frequency sequence f that step (7.1) utilizes mean filter that step (4.3) is obtained _n(n) carry out smothing filtering, the instantaneous frequency sequence f that obtains upgrading _n(n), filter length l _fBy the possible maximum symbol rate f of analyzed signal _DmaxSetting sample frequency f with signal _sObtain by following formula:

$l_f=m\frac{.f_s}{f_{d\max }},$ $m=\frac{1}{3}~\frac{1}{2};$

The instantaneous frequency sequence of the renewal that step (7.2) obtains according to step (7.1) recomputates the normalization instantaneous frequency rate of change sequence Δ f of signal _n(n): Δ f _n(n)=| f _n(n+2)-f _n(n) |;

The variance of the normalization instantaneous frequency rate of change sequence that step (7.3) obtains according to following formula calculation procedure (7.2)

${\mathrm{\&tau;}}_{\mathrm{\&Delta;}{f}_n}=\sqrt{\frac{1}{N-4}\underset{n=1}{\overset{N-3}{\mathrm{\&Sigma;}}}({\mathrm{\&Delta;f}}_n\left(n\right)-m_{{\mathrm{\&Delta;f}}_n}{)}^2},$ $m_{{\mathrm{\&Delta;f}}_n}=\frac{1}{N-3}\underset{n=1}{\overset{N-3}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;f}}_n\left(n\right);$

Step (7.4) search Δ f _n(n) sequence, if $\left|{\mathrm{\&Delta;f}}_n\left(n\right)\right|<{\mathrm{\&tau;}}_{{\mathrm{\&Delta;f}}_n},$ F then _n(n+1) belong to symbol stable region Г _Si, i=1,2 ..., N _s, N wherein _sBe the number of symbol stable region, if $\left|{\mathrm{\&Delta;f}}_n\left(n\right)\right|>{\mathrm{\&tau;}}_{{\mathrm{\&Delta;f}}_n}$ F then _n(n+1) belong to the interval Г of symbol saltus step _Ti, i=1,2 ..., N _t, N wherein _tNumber for symbol saltus step interval;

Step (7.5) changes the instantaneous frequency value of each symbol stable region in the step (7.4) average of all instantaneous frequency values of this interval into, gives up the instantaneous frequency value in all symbol saltus step intervals, constitutes new instantaneous frequency sequence f _n(n);

Step (8) basis " entropy principle " is the histogram of the described instantaneous frequency sequential value of analytical procedure (7.5) according to the following steps:

Step (8.1) is according to normalization instantaneous frequency sequence f _n(n) set up distribution histogram With this as f _nThe estimation of probability density function (n):

$\hat{p}\left(f_i\right)=\frac{\mathrm{count}\{f_n\left(n\right),f_i\&le;f_n\left(n\right)<f_i+\mathrm{\&Delta;f}\}}{N_{f_n}},$ f _i＝min{f _n(n)}+iΔf，i＝0，1，……，

Wherein, Be f _n(n) length, $\mathrm{\&Delta;f}=(\max \left\{f_n\left(n\right)\right\}-\min \left\{f_n\left(n\right)\right\})/N_{\hat{p}},$ Be the length of distribution histogram, $N_{\hat{p}}=512~4096,$ Count{x, y} represent the to satisfy condition number of x value of y;

Step (8.2) is determined thresholding τ _pThe hunting zone be $0~\max \left\{\hat{p}\left(f_i\right)\right\}/2,$ Step-size in search ${\mathrm{\&Delta;\&tau;}}_p=\max \left\{\hat{p}\left(f_i\right)\right\}/\left(2K\right),$ K＝10～100；

Step (8.3) is to each τ in the hunting zone _pValue is carried out following steps:

Step (8.3.1) search is satisfied $\hat{p}\left(f_i\right)>{\mathrm{\&tau;}}_p$ F _i, continuous f _iBe classified as one group;

Step (8.3.2) is calculated the probability density function of each group and is estimated Sum obtains the estimated value of the distribution probability of each code element of MFSK signal

Step (8.3.3) is calculated should τ _pEntropy e:

$e=-\underset{k=1}{\overset{\hat{M}}{\mathrm{\&Sigma;}}}{\hat{P}}_k\log \left({\hat{P}}_k\right),$ Be continuous f _iThe number of the group that value is formed;

Step (8.4) is got the τ of corresponding entropy maximum _p, with corresponding f _iThe value of the number of group As the initial estimate of code element state number M, every group minimum instantaneous frequency f _MinkWith the maximum instantaneous frequency f _MaxkFrequency range (f as each code element state _Mink, f _Maxk), k=1,2 ...,

Step (9) basis " equiprobability distribution principle ", the exact value of accurate estimating code element state number M according to the following steps:

Step (9.1) is set thresholding $0<{\mathrm{\&tau;}}_{P_{\min }}<1;$

Step (9.2) is removed ${\hat{P}}_k<{\mathrm{\&tau;}}_{P_{\min }}\&times;\max \{{\hat{P}}_k,k=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;,\hat{M}\}$ State, herein, Be a span at 0～1 thresholding, $\max \{{\hat{P}}_k,k=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;,\hat{M}\}$ Be

In maximum, the two multiplies each other and obtains a value that is used to screen rational code element state;

Step (9.3) recomputates the estimated value of code element state number M

Frequency range (f with each code element state _Mink, f _Maxk), k=1,2 ...,

Step (10) is according to accurately estimating code element state number M and each state centre frequency such as " frequency difference principle ", and its step is as follows:

Step (10.1) is calculated the centre frequency of each code element state

${\hat{f}}_{\mathrm{mk}}=\frac{\underset{i=i_{\min k}}{\overset{i_{\max k}}{\mathrm{\&Sigma;}}}\hat{p}\left(f_i\right)f_i}{\underset{i=i_{\min k}}{\overset{i_{\max k}}{\mathrm{\&Sigma;}}}\hat{p}\left(f_i\right)}$

Wherein, i _MinkAnd i _MaxkBe the minimum of every set of symbols state, the label of maximum instantaneous frequency,

It is the estimation of the probability density function of the instantaneous frequency sequence that calculates of step (8.1);

Step (10.2) is calculated the frequency difference between each state

${\hat{f}}_{\mathrm{dk},l}=|{\hat{f}}_{\mathrm{mk}}-{\hat{f}}_{\mathrm{ml}}|;$

Step (10.3) is calculated the weight w of each frequency difference _Kl:

Weight equals to obtain the smaller value that the probability distribution of two code element states of this frequency difference is estimated, the probability distribution that deducts all the code element states between these two code element states is estimated sum, if difference less than 0, then weight is 0:

$w_{\mathrm{kl}}=\left\{\begin{array}{cc}\max \{0,\min \{{\hat{P}}_k,{\hat{P}}_l\}-\underset{n=k+1}{\overset{l-1}{\mathrm{\&Sigma;}}}{\hat{P}}_n\}& ,k<l-1\\ \max \{0,\min \{{\hat{P}}_k,{\hat{P}}_l\}\}& ,k=l-1\end{array}\right.;$

Step (10.4) obtains the estimation of frequency difference to each frequency difference weighted average:

${\hat{f}}_d=\frac{\underset{k,l}{\mathrm{\&Sigma;}}{w}_{\mathrm{kl}}{\hat{f}}_{\mathrm{dk},l}}{\underset{k,l}{\mathrm{\&Sigma;}}{w}_{\mathrm{kl}}}$

Step (10.5) is set thresholding τ _Min, τ _Max: ${\mathrm{\&tau;}}_{\min }=0.4{\hat{f}}_d~0.8{\hat{f}}_d,$ ${\mathrm{\&tau;}}_{\max }=1.2{\hat{f}}_d~1.6{\hat{f}}_d;$

Step (10.6) for ${\hat{f}}_{\mathrm{dk},k+1}<{\mathrm{\&tau;}}_{\min }$ Adjacent states, merge code element state k and k+1, wherein,

Adjacent code element state k that obtains in the expression step (10.2) and the frequency difference between the k+1;

Step (10.7) for ${\hat{f}}_{\mathrm{dk},k+1}>{\mathrm{\&tau;}}_{\max }$ Adjacent states, between code element state k and k+1, insert

Individual code element state, wherein, Expression is greater than the smallest positive integral of " ", Adjacent code element state k that obtains in the expression step (10.2) and the frequency difference between the k+1;

Step (10.8) is upgraded the estimation of code element state number

Centre frequency with each code element state

Step (11) subsequent treatment:

Be calculated as follows the final estimation that obtains MFSK signal order of modulation

${\hat{M}}_0=2^{\left[{\log }_2^{\hat{M}}\right]},$

Wherein [] computing is represented to get from " " nearest integer.

Images