CN104536022B - Low complexity frequency detector designing method for GNSS weak signal tracking - Google Patents

Abstract

The invention belongs to the technical field of satellite navigation signal processing, and particularly relates to a low complexity frequency detector designing method which is used for decreasing the calculated amount needed in prolonging the coherent accumulation time in the weak satellite navigation signal tracking process. The low complexity frequency detector designing method for the GNSS weak signal tracking comprises the specific steps that a frequency detector is initialized, and a maximum likelihood combined estimator for frequency offset and telegraph text is designed; a possible value range of tracking frequency offset is set, and uniformly distributed sampling points are selected from the value range; a function used for conducting quantification processing to an azimuth angle of a vector quantity on a two-dimensional plane under a polar coordinate system is set; according to ft (t=1, 2, ..L), Et is obtained through calculation, Et=Eacc (ft), and an initial estimated value of tracking frequency offset is obtained accordingly; the estimated value of tracking frequency offset is further modified through a three-point interpolation method. According to the low complexity frequency detector designing method, partial coherent accumulative results are used for representing telegraph text combination braches, and through the quantification processing on coherent values under the polar coordinate system, prediction and estimation, the telegraph text branches with smaller coherent accumulative energy values are removed, and therefore the calculated amount can be increased along with the to-be-estimated telegraph text bit number in a linear relationship.

Description

The low complex degree frequency discriminator method for designing that a kind of GNSS weak signal is followed the tracks of

Technical field

The invention belongs to satellite navigation signals processing technology field and in particular to a kind of GNSS weak signal follow the tracks of low complexity Degree frequency discriminator method for designing, reduces in processing in weak satellite navigation signals tracking and extends based on needed for the coherent accumulation time Calculation amount.

Background technology

With going deep into further of satellite navigation application, receiver is faced with more and more harsh applied environment.The whole world is led Boat satellite system (Global Navigation Satellite System, guide number SS) signal blocks ring in urban canyons etc. Weak more than ten or even tens decibels than in open region under border, and improve tracking sensitivity and can improve receiver in weak signal conditions Under position success rate and strengthen robustness.Frequency discriminator is the key link maintaining signal(-) carrier frequency to follow the tracks of, and its performance is to carrying High tracking sensitivity has material impact.Frequency discriminator is used for estimating signal from the on-time arm correlator sequence of correlator output Carrier frequency deducts the tracking frequency offset of remnants after local carrier frequency.As a frequency estimator, the frequency for obtaining high is estimated Meter precision, should extend the length of the used data segment of frequency discriminator single Frequency Estimation.In the case of can not obtaining text auxiliary, Frequency discriminator needs to realize coherent accumulation by estimating text.Although introduce in modernization GNSS signal not modulating leading of text Frequency component, can be used alone pilot data and estimates frequency deviation, but pilot tone data component is used in combination and can obtain more preferable frequency discrimination Performance.Therefore on the premise of given tracking frequency offset, according to the correlation (N is integer) corresponding to adjacent N number of text, first estimate The text bit going out, then is multiplied to peel off text modulation therein with estimating the text obtaining and correlation, then to stripping N number of correlation after text is added up, and can effectively improve frequency and differentiate required signal to noise ratio in operation.Due to text Estimate there may be error code, a kind of this cumulative simply approximate coherent accumulation, be hereinafter all referred to as phase for simplifying description Dry cumulative.Realizing approximate coherent accumulation is the important channel improving GNSS receiver weak signal tracking sensitivity.

The text algorithm for estimating that traditional weak signal is followed the tracks of is directly in the value domain search of text combination, its amount of calculation and 2^NBecome Linear relationship, but because it is mostly directed to the gps signal design that text speed is 50bps, the value of N is less.For text speed Rate is the gps signal of 50bps, and N takes 5～8 can obtain signal to noise ratio after sufficiently high coherent accumulation, text group now to be searched Close number less.No. two satellite navigation systems of the Big Dipper of China are (referred to as：BDS geo-synchronous orbit satellite) is (referred to as：GEO) institute The text speed broadcasting signal is 500bps, compared to 50bps signal, reaches identical coherent integration time, text to be searched Bit number increases by 10 times.For example：Realize the coherent integration of 0.12s, using traditional text algorithm for estimating, gps signal only need to be searched Rope 2⁵=32 text combinations, and the GEO signal demand for BDS searches for 2⁵⁹Individual text combination, i.e. amount of calculation increase by 2⁵⁴Times, that is, Amount of calculation is increased with exponential law, leads to algorithm to be difficult to real-time implementation.Therefore for the GNSS signal with higher text speed, Traditional text algorithm for estimating is difficult to the coherent accumulation of real-time implementation long period, and this is accomplished by a kind of computation complexity that reduces Text method of estimation.

