CN104536022A - 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, be specifically related to the low complex degree frequency discriminator method for designing that a kind of GNSS weak signal is followed the tracks of, for following the tracks of the calculated amount reducing in process and extend needed for the coherent accumulation time at weak satellite navigation signals.

Background technology

Along with going deep into further of satellite navigation application, receiver is faced with more and more harsh applied environment.GLONASS (Global Navigation Satellite System) (Global Navigation Satellite System, guide number SS) signal under urban canyons etc. blocks environment than in open region weak tens even tens decibels, and improve tracking sensitivity and can improve receiver position success rate in weak signal conditions and strengthen robustness.Frequency discriminator is the key link maintaining signal(-) carrier frequency tracking, and its performance has material impact to raising tracking sensitivity.The tracking frequency offset that frequency discriminator is remaining after deducting local carrier frequency for the carrier frequency estimating signal the on-time arm correlator sequence from correlator output.As a frequency estimator, for obtaining high frequency estimation accuracy, the length of frequency discriminator single Frequency Estimation institute usage data section should be extended.When not obtaining text and be auxiliary, frequency discriminator needs by estimating that text realizes coherent accumulation.Although introduce the weight of pilot frequency not modulating text in modernization GNSS signal, can be used alone pilot data and estimate frequency deviation, conbined usage pilot tone and data component can obtain better frequency discrimination performance.Therefore under the prerequisite of given tracking frequency offset, according to the correlation (N is integer) corresponding to adjacent N number of text, the text bit first estimated, then with estimating that the text obtained is multiplied by with correlation the text peeled off wherein mutually and modulates, then the N number of correlation peeled off after text is added up, effectively can improve frequency and differentiate signal to noise ratio (S/N ratio) required in operation.Because text estimates to there is error code, this cumulative just a kind of approximate coherent accumulation, for simplified characterization is hereinafter all referred to as coherent accumulation.Realizing approximate coherent accumulation is the important channel of improving GNSS receiver weak signal tracking sensitivity.

The text algorithm for estimating that tradition weak signal is followed the tracks of directly is searched in the codomain of text combination, its calculated amount and 2 ⁿlinear, but be mostly the gps signal design of 50bps for text speed due to it, the value of N is less.Be the gps signal of 50bps for text speed, N gets 5 ~ 8 can obtain signal to noise ratio (S/N ratio) after sufficiently high coherent accumulation, and text number of combinations now to be searched is less.The Big Dipper No. two satellite navigation systems of China (be called for short: geo-synchronous orbit satellite BDS) (be called for short: GEO) broadcast signal text speed be 500bps, compared to 50bps signal, reach identical coherent integration time, text bit number to be searched increases by 10 times.Such as: the coherent integration realizing 0.12s, adopt traditional text algorithm for estimating, only need search for 2 to gps signal ⁵=32 text combinations, and 2 are searched for for the GEO signal demand of BDS ⁵⁹individual text combination, namely calculated amount increases by 2 ⁵⁴doubly, namely calculated amount increases with exponential law, causes 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 just needs a kind of text method of estimation reducing computation complexity.

Summary of the invention

When higher GNSS signal for text speed extends the coherent accumulation time by traditional text method of estimation, the problem that calculated amount increases with exponential law.The invention provides the low complex degree frequency discriminator method for designing that a kind of GNSS weak signal is followed the tracks of, the method to be peeled off based on the text in the high-performance frequency discriminator of maximal possibility estimation to realize one fast and relevant to be accumulated as target, owing to directly obtaining coherent accumulation values, therefore do not export the estimated value to text bit, concrete technical scheme steps is as follows:

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

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

After being located at given tracking frequency offset f, only searching for text and combine the maximum coherence accumulated energies value obtained

For E _acc(f), that is:

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

Wherein, k=1,2 ..., N, D=[d ₁, d ₂..., d _n], d _k{+1 ,-1} is the vector be made up of the text that each bit is to be estimated successively to ∈, and D represents the combined situation of N number of text bit, combines hereinafter referred to as text; F is tracking frequency offset; for the estimated result of tracking frequency offset, D _mlfor the estimated result of text combination; When arg () represents that objective function gets maximal value, the value of corresponding function argument.

(S2) the span [-f of tracking frequency offset f is set _max, f _max], and get L sampled point: f with equal interval delta f wherein ₁< f ₂< ... < f _l-1< f _l, L is integer;

