A refined frequency estimation method based on ESPRIT algorithm

Abstract

The present invention relates to a refined frequency estimation method based on an ESPRIT algorithm, comprising the following steps: S1 searching for M0 peak points of the FFT spectrum function amplitude of a signal to be estimated, obtaining a frequency index at each peak point, and obtaining an initial frequency estimation value of each peak point of the signal to be estimated; S2 determining the current refined frequency interval of the ESPRIT algorithm for each peak point; S3 performing complex modulation on the signal to be estimated within the current refined frequency interval to obtain a frequency shift signal; S4 performing digital low-pass filtering on the frequency shift signal to obtain a frequency band sequence to be analyzed; S5 estimating the number of real signals within the current refined frequency interval of each peak point to obtain an optimal refined frequency interval for each peak point; S6 using the ESPRIT algorithm to estimate the frequency of each signal within the optimal refined frequency interval of each peak point to obtain a final frequency estimation value of each signal at each peak point.

Description

Technical Field

[0001] The present invention relates to the field of frequency estimation, and more specifically to a refined frequency estimation method based on ESPRIT algorithm. Background Art

[0002] Signal frequency estimation is a common problem in engineering applications, and many scenarios require accurate frequency estimation of signals. For example, in a frequency modulated continuous wave (FMCW) radar ranging system, the relevant distance information can be obtained based on the difference frequency signal between the transmitted wave and the reflected wave, that is, the estimated distance is obtained by estimating the frequency of the difference frequency signal. Therefore, the accuracy of the frequency estimation of the difference frequency signal directly affects the accuracy of the estimated distance, and the measurement accuracy of the difference frequency signal directly determines the accuracy of the ranging.

[0003] In order to improve the accuracy of frequency estimation, many frequency estimation algorithms have been proposed. Among them, the most common is the fast Fourier transform (FFT) method, but it has fence effect and spectrum leakage problems, which greatly affect the accuracy of frequency estimation. In order to overcome the shortcomings of the FFT method, the prior art uses a super-resolution algorithm called ESPRIT (signal parameter estimation based on rotational invariance technology) to estimate the frequency, which uses the eigendecomposition (EVD) or singular value decomposition (SVD) of the correlation matrix of the received signal to separate the signal subspace and the noise subspace, and then uses the orthogonal relationship between the signal subspace and the noise subspace to accurately estimate the parameters of the required signal, including the estimation of the number of signals.

[0004] The existing super-resolution frequency estimation method based on ESPRIT, on the one hand, needs to calculate the eigendecomposition of high-dimensional matrices, so its computational complexity is very high; on the other hand, since the frequency refinement interval, frequency refinement multiple and the number of signals in the refined frequency interval are interrelated parameters, and the number of signals in the refined frequency interval is unknown, the ESPRIT algorithm is difficult to use for refined frequency estimation, which will result in low frequency estimation accuracy and frequency resolution. Summary of the invention

[0005] The object of the present invention is to provide a refined frequency estimation method based on the ESPRIT algorithm, which uses the ESPRIT algorithm for refined frequency estimation to improve frequency estimation accuracy and frequency resolution while reducing computational complexity.

[0006] The present invention provides a refined frequency estimation method based on the ESPRIT algorithm, comprising the following steps:

[0007] S1: Perform N-point FFT transformation on the signal to be estimated to obtain the FFT spectrum function, search for M0 peak points of the FFT spectrum function amplitude, obtain the frequency index at each peak point, and obtain the initial frequency estimate value of each peak point in the FFT amplitude spectrum of the signal to be estimated;

[0008] S2: Determine the current refined frequency interval of the ESPRIT algorithm for each peak point;

[0009] S3: performing complex modulation on the signal to be estimated within the current refined frequency interval to obtain a frequency shift signal;

[0010] S4: performing digital low-pass filtering on the frequency-shifted signal to obtain a frequency band sequence to be analyzed;

[0011] S5: estimating the number of real signals within the current refined frequency interval of each peak point, and obtaining the optimal refined frequency interval of each peak point;

[0012] S6: Using the ESPRIT algorithm within the optimal refined frequency interval of each peak point, the difference frequency of the difference frequency signal having a number equal to the number of real signals at each peak point is estimated to obtain the final frequency estimation value of each signal at each peak point.

[0013] Furthermore, step S2 further comprises:

[0014] S21: Determine the current signal number of each peak point;

[0015] S22: Calculating the number of rows and refinement multiples of the data matrix according to the current signal number of each peak point;

[0016] S23: Calculate and obtain the current refined frequency range of each peak point using the FFT peak point and the refinement multiple.