Content of the invention

When extending the coherent accumulation time for the higher GNSS signal of text speed with traditional text method of estimation, meter The problem that calculation amount is increased with exponential law.The invention provides the low complex degree frequency discriminator design side that a kind of GNSS weak signal is followed the tracks of Method, the method is accumulated as based on the text stripping in the high-performance frequency discriminator of maximal possibility estimation and correlation with quickly realizing one Target, the coherent accumulation values due to direct access, therefore do not export the estimate to text bit, concrete technical scheme step is such as Under：

(S1) initialize frequency discriminator, the maximum likelihood joint estimator of design frequency deviation and text is as follows：

After being located at given tracking frequency offset f, search only for the obtained maximum coherence accumulated energies value of text combination

For E_acc(f), that is,：

Wherein, k=1,2 ..., N, D=[d₁, d₂..., d_N], d_k∈ {+1, -1 } is successively by the electricity that each bit is to be estimated The vector that literary composition is constituted, D represents the combined situation of N number of text bit, hereinafter referred to as text combination；F is tracking frequency offset；For The estimated result of tracking frequency offset, D_mlEstimated result for text combination；When arg () represents that object function takes maximum, corresponding Function argument value.

(S2) set the span [- f of tracking frequency offset f_max, f_max], and Δ f takes L to adopt at equal intervals wherein Sampling point：f₁＜ f₂＜ ... ＜ f_L-1＜ f_L, L is integer；

(S3) by planeCarry out quantification treatment according to azimuth, be quantified as the equally distributed ray of M bar, M is integer； WithFor interval, M discrete sampling point is taken on continuum [0,2 π]：For appoint Meaning vector to be quantified { S_x, S_y, if its polar coordinate representation is：

FromIn find and azimuthImmediate phase placeAs the side after quantifying Parallactic angle, that is, the cumulative vector after quantifyingFor：

(S4) the sampled point f to each tracking frequency offset value_t, t=1,2 ..., L, it is calculated E_t=E_acc(f_t), and Therefrom find out the maximum E_n=max { E_t, and with corresponding frequency f_nAs initial tracking frequency offset estimate, n is integer；

(S5) revise tracking frequency offset estimate further with three point interpolation method：If E₁Or E_LDuring for maximum, do not enter at row interpolation Reason, directly by offset estimation value amplitude limit be f₁Or f_L, otherwise, the coefficient of first computational representation peakdeviation amountAgain the output of frequency discriminator is modified to：

Further, sampled value f according to certain tracking frequency offset in described step (S4)_t, t=1,2 ..., L, solve E_t's Detailed process is：

(S41) sampled value f according to tracking frequency offset_t, phase place is carried out to the original correlation x [k] of correlator output Obtain x_r[k], computing formula is as follows：

IP in above formula_r[k] and QP_r[k] is respectively x_rThe real part of [k] and imaginary part；Formula (6) is substituted into formula (3) can get：

By x_r[k] regards bivector as, then d_kTake and can regard as x when -1_r[k] rotates 180 degree, d_k1 correspondence is taken not revolve Turn；

(S42) define planeIn subsetAnd its recurrence relation：By d₁It is fixed as 1, m upset vector before note With forI.e.：

AllThe collection that possible value is constituted is combined intoIt is apparent from,In an only element x_r[1], and Recurrence relation is as follows：

(S43) recurrence calculationIn have what maximum coherence accumulated energies value element was constituted in each quantized directions SubsetInitializationThen m+1=2, the step of 3 ..., N is：

1) arrangeFor empty set；

2) travel throughIn each element y, after the recurrence relation according to formula (9) calculates cumulative correlation vector New two producingIn element, wherein：

