CN104535959A - Signal frequency and DOA joint measurement method and device under spatial-temporal sub-nyquist sampling - Google Patents

Abstract

The invention discloses a signal frequency and DOA joint measurement method and device under spatial-temporal sub-nyquist sampling. The method comprises the steps that DFT is carried out on L paths of signal samples, frequency and phase position correction is carried out on each path of DFT spectrum peak through a Candan interpolation estimator, and L pairs of normalization frequencies and phase position estimation values are obtained; by adoption of a closed type CRT algorithm, a signal frequency value is obtained through the L pairs of normalization frequencies; a receiving signal phase position difference estimated value of sensor array elements 1 and other array elements is obtained, then a phase position remainder is formed, and finally the spatial signal incidence angle is calculated. The device comprises the sensor array elements, A/D samplers, a DSP, an output driver and a display circuit of the output driver, wherein the spatial far-field narrowband signals reach the sensor array elements by a certain incidence angle theta, and array receiving signals are obtained; the A/D samplers on the sensor array elements sample the signals and input the obtained sample sequences into the DSP in parallel, frequency estimation and DOA estimation of the incident signals are obtained after processing by an internal algorithm of the DSP, and finally the estimation results are displayed through the output driver and the display circuit of the output driver.

Description

Signal frequency and DOA union measuring method and device under space-time lack sampling

Technical field

The present invention relates to digital processing field, particularly relate to one and utilize sensor array to carry out time domain, space lack sampling to space incident signal, by carrying out method and the device of high-precision frequency and DOA combined measurement to sample analysis.

Background technology

Array parameter is estimated as the important component part of Array Signal Processing, and rapidly, its application has related to the numerous areas such as sonar, radar, communication and biomedicine for very active in recent years and development.Wherein, the Frequency Estimation of space incident signal and direction of arrival (Direction of Arrival, DOA) estimate it is hot research content.And more go deep into extensive along with what apply, the requirement of people to parameter estimation is more and more higher, such as require that arithmetic accuracy is higher, fault-tolerance and robustness stronger, calculate simpler, hardware layout cost is lower, these require the development also having promoted array parameter estimation technique.Wherein, how under the condition of space-time lack sampling, the attention in obtaining academia in recent years is accurately measured to signal frequency and DOA parameter.

Find through literature survey, in classical way, such prerequisite is needed to the frequency of space incident signal, the isoparametric measurement of DOA: all will meet nyquist sampling condition to signal in the sampling of time domain and spatial domain, namely require that sampling rate is greater than the twice of incoming signal frequency, and the arrangement pitch of adjacent sensors array element is not more than the half-wavelength of incoming signal.But above condition is difficult to meet in some actual environments.Such as, document [1] is pointed out, the operating frequency range of military airborne radar investigation receiver is 2-18GHz, and it is unpractical for directly carrying out time domain nyquist sampling to the signal be distributed in broadband like this.Although there is the analog to digital converter of sampling rate at the GHz order of magnitude, its price is very expensive, and is also difficult to reach so large data throughput with the back-end digital processor of system support.Depend merely on its effect of data acquisition performance improving hardware device very limited, only propose new spectrum analysis theory method in signal transacting field and could fundamentally solve this kind of problem.Equally, the layout of so wide working band to sensor array it is also proposed rigors.At classics without in fuzzy DOA measuring method, require that the spacing of adjacent sensors array element on array is not more than the half-wavelength of incoming signal.When signal wavelength is less (half-wavelength as corresponding in 18GHz signal is 0.83cm), the difficulty that so short spacing not only can cause sensor array element installing, also will certainly strengthen the mutual coupling of low frequency signal between different sensors array element.Therefore, widen sensor array element distance, realize the sparse layout of array, under the condition that space lack sampling is carried out to signal, realize high precision DOA measure significant in actual applications.

But space-time lack sampling can cause frequency, the phase ambiguity of signal, for this problem, ancient Chinese remainder theorem (Chinese Remainder Theorem, CRT) [2] can be introduced this field and is solved.From CRT, if a certain positive integer is less than the lowest common multiple of one group of relatively prime between two integer, then this positive integer uniquely can be determined by the remainder after its modulo operation.But CRT is very responsive to Residue error, if remainder exists minimum error, larger reconstructed error can be caused.At radio interference location (Radio Interferometric Positioning, RIPS [3]), inevitably there is noise in the system such as radar [4], this must bring Residue error to CRT, and therefore the range of application of traditional C RT is very limited.Document [5] proposes a kind of CRT algorithm of robust, and be no more than the condition of certain limit at Residue error under, this algorithm can more accurately reconstruct unknown integer, but this algorithm just can must realize by two-dimensional search, and operand is large.Document improves the algorithm of document [5] in [6,7], and by two-dimensional search being reduced to linear search thus reducing operand, but when the number of institute's delivery is more or numerical value is larger, its search complexity is still very high.For fundamentally avoiding heavy search procedure, document [8] proposes a kind of CRT enclosed method for solving of robust, greatly reduces operand, is conducive to real-time implementation in engineering, thus for the invention provides suitable theoretical tool.