[0017] Furthermore, for the peak point j, the number of rows and the refinement multiple of the data matrix satisfy the following relationship:

[0018]

[0019] in, is the current signal number, K j is the number of rows in the data matrix, D j To refine the multiple, Indicates not greater than N / K j where N is the maximum integer of the FFT transform.

[0020] Furthermore, the current refined frequency interval (f j,1 ,f j,2 ) satisfies the following relationship:

[0021]

[0022] in, is the refinement interval of peak point j The frequency point where the middle value of j is the frequency index of peak point j, Δf j is the half bandwidth, Δf j =W j / 2,W j To refine the bandwidth, W j =f s / D j .

[0023] Furthermore, the frequency shift signal satisfies the following relationship:

[0024]

[0025] Where x(n) is the frequency shift signal, d(n) is the band estimation signal, q is the imaginary unit, and f s =NΔf FFT is the sampling frequency, Δf FFT is the spectral line spacing, and L0 is the center shift of the frequency.

[0026] Further, step S5 includes:

[0027] S51: Group the frequency band sequences to be analyzed to construct the data matrix Z j ;

[0028] S52: Calculate the correlation matrix of the data matrix;

[0029] S53: Obtaining the eigenvalues ​​of the correlation matrix;

[0030] S54: Estimate the estimated value of the number of signals at a given false alarm probability based on the eigenvalue, compare the current number of signals with the estimated value, if the two are equal and the minimum eigenvalue of the correlation matrix is ​​much smaller than other eigenvalues, the current number of signals at this time is the true number of signals, and the current refined frequency interval at this time is the optimal refined frequency interval, and directly execute step S6; if the current number of signals is equal to the estimated value, but the minimum eigenvalue of the correlation matrix is ​​not much different from the second smallest eigenvalue, add one to the current number of signals and re-execute steps S2-S5; if the current number of signals is less than the estimated value, replace the current number of signals with the estimated value of the number of signals at the false alarm probability, and re-execute steps S2-S5.

[0031] Furthermore, the number of rows, refinement multiple and signal number of the data matrix satisfy the following relationship:

[0032]

[0033] Among them, K j is the number of rows in the data matrix, D j is the refinement factor, C j is the peak point j in the current refined frequency interval The number of real signals in .

[0034] Furthermore, the data matrix Z j The number of columns is D j -1.

[0035] Furthermore, the correlation matrix satisfies the following relationship:

[0036]

[0037] Among them, Z j is the data matrix.

[0038] Further, step S6 includes:

[0039] S61: Calculate data correlation matrix The eigenvalue decomposition of The eigenvalues ​​are signal eigenvalues, where is an estimate of the true number of signals;

[0040] S62: From the data correlation matrix Corresponding to The eigenvector of the signal eigenvalues ​​determines the signal subspace E S , then, the signal subspace E S Calculate submatrix and in By E S The first K j -1 line, By E S After K j -1 line;

[0041] S63: Compute Matrix Bundle The generalized eigenvalue decomposition of Then, by The i-th diagonal element of The phase of the peak point j is obtained within the optimal refinement frequency interval The delay estimate of the difference frequency signal Where ζ=B / T c is the slope of the linear frequency modulation;

[0042] S64: Calculate refined difference frequency estimate The final frequency estimates for each signal are then obtained:

[0043]

[0044] Among them, f e,j is the center frequency of the optimal refined frequency interval of peak point j.

[0045] The refined frequency estimation method based on the ESPRIT algorithm of the present invention adopts an iterative method to calculate the optimal parameter solution according to the constraint relationship between parameters such as the refined frequency interval, the refinement multiple, the number of signals, the number of rows and columns of the data matrix when the number of signals is unknown, thereby realizing frequency estimation based on the ESPRIT algorithm within the optimal refined frequency interval, and the frequency estimation accuracy and frequency resolution are high, and the computational complexity is low. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] Figure 1A is an FFT amplitude spectrum of a discrete received signal according to an embodiment of the present invention;

[0047] Figure 1B for Figure 1A A local area map of

[0048] Figure 2 4 is a flow chart of a refined frequency estimation method based on the ESPRIT algorithm according to an embodiment of the present invention. DETAILED DESCRIPTION

[0049] The preferred embodiments of the present invention are given below in conjunction with the accompanying drawings and described in detail.

[0050] The ranging principle of FMCW radar is as follows:

[0051] Assume that the transmission signal s of the FMCW radar is tx (t) is:

[0052]

