A phase unwrapping algorithm and device based on cycle extension

Abstract

A phase unwrapping algorithm and device based on periodic extension, the method comprising: step 1, calculating the cross-spectral phase, finding the cross-spectral phase of two array element signals, the value being the phase difference wrapped between π and π; step 2, periodic extension, adding or subtracting integer multiples of 2π from the phase spectrum of one period obtained in step 1 to obtain different periods, forming a periodic extension phase matrix; step 3, selecting phase points, starting from the middle row of the first column, using the difference between the current position and the next column as the quality basis; if the difference is less than a threshold value, the phase point is directly selected; otherwise, the phase value predicted by linear fitting of the selected phase point and the difference between the next column are used as the quality basis, a secondary judgment is performed, and a straight line, namely the unwrapped phase, is obtained by traversing the entire matrix; step 4, linear fitting, using the least squares method to perform linear fitting on the selected phase points to obtain the linear slope, which is divided by 2π and the frequency resolution to obtain the estimated delay value. The present invention can fully utilize the phase information of the entire bandwidth of the signal, thereby improving the accuracy of delay estimation.

Description

Technical Field

[0001] The present invention belongs to the field of time delay estimation algorithms in signal processing, and in particular relates to a phase unwrapping algorithm and device suitable for ultra-wideband signals under ultra-short baseline array conditions. Background Art

[0002] In the field of passive positioning, high-precision delay estimation is an important research topic. Delay estimation can achieve ranging and direction finding of the target. The accuracy of delay estimation has a crucial impact on the performance of the entire positioning system. Currently, commonly used delay estimation methods are mainly based on optimal signal matching, such as cross-correlation method, generalized cross-correlation method, and high-order cumulant method. However, the performance of these methods is limited by the sampling rate of the signal, and may not meet the high precision requirements at low sampling rates. While upsampling methods can improve accuracy, they increase the computational complexity of the system and are easily affected by noise. Phase-based methods are another commonly used delay estimation method. They make full use of the amplitude and phase information of the signal to perform delay estimation, thereby improving estimation accuracy.

[0003] However, when the baseline length exceeds half the signal wavelength, phase ambiguity arises. The measured phase (i.e., ambiguous phase, which varies only between -π and π) cannot directly reflect the true phase (which may have multiple cycles). For broadband signals, the noise-free phase difference between two receiver signals is a straight line as a function of frequency, and the slope corresponds to the physical quantity of the delay value. The straight line of the broadband signal phase difference with respect to frequency becomes blurred again, which is called phase wrapping. The frequency phase wrapping effect greatly reduces bandwidth utilization, resulting in low delay estimation accuracy. The process of restoring the phase-wrapped sawtooth wave of the broadband signal to a straight line is called phase unwrapping or phase unwrapping.

[0004] When using cross-spectral phase for delay estimation, the higher the bandwidth utilization, the higher the delay estimation accuracy. Solutions to the phase wrapping problem mainly include spatial diversity algorithms and frequency diversity algorithms. The spatial diversity algorithm solves the phase wrapping problem by setting baselines of different lengths. This algorithm can reduce the impact of phase wrapping to a certain extent, but it increases the complexity of the system. The frequency diversity algorithm eliminates the phase wrapping problem through the relationship between the phases of different frequencies. Currently, the main algorithms used for phase unwrapping in delay estimation include classical methods and differential methods. However, these methods accumulate errors at each mutation point and cannot reduce the impact of noise on the wrapped phase spectrum, resulting in low estimation accuracy under low signal-to-noise ratio conditions. Summary of the Invention

[0005] The present invention aims to overcome the above-mentioned shortcomings of the prior art and provides a phase unwrapping algorithm and device based on cycle extension.

[0006] A first aspect of the present invention relates to a phase unwrapping algorithm based on cycle extension, comprising the following steps:

[0007] Step 1: Calculate the cross-spectral phase;

[0008] The cross-spectrum of the two array element signals—the phase difference between the two signals at the same frequency—is calculated. The phase spectrum is obtained, with the index corresponding to the frequency and the value being the phase difference. Affected by noise, it resembles a sawtooth wave, appearing as a row vector wrapped between -π and π. Phase wrapping can occur when the array element spacing exceeds half the signal wavelength.

[0009] Step 2: Cycle extension;