(S3) by plane carry out quantification treatment according to position angle, be quantified as the equally distributed ray of M bar, M is integer; With for interval, continuum [0,2 π] gets M discrete sampling point: for any vector { S to be quantified _x, S _y, if its polar coordinate representation is:

From in find and position angle immediate phase place as the position angle after quantification, the cumulative vector after namely quantizing for:

(S4) to the sampled point f of each tracking frequency offset value _t, t=1,2 ..., L, calculates E _t=E _acc(f _t), and therefrom find out the maximum E _n=max{E _t, and with the frequency f of correspondence _nas initial tracking frequency offset estimated value, n is integer;

(S5) tracking frequency offset estimated value is revised further by three point interpolation method: if E ₁or E _lfor time maximum, not carrying out interpolation processing, is directly f by frequency deviation estimated value amplitude limit ₁or f _l, otherwise, the coefficient of first computational representation peakdeviation amount again the output of frequency discriminator is modified to: $\hat{f}=f_n+2\mathrm{\&alpha;\&Delta;f}.$

Further, in described step (S4) according to the sampled value f of certain tracking frequency offset _t, t=1,2 ..., L, solves E _tdetailed process be:

(S41) according to the sampled value f of tracking frequency offset _t, phase rotating is carried out to the original correlation x [k] of correlator output and obtains x _r[k], computing formula is as follows:

$x_r\left[k\right]={\mathrm{IP}}_r\left[k\right]+{\mathrm{QP}}_r\left[k\right]\&CenterDot;1j=x\left[k\right]e^{-j2\mathrm{\&pi;}{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 obtain:

$E_t=\underset{D}{\max }\{S_i{\left(D\right)}^2+S_q{\left(D\right)}^2\}$

$S_i\left(D\right)=\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{IP}}_r\left[k\right]d_k,S_q\left(D\right)=\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{QP}}_r\left[k\right]d_k---\left(7\right)$

By x _r[k] regards bivector as, then d _kcan regard as x when getting-1 _r[k] revolves turnback, d _kget 1 corresponding non rotating;

(S42) plane is defined in subset and recurrence relation: by d ₁be fixed as 1, before note, m the vectorial sum of upset is that is:

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

All the set that possible value is formed is yi Zhi, in only have an element x _rand recurrence relation is as follows [1]:

(S43) recurrence calculation in there is the subset that maximum coherence accumulated energies value element forms in each quantized directions initialization then m+1=2,3 ..., the step of N is:

1) arrange for empty set;

2) travel through in each element y, calculate two that a cumulative correlation vector newly produces afterwards according to the recurrence relation of formula (9) in element, wherein:

Y ₊corresponding to d _m+1the situation of=1, i.e. y ₊=y+x _r[m+1];

Y _-corresponding to d _m+1the situation of=-1, i.e. y _-=y-x _r[m+1];

3) y is judged successively ₊or y _-add in which way in, specifically with y ₊for example illustrates determination methods, if y ₊according to the position angle after formula (4), formula (5) quantize be if in no quantization back bearing be element, then by y ₊add to in; Otherwise establish in an existing element z quantize back bearing and be also then compare y ₊with the modulus value size of z, if | y ₊| > | z|, then will in z y ₊replace, otherwise keep constant;

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

(S44) find out in there is the element y of maximum modulus value _m, then calculate E _t, namely get E _t=| y _m| ².

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

The present invention possesses following beneficial effect: the present invention effectively overcomes traditional text method of estimation when the coherent accumulation time for extending frequency discriminator, the problem that calculated amount and text bit number to be estimated increase with exponential law, provide a kind of calculated amount and the linear text method of estimation of text bit number to be estimated, frequency discriminator design provided by the invention simultaneously can the effective maximum likelihood tracking frequency offset estimator of optimum on approximation theory.The present invention can be used for improving with less cost the sensitivity that GNSS receiver tracking has the GNSS signal of higher text speed thus.

Accompanying drawing explanation

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

Fig. 2 represents the realization flow figure of frequency discriminator in the present invention;

Fig. 3 represents that tracking frequency offset estimator realizes block diagram;

Fig. 4 represents in the present invention the process flow diagram calculating maximum coherence accumulated energies value;

Fig. 5 represents the geometrical principle schematic diagram getting rid of text combination;

Fig. 6 is carrier-to-noise ratio 20dB-Hz, without two kinds of frequency discriminator tracking performance comparison diagrams during frequency change rate;

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

Embodiment