[0053] In the formula, A tx is the amplitude of the transmitted signal, f0 is the initial frequency, ζ=B / T c is the slope of the linear frequency modulation, B is the radar bandwidth, T c is the frequency sweep period, and the initial phase of the transmitted signal is φ0. Assuming that the distance between the target and the radar is R, the time delay of the received echo signal is:

[0054]

[0055] Where R is the distance between the target and the radar, and c is the speed of light.

[0056] Due to the time delay, the radar received signal s rxIt can be expressed as:

[0057]

[0058] In the formula, A rx is the amplitude of the received signal, φ rx is the phase of the received signal. Then, the received signal is multiplied by the transmitted signal and demodulated, and the resulting signal is called the in-phase component. for:

[0059]

[0060] After mixing, the in-phase component Will pass through the low-pass filter. At this time, the low-pass filtered in-phase component I(t) can be expressed as:

[0061] I(t)=A b cos(πζτ 2 -2πζtτ-2πf0τ+φ2) (5)

[0062] Similarly, the orthogonal phase component Q(t) after low-pass filtering can be expressed as:

[0063] Q(t)=-A b sin(πζτ 2 -2πζtτ-2πf0τ+φ2) (6)

[0064] Then, the low-pass filtered in-phase component I(t) and the low-pass filtered quadrature-phase component Q(t) are combined into a complex signal, namely the difference frequency signal s(t):

[0065] s(t)=I(t)+qQ(t)=A b exp(q(-πζτ 2 +2πζtτ+2πf0τ-φ2))=A b e qψ(t) (7)

[0066] In the formula, q is the imaginary unit, A b is the amplitude of the difference frequency signal, ψ(t)=-πζτ 2 +2πζtτ+2πf0τ-φ2 is the phase of the difference frequency signal. The instantaneous frequency of the signal is f τ is the derivative of formula (7):

[0067]

[0068] According to τ=2R / c, the instantaneous frequency f τ The relationship between the distance R between the target and the radar is:

[0069]

[0070] It can be seen that the distance between the target and the radar is proportional to the beat frequency. Therefore, the accurate distance estimation between the target and the radar depends on the estimation of the beat frequency.

[0071] If the transmitted signal s tx (t) is reflected by M targets at different distances, then the FMCW radar will receive the difference frequency signals of M targets. Assume that the difference frequency signal of the i-th target is Among them A b,i is the amplitude of the i-th target difference frequency signal, is the phase of the i-th target difference frequency signal, then the sum of the M difference frequency signals, that is, the received difference frequency signal d(t) of the FMCW radar can be expressed as:

[0072]

[0073] Where w(t) is additive white Gaussian noise (AWGN), with a mean of 0 and a variance of σ 2 .

[0074] Analog-to-digital conversion via Nyquist sampling, where the required sampling frequency is f s =1 / T s , the discrete received signal d(n) of the sampled de-linear frequency modulated output d(t) can be expressed as:

[0075]

[0076] Among them, T s is the sampling period, and N is the number of sample points.

[0077] Perform an N-point FFT transformation on the estimated signal d(n) to obtain the FFT spectrum function X0(k) of d(n). However, since the spectrum resolution of FFT is Δf FFT =f s / N, FFT will not be able to distinguish the difference frequency interval smaller than Δf FFT Therefore, each peak of the FFT amplitude spectrum may contain multiple signals. Assume that there are M0 peaks in the amplitude spectrum of the FFT spectrum function, and assume that there are C j signals, the total number of signals For ease of understanding, an example is used to illustrate the FFT amplitude spectrum of the received signal d(n). Assume that the scanning bandwidth of the FMCW radar is B = 999.47055MHz, the sampling rate is fs = 92.7835kHz, the data length is N = 1024, the number of received signals is M = 6, and the frequencies are f1 = 4600Hz, f2 = 5000Hz, f3 = 5070Hz, f4 = 5540Hz, f5 = 5640Hz, and f6 = 6000Hz. Figure 1A and 1B The FFT amplitude spectrum of the discrete digital received signal d(n) is given, where Figure 1B for Figure 1A A local area map of . Figure 1B As shown in the figure, there are 4 peak points in the amplitude spectrum, of which peak 1 contains 1 signal (f1=4600Hz), peak 2 contains 2 signals (f2=5000Hz and f3=5070Hz), peak 3 contains 2 signals (f4=5540Hz and f5=5640Hz), and peak 4 contains 1 signal (f6=6000Hz). However, in practical applications, the number of signals in each peak point is not known in advance. In order to distinguish multiple signals in the same peak, an algorithm with higher resolution is required.