[0010] The intertwined cross-spectral phase is still a periodic function. In the horizontal frequency direction, the period depends on the path difference, which is related to the array element spacing, the incident direction, and the speed of sound. However, in the vertical phase difference direction, the period is 2π. Therefore, in the vertical phase difference direction, adding or subtracting integer multiples of 2π yields different row vectors to form a periodic extension phase matrix. The maximum multiple depends on the bandwidth and array element spacing, and is taken as Bd / c, where B is the bandwidth, c is the speed, and d is the array element spacing. The next step is to select phase points from the phase period extension matrix across rows (and periods).

[0011] Step 3: Select phase point;

[0012] Phase unwrapping is achieved by correctly finding a line that meets the requirements across rows (and periods) in the periodic extension phase matrix. If the signal bandwidth starts at zero frequency, the phase is correct. If a frequency interval is selected, the slope of the line is known, meaning that the group delay equals the phase delay, and the correct phase value can still be recovered.

[0013] Specifically, it can be divided into the following sub-steps:

[0014] S31 selects a starting point;

[0015] The winding phase is a sawtooth wave. Assume that the expression of the sawtooth wave in one cycle is: Where a = 2πτ, b is a constant, and the period is From the Fourier transform of the sawtooth wave, we can get: Where f0 is the frequency corresponding to the maximum amplitude of the spectrum. First, perform FFT on the original single-line cross-spectrum phase. The number of FFT points is equal to the signal length, which is Nfft. Find the corresponding position with the maximum amplitude, and multiply the imaginary part of the corresponding Fourier transform value by 2πT / Nfft to obtain the specific delay value τ (T is the signal sampling time). Divide this by Δf (the FFT frequency resolution in step 1) to get the number of frequency cycles of the cross-spectral phase. Divide the cross-spectral phase into frames by period, with a repetition rate of 25%. Use the second-order difference of the phase points within each period as the quality criterion. The smaller the second-order difference, the better the quality. Select the columns of the best-quality period starting points as the starting points of the entire phase matrix. Because the slope is uncertain, the starting rows of the entire phase matrix are the middle rows of the matrix. This step ensures that even at low signal-to-noise ratios, the best-quality subband within the broadband can be used as the starting point.

[0016] S32 quality judgment;

[0017] Sub-step 1 determines the starting position, which is the current position in this sub-step. The difference between the current position and the first-order phase of the next column is used as the quality basis. The smaller the difference, the better the quality. The signal power and noise power are estimated by power spectrum density. The specific threshold value is the theoretical phase difference 2πΔfd / c plus the phase error tan affected by noise. -1 (1 / SNR), where SNR is the signal-to-noise ratio (SNR) calculated from the power spectrum of the signal and noise, expressed in true values. First, select the phase point with the best quality in the next column and compare the phase difference with the threshold. If the phase difference is less than the threshold, select that phase point and move the current position to the best quality position in the next column and repeat step S32. Otherwise, proceed to the next step.

[0018] S33 forecast revision;

[0019] When the phase difference exceeds the threshold, the phase value of the low-quality point is predicted by least squares fitting using the selected high-quality points. The predicted value and the first-order phase difference value of the next column are used as the quality basis. The best quality phase point in the next column is selected and the phase difference is compared with the threshold value. If it is less than the threshold, the current position is moved to the best quality position in the next column. Otherwise, the column is discarded and the next column is judged, and then step S32 is executed.

[0020] S34 traverses the matrix;

[0021] Because step 3 starts judging from the starting position to the right, the left of the starting position has not yet been judged, so the starting position is reset as the current position, and the traversal direction is changed to the left. The judgment rules are the same as step S32 and step S33. After all the traversals, the corresponding value of the current position of each column is the final unwrapping phase value.

[0022] Step 4: Fold in half;

[0023] In order to avoid errors caused by noise, the periodic extension matrix is ​​folded left and right, and step 3 is repeated. The number of selected phase points is used as the basis for quality judgment, and the one with more phase points is selected as the final unwrapping phase line.

[0024] Step 5: Linear fitting;

[0025] The phase points are fitted by the least squares fitting method to obtain the mutual spectrum slope k and the delay value to be estimated. where Δf is the FFT frequency resolution in step 1.