y₊Corresponding to d_m+1=1 situation, i.e. y₊=y+x_r[m+1]；

y_-Corresponding to d_m+1=-1 situation, i.e. y_-=y-x_r[m+1]；

3) judge y successively₊Or y_-It is added in which wayIn, specifically with y₊As a example illustrate determination methods, if y₊According to the azimuth after formula (4), formula (5) quantization it isIfIn no quantization back bearing be's Element, then by y₊It is added toIn；Otherwise setIn an existing element z quantify back bearing and be also Then compare y₊With the modulus value size of z, if | y₊| ＞ | z |, then willIn z y₊Replace, otherwise keepConstant；

According to step 1), 2), 3) method, obtain after iteration N-1 time

(S44) find outIn there is the element y of maximum modulus value_m, then calculate E_t, that is, take E_t=| y_m|².

Further, M value is 2⁹、2¹⁰Etc. less number.

The present invention possesses following beneficial effect：The present invention effectively overcomes traditional text method of estimation for extending mirror During the coherent accumulation time of frequency device, the problem that amount of calculation and text bit number to be estimated are increased with exponential law, there is provided Plant the amount of calculation text method of estimation linear with text bit number to be estimated, the frequency discriminator design that the present invention provides simultaneously Maximum likelihood tracking frequency offset estimator that can be optimum effectively on approximation theory.Can be used for improving with less cost by this present invention GNSS receiver follows the tracks of the sensitivity of the GNSS signal with higher text speed.

Brief description

Fig. 1 represents that weak signal conditions lower frequency locked loop realizes block diagram；

Fig. 2 represents the flowchart of frequency discriminator in the present invention；

Fig. 3 represents that block diagram realized by tracking frequency offset estimator；

Fig. 4 represents the flow chart calculating maximum coherence accumulated energies value in the present invention；

Fig. 5 represents the geometrical principle schematic diagram of exclusion text combination；

Fig. 6 be carrier-to-noise ratio 20dB-Hz, no frequency change rate when two kinds of frequency discriminator tracking performance comparison diagrams；

Fig. 7 is two kinds of frequency discriminator tracking performance comparison diagrams when carrier-to-noise ratio 23dB-Hz, frequency change rate 3Hz/s.

Specific embodiment

Below, in conjunction with the drawings and specific embodiments, the invention will be further described.

As shown in figure 1, the frequency discriminator being applied to the frequency locked loop (FLL) of the GEO satellite signal of BDS with the present invention is set It is calculated as example to illustrate.Receipt signal r (t) sends into correlator after code and carrier wave are peeled off, and the integrating range of correlator is just It is well a text code element, that is, the time of integration is T_c=2ms.Differentiate performance, the original correlation of correlator output for improving frequency Value sequence often N number of cache for send into after one group can Combined estimator text and frequency deviation frequency discriminator, frequency identification result passes through loop Control the digital controlled oscillator (NCO) producing local carrier, loop updates and is spaced apart T=NT after wave filter_c.In relatively low dynamic condition It is believed that down tracking frequency offset immobilizes in a loop update cycle, ignores code tracking error, with k=1,2 ..., N are successively Represent correlation sequence number, then one group of used correlation of frequency discriminator single work can be by following model representation：

In above formula：K=1,2 ..., N, C/N₀For the carrier-to-noise ratio of signal, A is the signal width later to noise energy normalization Degree, η_i、η_qFor the noise section in correlation, obey standardized normal distribution, f_errFor the tracing deviation of carrier frequency, referred to as follow the tracks of Frequency deviation, φ is the phase deviation of carrier track,For text symbol, value is that+1 or -1, IP [k] is respectively with QP [k] The in-phase component of on-time arm correlator and quadrature component, can be written as plural form x [k]=IP [k]+1j*QP [k], and j represents imaginary number list Position.

In FLL loop, frequency discriminator is used for estimating tracking frequency offset f according to x [k]_err, improve the key of tracking sensitivity It is to improve to f_errEstimated accuracy.Model according to formula (1) and theory analysis, for estimating f_errOptimal estimation device It is the following maximum likelihood estimator module with regard to text and tracking frequency offset：