Below, by reference to the accompanying drawings and embodiment, 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 designed to example and is described.Received signal strength r (t) sends into correlator after code and carrier wave are peeled off, and the integrating range of correlator is just in time a text code element, and namely integral time is T _c=2ms.Performance is differentiated for improving frequency, the every N number of buffer memory of original sequence of correlation values that correlator exports is that send into after one group can the frequency discriminator of Combined estimator text and frequency deviation, frequency identification result is by controlling the digital controlled oscillator (NCO) producing local carrier after loop filter, loop upgrades and is spaced apart T=NT _c.Can think that under lower dynamic condition tracking frequency offset immobilizes within a loop update cycle, ignore code tracking error, with k=1,2 ..., N represents correlation sequence number successively, then one group of correlation using of frequency discriminator single work can by following model representation:

$\begin{array}{cc}\begin{array}{c}\mathrm{IP}\left[k\right]=A\stackrel{\&OverBar;}{d_k}\cos (2\mathrm{\&pi;}{f}_{\mathrm{err}}\mathrm{Tk}+\mathrm{\&phi;})+{\mathrm{\&eta;}}_i\left[k\right]\\ \mathrm{QP}\left[k\right]=A\stackrel{\&OverBar;}{d_k}\sin (2\mathrm{\&pi;}{f}_{\mathrm{err}}\mathrm{Tk}+\mathrm{\&phi;})+{\mathrm{\&eta;}}_q\left[k\right]\end{array}& A=\sqrt{\frac{2{\mathrm{CT}}_C}{N_0}}\&CenterDot;\frac{\sin \left(\mathrm{\&pi;}{f}_{\mathrm{err}}{T}_C\right)}{\left(\mathrm{\&pi;}{f}_{\mathrm{err}}{T}_C\right)}---\left(1\right)\end{array}$

In above formula: k=1,2 ..., N, C/N ₀for the carrier-to-noise ratio of signal, A is the signal amplitude later to noise energy normalization, η _i, η _qfor the noise section in correlation, obey standardized normal distribution, f _errfor the tracing deviation of carrier frequency, be called for short tracking frequency offset, φ is the phase deviation of carrier track, for text symbol, value is in-phase component and the quadrature component that+1 or-1, IP [k] and QP [k] are respectively on-time arm correlator, and can be written as plural form x [k]=IP [k]+1j*QP [k], j represents imaginary unit.

In FLL loop, frequency discriminator is used for estimating tracking frequency offset f according to x [k] _err, the key improving tracking sensitivity is to improve f _errestimated accuracy.According to model and the theoretical analysis of formula (1), for estimating f _erroptimal estimation device be maximum likelihood estimator module as follows about text and tracking frequency offset:

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

In above formula: D=[d ₁, d ₂..., d _n], d _k{+1 ,-1} is the vector be made up of the text that each bit is to be estimated successively to ∈, represents the combined situation of N number of text bit, combines hereinafter referred to as text; F is tracking frequency offset; for the estimated result of tracking frequency offset, D _mlfor the estimated result of text combination.

The objective function value of calculation optimization problem formula (2) can be regarded as and first carry out phase rotating to original correlation, then peels off text modulation wherein, finally asks the energy value of coherent accumulation results.The impact that frequency deviation brings to coherent accumulation is overcome by phase rotating.After being located at given tracking frequency offset f, only search text combines the maximum coherence accumulated energies value obtained is E _acc(f), that is:

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

The present invention proposes a kind of for the approximate frequency discriminator design estimating tracking frequency offset according to formula (2), and its technical scheme is: utilize E _accf feature that () is continuous function, chooses some sampled points, and calculate E on these sampled point in the span that frequency deviation is possible _accf the value of (), then adopts three point interpolation method to estimate E _accf the peak of (), obtains the estimated value of tracking frequency offset thus.

As shown in Figure 2, frequency discriminator of the present invention realizes block diagram, specifically need perform following step S2 ~ S5:

(S2) possible the span [-f of tracking frequency offset is set _max, f _max], and get L sampled point: f with equal interval delta f wherein ₁< f ₂< ... < f _l-1< f _l, L is integer.