[0026] Furthermore, step 1 specifically includes: assuming x(t) and y(t) are the time domain signals of two array elements, and their Fourier transforms are X(f) and Y(f) respectively, then their cross spectrum P xy (f) can be expressed as: P xy (f) = X(f)Y * (f). Where Y * (f) represents the complex conjugate of Y(f). The phase in the specified frequency range is taken out to obtain the cross-spectral phase for the next step. The 3dB bandwidth or 10dB bandwidth of the signal can be selected.

[0027] The second aspect of the present invention provides a phase unwrapping device based on periodic extension, comprising a memory and one or more processors, wherein the memory stores executable code, and when the one or more processors execute the executable code, they are used to implement a phase unwrapping method based on periodic extension of the present invention.

[0028] A computer-readable storage medium of the present invention stores a program, which, when executed by a processor, implements a phase unwrapping method based on cycle extension of the present invention.

[0029] The present invention provides a phase unwrapping algorithm based on periodic extension, which can be used for time delay estimation in real-time positioning systems. Compared with traditional algorithms, the present algorithm has the advantages of simple implementation, high precision, and low formation dependence. The algorithm uses the unwrapped phase spectrum to perform time delay estimation, which can make full use of the phase information of the entire bandwidth of the signal and can effectively eliminate low-quality points, thereby improving the accuracy of time delay estimation. The algorithm does not require a special formation and is still applicable even when the baseline length is much larger than the signal wavelength. Therefore, it is particularly suitable for real-time positioning systems of ultra-short baseline arrays. The phase unwrapping algorithm based on periodic extension proposed in the present invention has high practical value and broad application prospects, and can be widely used in fields such as real-time positioning systems.

[0030] The advantages of this invention are that it fully utilizes the periodicity of the wrapped phase. Unlike existing unwrapping methods that select phase points within a single cycle, this invention directly obtains a period-extended phase matrix containing correct phase information through period extension. Subsequently, through steps such as quality assessment and prediction correction, the influence of noise is effectively reduced, thereby accurately unwrapping the phase values. This invention can fully utilize phase information across the entire signal bandwidth, thereby improving the accuracy of delay estimation. BRIEF DESCRIPTION OF THE DRAWINGS

[0031] Figure 1 It is the algorithm flow chart of the present invention.

[0032] Figure 2 Schematic diagram of phase winding of the present invention.

[0033] Figure 3 This is a period extension diagram of the present invention. The same color represents one period, and the phase values ​​of different colors differ by only an integer multiple of 2π.

[0034] Figure 4 This is a diagram of the unwrapping results of the periodic extension unwrapping algorithm of the present invention, and the bold straight line is the result of phase unwrapping.

[0035] Figure 5 Schematic diagram of the planar cross array structure of the present invention.

[0036] Figure 6 This is the cross-spectral entanglement phase diagram of array elements 1 and 2 obtained by calculation according to the present invention.

[0037] Figure 7 This is a comparison chart of the unwinding results of the array elements 1 and 2 of the present invention using the classical method and the method of the present invention.

[0038] Figure 8 This is a comparison chart of the unwinding results of array elements 3 and 4 of the present invention using the classical method and the method of the present invention. DETAILED DESCRIPTION

[0039] The technical solution of the present invention is further described below with reference to the accompanying drawings.

[0040] Example 1

[0041] This embodiment describes how the present invention is applied to planar cross array arrival angle estimation.

[0042] Planar cross array structure such as Figure 5 As shown, array element 1 and array element 2 are located at the positive and negative half axes of the x-axis respectively, array element 3 and array element 4 are located at the positive and negative half axes of the y-axis respectively, and the diameter of the array circumscribed circle d = 0.5m.

[0043] Assume the target signal is a broadband LFM signal between 50kHz and 100kHz, with a duration of 0.1s, a signal-to-noise ratio of 5dB, and an azimuth angle (angle with respect to the positive x-axis) of 20 degrees. In this case, the half-wavelength is much smaller than the array element spacing. The target is 500m from the array center, meeting the far-field condition. The sampling rate is 300kHz.

[0044] A phase unwrapping algorithm based on cycle extension of the present invention comprises the following steps:

[0045] Step 1: First, perform windowed FFT of the signals received by element 1 and element 2 with equal length points. Then, multiply the FFT result of element 1 by the conjugate of the FFT result of element 2, and take the phase spectrum corresponding to the frequency range of 50kHz to 100kHz to obtain the cross-spectral phase of the signals collected by element 1 and element 2. Figure 6 shown.