[0078] In order to distinguish multiple signals in the same peak, such as Figure 2 As shown, an embodiment of the present invention provides a refined frequency estimation method based on the ESPRIT algorithm. For each peak point j (1≤j≤M0), the method includes the following steps:

[0079] S1: Use the FFT spectrum function of the signal to be estimated to search for the peak point of the FFT spectrum function amplitude, and find the frequency index m′ of the peak point j in the amplitude spectrum of the FFT spectrum function j , and get the initial estimate of the frequency of the peak point j

[0080] In this embodiment, the signal to be estimated is a discrete received signal d(n) obtained by sampling the difference frequency received signal d(t) of the FMCW radar, as shown in formula (11). The initial estimated value of the frequency of the peak point j in the FFT amplitude spectrum of the signal to be estimated is:

[0081]

[0082] Where Δf FFT =f s / N is the frequency resolution of FFT.

[0083] S2: Determine the current refinement frequency interval (f j,1 ,f j,2 ),for example Figure 1BFor peak point 2, the current refined frequency interval is (f 2,1 ,f 2,2 ).

[0084] S2 specifically includes:

[0085] S21: Assume that the number of real signals contained in peak point j is C j , but in fact the real signal number C in the peak point j j is unknown, so the real signal number C j That is, the number of signals that need to be estimated. Assume that the current refined frequency interval (f j,1 ,f j,2 ) is defined as the initial value of the number of signals to be estimated in C j,max , C j,max It can be set according to the required spectral resolution or the spectral resolution of the ESPRIT algorithm. Assume that the required spectral resolution is Δf require , then we can set C j,max =(f j,2 -f j,1 ) / Δf require , where at the beginning f j,1 and f j,2 Can be set to and 1≤η≤2. Moreover, the current refined frequency interval (f j,1 ,f j,2 ) Set to the initial value C of the number of signals to be estimated j,max ,Right now

[0086] S22: According to the number of signals in peak point j Calculate the number of rows K in the data matrix j and refinement factor D j :

[0087] in It represents rounding down, that is, it represents the largest integer not greater than x. N is the number of sample points of the discrete received signal d(n), which is also the length of the FFT transformation.

[0088] S23: Using FFT peak points and refinement factor D j Calculate the current refined frequency interval of the peak point j;

[0089] First, calculate the refinement bandwidth W j =f s / D j and half bandwidth Δf j =W j / 2; then calculate the refinement interval At this time, the center frequency f of the current refined frequency interval e,j Set to:

[0090]

[0091] Therefore, the current refinement frequency interval is determined as:

[0092]

[0093] S3: performing complex modulation on the signal to be estimated within the current refined frequency interval to obtain a frequency shift signal;

[0094] The center frequency that needs to be analyzed in detail is f e,j and bandwidth is W j =f j,2 -f j,1 Frequency band of the region (f j,1 f 2,2 ), through frequency shift, the center of the observed frequency band is moved to the zero frequency position of the frequency domain coordinate, and becomes a low-frequency signal. To this end, d(n) is After complex modulation, the frequency shift signal is:

[0095]

[0096] In the formula, q is the imaginary unit, f s =NΔf FFT is the sampling frequency, Δf FFT is the line spacing, L 0,j =f e,j / Δf FFT is the center shift of the frequency, which corresponds to the center frequency f in the global spectrum e,j spectral line number, then f e,j =L 0,j Δf FFT It can be seen that complex modulation makes the frequency component f of d(n) e,j Move to x j (n), that is, the L0th spectral line in X0(k) moves to the zero frequency point of X0(k). j (k), where X0(k) is the FFT of d(n), and X j (k) is x j (n). In order to obtain X j (k) The partial refined spectrum near the zero point, whose sampling rate can be reduced to f s / D j , D j is the refinement multiple. The refinement multiple can be determined as D j =f s / Wj .

[0097] S4: performing digital low-pass filtering on the frequency-shifted signal to obtain a frequency band sequence to be analyzed;

[0098] In order to ensure that the resampled signal does not have spectrum aliasing during spectrum analysis, anti-aliasing filtering is required to filter out the frequency band signal that needs to be analyzed. j When the cutoff frequency of the digital low-pass filter is f c,j =f s / (2D j ). When D j =f s / W j When f c,j =0.5W j After digital low-pass filtering, the frequency band sequence y to be analyzed is obtained j (n), n=0,1,…,N-1.

[0099] S5: Estimate and update the number of real signals in the current refined frequency interval, and obtain the optimal refined frequency interval of the peak point j. The specific steps are as follows;