In above formula：D=[d₁, d₂..., d_N], d_k∈ {+1, -1 } be successively from each bit text to be estimated constitute to Amount, represents the combined situation of N number of text bit, hereinafter referred to as text combination；F is tracking frequency offset；For tracking frequency offset Estimated result, D_mlEstimated result for text combination.

The object function value of calculation optimization problem formula (2) can be regarded as and first carry out phase place to original correlation, Then peel off text modulation therein, finally seek the energy value of coherent accumulation results.Frequency deviation is overcome to phase by phase place The dry cumulative impact bringing.After being located at given tracking frequency offset f, search only for the obtained maximum coherence accumulated energies value of text combination For E_acc(f), that is,：

The present invention proposes a kind of frequency discriminator design for approximately estimating tracking frequency offset according to formula (2), and its technical scheme is： Using E_accF () is the feature of continuous function, choose some sampled points in the possible span of frequency deviation, and in these samplings E is calculated on point_accF the value of (), then to estimate E using three point interpolation method_accF the peak of (), thus obtains tracking frequency offset Estimate.

As shown in Figure 2, block diagram realized by the frequency discriminator of the present invention, specifically need to execute steps S2～S5：

(S2) set possible the span [- f of tracking frequency offset_max, f_max], and Δ f takes L at equal intervals wherein Sampled point：f₁＜ f₂＜ ... ＜ f_L-1＜ f_L, L is integer.

(S3) set to planeThe azimuth of middle vector carries out the function of quantification treatment, willIt is equal that plane is quantified as M bar The ray of even distribution, M is integer.WithFor interval, M discrete sampling point is taken on continuum [0,2 π]：For vector { S arbitrarily to be quantified_x, S_y, if its polar coordinate representation is：

Then quantification treatment is：FromIn find withImmediate phase placeAs quantization Azimuth afterwards, that is, the cumulative vector after quantifyingFor：

M is generally taken to be 2⁹、2¹⁰Left and right, correspondinglyForSufficiently high quantified precision now can be ensured.

(S4) the sampled point f to the possible value of each tracking frequency offset_t, t=1,2 ..., L, it is calculated E_t=E_acc (f_t), and therefrom find out the maximum E_n=max { E_t, and with f_nAs initial tracking frequency offset estimate.

(S5) revise tracking frequency offset estimate further with three point interpolation method：If E₁Or E_LDuring for maximum, do not enter at row interpolation Reason, directly by offset estimation value amplitude limit be f₁Or f_L, the otherwise coefficient of first computational representation peakdeviation amountAgain the output of frequency discriminator is modified to：

The present invention reduces it is critical only that of computation complexity and calculates E_accWhen (f), quantified using the azimuth under polar coordinate system To realize the minimizing of amount of calculation with dynamic programming algorithm thought.

As shown in Figure 4, for calculating energy value E after maximum coherence adds up_tFlow process, specifically include steps S41 ～S44：

(S41) sampled value f according to certain tracking frequency offset concrete_t, the original correlation x [k] of correlator output is carried out Phase place obtains x_r[k], computing formula is as follows：

IP in above formula_r[k] and QP_r[k] is respectively x_rThe real part of [k] and imaginary part.Formula (6) is substituted into formula (3) can get：

By x_r[k] regards bivector as, then d_kTake and can regard as x when -1_r[k] rotates 180 degree, d_k1 correspondence is taken not revolve Turn.Therefore calculate E_tGiven N number of direction can be regarded as can the bivector of 180 degree upset, and finding out a kind of upset combination makes Obtain itself and vectorial { S_i, S_qThere is maximum energy value.

(S42) define planeIn subsetAnd its recurrence relation.Due to by all of d_kAll negate (i.e. all of d_kSymbol negates simultaneously, and 1 is changed into -1, and -1 is changed into 1) do not affect E_tValue, therefore, by d₁It is fixed as 1.M turning vector before note Sum isI.e.：

AllThe collection that possible value is constituted is combined intoIt is apparent from,In an only element x_r[1], and Meet following recurrence relation：

(S43) recurrence calculationIn have what maximum coherence accumulated energies value element was constituted in each quantized directions SubsetInitializationM+1=2, the step of 3 ..., N is：

1) arrangeFor empty set.

2) travel throughIn each element y, after the recurrence relation according to formula (9) calculates cumulative correlation vector New two producingIn element, wherein：

y₊Corresponding to d_m+1=1 situation, i.e. y₊=y+x_r[m+1]；

y_-Corresponding to d_m+1=-1 situation, i.e. y_-=y-x_r[m+1].

3) judge y successively₊Or y_-It is added in which wayIn, with y₊As a example illustrate determination methods (y_-Sentence Disconnected situation is identical).If y₊According to the azimuth after formula (4), formula (5) quantization it isIfIn no quantization back side Angle isElement, then by y₊It is added toIn；Otherwise setIn an existing element z quantify back side Angle is alsoThen compare y₊With the modulus value size of z, if | y₊| ＞ | z |, then willIn z y₊Replace, otherwise protect HoldConstant.