Have benefited from the improvement to classical CRT and optimization [5-8], at signal frequency estimation, DOA, CRT estimates that fields such as [9] obtains application in recent years.In signal frequency estimation, if document [10-12] is by carrying out the sampling of multi-path low speed rate to signal, achieve the Accurate Reconstruction of signal frequency under time domain lack sampling in conjunction with CRT theorem.But in these methods, CRT remainder values is directly obtained by DFT spectrum peak search, when the leakage of existence spectrum can bring error to remainder thus affect estimated accuracy.DFT counts and is set to sampling rate value corresponding to each road (frequency resolution of corresponding DFT is 1Hz) by document [10-12], and supposition measured signal frequency values is integer (being namely the integral multiple of DFT frequency resolution just), thus avoid the observation Residue error brought because spectrum is leaked.Obviously integer-valued hypothesis is done to signal frequency and limit its range of application, and due to the value (usually larger) that is set to sampling rate and bring larger operand of counting by DFT; In addition, owing to using the CRT method based on search in document [10-12], when measured signal frequency is higher, the large and length consuming time of computational complexity, is unfavorable for real-time measurement.In DOA estimation, document [9] realizes the sparse layout of array element by the form extracting the grouping of sensor array element from imaginary even linear array, utilizes the phase differential structure remainder of sensor array tuple Received signal strength, carrys out the incident angle of estimated signal in conjunction with CRT.But, do not provide concrete sampling plan and the measuring method of signal phase in document [9]; Document [13] proposes CRT to be applied to signal frequency under space-time lack sampling and DOA Combined estimator, the method first utilizes separately certain sensor array element to carry out a point frequency multi-channel sampling to carry out estimated signal frequency, and then utilizes array grouping multiple repairing weld to estimate DOA.Because its sampling of the method in document [13] and parameter estimation procedure complete stage by stage, consuming time longer and sample utilization factor is low, its actual application value is restricted.At present, under space-time lack sampling, real-time combined measurement is carried out to the frequency of signal and DOA, still belong to technological gap.

Summary of the invention

The invention provides signal frequency and DOA union measuring method and device under a kind of space-time lack sampling, present invention achieves and the high precision of signal frequency to be estimated under time domain lack sampling and DOA high precision under spatial domain lack sampling is estimated, described below:

Signal frequency and DOA union measuring method under a kind of space-time lack sampling, described measuring method comprises the following steps:

Arrange the non-uniform linear arrays containing L sensor array element, with sensor array element 1 for reference array element, the spacing of definition sensor array element i and sensor array element 1 is d _{i, 1};

Interval T at one time ₀in, L sensor array element is respectively with sampling rate F ₁~ F _lparallel lack sampling is done to space incident signal, wherein, F ₁~ F _lthere is positive common divisor M _f;

DFT is done to L road sample of signal, with Candan interpolation estimator, frequency, phase correction is done to each road DFT spectrum peak, obtain L to normalized frequency and phase estimation value;

Adopt enclosed CRT algorithm, by L, signal frequency value is obtained to normalized frequency

Obtain sensor array element 1 estimated value poor with the phase of received signal of other array element, and form phase place remainder further, calculate spacing wave incident angle.

Described employing enclosed CRT algorithm, obtains signal frequency value by L to normalized frequency step be specially:

? as frequency remainder, using the value of each sensor array element sampling rate as the mould needed for CRT, remainder group prime number group Γ ₁... Γ _i... Γ _land M _fsubstitute into enclosed CRT and try to achieve signal frequency value

The described phase of received signal difference estimated value obtaining sensor array element 1 and other array element, and form phase place remainder further, calculate spacing wave incident angle and be specially:

Obtain sensor array element 1 estimated value poor with the phase of received signal of other array element and form phase place remainder further; By phase place remainder group, prime number group η ₁~ η _l-1and M _θthe CRT substituting into closed form reconstructs nonnegative constant N ₀, and then by N ₀calculate spacing wave incident angle

Under described space-time lack sampling, signal frequency and DOA combined measurement device comprise: sensor array element, A/D sampling thief, DSP and output driving and display circuit thereof,

Space far-field narrow band signal s (t) arrives each sensor array element with a certain incidence angle θ, obtains array received signal; A/D sampling thief in each sensor array element is respectively with F ₁, F ₂..., F _lspeed signal parallel is sampled, be input to DSP device by parallel for the sample sequence obtained, through the internal algorithm process of DSP, the Frequency Estimation and the DOA that obtain incoming signal estimate, finally drive by output and display circuit display estimated result.