[0046] Step 2: Then, the phase spectrum obtained in the first step is periodically extended in the phase direction, that is, the original phase spectrum is added or subtracted by integer multiples of 2π to obtain a phase period extension matrix.

[0047] Step 3: Next, phase point selection is performed. First, the original single-row cross-spectral phase is subjected to an FFT to calculate the number of cycles. The cross-spectral phase is then framed by cycle with a repetition rate of 25%. Finally, the second-order difference of the phase points within each cycle is used as the quality criterion. The smaller the second-order difference, the better the quality. The index of the starting point of the cycle with the best quality is selected as the column index of the starting point, and the middle row of the entire phase matrix is ​​used as the row index of the starting point. The phase difference between the current position and the next column is used as the quality criterion. The smaller the difference, the better the quality. The threshold value is calculated to be 0.3272. The best-quality phase value in the next column is selected and compared with the threshold. If it is less than the threshold, the phase point is selected. Otherwise, the phase difference between the current position and the next column is used as the quality criterion using a linear fit of the selected phase points. A second evaluation is performed. If it is less than the threshold, the phase point is selected. Otherwise, the phase point is directly evaluated in the next column until the entire matrix is ​​traversed to the right. At this point, the phase to the left of the starting position has not yet been determined. The starting position is reset as the current position, and the traversal direction is changed to the left. The quality evaluation and prediction correction are repeated to complete the phase point selection and obtain the unwrapped phase line.

[0048] Step 4: Fold the periodic extension matrix left and right, repeat step 3, select the same number of phase points, and use the unwrapping phase line selected in step 3 as the final unwrapping phase line.

[0049] Step 5: Use the least square fitting method to fit the unwrapped phase line to obtain the mutual spectrum slope k. The estimated delay value can be obtained by The delay estimation error of array elements 1 and 2 is 4.6915×10 -7 s, while the error of cross-correlation delay estimation is 2.1323×10 -4 s.

[0050] Follow the above steps again to obtain the phase unwrapping results and delay estimation value τ of array element 3 and array element 4 43 ,like Figure 8 As shown. The azimuth angle of the plane cross array can be estimated by Directly obtain, the estimated angle error is 0.2553 degrees.

[0051] In summary, the present invention is applied to planar cross array arrival angle estimation, and the delay estimation value obtained by detangling has a small error, which is better than the cross-correlation and classical detangling methods.

[0052] Example 2

[0053] This embodiment relates to a phase dewrapping device based on periodic extension, comprising a memory and one or more processors, wherein the memory stores executable code, and when the one or more processors execute the executable code, they are used to implement a phase dewrapping method based on periodic extension of embodiment 1.

[0054] Example 3

[0055] This embodiment relates to a computer-readable storage medium having a program stored thereon. When the program is executed by a processor, a phase unwrapping method based on period extension in embodiment 1 is implemented.

[0056] The foregoing is merely an embodiment of the present invention and is not intended to limit the present invention. It will be apparent to those skilled in the art that various modifications and variations of the present invention are possible. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention are intended to be included within the scope of the claims of the present invention.

Claims