Obtain after iteration N-1 time

(S44) find outIn there is the element y of maximum modulus value_m, then calculate E_t, that is, take E_t=| y_m|².

For fully understanding the present invention, its relative theory is briefly described as follows：

Say from the strict sense, Ying CongThe element that middle search has maximum modulus value just can obtain E_t, but adopt this kind of side Method will not reduce amount of calculation.Because the noise section in correlation is uncorrelated, accumulation result can be confined to of two dimensional surface In smaller area.Increase with mIn element number exponentially regular increase, for larger m,Middle element Distribution will assume the feature of a large amount of gatherings.Exactly make use of These characteristics, the present invention allow certain error in the case of, fromIn select a small amount of most probable and continue the decision-making that the cumulative element obtaining maximum energy value participates in next stage, substantial amounts of text Combination branch is excluded in advance.CalculateIterative process both embodied the thought of above-mentioned minimizing amount of calculation, exclude InDuring element in addition, employ following geometrical principle and approximate processing：

Geometrical principle：ForIn two azimuth identicals vectorsAs long asThen fromContinue the cumulative vector that can not possibly obtain and there is maximum modulus value.Above-mentioned principle provable as follows：As shown in Figure 5, if fromSet out and continue the cumulative vector having obtained ceiling capacity, that is,：After subsequently the cumulative newly-increased vector of N-m step is υ,Tool There is maximum modulus value.Due toIt is also a vector that may add up out, thereforeAnother Aspect, due toIt is also the vector obtaining that can add up, thereforeThus obtainIn r Amass as nonnegative number, that is, angle therebetween is less than or equal to 90 degree.In the case, from accompanying drawing 5 it will be evident thatThis is just derived contradiction.

Approximate processing：For twoIn vector, as long as its quantization after azimuth identical, then be approximately considered its In same direction, willIn sector region come approximately with ray in the middle of it.

According toCalculating process understand, be quantized in the sector region of ray at each and at most only haveIn An element, thereforeMiddle element number is less than M.Obviously calculateThe amount of calculation of each iterative step with Middle element number is directly proportional, and therefore the computation complexity in the present invention is less than O (M × N).

To verify that by numerical simulation the present invention is directed to the practical application effect of the GEO signal of BDS below.

Emulation is carried out for the GEO satellite of BDS, and specific parameter setting is as follows：

Participate in text times N=60 that single frequency differentiates, that is, the coherent accumulation time is T=0.12s；

Azimuth quantified precision parameter is taken as M=2¹⁰；

The span of tracking frequency offset is set to [- 10Hz, 10Hz], and thinks that 2Hz is that interval takes 11 sampled points wherein.

With text algorithm for estimating and the frequency discriminator of C++ programming realization this paper during emulation, using 64 bit CPUs and dominant frequency For running on the server of 2.4GHz, with run time come the amount of calculation of quantization means algorithm.Obtained by counting run time, Single text is estimated (to calculate an E_t), 14ms is about using the inventive method required time, and traditional exhaustion Time needed for searching method is more than 1500s, and that is, context of methods is used for amount of calculation during search 60 bit text is traditional algorithm Less than 0.00001 times so that for 500bps signal be up to the coherent accumulation of 0.12s may real-time implementation.Calculate herein Method needs storage at most may have the set of M elementTherefore required amount of storage is linear with M, and tradition is calculated Method only needs to storage and currently has stepped through the person that has ceiling capacity in text combination, and required storage is few.Thus reflect dynamic Planning algorithm changes the feature of time complexity with space complexity.

