Disclosure of Invention
The invention aims to provide a high-precision iteration estimation method for the frequency of a sine wave coherent pulse train signal based on the combination of an FFT algorithm and phase straight line fitting, which can be used for high-precision frequency estimation of the sine wave coherent pulse train signal with arbitrarily changed repetition frequency.
The technical solution for realizing the purpose of the invention is as follows: a sine wave coherent pulse train signal frequency estimation method based on phase straight line fitting comprises the following steps:
step 1, for input complex analysis digitalFFT operation is carried out on a plurality of pulses of the signal, the frequency value corresponding to the strongest spectral line is taken as the primary estimated value of each pulse frequency, and the average value of the primary estimated values of each pulse frequency is taken as f e ;
Step 2, using the average value f of the primary estimation values of each pulse frequency e Performing down-conversion processing on the input complex analysis digital signal;
step 3, extracting phase angle information of each pulse from the complex analysis digital signal after the down-conversion treatment;
step 4, preprocessing the obtained phase angle data, and adjusting an initial phase angle;
step 5, extracting the phase mean value and the center moment of each pulse;
step 6, performing linear fitting operation on the phase average value of each pulse, and calculating a frequency estimation error delta f according to the result of linear approximation e ;
Step 7, Δf e Comparing with a set threshold: when Deltaf e When the value is smaller than the set threshold value, directly outputting f e As the final estimated value of the frequency, otherwise, go to step 8;
step 8, frequency estimation value f e Adjusted to f e =f e +Δf e By Δf e And (3) performing down-conversion treatment on the input complex analysis digital signal subjected to down-conversion and initial phase adjustment, and then turning to step (3).
Further, in step 1, the plurality of pulses of the input complex analysis digital signal are respectively subjected to FFT operation, the frequency value corresponding to the strongest spectral line is used as the initial estimation value of each pulse frequency, and the average value of the initial estimation values of each pulse frequency is denoted as f e The method is characterized by comprising the following steps:
recording deviceThe method comprises the steps of analyzing sine wave coherent pulse train signals for noise-containing complex to be processed, wherein T is more than or equal to 0 and less than or equal to T; w (t) is observed white noise with a mean value of zero; f (f) 0 、θ、T、T i 、τ i And P is the signal frequency, initial phase, duration, and center of the ith pulse, respectivelyTime, duration of the ith pulse and number of pulses; />Is of duration tau i The time center is T i Is a rectangular pulse of (2); f (f) 0 The measurement is to be estimated; at the same time, signal s is recorded 0 The sampling frequency and the sampling point number of (t) are f respectively s N, the signal sequence obtained after sampling is s 0 (n),n=1,2,3,…,N-1;
Performing M-point FFT operation on the ith pulse signal, wherein M is greater than the sampling point number of the ith pulse signal, and the frequency value corresponding to the strongest spectral lineAll +.>Then, the average value of the primary estimation values of the pulse frequencies is denoted as f e Then->
Further, step 2 is described as using the average value f of the primary estimation values of each pulse frequency e The down-conversion processing is carried out on the input complex analysis digital signal, and the down-conversion processing is concretely as follows:
by f e For a sequence of sampled signals s 0 (n) down-conversion processing, in which the portions other than the pulse durations are set to zero in the subsequent signal processing, i.eObtain a down-converted sample signal sequence +.>n=1,2,3,…,N-1。
Further, in step 3, the phase angle information of each pulse is extracted from the complex analysis digital signal after the down-conversion processing, which specifically includes the following steps:
computing a data sequencePhase angle information θ of each pulse signal in (a) i (n),n=0,1,2,…,Q i -1, i=1, 2,3, …, P; wherein Q is i The number of sampling points for the ith pulse; the sampling data sequence of the ith pulse is marked as +.>Then->-π≤θ i (n)≤π,n=0,1,2,…,Q i -1。
Further, the step 4 is to preprocess the obtained phase angle data, adjust the initial phase angle, and specifically includes the following steps:
1) Firstly, judging whether the percentage of the phase folding phenomenon occurring in the pulse exceeds a threshold value or not according to the extracted pulse phase data sequences: if the phase folding phenomenon occurs in the pulse phase data of more than 75% of the pulses, the signal initial phase angle is adjusted, the phase data sequence of each pulse is re-extracted, and the specific adjustment method is as follows: order then=0, 1,2, …, N-1; then, extracting the phase data sequence of each pulse; if the pulse number of the phase folding phenomenon of the phase data in the pulse is less than 75%, performing the phase folding operation on the phase data in the pulse with the phase folding phenomenon;
2) The method for eliminating the phase data of which the phase estimation value does not meet the requirement in each pulse comprises the following specific steps: first calculate psi i (n)=θ i (n+1)-θ i (n), n=0, 1,2, …, Q-1, and then solving for the mean value and standard deviation of ψ (n), when the difference between ψ (n) and its mean value is greater than 3 times its standard deviation, either ψ (n-1) or ψ (n+1) is substituted for ψ (n).
Further, the step 5 of extracting the phase mean value and the center moment of each pulse is specifically as follows;
calculating the central moment t of each pulse after phase preprocessing ci And the average value theta of each pulse phase data sequence i Obtaining a new phase data sequence theta i (t ci ),i=1,2,…,P。
Further, in step 6, a line fitting operation is performed on the phase average value of each pulse, and then a frequency estimation error Δf is calculated according to the result of the line fitting e The method is characterized by comprising the following steps:
for phase data sequence theta i (t ci ) Performing linear fitting operation, wherein i=p, p+1, p+2, …, L-1, p and L are index numbers of a start pulse and a final pulse corresponding to a phase data segment with no phase folding at maximum respectively;
assume that the resulting fitted line is Θ (t ci )=2πΔf e t ci +θ 0 I=p, p+1, p+2, …, L-1, Δf e I.e. the frequency estimation error.
Further, step 8 is described as Δf e And (3) performing down-conversion treatment on the down-converted and initially-phase-adjusted input complex analysis digital signal again, wherein the formula is as follows:
where n=0, 1,2, …, N-1.
Compared with the prior art, the invention has the remarkable advantages that: (1) The necessary preprocessing operation is carried out on the intra-pulse phase data and inter-pulse phase data, so that the frequency estimation accuracy is better; (2) When phase straight line fitting is performed by utilizing intra-pulse phase data, no constraint condition exists on inter-pulse intervals; (3) The method can be used for high-precision frequency estimation of fixed-repetition-frequency sine-wave coherent pulse train signals, and also can be used for high-precision frequency estimation of random-repetition-frequency sine-wave coherent pulse train signals, and the application range is wider.
Detailed Description
The invention relates to a sine wave coherent pulse train signal frequency estimation method based on phase straight line fitting, which comprises the following steps:
step 1, FFT operation is respectively carried out on a plurality of pulses of an input complex analysis digital signal, a frequency value corresponding to the strongest spectral line is used as an initial estimation value of each pulse frequency, and an average value of the initial estimation values of each pulse frequency is recorded as f e ;
Step 2, using the average value f of the primary estimation values of each pulse frequency e Performing down-conversion processing on the input complex analysis digital signal;
step 3, extracting phase angle information of each pulse from the complex analysis digital signal after the down-conversion treatment;
step 4, preprocessing the obtained phase angle data, and adjusting an initial phase angle;
step 5, extracting the phase mean value and the center moment of each pulse;
step 6, performing linear fitting operation on the phase average value of each pulse, and calculating a frequency estimation error delta f according to the result of linear approximation e ;
Step 7, Δf e Comparing with a set threshold (the threshold is determined by the required estimation accuracy): when Deltaf e When the value is smaller than the set threshold value, directly outputting f e As the final estimated value of the frequency, otherwise, go to step 8;
step 8, frequency estimation value f e Adjusted to f e =f e +Δf e By Δf e And (3) performing down-conversion treatment on the input complex analysis digital signal subjected to down-conversion and initial phase adjustment, and then turning to step (3).
Further, in step 1, the plurality of pulses of the input complex analysis digital signal are respectively subjected to FFT operation, the frequency value corresponding to the strongest spectral line is used as the initial estimation value of each pulse frequency, and the average value of the initial estimation values of each pulse frequency is denoted as f e The method is characterized by comprising the following steps:
recording deviceThe method comprises the steps of analyzing sine wave coherent pulse train signals for noise-containing complex to be processed, wherein T is more than or equal to 0 and less than or equal to T; w (t) is observed white noise with a mean value of zero; f (f) 0 、θ、T、T i 、τ i And P is the signal frequency, initial phase, duration, central time of the ith pulse, duration of the ith pulse and pulse number respectively; i=1, 2,3, …, P;is of duration tau i The time center is T i Is a rectangular pulse of (2); f (f) 0 The measurement is to be estimated; at the same time, signal s is recorded 0 The sampling frequency and the sampling point number of (t) are f respectively s N, the signal sequence obtained after sampling is s 0 (n),n=1,2,3,…,N-1;
Performing M-point FFT operation on the ith pulse signal (M should be greater than the sampling point of the ith pulse signal, but to increase the frequency estimation accuracy as much as possible, and to avoid the inaccuracy of the estimation of each pulse starting and stopping time under the condition of low input signal-to-noise ratio and the introduction of pure receiving noise signal, more conservative pulse starting should be adopted)Time-of-day estimation to ensure that no pure received noise signal portion is introduced), the frequency value corresponding to the strongest spectral lineAll of the pulse frequencies are obtained as initial estimation valuesThen, the average value of the primary estimation values of the pulse frequencies is denoted as f e Then->When the input signal-to-noise ratio is high, for example, greater than 10dB, only the frequency estimation value of a small number of pulses is calculated to reduce the operation amount, and then the average value is obtained.
Further, step 2 is described as using the average value f of the primary estimation values of each pulse frequency e The down-conversion processing is carried out on the input complex analysis digital signal, and the down-conversion processing is concretely as follows:
by f e For a sequence of sampled signals s 0 (n) down-conversion processing, in which the portions other than the pulse durations are set to zero in the subsequent signal processing, i.eObtain a down-converted sample signal sequence +.>n=1,2,3,…,N-1。
Further, in step 3, the phase angle information of each pulse is extracted from the complex analysis digital signal after the down-conversion processing, which specifically includes the following steps:
computing a data sequencePhase angle information θ of each pulse signal in (a) i (n),n=0,1,2,…,Q i -1, i=1, 2,3, …, P; wherein Q is i The number of sampling points for the ith pulse; record the ithThe sampled data sequence of each pulse isThen->-π≤θ i (n)≤π,n=0,1,2,…,Q i -1。
Further, the step 4 is to preprocess the obtained phase angle data, adjust the initial phase angle, and specifically includes the following steps:
1) Firstly, judging whether the percentage of the phase folding phenomenon occurring in the pulse exceeds a threshold value or not according to the extracted pulse phase data sequences: if the phase folding phenomenon occurs in the pulse phase data of more than 75% of the pulses, the signal initial phase angle is adjusted, the phase data sequence of each pulse is re-extracted, and the specific adjustment method is as follows: order then=0, 1,2, …, N-1; then, extracting the phase data sequence of each pulse; if the pulse number of the phase folding phenomenon of the phase data in the pulse is less than 75%, performing the phase folding operation on the phase data in the pulse with the phase folding phenomenon;
2) The method for eliminating the phase data of which the phase estimation value does not meet the requirement in each pulse comprises the following specific steps: first calculate psi i (n)=θ i (n+1)-θ i (n), n=0, 1,2, …, Q-1, and then solving for the mean value and standard deviation of ψ (n), when the difference between ψ (n) and its mean value is greater than 3 times its standard deviation, either ψ (n-1) or ψ (n+1) is substituted for ψ (n).
Further, the step 5 of extracting the phase mean value and the center moment of each pulse is specifically as follows;
calculating the central moment t of each pulse after phase preprocessing ci And the average value theta of each pulse phase data sequence i Obtaining a new phase data sequence theta i (t ci ),i=1,2,…,P。
Further, in step 6, a straight line fitting operation is performed on the phase average value of each pulse, and thenCalculating frequency estimation error Deltaf according to the result of straight line similarity e The method is characterized by comprising the following steps:
for phase data sequence theta i (t ci ) Performing linear fitting operation, wherein i=p, p+1, p+2, …, L-1, p and L are index numbers of a start pulse and a final pulse corresponding to a phase data segment with no phase folding at maximum respectively;
assume that the resulting fitted line is Θ (t ci )=2πΔf e t ci +θ 0 I=p, p+1, p+2, …, L-1, Δf e I.e. the frequency estimation error.
The phase data sequence theta obtained in the step 6 can also be used for the following steps before the step 6 is carried out i (t ci ) (i=1, 2, …, P) performs a linearization operation, i.e., an unfolding operation, so that more pulse phase data participates in the line fitting operation.
Further, step 8 is described as Δf e And (3) performing down-conversion treatment on the down-converted and initially-phase-adjusted input complex analysis digital signal again, wherein the formula is as follows:
where n=0, 1,2, …, N-1.
Compared with a two-stage FFT algorithm, the frequency estimation accuracy is better although the calculation complexity is increased, and the limitation of the two-stage FFT algorithm is overcome, so that the method can be used for high-accuracy frequency estimation of a fixed-repetition-frequency sine-wave coherent pulse train signal and high-accuracy frequency estimation of a sine-wave coherent pulse train signal with random change of the repetition frequency, and the application range is wider.
The invention is described in further detail below with reference to the accompanying drawings and specific examples.
Examples
Fig. 1 shows specific processing steps of a single sine wave coherent burst signal frequency estimation method employed in the present invention. The noise-containing complex analysis sine type coherent pulse train signal to be processed in the invention can be expressed as a form shown in (1).
In the formula (1), w (t) is observed white noise with a mean value of zero; f (f) 0 、θ、T、T i 、τ i And P is the signal frequency, initial phase, duration, central time of the ith pulse, duration of the ith pulse and pulse number respectively;is of duration tau i The time center is T i Is a rectangular pulse of (2); f (f) 0 Is the quantity to be estimated. At the same time, signal s is recorded 0 The sampling frequency and the sampling point number of (t) are f respectively s N, the signal sequence obtained after sampling is s 0 (N) (n=1, 2,3, …, N-1). The phase of the signal sequence can be expressed in the form shown in formula (2) (although the invention processes the sampled digital signal, the correspondence between signal phase and time is represented by the continuous time variable t, which does not affect the description of the invention).
In the formula (2), the amino acid sequence of the compound,is the phase noise generated by w (t). When w (t) is random noise with mean value zero, ++>Also random noise with zero mean. Theoretically, ignore +.>Phi (t) varies linearly with time t, both during the respective pulse duration and during the whole observation time. But in the opposite directionWhen phi (t) is calculated numerically, phi (t) can only take values within the (-pi, pi) range, so that phi (t) is in fact periodically and linearly changed in the whole observation time unless the sine wave coherent pulse train signal to be processed is changed in only one period in the whole observation time. When->Cannot be ignored unless +.>Is strong enough to drown out phi (t), otherwise such a law of variation of phi (t) should still be observed.
Performing M-point FFT operation on the ith (i=1, 2,3, …, P) pulse signal (M should be greater than the sampling point of the ith pulse signal, but should increase the value of M as much as possible to improve the frequency estimation accuracy; to avoid the inaccurate estimation of the starting and ending moments of each pulse under the condition of lower input signal-to-noise ratio and the introduction of pure receiving noise signals, more conservative estimation of the starting and ending moments of the pulse should be adopted to ensure that the part of the pure receiving noise signals is not introduced), obtaining the frequency value corresponding to the strongest spectral line, and recording asAll are obtainedThen, obtain the coarse frequency estimation value +.>
By f e For signal sequence s 0 (N) (n=1, 2,3, …, N-1) after the down-conversion processing, the extracted signal phase information can be expressed in the form of expression (3).
In formula (3), Δf=f 0 -f e . Due to frequency estimationThe gauge error Δf is typically much smaller than f 0 So, the extracted signal phase information may no longer be periodically varying, at least the period of variation will increase significantly.
Subject to phase noiseAnd the influence of the initial value theta of the phase on the signal sequence s 0 (N) (n=1, 2,3, …, N-1) there is a possibility that the intra-pulse phase information extracted after the down-conversion processing has out-of-tolerance or phase-folded data. The necessary preprocessing of the extracted phase data is a necessary condition for obtaining a high-precision frequency estimate.
FIG. 2 shows the embodiment in using f e For signal sequence s 0 (N) (n=1, 2,3, …, N-1) after the down-conversion processing (the down-converted signal sequence is recorded) the extracted pulse phase data sequences. As can be seen from fig. 2, the extracted pulse phase data sequences all have a phase folding phenomenon. In this case, the down-converted signal sequence is multiplied by a factor e jπ Then, the phase data sequence of each pulse is re-extracted. FIG. 3 shows the embodiment in using f e For signal sequence s 0 (N) (n=1, 2,3, …, N-1) each pulse phase data sequence extracted after the down-conversion processing. As can be seen from fig. 3, the phase data sequences of the second half pulse all have phase folding phenomenon, while the phase data sequences of the first half pulse have no phase folding phenomenon, but have a problem that the individual phase estimation values are significantly out of tolerance. In this case, for the problem that only the individual phase estimates are significantly out of tolerance within the pulse, the estimate may be eliminated from participating in the subsequent calculation process. For the data with phase folding and individual out-of-tolerance in pulse, the phase folding should be firstly released, and then the individual out-of-tolerance data should be removed.
After preprocessing the extracted intra-pulse phase data, calculating an intra-pulse phase average value of each pulse to obtain an intra-pulse phase average value data sequence
In the formula (4), the amino acid sequence of the compound,is the noise part in the phase mean, +.>FIG. 4 is an intra-pulse phase average sequence extracted after preprocessing the intra-pulse phase data sequence shown in FIG. 3 according to the present invention. As can be seen from fig. 4, this phase-averaged sequence has a phase folding phenomenon.
In straight line fitting the pulse phase averaging sequence shown in fig. 4, two methods can be employed: firstly, performing an unfolding phase folding operation on an intra-pulse phase mean value sequence with a phase folding phenomenon, and then performing straight line fitting; and secondly, taking the intra-pulse phase average value sequence with the maximum phase folding, and performing straight line fitting on the first 11 intra-pulse phase average value sequences in fig. 4.
Performing straight line fitting on the intra-pulse phase average value sequence to obtain an estimated value delta f of the frequency difference e After that, if Δf e Less than or equal to a given threshold value, the frequency estimate f is output e The frequency estimation process ends. When Deltaf e When the frequency is greater than a given threshold value, the frequency estimation value is adjusted according to the formula (5), and delta f is utilized e For down-converted signal sequence(n=0, 1,2, …, N-1) performing down-conversion again, and then sequentially performing extraction, pretreatment, intra-pulse phase average value sequence extraction and frequency difference re-iteration estimation on each pulse phase data sequence until the frequency difference estimated value is smaller than a given threshold. Generally, the frequency estimation process can be ended through two iterative operations of the frequency error estimation value.
Fig. 5 and 6 are graphs showing the frequency offset estimation Δf obtained by performing a phase straight line fitting using the phase average data of the first 11 pulses shown in fig. 4 e (accordingly, f e =f e +Δf e ) The estimated value is used for down-converting the signal sequence(n=0, 1,2, …, N-1) the pulse phase data sequences extracted after the re-down conversion and intra-pulse phase mean value sequences.
FIG. 7 shows comparison diagrams of estimation performance simulation results of the method and the secondary FFT algorithm under several input signal-to-noise ratio conditions (signal sampling frequency, coherent pulse train duration, pulse repetition frequency and pulse width are respectively 100MHz, 1ms, 20kHz and 2 μs; the frequency of a signal to be estimated is uniformly valued within the range of 5 MHz+/-1.05 kHz; the number of simulation times under the condition of single input signal-to-noise ratio is 500), the same frequency coarse estimation value is used in the method and the secondary FFT algorithm, the digital resolution of the FFT algorithm used is 7.3Hz when the frequency error estimation is carried out in the secondary FFT algorithm, and the starting and ending time of each pulse is assumed to be known when the simulation calculation is carried out. As can be seen from fig. 7, in a given input signal-to-noise ratio range, the estimation performance of the frequency estimation method provided by the present invention is better than that of the two-stage FFT algorithm; as the input signal-to-noise ratio decreases, the estimation performance of the two frequency estimation methods gets closer and closer.
The invention carries out necessary preprocessing operation on the intra-pulse and inter-pulse phase data, and has better frequency estimation precision; when phase straight line fitting is performed by utilizing intra-pulse phase data, no constraint condition exists on inter-pulse intervals; the method can be used for high-precision frequency estimation of fixed-repetition-frequency sine-wave coherent pulse train signals, and also can be used for high-precision frequency estimation of random-repetition-frequency sine-wave coherent pulse train signals, and the application range is wider.