(S3) setting is to plane carry out the function of quantification treatment at the position angle of middle vector, will plane is quantified as the equally distributed ray of M bar, M is integer.With for interval, continuum [0,2 π] gets M discrete sampling point: for any vector { S to be quantified _x, S _y, if its polar coordinate representation is:

Then quantification treatment is: from in find with immediate phase place as the position angle after quantification, the cumulative vector after namely quantizing for:

Usually getting M is 2 ⁹, 2 ¹⁰left and right, correspondingly for now can ensure sufficiently high quantified precision.

(S4) to the sampled point f of each tracking frequency offset possibility value _t, t=1,2 ..., L, calculates 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 estimated value.

(S5) tracking frequency offset estimated value is revised further by three point interpolation method: if E ₁or E _lfor time maximum, not carrying out interpolation processing, is directly f by frequency deviation estimated value amplitude limit ₁or f _l, otherwise the coefficient of first computational representation peakdeviation amount $\mathrm{\&alpha;}=\frac{E_{n+1}-E_{n-1}}{E_{n+1}+E_{n-1}},$ Again the output of frequency discriminator is modified to: $\hat{f}=f_n+2\mathrm{\&alpha;\&Delta;f}.$

The key that the present invention reduces computation complexity is to calculate E _acctime (f), the quantification of the position angle under polar coordinate system and dynamic programming algorithm thought is adopted to realize the minimizing of calculated amount.

As shown in Figure 4, for calculating the cumulative rear energy value E of maximum coherence _tflow process, specifically comprise following step S41 ~ S44:

(S41) according to the sampled value f of certain tracking frequency offset concrete _t, phase rotating is carried out to the original correlation x [k] of correlator output and obtains x _r[k], computing formula is as follows:

$x_r\left[k\right]={\mathrm{IP}}_r\left[k\right]+{\mathrm{QP}}_r\left[k\right]\&CenterDot;1j=x\left[k\right]e^{-j2\mathrm{\&pi;}{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 obtain:

$E_t=\underset{D}{\max }\{S_i{\left(D\right)}^2+S_q{\left(D\right)}^2\}$

$S_i\left(D\right)=\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{IP}}_r\left[k\right]d_k,S_q\left(D\right)=\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{QP}}_r\left[k\right]d_k---\left(7\right)$

By x _r[k] regards bivector as, then d _kcan regard as x when getting-1 _r[k] revolves turnback, d _kget 1 corresponding non rotating.Therefore E is calculated _tgiven N number of direction can be regarded as and can have 180 degree of bivectors overturn, find out a kind of upset combination and make itself and vector { S _i, S _qthere is maximum energy value.

(S42) plane is defined in subset and recurrence relation.Due to by all d _kall negates (namely all d _ksymbol negate simultaneously, 1 becomes-1 ,-1 becomes 1) do not affect E _tvalue, therefore, by d ₁be fixed as 1.Before note, m turning vectorial sum is that is:

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

All the set that possible value is formed is yi Zhi, in only have an element x _rand meet following recurrence relation [1]:

(S43) recurrence calculation in there is the subset that maximum coherence accumulated energies value element forms in each quantized directions initialization step be:

1) arrange for empty set.

2) travel through in each element y, calculate two that a cumulative correlation vector newly produces afterwards according to the recurrence relation of formula (9) in element, wherein:

Y ₊corresponding to d _m+1the situation of=1, i.e. y ₊=y+x _r[m+1];

Y _-corresponding to d _m+1the situation of=-1, i.e. y _-=y-x _r[m+1].

3) y is judged successively ₊or y _-add in which way in, with y ₊for example illustrates determination methods (y _-estimate of situation identical).If y ₊according to the position angle after formula (4), formula (5) quantize be if in no quantization back bearing be element, then by y ₊add to in; Otherwise establish in an existing element z quantize back bearing and be also then compare y ₊with the modulus value size of z, if | y ₊| > | z|, then will in z y ₊replace, otherwise keep constant.

Obtain after iteration N-1 time

(S44) find out in there is the element y of maximum modulus value _m, then calculate E _t, namely get E _t=| y _m| ².

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

Say from the strict sense, Ying Cong the element that middle search has maximum modulus value just can obtain E _t, but adopt this kind of method to reduce calculated amount.Because the noise section in correlation is uncorrelated, accumulation result can be confined to one of two dimensional surface comparatively in zonule.Along with the increase of m in element number exponentially regular increase, for larger m, the distribution of middle element will present a large amount of feature of assembling.Make use of These characteristics just, the present invention when allow certain error, from in select a small amount of most probable and continue the decision-making that the cumulative element obtaining maximum energy value participates in next stage, a large amount of text combination branches is excluded in advance.Calculate iterative process both embodied the thought of above-mentioned minimizing calculated amount, exclude in during element in addition, have employed following geometrical principle and approximate processing:

Geometrical principle: for in the identical vector in two position angles as long as then from continue to add up and can not obtain the vector with maximum modulus value.Above-mentioned principle provable as follows: as shown in Figure 5, if from set out and continue the cumulative vector obtaining ceiling capacity, that is: the cumulative newly-increased vector of follow-up N-m step is after v, there is maximum modulus value.Due to also be a vector that possible add up out, therefore on the other hand, due to also be the vector obtained that can add up, therefore obtain thus be nonnegative number with the inner product of r, namely angle is therebetween less than or equal to 90 degree.In the case, can obviously find out from accompanying drawing 5, this is just derived contradiction.

Approximate processing: for two in vector, as long as its quantize after position angle identical, then approximately think that it is in the same direction, is about to in sector region be similar to ray in the middle of it.

According to computation process known, only have at most in the sector region that each is quantized into ray in an element, therefore middle element number is no more than M.Obvious calculating each iterative step calculated amount with middle element number is directly proportional, and the computation complexity therefore in the present invention is no more than O (M × N).

The practical application effect of the GEO signal that the present invention is directed to BDS will be verified below by numerical simulation.

The GEO satellite emulated for BDS carries out, and concrete optimum configurations is as follows:

Participate in text times N=60 that single frequency is differentiated, namely the coherent accumulation time is T=0.12s;

Position angle 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 11 sampled points are got at interval wherein.

With C++ programming realization text algorithm for estimating herein and frequency discriminator during emulation, at employing 64 bit CPU and dominant frequency be 2.4GHz server on run, by the calculated amount of quantization means algorithm working time.Obtained working time by statistics, single text is estimated (namely to calculate an E _t), the inventive method required time is adopted to be about 14ms, and traditional time needed for exhaustive search method is greater than 1500s, namely when context of methods is for searching for 60 bit text, calculated amount is less than 0.00001 times of traditional algorithm, thus make the signal for 500bps reach the coherent accumulation of 0.12s may real-time implementation.Algorithm needs to store the set that may have at most M element herein therefore required memory space is linear with M, and traditional algorithm only needs to store the current person that has ceiling capacity of having traveled through during text combines, and required storage is few.Reflect dynamic programming algorithm changes time complexity feature with space complexity thus.

The present invention realizes the tracking performance that long-time coherent accumulation brings and promotes, be reflected in below to whole FLL loop in the emulation of Static and dynamic two kinds of situations.Adopt conventional second-order F LL loop during emulation, adopt frequency discriminator of the present invention and traditional frequency discriminator (single frequency discrimination institute usage data length is also 0.12s) to obtain two track loop respectively.Doppler frequency initial value is-180Hz, emulated data total length 120s, and loop is started working from tenacious tracking state, and namely initial velocity and acceleration tracking error are 0Hz.Only consider carrier track, suppose nonexistent code tracing deviation.With loop noise bandwidth B _lthe product B of interval T is upgraded with loop _lt describes the principal feature of track loop.Be the quiescent conditions of 0 for doppler changing rate, adopt less loop bandwidth, according to B _lt=0.01 designs loop.Be 3Hz/s current intelligence for doppler changing rate, adopt larger loop bandwidth, according to B _lt=0.025 designs loop.Under quiescent conditions, for 20dB-Hz signal, the tracking error of two kinds of loops as shown in Figure 6; Under current intelligence, for 23dB-Hz signal, the tracking error of two kinds of loops as shown in Figure 7.Under quiescent conditions, when adopting new frequency discriminator, the standard deviation of carrier frequency tracking error is 2.18Hz, and when adopting traditional frequency discriminator, the standard deviation of carrier frequency tracking error is 14.54Hz.Under current intelligence, when adopting new frequency discriminator, the standard deviation of carrier frequency tracking error is 1.27Hz, and when adopting traditional frequency discriminator, the standard deviation of carrier frequency tracking error is 19.96Hz.In weak signal conditions, no matter with or without signal dynamics, the frequency tracking error of track algorithm is less than 5Hz herein, and the frequency tracking error of traditional track algorithm can reach 30 ~ 40Hz.Further emulation shows, be reach the tracking accuracy suitable with frequency discriminator of the present invention for traditional frequency discriminator, be respectively 27dB-Hz and 30dB-Hz for the carrier-to-noise ratio needed for dynamic and static state signal, namely frequency discriminator of the present invention makes the tracking sensitivity of FLL track loop improve about 7dB.