The present invention realizes the tracking performance lifting that long-time coherent accumulation brought, be reflected in following to whole FLL ring Road is in the emulation of static and dynamic two kinds of situations.Using conventional second-order F LL loop during emulation, it is respectively adopted the mirror of the present invention Frequency device and traditional frequency discriminator (the used data length of single frequency discrimination is also 0.12s) obtain two track loop.Doppler frequency Initial value is -180Hz, emulates total length of data 120s, and loop is started working from tenacious tracking state, i.e. initial velocity and acceleration Degree tracking error is 0Hz.Only consider carrier track it is assumed that nonexistent code tracing deviation.With loop noise bandwidth B_LWith loop Update the product B of interval T_LT is describing the main feature of track loop.The quiescent conditions being 0 for doppler changing rate, adopt Less loop bandwidth, according to B_LT=0.01 designs loop.It is 3Hz/s current intelligence for doppler changing rate, using larger Loop bandwidth, according to B_LT=0.025 designs loop.Under quiescent conditions, for 20dB-Hz signal, the tracking of two kinds of loops is by mistake Difference is as shown in Figure 6；Under current intelligence, for 23dB-Hz signal, the tracking error of two kinds of loops is as shown in Figure 7.Static feelings Under condition, the standard deviation using carrier frequency tracking error during new frequency discriminator is 2.18Hz, using carrier wave during traditional frequency discriminator The standard deviation of frequency tracking error is 14.54Hz.Under current intelligence, using carrier frequency tracking error during new frequency discriminator Standard deviation is 1.27Hz, and the standard deviation using carrier frequency tracking error during traditional frequency discriminator is 19.96Hz.In weak signal conditions Under, no matter having or not signal dynamics, the frequency tracking error of this paper track algorithm is less than 5Hz, and the frequency of traditional track algorithm with Track error is up to 30～40Hz.Further emulation shows, for traditional frequency discriminator be reach suitable with frequency discriminator of the present invention Tracking accuracy, is respectively 27dB-Hz and 30dB-Hz, i.e. frequency discriminator of the present invention for the carrier-to-noise ratio dynamically and needed for stationary singnal The tracking sensitivity making FLL track loop improves about 7dB.

The above is only the preferred embodiment of the present invention, protection scope of the present invention is not limited merely to above-described embodiment, All technical schemes belonging under thinking of the present invention belong to protection scope of the present invention.It should be pointed out that for the art For those of ordinary skill, some improvements and modifications without departing from the principles of the present invention, should be regarded as the protection of the present invention Scope.

Claims (3)

1. the low complex degree frequency discriminator method for designing that a kind of GNSS weak signal is followed the tracks of is it is characterised in that comprise the following steps：

(S1) initialize frequency discriminator, the maximum likelihood joint estimator of design frequency deviation and text is as follows：

$\{{\hat{f}}_{ml},D_{ml}\}=\arg \left\{\underset{f}{max}\underset{D}{max}\right\{|\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{d}_{k}x\[k\]e^{-j2\mathrm{\π}fk}{|}^2\}\}---\left(2\right)$

After given tracking frequency offset f, searching only for the obtained maximum coherence accumulated energies value of text combination is E_acc(f), that is,：

$E_{acc}\left(f\right)=\underset{D}{max}\left\{\right|\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{d}_{k}x\[k\]e^{-j2\mathrm{\π}fk}{|}^2\}---\left(3\right)$

Wherein, k=1,2 ..., N, D=[d₁, d₂..., d_N], d_k∈ {+1, -1 } is to be made up of each bit text to be estimated successively Vector, D represents the combined situation of N number of text bit, hereinafter referred to as text combination；F is tracking frequency offset；For following the tracks of frequency Inclined estimated result, D_mlEstimated result for text combination；

(S2) set the span [- f of tracking frequency offset f_max, f_max], and Δ f takes L sampled point at equal intervals wherein： f₁＜ f₂＜ ... ＜ f_L-1＜ f_L, L is integer；