[0100] Step S51: Group the frequency band sequences to be analyzed to construct a data matrix Z j Specifically, the frequency band sequence y of N points to be analyzed after low-pass filtering is j (n) divided into K j groups, the number of data in each group is L j =D j -1. At this time, the number of data groups divided is Assume that the i-th group of data can be expressed as:

[0101] z i =[y((i-1)D j )y((i-1)D j +1)…y((i-1)D j +D j -2)] (16)

[0102] Then, construct K according to the grouping results j ×(D j -1)-dimensional data matrix Z j :

[0103]

[0104] In the formula, the data matrix Z j The dimension is K j ×(D j-1). At this time, the steering vector a(τ) of the difference frequency signal with a time delay of τ is:

[0105] a(τ)=[1,exp(j2πζτD j T s ),…,exp(j2πζτD j T s (K j -1))] T (18)

[0106] In order to ensure that the refined ESPRIT algorithm can properly distinguish each signal, the data matrix Z is required j The number of rows K j Greater than the number of signals. When the actual number of signals C j When known, At this time, the spectrum refinement factor Number of rows in the matrix With the number of columns and the real signal number C j The following conditions must be met:

[0107]

[0108] However, the actual signal number C j is unknown, and can only be assumed to be the current number of signals To construct the data matrix Z j At this time, the spectrum refinement factor D j , the number of rows of the matrix K j and the number of columns (D j -1) and the assumed current number of signals The following conditions must be met:

[0109]

[0110] If the current signal number set above Much larger than the real signal number C j , then the number of rows in the data matrix is ​​K j will be much larger than the actual required value This will greatly affect the estimation accuracy of the ESPRIT algorithm. On the other hand, if the current signal number Less than the real signal number C j , then the ESPRIT algorithm cannot estimate the frequency of the signal. Therefore, it is necessary to j and refinement factor D j Further optimization is performed by iteratively estimating the number of true signals C j , and assume that the current number of signals Updates are made to obtain the optimal signal number estimation value and the optimal refined frequency interval.

[0111] S52: To estimate the number of signals, calculate the data matrix Z j The correlation matrix (At the beginning and However, as the following algorithm iterates, will be iterated and updated by subsequent estimates of the signal, and similarly K j will also be updated accordingly):

[0112]

[0113] S53: Then the correlation matrix Perform eigendecomposition to obtain the correlation matrix The eigenvalue of . Assume that its eigenvalue is

[0114] S54: Based on eigenvalues The conventional AIC criterion and MDL criterion are used to estimate the number of signals. Or use the RMT algorithm based on random matrix theory in the literature Kritchman, S., Nadler, B.: Non-parametric detection of the number of signals: hypothesis testing and random matrix theory. IEEE Trans. Signal Process. 57 (10), 3930–3941 (2009) to estimate the estimated value of the number of signals when the false alarm probability is given Based on the assumed current number of signals Estimated value of the number of signals when and false alarm probability The relationship between them is as follows:

[0115] 1) If and Indicates Eigenvalues is indeed the characteristic value corresponding to the noise. Therefore, the current signal value is That is, the number of real signals C in the current refined frequency interval j The best estimate of Moreover, the current refined frequency interval is the optimal refined frequency interval, and the ESPRIT algorithm can be used to estimate The frequency of the signal; that is, executing step S6;

[0116] 2) If And the first Eigenvalues and Eigenvalues When they are about the same size (for example: the ratio of the two ), it indicates that The eigenvalues ​​may correspond to the eigenvalues ​​of the signal. Therefore, the current signal number May be smaller than the actual signal number C j , then the ESPRIT algorithm cannot be used for frequency estimation. At this point, it is necessary to further estimate the actual signal number C j To do this, the current signal number Updated to Right now Then go to step S22 to refine the multiple D j Sum signal number Iterative estimation and update of ;

[0117] 3) If Indicates the current number of signals Greater than the real signal number C j , the spectrum refinement multiple at this time has not reached the optimal value, which will reduce the frequency estimation accuracy and resolution. In order to optimize the refinement frequency interval and refinement multiple D j , it is necessary to further optimize and refine the frequency range and Updated to an estimate of the number of signals Right now Then go to step S22 to refine the multiple D j Sum signal number Iterative estimation and update of .

[0118] S6: Using the data matrix Z j The correlation matrix The ESPRIT algorithm is used to perform the optimal refinement frequency interval of the peak point j. The difference frequency of the difference frequency signals is estimated to obtain the final frequency estimation value of each signal;

[0119] S6 specifically includes:

[0120] S61: Calculate the data matrix Z j The correlation matrix The eigenvalue decomposition (EVD) of The characteristic values ​​are taken as signal characteristic values;

[0121] S62: From the data correlation matrix Corresponding to The eigenvector of the signal eigenvalues ​​determines the signal subspace E S , then, the signal subspace ES Calculate submatrix and in By E S The first (K-1) rows of By E S The last (K-1) rows of ;

[0122] S63: Compute Matrix Bundle The generalized eigenvalue decomposition of Then, by The i-th diagonal element of The phase of the peak point j is obtained within the optimal refinement frequency interval The delay estimate of the difference frequency signal where φ j，i for The i-th diagonal element of c is the slope of the linear frequency modulation;

[0123] S64: Calculate the optimal refinement frequency range of the peak point j Refined beat frequency estimate of a beat frequency signal Then the center frequency f of the optimal refined frequency interval of peak point j is e，j Refined difference frequency estimate Combined, we get The final frequency estimate of the difference frequency signal is:

[0124]

[0125] Among them, f e，j is the center frequency of the optimal refined frequency interval of peak point j, and its calculation formula is similar to the above formula (13), f e，j It is the mean of the upper and lower limits of the optimal refinement frequency interval.

[0126] For the M0 peak points of the FFT spectrum function X0(k) of d(n), the above process S2 to S6 is used to obtain the real signal number C in the M0 peak points. j Estimated value of and the corresponding frequency estimates The total number of signals is estimated

[0127] After obtaining the frequency estimation value of each signal, the distance estimation value of each reflecting target can be calculated based on the frequency estimation value.

[0128] It should be noted that, although the FMCW radar ranging is used as an example in this embodiment to illustrate the specific steps of the frequency estimation method of the present invention, it is not limited to being applied to FMCW radar ranging, but can also be extended to other frequency estimation applications.

[0129] The refined frequency estimation method based on the ESPRIT algorithm provided in an embodiment of the present invention adopts an iterative method to calculate the optimal parameter solution according to the constraint relationship between parameters such as the refined frequency interval, the refinement multiple, the number of signals, the number of rows and columns of the data matrix when the number of signals is unknown, thereby realizing frequency estimation based on the ESPRIT algorithm within the optimal refined frequency interval, and the frequency estimation accuracy and frequency resolution are high, and the computational complexity is low.

[0130] In order to verify the above technical effects, the method of the present invention is simulated and evaluated, and compared with the existing frequency estimation methods. The compared methods include conventional FFT method, RIFE method, FFT-based spectrum refinement estimation method (FFT CZT), cursor FFT method (zoom FFT), and conventional super-resolution method based on ESPRIT.

[0131] The parameters used in the simulation process are: the scanning bandwidth of the FMCW radar B = 999.47055MHz, the sampling rate fs = 92.7835kHz, and the data length N = 1024. Assume that the frequency estimated by the FFT method is f0, and the refined frequency intervals of the FFT CZT method and the zoom FFT method are (f0-110, f0+110), where the refinement factor of the FFT CZT method is 512 and the downsampling rate of the zoom FFT method is 210. The refined frequency interval of the proposed refined ESPRIT method is (f0-140, f0+140). In addition, the ESPRIT algorithm calculates the eigendecomposition of a 256×256 dimensional matrix. Regarding the ESPRIT frequency refinement estimation method, since the observed bandwidth bw = 280Hz, the frequency refinement factor Therefore, a 3×331-dimensional data matrix Z is constructed. Therefore, the ESPRIT frequency refinement estimation method calculates the eigendecomposition (EVD) of the 3×3-dimensional matrix or the singular value decomposition (SVD) of the 3×331-dimensional data matrix Z. The simulation results are statistical values ​​of 20,000 results.

[0132] First, consider that there is only one difference frequency signal in FMCW, and assume that its frequency is f=5100Hz. Table 1 and Table 2 respectively give the mean and standard deviation of the frequency estimation values ​​of various methods. It can be seen from Table 1 and Table 2 that the ESPRIT refined frequency estimation method of the present invention is superior to the conventional FFT method, RIFE method and cursor FFT method (zoom FFT). The standard deviation of the ESPRIT refined frequency estimation method of the present invention is slightly larger than the FFT CZT method, but its complexity is lower than the FFT CZT method. In addition, when the signal-to-noise ratio is low, the ESPRIT refined frequency estimation method of the present invention is far superior to the conventional ESPRIT method; when the signal-to-noise ratio is high, the standard deviation of the ESPRIT refined frequency estimation method of the present invention is slightly lower than the conventional ESPRIT method.