Below be only the preferred embodiment of the present invention, protection scope of the present invention be not only confined to above-described embodiment, all technical schemes belonged under thinking of the present invention all belong to protection scope of the present invention.It should be pointed out that for those skilled in the art, some improvements and modifications without departing from the principles of the present invention, should be considered as protection scope of the present invention.

Claims (3)

1. a low complex degree frequency discriminator method for designing for GNSS weak signal tracking, is characterized in that, comprise the following steps:

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

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

After given tracking frequency offset f, only search text combines the maximum coherence accumulated energies value obtained is E _acc(f), that is:

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

Wherein, k=1,2 ..., N, D=[d ₁, d ₂..., d _n], d _k{+1 ,-1} is the vector be made up of the text that each bit is to be estimated successively to ∈, and D represents the combined situation of N number of text bit, combines hereinafter referred to as text; F is tracking frequency offset; for the estimated result of tracking frequency offset, D _mlfor the estimated result of text combination;

(S2) the span [-f of tracking frequency offset f is set _max, f _max], and get L sampled point: f with equal interval delta f wherein ₁< f ₂< ... < f _l-1< f _l, L is integer;

(S3) by plane carry out quantification treatment according to position angle, be quantified as the equally distributed ray of M bar, M is integer; With for interval, continuum [0,2 π] gets M discrete sampling point: for any vector { S to be quantified _x, S _y, if its polar coordinate representation is:

From r=1,2 ..., M, in find and position angle immediate phase place as the position angle after quantification, the cumulative vector after namely quantizing for:

(S4) to the sampled point f of each tracking frequency offset value _t, t=1,2 ..., L, calculates E _t=E _accand therefrom find out the maximum E (ft) _n=max{E _t, and with the frequency f of correspondence _nas initial tracking frequency offset estimated value, n is integer;

(S5) tracking frequency offset estimated value is revised further by three point interpolation method: if E ₁or E _lfor time maximum, not carrying out interpolation processing, is directly f by frequency deviation estimated value amplitude limit ₁or f _l, otherwise, the coefficient of first computational representation peakdeviation amount again the output of frequency discriminator is modified to: $\hat{f}=f_n+2\mathrm{\&alpha;\&Delta;f}.$

2. the low complex degree frequency discriminator method for designing of a kind of GNSS weak signal tracking as claimed in claim 1, is characterized in that, according to the sampled value f of tracking frequency offset in described step (S4) _t, t=1,2 ..., L, solves E _tdetailed process be:

(S41) according to the sampled value f of tracking frequency offset _t, phase rotating is carried out to the original correlation x [k] of correlator output and obtains x _r[k], computing formula is as follows:

$x_r\left[k\right]={\mathrm{IP}}_r\left[k\right]+Q{P}_r\left[k\right]\&CenterDot;1j=x\left[k\right]e^{-j2\mathrm{\&pi;}{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 obtain:

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

By x _r[k] regards bivector as, then d _kcan regard as x when getting-1 _r[k] revolves turnback, d _kget 1 corresponding non rotating;

(S42) plane is defined in subset and recurrence relation: by d ₁be fixed as 1, before note, m the vectorial sum of upset is that is:

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

All the set that possible value is formed is yi Zhi, in only have an element x _rand recurrence relation is as follows [1]:

(S43) recurrence calculation in there is the subset that maximum coherence accumulated energies value element forms in each quantized directions initialization then m+1=2,3 ..., the step of N is:

1) arrange for empty set;

2) travel through in each element y, calculate two that a cumulative correlation vector newly produces afterwards according to the recurrence relation of formula (9) in element, wherein:

Y+ corresponds to d _m+1the situation of=1, i.e. y+=y+x _r[m+1];

Y-corresponds to d _m+1the situation of=-1, i.e. y-=y-x _r[m+1];

3) judge that y+ or y-adds in which way successively in, specifically for y+, determination methods is described, if y+ according to the position angle after formula (4), formula (5) quantification is if in no quantization back bearing be element, then y+ is added to in; Otherwise establish in an existing element z quantize back bearing and be also then compare the modulus value size of y+ and z, if | y+| > | z|, then will in z y+ replace, otherwise keep constant;

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

(S44) find out in there is the element y of maximum modulus value _m, then calculate E _t, namely get E _t=| y _m| ².

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