Signal frequency under the space-time lack sampling that the present invention proposes and DOA union measuring method, if be applied to Practical Project field, can produce following beneficial effect:

The first, under achieving time domain lack sampling, the high precision of signal frequency is estimated;

In traditional frequency measurement method, require that sampling rate is higher than the twice of measured signal frequency, otherwise will frequency ambiguity be caused thus correct frequency estimation cannot be obtained.And the present invention adopts the scheme of multi-path low speed rate time domain lack sampling, achieve the measurement to high-frequency signal, considerably increase the measurement range of frequency.To the low sampling rate F in L road ₁~ F _l, the signal frequency range utilizing the present invention accurately to measure is (0, f _max], wherein f _maxequal F ₁~ F _llowest common multiple.

Such as, in experiment, the sampling rate of sensor array element is respectively F ₁=52400Hz, F ₂=54800Hz, F ₃=55600Hz, with the measurement range of classical DFT method then frequency be (0,27800], according to the present invention then maximum detection frequency be f _max=997853200Hz, improves 4 orders of magnitude than classic method.

The second, the DOA high precision achieved under spatial domain lack sampling is estimated, and is achieved the flexible arrangement of sensor array element;

In tradition without in fuzzy DOA measuring method, require that the spacing of adjacent sensors array element is less than the half-wavelength of incoming signal, so not only can cause the difficulty that sensor array element is being installed, also will certainly strengthen the mutual coupling of low frequency signal between different sensors array element, reduce the precision of each estimating DOA forsingals.The present invention proposes the sparse arrangement of concrete array, and can according to actual conditions flexible configuration sensor array element distance.

Such as, the incoming signal wavelength supposed in experiment is λ ₀=0.3426m, if according to tradition without fuzzy DOA estimation method, then adjacent sensors array element distance d need meet d < λ ₀/ 2=0.1713m, and according to the array that the present invention arranges, the minimum spacing of its adjacent sensors array element is 0.3m, breaches the restriction of half-wavelength.And by changing the value of systematic parameter, the flexible configuration of sensor array element distance can be realized.

Three, make full use of array sample, improve efficiency and the real-time of parametric joint measurement.

Its sampling of the parametric joint estimation technique of document [13] and measuring process are divided into two stages, namely first carry out frequeney division multiple (FDM) sampling by certain sensor array element independent and carry out measuring-signal frequency, and then utilize sensor array parallel sampling to carry out measuring-signal DOA, longer being unfavorable for consuming time of sampling measures in real time, and the data of two sample phase can not be multiplexing, sample utilization factor is low.The present invention directly utilizes sensor array to carry out parallel duplex sampling, then must obtain frequency and the DOA measured value of signal, being more conducive to engineer applied by analyzing a few sample collected.

Accompanying drawing explanation

Fig. 1 is frequency and the DOA combined measurement process flow diagram of signal under space-time lack sampling;

Fig. 2 is that sensor array element arranges schematic diagram;

Fig. 3 (a) is signal to noise ratio (S/N ratio) and signal frequency detection probability relation curve;

Fig. 3 (b) estimates square error relation curve for signal to noise ratio (S/N ratio) and signal frequency;

Fig. 4 is signal to noise ratio (S/N ratio) and signal DOA estimation property relationship curve;

Fig. 5 is hardware implementation figure of the present invention;

Fig. 6 is DSP internal processes flow graph.

Embodiment

For making the object, technical solutions and advantages of the present invention clearly, below embodiment of the present invention is described further in detail.

The present invention is in conjunction with the enclosed CRT algorithm of robust ^[8]frequency and the DOA combined measurement of signal under space-time lack sampling is carried out with DFT spectrum correction method.Owing to introducing enclosed CRT algorithm, avoid the search in calculating process thus to decrease computing consuming time.Compared with the method for document [13], present invention eliminates the step of first carrying out multiple repairing weld frequency measurement by a certain separated sensor array element, consuming time shorter and sample utilization factor is higher.Different from document [10-12], DFT when processing sample of signal in the present invention counts and elects number of samples as but not sampling rate value, greatly reduces calculated amount and the requirement to hardware.Meanwhile, be the Residue error that " fence effect " that overcome DFT causes, the present invention introduces Candan interpolation estimator ^[14,15]frequency spectrum and phase correction are carried out to DFT result, improves Parameter Estimation Precision.Therefore, the present invention has higher practical value and application prospect widely.

101: arrange the non-uniform linear arrays containing L sensor array element, with sensor array element 1 for reference array element, the spacing of definition sensor array element i and sensor array element 1 is d _{i, 1}and meet d _{i, 1}=K η ₁η ₂... η _l-1/ η _i-1, i=2 ..., L, wherein η ₁~ η _l-1for array relatively prime between two, K be according to actual conditions setting constant and be required to meet 0 < K < λ _min(λ _minfor system maximum detection frequency f _maxcorresponding signal wavelength).