(S3) by planeCarry out quantification treatment according to azimuth, be quantified as the equally distributed ray of M bar, M is integer；With For interval, M discrete sampling point is taken on continuum [0,2 π]：For the arbitrarily amount for the treatment of Change vectorial { S_x, S_y, if its polar coordinate representation is：

FromIn find and azimuthImmediate phase placeAs the azimuth after quantifying, Cumulative vector after quantifyingFor：

(S4) the sampled point f to each tracking frequency offset value_t, t=1,2 ..., L, it is calculated E_t=E_acc(f_t), and therefrom Find out the maximum E_n=max { E_t, and with corresponding frequency f_nAs initial tracking frequency offset estimate, n is integer；

(S5) revise tracking frequency offset estimate further with three point interpolation method：If E₁Or E_LDuring for maximum, do not carry out interpolation processing, Direct is f by offset estimation value amplitude limit₁Or f_L, otherwise, the coefficient of first computational representation peakdeviation amount Again the output of frequency discriminator is modified to：

2. a kind of GNSS weak signal is followed the tracks of as claimed in claim 1 low complex degree frequency discriminator method for designing it is characterised in that Sampled value f according to tracking frequency offset in described step (S4)_t, t=1,2 ..., L, solve E_tDetailed process be：

(S41) sampled value f according to tracking frequency offset_t, the original correlation x [k] of correlator output is carried out by phase place obtains x_r [k], computing formula is as follows：

$x_r\[k\]={\mathrm{IP}}_r\[k\]+{\mathrm{QP}}_r\[k\]\·1j=x\[k\]e^{-j2{\mathrm{\πf}}_{t}k}---\left(6\right)$

IP in above formula_r[k] and QP_r[k] is respectively x_rThe real part of [k] and imaginary part；Formula (6) is substituted into formula (3) can get：

$\begin{array}{c}{E}_t=\underset{D}{\max }\left\{{S_i}^2\left(D\right)+{S_q}^2\left(D\right)\right\}\\ S_i\left(D\right)=\underset{k=1}{\overset{N}{\Σ}}{\mathrm{IP}}_r\[k\]d_k,S_q\left(D\right)=\underset{k=1}{\overset{N}{\Σ}}{\mathrm{QP}}_r\[k\]d_k\end{array}---\left(7\right)$

By x_r[1] regard bivector as, then d_kTake and can regard as x when -1_r[k] rotates 180 degree, d_k1 correspondence is taken not rotate；

(S42) define planeIn subsetAnd its recurrence relation：By d₁It is fixed as 1, before note, m upset vector sum isI.e.：

$S_i^m=\underset{k=1}{\overset{m}{\Σ}}{\mathrm{IP}}_r\[k\]d_k,S_q^m=\underset{k=1}{\overset{m}{\Σ}}{\mathrm{QP}}_r\[k\]d_k,m=1,2,...,N---\left(8\right)$

AllThe collection that possible value is constituted is combined intoIt is apparent from,In an only element x_rAnd recursion [1], Relation is as follows：

(S43) recurrence calculationIn the subset that maximum coherence accumulated energies value element is constituted is had on each quantized directionsInitializationThen m+1=2, the step of 3 ..., N is：

1) arrangeFor empty set；

2) travel throughIn each element y, the recurrence relation according to formula (9) newly produces after calculating cumulative correlation vector Raw twoIn element, wherein：

y₊Corresponding to d_m+1=1 situation, i.e. y₊=y+x_r[m+1]；

y_-Corresponding to d_m+1=-1 situation, i.e. y_-=y-x_r[m+1]；

3) judge y successively₊Or y_-It is added in which wayIn, specifically with y₊As a example illustrate determination methods, if y₊Press Azimuth after illuminated (4), formula (5) quantify isIfIn no quantization back bearing beUnit Element, then by y₊It is added toIn；Otherwise setIn an existing element z quantify back bearing and be alsoThen Relatively y₊With the modulus value size of z, if | y₊| ＞ | z |, then willIn z y₊Replace, otherwise keepConstant；

According to step 1), 2), 3) method, obtain after iteration N-1 time

(S44) find outIn there is the element y of maximum modulus value_m, then calculate E_t, that is, take E_t=| y_m|².

3. a kind of GNSS weak signal is followed the tracks of as claimed in claim 1 low complex degree frequency discriminator method for designing it is characterised in that M value is 2¹⁰.