[0133] Table 1: Mean frequency estimation values ​​of different frequency estimation methods (f = 5100 Hz)

[0134] SNR(dB) -14 -12 -10 -8dB -6dB -4 0 4 FFT 5074.2 5074.1 5074.1 5074.1 5074.1 5074.1 5074.1 5074.1 Rife 5092.9 5.095.2 5097.5 5099.0 5099.8 5100.0 5100.0 5100.0 FFT CZT 5100.0 5100.0 5100.0 5100.0 5100.0 5100.0 5100.0 5100.0 Zoom FFT 5100.6 5100.3 5100.2 5100.2 5100.2 5100.2 5100.2 5100.2 ESPRIT method 5099.6 5100.0 5100.0 5100.0 5100.0 5100.0 5100.0 5100.0 The present invention 5100.0 5100.0 5100.0 5100.0 5100.0 5100.0 5100.0 5100.0

[0135] Table 2: Standard deviation of frequency estimation values ​​for different frequency estimation methods (f = 5100 Hz)

[0136] SNR(dB) -14 -12 -10 -8 -6 -4 0 4 FFT 2.3094 0 0 0 0 0 0 0 Rife 20.4957 16.7575 12.4272 8.1314 4.6913 2.5480 1.2629 0.7965 FFT CZT 5.6319 4.4578 3.5112 2.7845 2.2043 1.7579 1.1177 0.7104 zoom FFT 8.5369 6.5465 4.9968 3.8454 2.9840 2.3630 1.4896 0.9426 ESPRIT 70.4894 5.7679 4.2470 3.2486 2.4957 1.9723 1.2182 0.7654 The present invention 7.3773 5.5033 4.1758 3.2405 2.5451 2.0054 1.2608 0.7902

[0137] In order to verify the ability of the method of the present invention to distinguish multiple difference frequency signals, consider that there are two difference frequency signals in FMCW, and the two frequencies are relatively close, one of which has a frequency of f1=5000Hz, and the other has a frequency of f2=5090Hz. Tables 3 and 4 respectively give the mean and standard deviation of the frequency estimation values ​​of various methods. It can be seen from Tables 3 and 4 that, firstly, the traditional FFT method and Rife method cannot distinguish the two signals, and can only estimate the frequency of one of the signals, and there is a large deviation; although the standard deviation of the estimation results of the FFT CZT method and the zoom FFT method is smaller than that of the ESPRIT refined frequency estimation method of the present invention, the estimation deviations of these two methods are much larger than that of the ESPRIT refined frequency estimation method of the present invention; finally, although the estimation deviation of the traditional ESPRIT method is smaller than that of the ESPRIT refined frequency estimation method of the present invention, its standard deviation is much larger than that of the ESPRIT refined frequency estimation method of the present invention when the signal-to-noise ratio is low.

[0138] Table 3: Mean values ​​of the estimated values ​​of different frequency estimation methods (f1 = 5000 Hz, f2 = 5090 Hz)

[0139]

[0140] Table 4: Standard deviation of the estimated values ​​of different frequency estimation methods (f1 = 5000 Hz, f2 = 5090 Hz)

[0141]

[0142]

[0143] The above is only a preferred embodiment of the present invention, and is not intended to limit the scope of the present invention. The above embodiment of the present invention can also be modified in various ways. That is, all simple, equivalent changes and modifications made according to the claims and the description of the present invention fall within the scope of protection of the claims of the present invention. The contents not described in detail in the present invention are all conventional technical contents.

Claims

1. A refined frequency estimation method based on ESPRIT algorithm, characterized in that: The following steps are involved: S1: Perform N-point FFT transformation on the signal to be estimated to obtain the FFT spectrum function, where N is the length of the FFT transformation; search for M0 peak points of the FFT spectrum function amplitude, obtain the frequency index at each peak point, and obtain the initial frequency estimate value of each peak point in the FFT amplitude spectrum of the signal to be estimated; S2: Determine the current refined frequency interval of the ESPRIT algorithm for each peak point; S3: performing complex modulation on the signal to be estimated within the current refined frequency interval to obtain a frequency shift signal; S4: performing digital low-pass filtering on the frequency-shifted signal to obtain a frequency band sequence to be analyzed; S5: estimating the number of real signals within the current refined frequency interval of each peak point, and obtaining the optimal refined frequency interval of each peak point; S6: Using the ESPRIT algorithm within the optimal refined frequency interval of each peak point, the difference frequency of the difference frequency signal having a number equal to the number of real signals at each peak point is estimated to obtain the final frequency estimation value of each signal at each peak point.