To comprise the linear array of 3 sensor array elements, make K=0.9 λ _min, η ₁=3, η ₂=2, then d _2,1=K η ₂, d _3,1=K η ₁=2.7 λ _min.Obviously, apart nearest array element is sensor 2 and 3 in the array, its spacing d _3,2=0.9 λ _min> λ _min/ 2, breach in classic method and require that array element distance is not more than the restriction of signal half-wavelength.In actual applications, by changing prime number group η ₁~ η _l-1and the value of K can realize the flexible configuration of sensor array element distance.

102: interval T at one time ₀in, L sensor array element is respectively with sampling rate F ₁~ F _lparallel lack sampling is done to space incident signal and (requires F ₁~ F _lthere is positive common divisor M _f, and F ₁~ F _ldivided by M _fafter relatively prime between two);

In radar, sonar etc. much application, incoming signal meets as drag:

s(t)＝Aexp(j2πf ₀t)+ω(t) (1)

Wherein, A is non-zero constant, and ω (t) is additive noise, f ₀for incoming signal frequency.The array received comprising L sensor array element to signal be:

y(t)＝a·s(t) (2)

Wherein, a is direction vector, is jointly determined by array structure and incoming signal.The matrix form of formula (2) is

$\left[\begin{array}{c}{y}_1\left(t\right)\\ .\\ .\\ .\\ y_i\left(t\right)\\ .\\ .\\ .\\ y_L\left(t\right)\end{array}\right]=\left[\begin{array}{c}1\\ .\\ .\\ .\\ e^{-j2\mathrm{\&pi;}{d}_{i,1}\sin \mathrm{\&theta;}/\mathrm{\&lambda;}}\\ .\\ .\\ .\\ e^{-j2\mathrm{\&pi;}{d}_{L,1}\sin \mathrm{\&theta;}/\mathrm{\&lambda;}}\end{array}\right]s\left(t\right)---\left(3\right)$

Wherein θ is the space incident angle of signal, and λ is incoming signal wavelength, y _ifor the signal that sensor array element i receives.Interval T at one time ₀interior all the sensors array element is sampled to space incident signal, supposes that the sampling rate of sensor array element i is F _i, i=1 ... L, then the sample of signal that sensor array element i obtains is:

$y_i\left(n\right)=y_i\left(t\right){|}_{t=n/F_i}=A{e}^{-j2\mathrm{\&pi;}{d}_{i,1}\sin \mathrm{\&theta;}/\mathrm{\&lambda;}}[\exp (j2\mathrm{\&pi;}{f}_{0}n/F_i)+\mathrm{\&omega;}\left(n_i\right)],n=1,...,N_i---\left(4\right)$

Wherein, N _ifor the sample number that sensor array element i gathers.

103: DFT is done to L road sample of signal, with Candan interpolation estimator, frequency, phase correction are done to each road DFT spectrum peak, obtain L to normalized frequency and phase estimation value