1. A phase unwrapping algorithm based on cycle extension, comprising the following steps: Step 1: Calculate the cross-spectral phase; Calculate the cross-spectrum of the two array element signals, that is, the phase difference between the two signals at the same frequency. This gives the phase spectrum, where the index corresponds to the frequency and the value is the phase difference. Affected by noise, it resembles a sawtooth wave, appearing as a row vector wrapped between -π and π. When the array element spacing exceeds half the signal wavelength, phase wrapping can occur. Step 2: Cycle extension; The intertwined cross-spectral phase is still a periodic function. In the horizontal frequency direction, the period depends on the path difference, which is specifically related to the array element spacing, the incident direction, and the speed of sound. However, in the vertical phase difference direction, the period is 2π. Therefore, in the vertical phase difference direction, adding or subtracting integer multiples of 2π results in different row vectors to form a periodic extension phase matrix. The maximum multiple depends on the bandwidth and array element spacing, and is taken as Bd / c, where B is the bandwidth, c is the speed, and d is the array element spacing. The next step is to cross rows from the phase period extension matrix, that is, to traverse across periods to select phase points. Step 3: Select phase point; Phase unwrapping is achieved by finding a straight line that meets the conditions across rows in the periodic extension phase matrix. If the signal bandwidth starts from zero frequency, the correct phase is obtained. If a frequency interval is selected, the slope of the straight line is known, that is, the group delay is equal to the phase delay, and the correct phase value can still be recovered. Step 4: Fold in half; To avoid errors caused by noise, fold the periodic extension matrix left and right, repeat step 3, and use the number of selected phase points as the quality judgment basis, and select the one with more phase points as the final unwrapping phase line; Step 5: Linear fitting; The phase points are fitted by the least squares fitting method to obtain the mutual spectrum slope k and the delay value to be estimated. Where Δf is the FFT frequency resolution in step 1; Step 3 specifically includes the following sub-steps: S31 selects a starting point; The winding phase is a sawtooth wave. Assume that the expression of the sawtooth wave in one cycle is: Where a = 2πτ, b is a constant, and the period is From the Fourier transform of the sawtooth wave, we can get: Where f0 is the frequency corresponding to the maximum amplitude of the spectrum; first perform FFT on the original single-line cross-spectrum phase, the number of FFT points is equal to the signal length Nfft, find the corresponding position with the maximum amplitude, and multiply the imaginary part of the corresponding Fourier transform value by 2πT / Nfft to obtain the specific delay value τ, T is the signal sampling time length, Divide by Δf to get the number of cycles in the frequency direction of the cross-spectrum phase. Divide the cross-spectrum phase into frames according to the period, with a repetition rate of 25%. Use the second-order difference of the phase point in each period as the quality basis. The smaller the second-order difference, the better the quality. Select the column number of the period starting point with the best quality as the column number of the starting point of the entire phase matrix. Because the slope is uncertain, the number of rows of the starting point of the entire phase matrix is ​​the middle row of the matrix. This step ensures that the best quality subband within the broadband can be used to start selection even at low signal-to-noise ratios. S32 quality judgment; Determine the starting position, which is the current position in this sub-step. Use the first-order phase difference between the current position and the next column as the quality basis. The smaller the difference, the better the quality. Estimate the signal power and noise power through the power spectrum density. The specific threshold value is the theoretical phase difference 2πΔfd / c plus the phase error tan affected by noise. -1 (1 / SNR), where SNR is the signal-to-noise ratio calculated from the power spectrum of the signal and noise, and is expressed in true values. First, the phase point with the best quality in the next column is selected, and the phase difference is compared with the threshold value. If it is less than the threshold, the phase point is selected, and the current position is moved to the position with the best quality in the next column and sub-step S32 is repeated. Otherwise, the next step is performed. S33 forecast revision; When the phase difference exceeds the threshold, the phase value of the low-quality point is predicted by least squares fitting using the selected high-quality points. The predicted value and the first-order phase difference value of the next column are used as the quality basis. The best-quality phase point in the next column is selected and the phase difference is compared with the threshold value. If the phase difference is less than the threshold, the current position is moved to the best-quality position in the next column. Otherwise, the column is discarded and the next column is judged, and then sub-step S32 is executed; S34 traverses the matrix; Because step 3 starts judging from the starting position to the right, the left of the starting position has not yet been judged, so the starting position is reset as the current position, and the traversal direction is changed to the left. The judgment rules are the same as step S32 and step S33. After all the traversals, the corresponding value of the current position of each column is the final unwrapping phase value.

2. A phase unwrapping algorithm based on cycle extension as claimed in claim 1, characterized in that: Step 1 specifically includes: let x(t) and y(t) be the time domain signals of two array elements, and their Fourier transforms are X(f) and Y(f) respectively, then their cross spectrum P xy (f) can be expressed as: xy (f) = X(f)Y * (f); where Y * (f) represents the complex conjugate of Y(f). The phase in the specified frequency range is taken out to obtain the cross-spectral phase for the next step. The 3dB bandwidth or 10dB bandwidth of the signal can be selected.

3. A phase unwrapping device based on period extension, characterized in that: It comprises a memory and one or more processors, wherein the memory stores executable code, and when the one or more processors execute the executable code, they are used to implement a phase unwrapping method based on period extension according to any one of claims 1 to 2.

4. A computer-readable storage medium, characterized in that A program is stored thereon, and when the program is executed by a processor, a phase unwrapping method based on periodic extension according to any one of claims 1 to 2 is implemented.

Images