2. The refined frequency estimation method based on the ESPRIT algorithm according to claim 1, characterized in that: Step S2 further comprises: S21: Determine the current signal number of each peak point; S22: Calculating the number of rows and refinement multiples of the data matrix according to the current signal number of each peak point; S23: Calculate and obtain the current refined frequency range of each peak point using the FFT peak point and the refinement multiple.

3. The refined frequency estimation method based on the ESPRIT algorithm according to claim 2, characterized in that: For peak point j, the number of rows and refinement multiples of the data matrix satisfy the following relationship: in, is the current signal number, K j is the number of rows in the data matrix, D j To refine the multiple, Indicates not greater than N / K j The maximum integer.

4. The refined frequency estimation method based on the ESPRIT algorithm according to claim 3, characterized in that: Current refined frequency range (f j,1 ,f j,2 ) satisfies the following relationship: in, is the refinement interval of peak point j The frequency point where the middle value of j is the frequency index of peak point j, Δf j is the half bandwidth, Δf j =W j / 2,W j To refine the bandwidth, W j =f s / D j , f s is the sampling frequency.

5. The refined frequency estimation method based on the ESPRIT algorithm according to claim 4, characterized in that: The frequency shift signal satisfies the following relationship: Among them, x j (n) is the frequency shift signal, d(n) is the band estimation signal, q is the imaginary unit, f s =NΔf FFT is the sampling frequency, Δf FFT is the spectral line spacing, L0 is the center shift of the frequency, f e,j is the center frequency of the current refined frequency interval, n is the sequence number of the signal at different sampling points, n = 0, 1, ..., N-1.

6. The refined frequency estimation method based on the ESPRIT algorithm according to claim 1, characterized in that: Step S5 includes: S51: Group the frequency band sequences to be analyzed to construct the data matrix Z j ; S52: Calculate the correlation matrix of the data matrix; S53: Obtaining the eigenvalues ​​of the correlation matrix; S54: Estimate the estimated value of the number of signals at a given false alarm probability based on the eigenvalue, compare the current number of signals with the estimated value, if the two are equal and the minimum eigenvalue of the correlation matrix is ​​much smaller than other eigenvalues, the current number of signals at this time is the true number of signals, and the current refined frequency interval at this time is the optimal refined frequency interval, and directly execute step S6; if the current number of signals is equal to the estimated value, but the minimum eigenvalue of the correlation matrix is ​​not much different from the second smallest eigenvalue, add one to the current number of signals and re-execute steps S2-S5; if the current number of signals is less than the estimated value, replace the current number of signals with the estimated value of the number of signals at the false alarm probability, and re-execute steps S2-S5.

7. The refined frequency estimation method based on the ESPRIT algorithm according to claim 6, characterized in that: The number of rows, refinement multiples and signal numbers of the data matrix satisfy the following relationship: Among them, K j is the number of rows in the data matrix, D j is the refinement factor, C j is the peak point j in the current refined frequency interval The number of real signals in .

8. The refined frequency estimation method based on the ESPRIT algorithm according to claim 7, characterized in that: The data matrix Z j The number of columns is D j -1.

9. The refined frequency estimation method based on the ESPRIT algorithm according to claim 7, characterized in that: The correlation matrix satisfies the following relationship: Among them, Z j is the data matrix, is the conjugate transpose of the data matrix.

10. The refined frequency estimation method based on the ESPRIT algorithm according to claim 1, characterized in that: Step S6 includes: S61: Calculate data correlation matrix The eigenvalue decomposition of The eigenvalues ​​are signal eigenvalues, where is an estimate of the true number of signals; S62: From the data correlation matrix Corresponding to The eigenvector of the signal eigenvalues ​​determines the signal subspace E S , then, the signal subspace E S Calculate submatrix and in By E S The first K j -1 line, By E S After K j -1 line; among them, K j is the number of rows in the data matrix; S63: Compute Matrix Bundle The generalized eigenvalue decomposition of Then, by The i-th diagonal element φ j,i =exp(q2πζτ j,i D j T s ) The phase of the peak point j is obtained within the optimal refinement frequency interval The delay estimate of the difference frequency signal Where ζ=B / T c is the slope of the linear frequency modulation, B is the radar bandwidth, T c is the frequency sweep period, D j is the refinement factor, T s is the sampling period; S64: Calculate refined difference frequency estimate The final frequency estimates for each signal are then obtained: Among them, f e,j is the center frequency of the optimal refined frequency interval of peak point j.

Images