The amplitude A=2 of hypothesis space incoming signal s (t), frequency f ₀=6 × 10 ⁸hz.Appoint and get a certain sensor array element i, suppose that the signal first phase arriving this sensor array element is zero, if sensor array element sampling rate F _i=54321Hz, then frequency remainder should be f ₀modF _i=24555 (remainder operation is got in mod representative).If the sample number N that this sensor array element gathers _i=64, sample is counted equal F respectively _iand N _idFT (for F _ithe DFT of=54321 needs to carry out zero padding to sample sequence), the position of spectrum peak is obtained by spectrum peak search frequency remainder (the F as shown in table 1 of signal phase estimated value and correspondence _ithe frequency remainder of some DFT n _ithe frequency remainder of some DFT

Table 1 DFT spectrum peak position and phase information

By the DFT analysis of spectrum result that in comparison sheet 1, two kinds are counted different, find the contradiction existed between calculated amount and analysis precision.For the large (F that counts _i=54321 points) DFT, analysis precision is high, and frequency remainder and phase value error are almost 0, but computing length consuming time; And for the little (N that counts _i=64 points) DFT, computing is consuming time few, but analysis precision is poor, frequency remainder exist | the deviation of 24614-24555|=59Hz, phase estimation deviation reach-12.3592 degree.

Because CRT is responsive to Residue error, the DFT spectrum analysis errors under therefore needing introducing spectrum alignment technique to reduce small point situation, solves the contradiction between calculated amount and analysis precision.The present invention adopts Candan interpolation estimator to correct each road frequency, the frequency deviation value that recycling frequency correction obtains phase place is corrected.Idiographic flow is as follows

Step1 counts to the data sample that sensor array element i gathers and equals sample number N _idFT, obtain amplitude spectrum Y _i(k _i) and phase spectrum

Step2 utilizes place's peak value spectrum with its about two other spectral lines to frequency deviation value do to estimate by following formula, wherein real represents real part.

${\widehat{\mathrm{\&delta;}}}_i=\frac{\tan (\mathrm{\&pi;}/N_i)}{\mathrm{\&pi;}/N_i}\&CenterDot;\mathrm{teal}\left\{\frac{Y_i\left(k_i^{*}\right)-Y_i(k_i^{*}+1)}{2Y_i\left(k_i^{*}\right)-Y_i(k_i^{*}-1)-Y_i(k_i^{*}+1)}\right\}---\left(5\right)$

Step3 is according to the frequency deviation value obtained carry out frequency and phase correction, the normalized frequency after correction and phase estimation be respectively:

${\hat{f}}_i=(k_i^{*}+{\widehat{\mathrm{\&delta;}}}_i)/N_i---\left(6\right)$

Step4 is through above correction process, and can calculate the CRT remainder estimated value corresponding with Frequency Estimation is

${\hat{r}}_i=(k_i^{*}+{\widehat{\mathrm{\&delta;}}}_i)\&CenterDot;F_i/N_i---\left(8\right)$

Learnt by emulation, under identical condition N _ithe frequency remainder that some DFT obtains after overcorrect is 24555, and phase correcting value is 4.8298 × 10 ^-5degree, error is almost nil.

104: adopt the enclosed CRT algorithm that document [8] proposes, as frequency remainder, using the value of each sensor array element sampling rate as the mould needed for CRT, (L corresponding relatively prime integer is for Γ _i=F _i/ M _f, i=1...L).Remainder group prime number group Γ ₁... Γ _i... Γ _land M _fsubstitute into enclosed CRT and try to achieve signal frequency value algorithm idiographic flow is as follows:

Step1 is from given frequency remainder calculating parameter wherein

${\hat{q}}_{i,1}=\left[\frac{{\hat{r}}_i-{\hat{r}}_1}{M_f}\right],2\&le;i\&le;L---\left(9\right)$

Step2 calculates mould is except Γ _iremainder:

${\widehat{\mathrm{\&xi;}}}_{i,1}={\hat{q}}_{i,1}{\stackrel{\&OverBar;}{\mathrm{\&Gamma;}}}_{i,1}\mathrm{mod}{\mathrm{\&Gamma;}}_i,2\&le;i\&le;L---\left(10\right)$

Wherein, Γ ₁about Γ _imould inverse.

Step3 calculates fuzzy integer

${\hat{n}}_1=\underset{i=2}{\overset{L}{\mathrm{\&Sigma;}}}{\widehat{\mathrm{\&xi;}}}_{i,1}{b}_{i,1}\frac{{\mathrm{\&gamma;}}_1}{{\mathrm{\&Gamma;}}_i}\mathrm{mod}{\mathrm{\&gamma;}}_1---\left(11\right)$

Wherein, b _{i, 1}be about Γ _imould inverse, and γ _i=Γ/Γ _i, Γ is Γ ₁~ Γ _lproduct.

Step4 calculates other fuzzy integer

${\hat{n}}_i=\frac{{\hat{n}}_1{\mathrm{\&Gamma;}}_1-{\hat{q}}_{i,1}}{{\mathrm{\&Gamma;}}_i}---\left(12\right)$

Step5 calculates each road frequency estimation 1≤i≤L, gets its mean value as Frequency Estimation namely

${\hat{f}}_{0i}={\hat{n}}_i{M}_f{\mathrm{\&Gamma;}}_i+{\hat{r}}_i,1\&le;i\&le;L---\left(13\right)$

${\hat{f}}_0=\frac{1}{L}\underset{i=1}{\overset{L}{\mathrm{\&Sigma;}}}{\hat{f}}_{0i}---\left(14\right)$

From document [8], maximum detection frequency of the present invention is f _max=lcm (F ₁..., F _l), wherein lcm represents lowest common multiple.

105: obtain sensor array element 1 estimated value poor with the phase of received signal of other array element and form phase place remainder further.By phase place remainder group, prime number group η ₁~ η _l-1and M _θthe CRT substituting into closed form reconstructs nonnegative constant N ₀, and then by N ₀calculate spacing wave incident angle concrete steps are as follows:

Step1 gets any positive integer M _θ, definition for the spacing of sensor array element i and sensor array element 1 and phase of received signal difference and then phase place remainder is constructed

Step2 is by prime number group η ₁~ η _l-1, phase place remainder group and M _θsubstitute into enclosed CRT (its solution procedure is similar to Frequency Estimation process of the present invention, does not repeat) here and reconstruct nonnegative constant

Step3 by obtain signal wavelength estimated value (c is constant propagation velocity of electromagnetic wave), then calculates incident angle estimated value wherein parameter d ₀=KM _θη ₁η ₂... η _l-1.

The principle of above DOA estimating step and optimum configurations is as follows:

When adjacent sensors array element distance is greater than a half of incoming signal wavelength X, the phase differential of corresponding Received signal strength its estimated value scope be [-π, π], it is fuzzy observed reading with 2 π, namely

Wherein n is unknown fuzzy integer, and ε is phase measurement error.If θ is the incident angle of spacing wave, λ is signal wavelength, and d is the spacing of adjacent sensors, then

Parameter d is multiplied by formula (16) both sides ₀and carry out equation conversion,

Use N ₀represent d ₀sin θ/λ, obvious N ₀be a nonnegative constant, CRT can be utilized reconstruct.Defined parameters d ₀, sensor array element i and array element 1 spacing d _{i, 1}with phase of received signal difference be

d ₀＝KM _θη ₁η ₂...η _L-1

d _i,1＝Kη ₁...η _i-2η _i...η _L-1(18)

Wushu (18) substitutes into, even d=d _{i, 1}, n=n _{i, 1}, then (17) are rewritten as

$N_0=\frac{d_0\sin \mathrm{\&theta;}}{\mathrm{\&lambda;}}=n_{i-1}{M}_{\mathrm{\&theta;}}{\mathrm{\&eta;}}_{i-1}+{\hat{r}}_{\mathrm{\&theta;},i-1},i=2,...,L$

(19)

Wherein ε _{i, 1}for the measuring error of phase differential.Formula (19) meets the data model of CRT, by prime number group η ₁~ η _l-1, remainder group and parameter M _θsubstitute into enclosed CRT to obtain reconstructing nonnegative constant then binding signal wavelength estimated value the incident angle estimated value of signal can be obtained

$\widehat{\mathrm{\&theta;}}=\arcsin \frac{\widehat{\mathrm{\&lambda;}}{\hat{N}}_0}{d_0}---\left(20\right)$

By Chinese remainder theorem, nonnegative constant N ₀need meet the following conditions

$0\&le;N_0=\frac{d_0\sin \mathrm{\&theta;}}{\mathrm{\&lambda;}}<M_{\mathrm{\&theta;}}{\mathrm{\&eta;}}_1{\mathrm{\&eta;}}_2...{\mathrm{\&eta;}}_{L-1}---\left(21\right)$

By d ₀=KM _θη ₁η ₂... η _lsubstitution formula (21), has

$K<\frac{\mathrm{\&lambda;}}{\sin \mathrm{\&theta;}}---\left(22\right)$

Suppose λ _minfor the signal wavelength that system maximum detection frequency is corresponding, for ensureing that the measurement range of incident angle meets θ ∈ [-pi/2, pi/2], in practical application, K being set to and being less than λ _mina constant.In order to ensure parameter energy Accurate Reconstruction, CRT requires that Residue error is less than 1/4th of lowest common multiple, uses ε _{i, 1}represent the phase measurement error of sensor array element i and array element 1,

$\frac{M_{\mathrm{\&theta;}}}{4}>\frac{M_{\mathrm{\&theta;}}{\mathrm{\&eta;}}_{i-1}{\mathrm{\&epsiv;}}_{i,1}}{2\mathrm{\&pi;}}$

(23)

Namely ${\mathrm{\&epsiv;}}_{i,1}<\frac{\mathrm{\&pi;}}{2{\mathrm{\&eta;}}_{i-1}},i=2,...,L$

Obtain prime number group η thus ₁~ η _l-1the relation with DOA estimated accuracy of choosing.

Emulation experiment and result

According to Fig. 1 scheme, arrange the linear array comprising three sensor array elements.Sampling rate corresponding to each sensor array element is respectively F ₁=52400Hz, F ₂=54800Hz, F ₃=55600Hz, then the common divisor M of each sensor array element sampling rate _f=400, corresponding prime number group { Γ ₁, Γ ₂, Γ ₃}={ 131,137,139}.

The signal frequency range that system can be surveyed is (0, f _max], wherein f _max=lcm (F ₁..., F ₃)=997853200Hz, corresponding signal wavelength lambda _min=c ₀/ f _max=0.3006m.If K=0.3, meet K < λ _mincondition.If η ₁=3, η ₂=2, then corresponding sensor array element distance d _2,1=K η ₂=0.6m, d _3,1=K η ₁=0.9m (obviously the spacing of the adjacent array element of this array is all greater than half-wavelength).If space incident signal amplitude A=2, signal frequency is for being in (0, f _max] f in measurement range ₀=875697821.6Hz (obvious each road sampling rate, much smaller than measured signal frequency, belongs to time domain lack sampling situation), the additive noise of sensor array element Received signal strength is the white Gaussian noise of zero-mean.Suppose the sampling time T of array ₀=0.1s, then correspond to three sampling rate F within this period ₁~ F ₃interior number of samples is respectively N ₁=5250, N ₂=5480, N ₃=5560, thus respectively N is done to the sample that three sensor array elements are collected ₁, N ₂, N ₃the DFT of point, then utilizes Candan interpolation estimator to carry out spectrum and corrects, obtain frequency and phase place remainder.In DOA estimation procedure, setting parameter M _θbe 3.

In order to verify parameter estimation performance of the present invention under space-time lack sampling, follow Fig. 1 scheme 5 steps, under different signal to noise ratio (S/N ratio), carry out 1000 Monte-carlo experiment, Fig. 3 (a), Fig. 3 (b), Fig. 4 give signal frequency when changing with signal to noise ratio (S/N ratio) and DOA estimated performance relation curve.

The frequency detecting probability of Fig. 3 (a) is as undefined: when frequency detecting relative error exceeds 0.1%, be then judged as detecting successfully, otherwise be judged as detecting unsuccessfully.Can find out from Fig. 3 (a), when signal to noise ratio (S/N ratio) is higher than threshold value-22dB, reach the detection probability of success of 100%.Can find out from Fig. 3 (b), when signal to noise ratio (S/N ratio) is higher than threshold value-22dB, the square error of its frequency estimation is tending towards 0.

Detect square error curve from the DOA of Fig. 4 can find out, when signal to noise ratio (S/N ratio) is higher than threshold value-22dB, its DOA estimates that square error can be controlled within 5 °.

See Fig. 5, under space-time lack sampling, signal frequency and DOA combined measurement device comprise: sensor array element, A/D sampling thief, DSP and output driving and display circuit thereof, space far-field narrow band signal s (t) arrives each sensor array element with a certain incidence angle θ, obtain array received signal y (t)=as (t), wherein a is array direction vector.A/D sampling thief in each sensor array element is respectively with F ₁, F ₂..., F _lspeed signal parallel is sampled, will obtain that sample sequence is parallel is input to DSP (Digital Signal Processor, digital signal processor) device.Through the internal algorithm process of DSP, the Frequency Estimation and the DOA that obtain incoming signal estimate, finally drive by output and display circuit display estimated result.

Wherein, the DSP in Fig. 5 is core component, in whole estimation procedure, complete following function:

1, call core algorithm, complete the DFT process of each road sample of signal and correct the estimated value of frequency, phase place, the frequency and the DOA that complete spacing wave measure;

2, sampling rate F is adjusted in time according to actual needs ₁, F ₂..., F _l, make it meet actual needs;

3, measurement result is exported to driving and display module.

Wherein, the principal element determining the complexity of Fig. 5 system, degree of accuracy and degree of stability is the kernel estimation algorithm that DSP internal program memory stores.The internal processes flow process of DSP device is as shown in Figure 6:

The flow process of Fig. 6 is divided into following step:

First according to actual requirement, the frequency range of guestimate incoming signal, and according to real needs setpoint frequency measurement range and each road sampling rate;

CPU primary controller reads sampled data from I/O port, enters internal RAM;

DC processing.Flip-flop in measured signal can reduce measuring accuracy, therefore needs to eliminate direct current impact;

By the most crucial part that the processing procedure of Fig. 1 carries out frequency, DOA measurement is DSP algorithm, by will frequency and the DOA estimated value of incoming signal be obtained after this algorithm process;

Judge whether measurement result meets engineering demand, if do not meet, then reset sampling rate and frequency measurement scope according to result and actual demand;

If the measurement result obtained meets the demands, then export driving or display device to by DSP output bus.

The embodiment of the present invention is to the model of each device except doing specified otherwise, and the model of other devices does not limit, as long as can complete the device of above-mentioned functions.

List of references

[1]Tsui J.Digital techniques for wideband receivers[M].SciTech Publishing,2004.

[2]Arazi B.A generalization of the Chinese remainder theorem[J].Pacific Journal of Mathematics,1977,70(2):289-96.

[3]Maróti M, P,Dóra S,et al.Radio interferometric geolocation[C].Proceedings of the Proceedings of the 3rd international conference on Embedded networked sensor systems,ACM,2005:1-12.

[4]Ruegg M,Meier E,Nuesch D.Capabilities of dual-frequency millimeter wave SAR with monopulse processing for ground moving target indication[J].Geoscience and Remote Sensing,IEEE Transactions on,2007,45(3):539-53.

[5]Xia X-G,Wang G.Phase unwrapping and a robust Chinese remainder theorem[J].IEEE Signal Processing Letters,2007,14(4):247.

[6]Li G,Xu J,Peng Y-N,et al.An efficient implementation of a robust phase-unwrapping algorithm[J].Signal Processing Letters,IEEE,2007,14(6):393-6.

[7]Xia X L X-G.A fast robust Chinese remainder theorem based phase unwrapping algorithm[J].Signal Processing,2008.

[8]Wang W,Xia X-G.A closed-form robust Chinese remainder theorem and its performance analysis[J].Signal Processing,IEEE Transactions on,2010,58(11):5655-66.

[9] Liang Hong, Zhang Qi, Yang Changsheng. the sane Chinese remainder theorem of broad sense and the application [J] in space undersampled signal DOA estimates thereof. Northwestern Polytechnical University's journal, 2010,3): 409-14.

[10]Xia X-G,Liu K.A generalized Chinese remainder theorem for residue sets with errors and its application in frequency determination from multiple sensors with low sampling rates[J].Signal Processing Letters,IEEE,2005,12(11):768-71.

[11]Li X,Liang H,Xia X-G.A robust Chinese remainder theorem with its applications in frequency estimation from undersampled waveforms[J].Signal Processing,IEEE Transactions on,2009,57(11):4314-22.

[12] Liang Hong, Zhang Qi, Yang Changsheng. the Chinese remainder theorem that broad sense is sane and the application [J] in lack sampling rate lower frequency is estimated thereof. electronics and information journal, 2010,32 (8): 1802-5.

[13] Liang Hong, Zhang Heng. multiple goal frequency and orientation joint estimate new method [J] under lack sampling time empty. Northwestern Polytechnical University's journal, 2012,30 (5): 694-8.

[14]Candan C.A method for fine resolution frequency estimation from three DFT samples[J].Signal Processing Letters,IEEE,2011,18(6):351-4.

[15]Candan C.Statistical Analysis and Improvement of Fine Resolution Frequency Estimation Method From Three DFT Samples[J].2013.

It will be appreciated by those skilled in the art that accompanying drawing is the schematic diagram of a preferred embodiment, the invention described above embodiment sequence number, just to describing, does not represent the quality of embodiment.

The foregoing is only preferred embodiment of the present invention, not in order to limit the present invention, within the spirit and principles in the present invention all, any amendment done, equivalent replacement, improvement etc., all should be included within protection scope of the present invention.

Claims (4)

1. signal frequency and a DOA union measuring method under space-time lack sampling, it is characterized in that, described measuring method comprises the following steps:

Arrange the linear nonuniform noise containing L sensor array element, with sensor array element 1 for reference array element, the spacing of definition sensor array element i and sensor array element 1 is d _{i, 1};

Interval T at one time ₀in, L sensor array element is respectively with sampling rate F ₁~ F _lparallel lack sampling is done to space incident signal, wherein, F ₁~ F _lthere is positive common divisor M _f;

DFT is done to L road sample of signal, with Candan interpolation estimator, frequency, phase correction is done to each road DFT spectrum peak, obtain L to normalized frequency and phase estimation value;

Adopt enclosed CRT algorithm, by L, signal frequency value is obtained to normalized frequency

Obtain sensor array element 1 and the phase of received signal difference estimated value of other array element, form phase place remainder further, finally calculate spacing wave incident angle.

2. signal frequency and DOA union measuring method under a kind of space-time lack sampling according to claim 1, is characterized in that, described employing enclosed CRT algorithm, obtains signal frequency value by L to normalized frequency step be specially:

? as frequency remainder, using the value of each sensor array element sampling rate as the mould needed for CRT, remainder group prime number group Γ ₁Γ _iΓ _land M _fsubstitute into enclosed CRT and try to achieve signal frequency value

3. signal frequency and DOA union measuring method under a kind of space-time lack sampling according to claim 1, it is characterized in that, the described phase of received signal difference estimated value obtaining sensor array element 1 and other array element, and form phase place remainder further, calculate spacing wave incident angle and be specially:

Obtain sensor array element 1 estimated value poor with the phase of received signal of other array element and form phase place remainder further; By phase place remainder group, prime number group η ₁~ η _l-1and parameter M _θsubstitute into closed form CRT and reconstruct nonnegative constant N ₀, and then by N ₀calculate spacing wave incident angle

4. signal frequency and DOA union measuring method under a kind of space-time lack sampling according to claim 1, it is characterized in that, under described space-time lack sampling, signal frequency and DOA combined measurement device comprise: sensor array element, A/D sampling thief, DSP and output driving and display circuit thereof

Space far-field narrow band signal s (t) arrives each sensor array element with a certain incidence angle θ, obtains array received signal; A/D sampling thief in each sensor array element is respectively with F ₁, F ₂..., F _lspeed signal parallel is sampled, be input to DSP device by parallel for the sample sequence obtained, through the internal algorithm process of DSP, the Frequency Estimation and the DOA that obtain incoming signal estimate, finally drive by output and display circuit display estimated result.