Platform multi-harmonic self-noise interference parameter estimation method

Abstract

The invention discloses a platform multi-harmonic self-noise interference parameter estimation method, which can realize accurate and steady estimation of harmonic self-noise interference parameters with multi-fundamental frequency characteristics generated by equipment such as various motors in an unmanned platform. Comprising the following steps: carrying out background equalization and line spectrum extraction on a platform noise spectrum, and obtaining a corresponding power spectrum and local signal-to-noise ratio data; separating the harmonic sequence based on a frequency multiplication method; roughly estimating a harmonic fundamental frequency based on a greatest common divisor method; the fundamental frequency is accurately estimated based on weighted least squares, the fence effect of DFT can be overcome to a certain extent, and the influence of frequency resolution and high and low frequency changes on the harmonic fundamental frequency estimation accuracy is avoided; calculating a multi-frame fundamental frequency statistical distribution condition based on the histogram; and performing sequential harmonic fundamental frequency estimation based on a multi-frame fundamental frequency statistical distribution condition. Finally, multi-harmonic fundamental frequency stable and accurate estimation under precision constraint is realized, and an accurate harmonic interference parameter estimation method is provided for suppressing harmonic interference of various unmanned platforms.

Description

Technical Field

[0001] The present invention relates to the technical field of underwater acoustic signal processing, and particularly to a method for estimating the parameters of self-noise interference of multiple harmonics of a platform. Background Art

[0002] The detection of underwater pulse signals is usually carried out under non-cooperative conditions. Underwater acoustic signals are received and processed by a passive sonar to extract the required underwater environment information. This method can effectively monitor the underwater environment and provide a decision-making basis for underwater operations of one's own side on the premise of ensuring its own safety, thus showing broad application prospects in both military and civilian fields. However, under non-cooperative conditions, due to the lack of prior information about the signal, key parameters such as the form, existence, direction of arrival, arrival time, frequency, and amplitude of the signal are unknown, which makes it impossible to use traditional methods such as matched filtering to achieve the best detection performance of the signal. In addition, the complexity of the ocean environment noise and the self-noise of the sonar installation platform further exacerbate the challenges of pulse signal detection under non-cooperative conditions, resulting in a significant decline in detection performance. Especially when the interfering line spectrum shows amplitude fluctuations or has pulse characteristics, false alarms are easily triggered, thus increasing the difficulty of detection. Therefore, the interception and detection of pulse signals under non-cooperative conditions are still an important and unsolved problem in the field of underwater acoustic signal processing. The underwater platform for installing the sonar pulse signal acquisition and detection system is usually an unmanned underwater vehicle (UUV). The various devices of such platforms are highly concentrated, and the electromagnetic interference and acoustic interference problems between different devices are serious. The high-speed switching of the power switching devices of the propulsion motor and the large number of high-order harmonics brought by the periodic pulse width modulation signal, as well as the nonlinearity of circuits and devices, are the main noise sources of passive sonar electromagnetic interference. Among them, the interference with high-order harmonic characteristics has the characteristics of a continuous spectrum and has the greatest impact on pulse interception and detection. Especially in the case of a small signal-to-noise ratio, the signal will be masked in the noise of harmonic interference. Especially for different platforms, the harmonics generated may have multiple fundamental frequencies, making it more difficult to identify and suppress the harmonics.

[0003] Zeng Qingjun proposed a method based on non - linear polynomial fitting curve for continuous spectrum feature extraction. This method uses the least - squares method to estimate the parameters of the fitting curve, and then obtains the continuous spectrum profile. By segmenting and averaging the continuous spectrum profile, the feature of the average intensity of the continuous spectrum is obtained. Yin Jingwei proposed a method to extract the harmonic fundamental frequency from the purified DEMON line spectrum using the greatest common divisor method, so as to obtain the blade frequency of the propeller shaft. Bai Jingxian improved the greatest common divisor method on this basis, combined with the remainder threshold algorithm for fundamental frequency estimation, effectively solving the problem of large fundamental frequency estimation error in the traditional greatest common divisor method. These methods were initially mainly used to resist the interference of ship - radiated noise. However, the interference faced by the underwater unmanned platform (UUV) equipped with a pulse detection system under non - cooperative conditions is different from the ship - radiated noise interference. Its interference frequency is higher and has complex click - type interference characteristics. Summary of the Invention

[0004] Object of the Invention: Aiming at the problems and deficiencies in the prior art, the present invention provides a method for estimating the parameters of multi - harmonic self - noise interference of a platform. The interference line spectrum brought by the underwater unmanned platform becomes an important reason for the decline in the performance of conventional frequency - domain pulse signal interception detection and subsequent parameter estimation.

[0005] Technical Solution: To reduce the influence of the platform interference line spectrum, it is necessary to identify the platform interference line spectrum. Through experimental data confirmation, the platform interference that has a greater impact on signal interception detection has obvious harmonic characteristics. Therefore, the present invention proposes a method for estimating the parameters of multi - harmonic self - noise interference of a platform, which is used to provide a reference for subsequent suppression of platform interference and thus improve the interception detection performance of the system.

[0006] A method for estimating the parameters of multi - harmonic self - noise interference of a platform includes the following steps:

[0007] Step S1: Perform background equalization on the platform noise, calculate the corresponding power spectrum, and obtain M ordered platform interference line spectra {f m}, m = 1, 2,..., M after line spectrum extraction, and calculate the local signal - to - noise ratio.

[0008] Step S2: Since each harmonic of a harmonic signal is an integer multiple of the corresponding fundamental frequency, for a composite harmonic signal, assume that the fundamental frequencies at this time are f F1 , f F2 , f F3 ,... f Fj ,..., f FJ , f Fj is the j - th fundamental frequency, and J is the number of composite harmonic fundamental frequencies. For example, the frequencies of each harmonic corresponding to the fundamental frequency f F1 are f F1,i = i × f F1,i=1,2,…, i is the order of the corresponding fundamental frequency, so the harmonic fundamental frequency is the greatest common divisor of the corresponding harmonic frequencies of each order. Therefore, the line spectrum sequences with different fundamental frequencies are separated based on the frequency doubling method.

[0009] Step S201, preliminary search: setting the search harmonic range based on the possible power switch frequency range of the motor [f Min ,f Max ], and arrange the M platform interference line spectra extracted based on step S1 in ascending order of frequency.

[0010] Since the energy at the fundamental frequency of the harmonic signal is large, we first start from [f Min ,f Max ] range. Select the lowest frequency f 1 At the line spectrum, set the low frequency frequency drift error Δf L =3·fs / WL, mid-high frequency frequency drift error Δf H =5·fs / WL, WL is the number of window length points for calculating the power spectrum in step S1, Δf=f s / WL is the frequency resolution. The frequency boundary between low frequency and medium and high frequency is defined as F mid K is an integer, and f is determined by the following formula 1 Is it the base frequency?

[0011]

[0012] Traverse f F1,i For {f m Each line spectrum frequency f in m , determine f that meets the above formula m is the fundamental frequency f F1 The i-th frequency f F1,i And when, when the traversal is completed, N 1 ≥2, then judge f 1 is the fundamental frequency f F1 , is the fundamental frequency f F1 The corresponding harmonic line spectrum sequence, where N 1 For preliminary judgment, f F1 is the number of line spectra of the fundamental frequency.

[0013] Step S202, secondary search: continue searching for values ​​greater than f 1 The frequency and corresponding harmonics are searched to determine the second line spectrum frequency f 2 With f F1 relationship, that is, there are two situations If the situation S 1 Then execute the judgment:

[0014]

[0015] Traverse f F2,i For each line spectral frequency greater than f in {f m}, determine the f that satisfies the above formula 2 as the fundamental frequency f m and its i-th harmonic frequency f F2 . When the traversal ends and N F2,i ≥2, then determine that f 2 = f 2 and the fundamental frequency f F2 and the harmonic frequencies F2 are the harmonic line spectral sequence, where N is the number of line spectra initially determined with f 2 as the fundamental frequency. F2

[0016] If condition S 2 is true, then continue to search for the next line spectrum f 3 and continue to execute this step.

[0017] Step S203, Three - time search: The number of fundamental frequencies of the test platform harmonics is one or two. Therefore, after determining f F1 , f F2 , one more search is sufficient to meet the harmonic estimation requirements of most platforms. That is, search for the next line spectrum f Min , f Max within the frequency range [f 3 . If the line spectral frequency f 3 is a harmonic of f F1 , f F2 within the allowable error, then end the search. If the frequency f 3 is not a harmonic of f F1 , f F2 and there are still harmonics of f 3 among the subsequent line spectral frequencies, then determine that f 3 = f F3 and the frequencies with a harmonic relationship within the allowable error with it as f F3 correspond to the harmonic line spectral sequence

[0018] Step S3, Calculate the difference frequency: Based on Step S2, the harmonic line spectral sequences corresponding to different fundamental frequencies are obtained. For the convenience of subsequent step description, here take the line spectral sequence F1 corresponding to the fundamental frequency f for subsequent processing description. Ideally, the difference frequency between two harmonic lines is the fundamental frequency or an integer multiple of the fundamental frequency. By calculating the difference frequency, possible values of the greatest common divisor can be obtained. For the N 1 extracted ordered line spectra, calculate the difference frequencies F 1,i,j ​

[0019] F 1,i,j = abs(f F1,i - f F1,j ), i, j = 1, 2, …, N 1 ; i ≠ j

[0020] And discard the difference frequencies that fall outside the frequency range [f Min , f Max . Sort the remaining difference frequency data in ascending order to obtain the difference frequency array F 1,i,j .

[0021] Step S4. To accurately evaluate the calculated fundamental frequency and difference frequency quality, the quality factor is calculated in this step, and the fundamental frequency is found based on the greatest common divisor method. The following steps are analyzed based on the separated line spectrum sequence, i.e., the single fundamental frequency line spectrum sequence. The analysis steps for different fundamental frequencies can refer to the same analysis.

[0022] Step S401. Define the quality factor as the elements in the corresponding difference frequency array F 1,i,j that may contain equal elements or elements with a deviation less than Δf between difference frequencies, i.e., approximately equal, and count the number of the same difference frequencies within the error tolerance range Δf as the number of line spectra of F n , defined as s n , and define the new difference frequency array, i.e., the possible fundamental frequency candidate array, as {F n}, where n = 1, 2,... N s , that is, there are N s difference frequencies in this new difference frequency array.

[0023] Step S402. Define the quality factor q n , and initialize all q n to 0. Calculate Q i,n = f F1,i / F n , i = 1, 2,... N 1 , when |Q i,n - K| < Δ, i = 1, 2, …, N 1 ; n = 1, 2, … N s holds, the quality factor q n corresponding to F n is incremented by 1, where K = round(Q i,n ) is a positive integer. Δ is the deviation control amount. The F n with the largest q n is the greatest common divisor of the difference frequencies.

[0024] Step S403. To ensure the accuracy of the quality factor estimation for subsequent judgment, define the quality factor m n = sn ·q n , that is, define a new quality factor m n as the quality factor s n and the quality factor q n The product of, has the maximum quality factor m n The difference frequency F n is the corresponding line spectrum sequence The greatest common divisor of

[0025] Step S5. Since the DFT algorithm used in line spectrum estimation has a fence effect and the extracted line spectrum frequencies are located at discrete frequency points, the frequencies of the harmonic line spectra corresponding to the same fundamental frequency are not in an absolute fundamental frequency harmonic relationship. Especially in the high-frequency band, the error caused by the large amplification factor is greater. Therefore, based on the rough fundamental frequency and the corresponding quality factor obtained by the greatest common divisor in step S4, this step uses the weighted least squares (WLS) method for accurate fundamental frequency estimation.

[0026] Step S501. Explain the fundamental frequency estimation expression of the weighted least squares method. Based on the local line spectrum signal-to-noise ratio obtained in step S1, the received signal model of the harmonic sine signal in colored background noise can be expressed as

[0027]

[0028] where w 0 is the fundamental frequency, α k is the amplitude of the k-th harmonic signal, φ k is the initial phase of the k-th harmonic, j is the imaginary unit, v(n) is the colored noise with zero mean, and the power spectral density is σ 2 (w), and N is the data length. Let

[0029]

[0030] where θ′ is the observation vector and η′ is the parameter estimation vector; based on the local line spectrum signal-to-noise ratio obtained in step S1, the observation vector θ′ has obtained the initial parameter estimation through a certain line spectrum estimation method that ignores the harmonic structure where is the estimation of w k =kw 0 Assume there exists a full column rank matrix such that θ′ = S′η′. At the same time, based on the Gauss-Markov theorem, the weighted matrix W′ is the inverse matrix of the error autocorrelation matrix R v The weighted least squares estimation result of the parameter estimation vector η′ is solved​ Considering that the amplitude estimation is independent of other parameter estimations, the observed vector θ after removing the amplitude parameter and the parameter estimation vector η are set as

[0031]

[0032] Meanwhile, it is assumed that there still exists a matrix with full column rank such that θ = Sη. Assume that the initial parameter estimation of θ is To obtain the optimal estimation, a weighted residual quadratic form is constructed as the cost function Then the weighted least squares estimation result of the parameter estimation vector η In the case of unbiased estimation, the lower bound of the error autocorrelation matrix R v is given by the Cramér-Rao lower bound (CRLB) matrix, that is, R v = R CRLB . When the number of observed data N is large enough, the weighted matrix W = R CRLB -1 , the CRLB matrix can be approximated as a block diagonal form, where

[0033]

[0034] where is the local signal-to-noise ratio of the k-th sine signal. Since W K ′ is a block diagonal matrix, the WLS method makes no contribution to the amplitude estimation accuracy. Based on this, the estimation cost function J(η) is derived. We get

[0035]

[0036] where w 0 is the fundamental frequency, α k is the amplitude of the k-th harmonic signal, K is the highest order, φ k is the initial phase of the k-th harmonic, is the estimated phase of the k-th harmonic, N is the data length, find the partial derivative with respect to φ k Take the partial derivative and set the result to 0, we get

[0037]

[0038] And substitute the result into J(η) and take the partial derivative with respect to the final estimation target w 0 Take the partial derivative.

[0039]

[0040] where SNR k is the local signal-to-noise ratio of the line spectrum obtained for the k-th frequency point based on step S1

[0041] Step S502. Based on the fundamental frequency f obtained in step S4 F1 its corresponding frequencies, harmonic frequencies, orders, and the local signal-to-noise ratio of the line spectrum after background equalization at the k-th frequency point obtained in step S1

[0042] Step S5021. Calculate the order of the fundamental frequency. For the extracted line spectrum in which those satisfying |Q i,n -K| < Δ, i = 1, 2, …, N 1 ; n = 1, 2, … N s the required line spectra are updated and rearranged as The corresponding line spectra of other fundamental frequencies are also rearranged in the same way.

[0043] Step S503. Correct the rough fundamental frequency based on the previous calculations. For the line spectra within, calculate the and respectively, and take the result as the estimation result.

[0044] Step S504. Re-confirm the f estimated in step S503, that is, compare the line spectrum frequencies within F1,WLS with the integer multiples of f . F1,WLS

[0045]

[0046] where F mid is the cut-off frequency between low frequency and high frequency. The number of line spectra satisfying the above formula is the quality factor q of the fundamental frequency f F1,WLS . When q n > N n , it is determined that the fundamental frequency f e is accurate. N is determined by experience and generally needs to satisfy N F1,WLS ≥ 3. e e

[0047] Step S6. Based on the corrected fundamental frequency, perform sequential harmonic fundamental frequency estimation by statistical histogram to obtain a more robust harmonic group estimation effect.

[0048] Step S601. Loop and execute the above operations from step S2 to step S5 until the number of consecutive estimated frames is greater than the number of frames P required for the received signal to stabilize. Assume that a total of P frames of platform noise are estimated to obtain P-frame fundamental frequency estimation results {f F,WLS} = {f F1,WLS , f F2,WLS ,... f​​​Fp,WLS}. Take the maximum value f of the fundamental frequency of the result F,WLSMax and the minimum value f F,WLSMin , and divide the interval [f F,WLSMax , f F,WLSMin into P - 1 parts.

[0049] Step S602. Take the center frequency of the interval with the highest probability of the estimated fundamental frequency in the previous P frames as the fundamental frequency result f of the robust estimation of the P frames F,Final ; when the current frame number is greater than P, f F,WLS is the fundamental frequency result f of the robust estimation of the data of the previous frame of the current estimated frame F,Final , is the result based on the weighted least squares estimation of the current frame, and perform sequential statistics based on the histogram for robust fundamental frequency estimation. is the probability that the estimation result in the previous P frames of the current frame falls within the interval centered on the previous robust f F,WLS , is the number of frames in the P frames that fall within the interval centered on the previous robust f F,WLS , The probability that the estimated fundamental frequency result of the current estimated P frames falls within the interval centered on the current frame estimation , is the number of frames that fall within the interval centered on the current frame estimation , and perform the following formula to update the robust fundamental frequency f of the current estimated frame F,Final Update

[0050]

[0051] where p L and p H are the probabilities for evaluating the robustness of the fundamental frequency estimation respectively, and can be selected according to experience.

[0052] A computer device, which includes a memory, a processor, and a computer program stored on the memory and executable on the processor. When the processor executes the above computer program, it implements the steps of the above - mentioned platform multi - harmonic self - noise interference parameter estimation method.

[0053] A computer - readable storage medium, which stores a computer program for executing the above - mentioned platform multi - harmonic self - noise interference parameter estimation method.

[0054] The beneficial effects of the present invention are as follows: Compared with the prior art, the platform multi-harmonic self-noise interference parameter estimation method provided by the present invention starts from the frequency-domain characteristics of platform interference in view of the detection difficulties in the actual application scenario of underwater unmanned platform pulse detection under non-cooperative conditions. Under the condition of multi-fundamental frequency harmonic interference, it can autonomously achieve: (1) harmonic sequence separation based on the frequency doubling method; (2) accurate estimation of fundamental frequency harmonics based on the least squares method; (3) robust estimation of sequential harmonic fundamental frequencies based on histogram statistics. Furthermore, by removing the interference frequency band line spectra of the corresponding fundamental frequency and its harmonics in the received signal, the purpose of anti-platform interference is achieved. This method overcomes the defect that single-harmonic fundamental frequency estimation cannot adapt to different unmanned platforms. At the same time, the accurate estimation method based on the least weighted squares method can overcome the fence effect of the DFT to a certain extent, avoiding the influence of frequency resolution and high and low frequency changes on the accuracy of harmonic fundamental frequency estimation. BRIEF DESCRIPTION OF THE DRAWINGS

[0055] Figure 1 is the flowchart of the method of the present invention;

[0056] Figure 2 is the typical single-harmonic power spectrum;

[0057] Figure 3 is the typical composite-harmonic power spectrum;

[0058] Figure 4 is the schematic diagram of line spectrum detection in the embodiment;

[0059] Figure 5 is the estimated value of the fundamental frequency harmonic based on the weighted least squares method in the 45th frame of the embodiment;

[0060] Figure 6 is the robust estimated value of the composite-harmonic double fundamental frequencies in multiple frames based on the least squares method and histogram statistics in the embodiment. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0061] The following further clarifies the present invention in conjunction with specific embodiments. It should be understood that these embodiments are only used to illustrate the present invention and not to limit the scope of the present invention. After reading the present invention, various equivalent modifications of the present invention by those skilled in the art fall within the scope defined by the appended claims of this application.

[0062] A method for estimating platform multi-harmonic self-noise interference parameters. The overall process of the method is as Figure 1As shown. Based on the frequency characteristics of harmonics, the line spectra extracted from the underwater acoustic signals received by the underwater unmanned platform detection system are analyzed. First, the frequency doubling method is used to separate the harmonic line spectra with different fundamental frequencies. Then, the weighted least squares method is used to accurately correct the harmonic fundamental frequencies roughly estimated based on the greatest common divisor. Then, based on the harmonic fundamental frequencies of multiple frames of data estimated in the previous step, sequential correction and update are performed using histogram statistics, thereby realizing a robust estimation of the harmonic fundamental frequencies of the underwater unmanned platform interference. Furthermore, by removing the interference band line spectra of the corresponding fundamental frequencies and their multiples in the received signals, the purpose of anti-platform interference is achieved. The steps are as follows:

[0063] Step 1: After background equalization and line spectrum extraction of the platform noise spectrum, M ordered platform interference line spectra {f m} are obtained, where m = 1, 2, …, M, and the corresponding power spectra and local signal-to-noise ratios are calculated.

[0064] For the platform noise, the short-time power spectrum of the signal is obtained by using the method based on the short-time Fourier transform. The short-time power spectrum calculation formula

[0065]

[0066] where n is the discrete time series, w N (n) is the sliding window function, N is the window length, l is the sliding window step size, m is the number of sliding window slides, s(n) is the received signal, and S N (k, m) is the short-time spectrum, and P(k, m) is the short-time power spectrum of the corresponding data. And based on the short-time power spectrum P(k, m) of this data and the ratio of P(k, m) to the smoothed value obtained by median filtering P M of order P smooth (k, m) is used as the characterization value of the local signal-to-noise ratio SNR(k, m).

[0067] Step 2: Since each harmonic of the harmonic signal is an integer multiple of the corresponding fundamental frequency, the power spectra of the typical single harmonic and the typical composite harmonic are as Figure 2 and Figure 3 shown. Therefore, the line spectrum sequences with different fundamental frequencies are first separated based on the frequency doubling method.

[0068] Specifically, in this embodiment, Step 2 includes the following steps:

[0069] Step 2-1: Preliminary search: Set the search harmonic range [f Min , f Max based on the possible frequency range of the power switching device of the motor, and arrange the M platform interference line spectra extracted in Step S1 in ascending order of frequency.

[0070] Since the energy at the fundamental frequency of the harmonic signal is relatively large, first start from [fMin , f Max range. Select the lowest frequency f 1 for the line spectrum at this point, the low-frequency frequency drift error Δf L = 3·fs / WL, the medium- and high-frequency frequency drift error Δf H = 5·fs / WL, define the frequency boundary between low frequency and medium-high frequency as F mid . K is an integer, N M is the number of line spectra, judge

[0071]

[0072] If the above formula holds, then determine that f 1 is the fundamental frequency f F1 , f F1,i is f F1 corresponding harmonic line spectrum sequence where N 1 is the number of line spectra initially determined with f F1 as the fundamental frequency and N 1 ≥ 2.

[0073] Step 2-2, secondary search: Continue to search for the harmonics corresponding to the subsequent frequencies, and judge the relationship between the second line spectrum frequency f 2 and f F1 , that is, there are two cases S 1 In this case, execute

[0074]

[0075] If satisfied, then determine the harmonic line spectrum sequence

[0076] Step 2-3, tertiary search: The number of fundamental frequencies of the harmonics on the test platform is one or two. Therefore, after determining f F1 , f F2 , one more search can meet the harmonic estimation requirements of most platforms. That is, look for the next line spectrum f Min , f Max range. If the line spectrum frequency f 3 is a multiple frequency of f 3 within the error tolerance, end the search. F1 , f F2 in the error allowable range, end the search.

[0077] Step 3, obtain the difference frequency: Based on Step S2, the harmonic line spectrum sequences corresponding to two possible fundamental frequencies are obtained. For the convenience of subsequent step description, here take the line spectrum sequence corresponding to the fundamental frequency f 1 Subsequent processing instructions are provided. Ideally, the difference frequencies between two harmonic line spectra are the fundamental frequency or an integer multiple of the fundamental frequency. By calculating the difference frequencies, possible values of the greatest common divisor are obtained. For the N 1 extracted ordered line spectra, calculate the difference frequencies F 1,i,j

[0078] F 1,i,j = abs(f F1,i - f F1,j ), where i, j = 1, 2, …, N 1 ; i ≠ j

[0079] And discard the difference frequencies that fall outside the frequency range [f Min , f Max . Sort the remaining difference frequency data in ascending order to obtain the difference frequency array F 1,i,j .

[0080] Step 4: To accurately evaluate the calculated fundamental frequency and difference frequency quality, the quality factor is calculated in this step, and then the fundamental frequency is found based on the greatest common divisor method.

[0081] Specifically, in this embodiment, Step 4 includes the following steps:

[0082] Step 4-1: Define the quality factor as the elements in the corresponding difference frequency array F 1,i,j that may be equal or have a deviation less than Δf between the difference frequencies, i.e., are approximately equal, and count the number of line spectra with the same difference frequency within the deviation less than the error tolerance Δf as F n , which is defined as s n , and define a new difference frequency array as {F n}, where n = 1, 2,...N s , that is, there are N s difference frequencies in this new difference frequency array.

[0083] Step 4-2: Define the quality factor q n , and initialize all q n to 0. For each line spectrum frequency taken in M 1 , calculate the quotient Q n with a deviation from the integer multiple when divided by F n in {F i,n}. When

[0084] |Q i,n - K| < Δ, i = 1, 2, …, N 1 ; n = 1, 2, …N

[0085] holds, increment the quality factor corresponding to F n by 1, where K = round(Q i,n) is a positive integer. Δ is the deviation control quantity. This process is the process of finding the greatest common divisor.

[0086] Step 4-3, define the quality factor m n = s n ·q n , that is, the comprehensive quality factor s n and the quality factor q n . The difference frequency F n in the one with the largest quality factor m n is the greatest common divisor of the corresponding line spectrum sequence .

[0087] Step 5, based on the rough fundamental frequency and the corresponding quality factor obtained in step S4 based on the greatest common divisor, this step uses the weighted least squares (WLS) method for accurate estimation of the fundamental frequency.

[0088] Specifically, in this embodiment, step 5 includes the following steps:

[0089] Step 5-1, based on the fundamental frequency f F1 obtained in step S4 and its corresponding frequency, multiple frequency, order, and the local line spectrum signal-to-noise ratio after background equalization of the kth frequency point obtained in step S1

[0090] Step 5-1-1, calculate the order of the fundamental frequency. For the extracted line spectrum that satisfies |Q i,n - K| < Δ, i = 1, 2,..., N 1 ; n = 1, 2,..., N in step S402, update the line spectrum and reconstruct it as Reconstruct the corresponding line spectra of other fundamental frequencies in the same way.

[0091] Step 5-2, correct the rough fundamental frequency based on the previous calculation. For the line spectra within , find and respectively and take the result as the estimation result.

[0092] Step 5-3, reconfirm the f F1,WLS estimated in step S503, that is, compare the line spectrum frequencies within with the integer multiples of f F1,WLS .

[0093]

[0094] where F mid is the cut-off frequency between low frequency and high frequency, and the number of line spectra satisfying the above formula is the fundamental frequency fF1,WLS Quality factor q n When q n > N e is satisfied, it is determined that The fundamental frequency f F1,WLS is accurate. N e is determined by experience and generally needs to satisfy N e ≥ 3.

[0095] Step 6: Update the center frequency of the harmonic group obtained based on the histogram statistics according to the corrected fundamental frequency, so as to obtain a more robust harmonic group estimation effect.

[0096] Specifically, in this embodiment, Step 6 includes the following steps:

[0097] Step 6-1: Loop and execute the above operations from Step S2 to Step S5. Assume that a total of P frames of platform noise are estimated to obtain the fundamental frequency estimation results {f F,WLS} = {f F1,WLS , f F2 , WLS,... f Fp,WLS}. Take the maximum fundamental frequency f F,WLSMax and f F,WLSMin of the results, and divide the interval [f F,WLSMax , f F,WLSMin into P - 1 parts.

[0098] 1. Step 6-2: Take the maximum fundamental frequency f F,WLSMax and f F,WLSMin of the results of the first P frames, and divide the interval [f F,WLSMax , f F,WLSMin into P - 1 parts. is the probability that the P frames fall within the previous robust f F,WLSMax interval, is the corresponding number in the P frames, is the probability that the current estimated P frames fall within the current frame estimated f F,WLS interval, is the corresponding number, f F,Final is the finally estimated robust fundamental frequency.

[0099]

[0100] where p L and p H are the probabilities for evaluating the robustness of the fundamental frequency estimation respectively, and can be selected according to experience.

[0101] Simulation example 1:

[0102] This example constructs a platform multi - fundamental - frequency harmonic interference noise containing pulse signals and white noise based on MATLAB simulation. The pulse signal is set with a frequency f = 22 kHz, a pulse width τ = 9 ms, a period T = 1 s, a signal length of 4T = 4 s, a signal - to - noise ratio of - 10 dB, and the fundamental frequencies of the composite harmonic interference are set to 20 kHz and 25 kHz respectively, with a duty cycle of 50% and a signal - to - interference ratio of - 7 dB. By Figure 4 It can be seen that a large number of harmonic line spectra are included in the detected line spectra, and the real signal line spectra may be submerged. Figure 5 It can be seen that the harmonic fundamental frequencies estimated based on the weighted least - squares method are around 20 kHz and 25 kHz. Figure 6 The error of the robust estimation result of the sequential harmonic fundamental frequency based on histogram statistics is within 5 Hz, which can meet the application requirements of anti - platform interference based on the sub - harmonic fundamental frequency.

[0103] From the results of the embodiment, it can be seen that under the conditions of a signal - to - noise ratio of - 10 dB and a signal - to - interference ratio of - 7 dB, accurate and robust estimation of the fundamental frequencies of the composite harmonic interference can be achieved.

[0104] Obviously, those skilled in the art should understand that each step of the method for estimating the platform multi - harmonic self - noise interference parameters in the above - mentioned embodiments of the present invention can be implemented by a general - purpose computing device. They can be concentrated on a single computing device or distributed on a network composed of multiple computing devices. Optionally, they can be implemented by program codes executable by the computing device. Thus, they can be stored in a storage device and executed by the computing device. And in some cases, the steps shown or described can be executed in a sequence different from here, or they can be made into individual integrated circuit modules respectively, or multiple modules or steps among them can be made into a single integrated circuit module to be implemented. In this way, the embodiments of the present invention are not limited to any specific combination of hardware and software. Where the present invention is not elaborated, they are all well - known technologies to those skilled in the art.

Claims

1. A platform multi-harmonic self-noise interference parameter estimation method, characterized in that: The steps include: Step S1: After background equalization and line spectrum extraction of the platform noise spectrum, M ordered platform interference line spectra are obtained as {f m }, m = 1, 2, ..., M, and calculate the corresponding power spectrum and local signal-to-noise ratio; Step S2, firstly separating the line spectrum sequences with different fundamental frequencies based on the frequency doubling method; Step S3, obtaining harmonic line spectrum sequences corresponding to different fundamental frequencies based on step S2, and obtaining difference frequencies based on the harmonic line spectrum sequences; Step S4, in order to accurately evaluate the quality of the calculated base frequency and difference frequency, the quality factor is calculated, so as to roughly estimate the base frequency based on the greatest common divisor method; Step S5, based on the rough fundamental frequency and the corresponding quality factor obtained based on the greatest common divisor in step S4, the weighted least square method is used to accurately estimate the fundamental frequency; Step S6: perform sequential harmonic fundamental frequency estimation based on histogram statistics according to the corrected fundamental frequency, thereby obtaining a more robust harmonic fundamental frequency estimation effect.

2. The platform multi-harmonic self-noise interference parameter estimation method according to claim 1 is characterized in that: In step S2, suppose the current base frequency is f F1 ,f F2 ,f F3 ,...f Fj ,...,f FJ , f Fj is the jth fundamental frequency, J is the number of composite harmonic fundamental frequencies, for example f F1 The frequencies of the harmonics corresponding to the fundamental frequency are f F1,i =i×f F1 , i = 1, 2, ..., i is the order of the frequency doubling of the corresponding fundamental frequency, and the harmonic fundamental frequency is the greatest common divisor of the corresponding harmonic frequencies of each order; firstly, the line spectrum sequences with different fundamental frequencies are separated based on the frequency doubling method; Step S201, preliminary search: set the search harmonic range [f Min ,f Max ], arranging the M platform interference line spectra extracted based on step S1 in ascending order of frequency; First, [f Min ,f Max ] range, select the line spectrum at the lowest frequency f1, the low-frequency frequency drift error Δf L =3·fs / WL, setting the mid-high frequency frequency drift error Δf H =5·fs / WL, the frequency boundary between low frequency and medium and high frequency is defined as F mid ; WL is the number of window length points for calculating the power spectrum in step S1, Δf = f s / WL is the frequency resolution; K is an integer, N M is the number of line spectra, and the following formula is used to determine whether f1 is the fundamental frequency: Traverse f F1,i For {f m Each line spectrum frequency f in m , determine f that meets the above formula m is the fundamental frequency f F1 The i-th frequency f F1,i And when the traversal is completed and N1≥2 is satisfied, f1 is determined to be the fundamental frequency f F1 , is the fundamental frequency f F1 The corresponding harmonic line spectrum sequence, where N1 is the initial judgment based on f F1 is the number of line spectra of the fundamental frequency; Step S202, secondary search: continue to search for frequencies greater than f1 and corresponding harmonics, and determine whether the second line spectrum frequency f2 is F1 relationship, that is, there are two situations If S1 is the case, the judgment is executed: Traverse f F2,i For {f m }, for each line spectrum frequency greater than f2, determine the f that meets the above formula m is the fundamental frequency f F2 The i-th frequency f F2,i When N2≥2 is satisfied after the traversal, it is determined that f2=f F2 and fundamental frequency f F2 Multiplication frequency is a harmonic line spectrum sequence, where N2 is the initial determination of f F2 is the number of line spectra of the fundamental frequency. If the situation is S2, the next line spectrum f3 will be searched and the step will be continued; Step S203, three searches: after determining f F1 ,f F2 After that, at the frequency [f Min ,f Max ] range to find the next line spectrum, if the line spectrum frequency f3 is f F1 ,f F2 If the frequency f3 is not within the allowable error range, the search ends. F1 ,f F2 If there is a frequency multiple of f3 in the measured line spectrum frequency, then f3 is determined to be the fundamental frequency f F3 The frequency with which the error is within the allowable range is f F3 The corresponding harmonic line spectrum sequence 3. The platform multi-harmonic self-noise interference parameter estimation method according to claim 1 is characterized in that: In step S3, the difference frequency is obtained: based on step S2, the harmonic line spectrum sequences corresponding to different fundamental frequencies are obtained, and the line spectrum sequence corresponding to the fundamental frequency f1 is taken. Follow-up processing instructions: for the extracted N1 ordered line spectra, calculate the difference frequency F between them 1,i,j F 1,i,j =abs(f F1,i -f F1,j ),i,j=1,2,…,N1;i≠j And falls within the frequency range [f Min ,f Max ] are discarded, and the retained difference frequency data are sorted in order from small to large to obtain the difference frequency array F 1,i,j .

4. The platform multi-harmonic self-noise interference parameter estimation method according to claim 1, characterized in that: In step S4, in order to accurately evaluate the quality of the calculated base frequency and difference frequency, the quality factor is calculated in this step, so as to find the base frequency based on the greatest common divisor method; Step S401, define the quality factor as the corresponding difference frequency array F 1,i,j There may be elements that are equal or whose difference frequencies have a deviation less than Δf, that is, they are approximately equal, and the same difference frequencies within the error tolerance range Δf are counted as F n The number of line spectra is defined as s n , and define a new difference frequency array as {F n }, where n = 1, 2, ... N s , that is, the new difference frequency array has N s difference frequency; Step S402: define the quality factor q n , all q n Initialized to 0, Calculate the frequency of each line spectrum in i,n =f F1,i / F n ,i=1,2,...N1,when |Q i,n -K|<Δ,i=1,2,…,N1; n=1,2,…N s When established, F n The corresponding quality factor q n Add 1, where K = round(Q i,n ) is a positive integer, Δ is the deviation control amount, q n The largest F n That is the greatest common divisor of the difference frequencies; Step S403: define the quality factor m n =s n ·q n , that is, define a new quality factor m n is the quality factor s n and quality factor q n The product of has the maximum quality factor m n The difference frequency F n The corresponding line spectrum sequence The greatest common divisor of .

5. The platform multi-harmonic self-noise interference parameter estimation method according to claim 1 is characterized in that: In step S5, based on the rough fundamental frequency and the corresponding quality factor obtained based on the greatest common divisor in step S4, this step uses the weighted least squares method to accurately estimate the fundamental frequency; Step S501, illustrating the fundamental frequency estimation expression of the weighted least squares method. Based on the line spectrum local signal-to-noise ratio obtained in step S1, the received signal model of the harmonic sinusoidal signal in the colored background noise can be expressed as Where w0 is the fundamental frequency, α k is the amplitude of the kth harmonic signal, j is the imaginary unit, φ k is the initial phase of the kth harmonic, v(n) is the zero-mean colored noise, and the power spectral density is σ 2 (w), N is the data length; let Where θ′ is the observation vector and η′ is the parameter estimation vector. Based on the line spectrum local signal-to-noise ratio obtained in step S1, θ′ as the initial parameter estimation has been obtained by a line spectrum estimation method that ignores the harmonic structure, that is, in w k = estimate of kw0; suppose there exists a matrix with full column rank So that θ′=S′η′; at the same time, based on the Gauss-Markov theorem, the weighting matrix W′ is the error autocorrelation matrix R v The inverse matrix WLS estimation of parameter estimation vector η′ Solved The observation vector θ and parameter estimation vector η after removing the amplitude parameter are At the same time, suppose there still exists a matrix with full column rank Let θ = Sη, to obtain the optimal estimate, construct the weighted residual quadratic form as the cost function The weighted least squares estimate of the parameter estimation vector η is In the case of unbiased estimation, the error autocorrelation matrix R v The lower limit of is given by the Cramer-Rao lower bound matrix, namely R v =R CRLB ; When the observed data N is large enough, the weighted matrix W = R CRLB -1 , the CRLB matrix can be approximated as a block diagonal form, in is the local signal-to-noise ratio of the kth sinusoidal signal, since W K ′ is a block diagonal matrix, and the WLS method has no contribution to the amplitude estimation accuracy; on this basis, the estimation cost function J(η) is derived; Where w0 is the fundamental frequency, α k is the amplitude of the kth harmonic signal, K is the highest order, φ k is the initial phase of the kth harmonic, is the estimated phase of the kth harmonic, N is the data length, find the value about φ k Take the partial derivative and set the result to 0, and we get Substitute the result into J(η) and find the partial derivative of the final estimated target w0; Where SNR k is the local signal-to-noise ratio of the line spectrum obtained based on the kth frequency point in step S1 Step S502: based on the fundamental frequency f obtained in step S4 F1 The corresponding frequency, frequency multiple, order and the line spectrum local signal-to-noise ratio after background equalization based on the kth frequency point obtained in step S1 Calculate the order of the fundamental frequency: For the extracted line spectrum In step S402, i,n -K|<Δ,I=1,2,…,N1; n=1,2,…N S The required line score was updated and reorganized as The corresponding line spectra of the other fundamental frequencies are rearranged in the same way; Step S503: Find the line spectra in and And take The results are given as estimated results; Step S504: The estimated value f in step S503 is F1,WLS Reconfirmation will be Internal spectrum frequency and f F1,WLS Integer multiples of . Among them, F mid is the dividing frequency between low frequency and high frequency. The number of line spectra that satisfy the above formula is the fundamental frequency f F1,WLS The quality factor q n ; When q is satisfied n >N e ,N e ≥3, then The fundamental frequency f F1,WLS precise.

6. The platform multi-harmonic self-noise interference parameter estimation method according to claim 1, characterized in that: Step S6, performing sequential harmonic fundamental frequency estimation based on the harmonic group center frequency obtained based on the histogram statistics according to the corrected fundamental frequency, thereby obtaining a more robust harmonic group estimation effect; Step S601, loop through steps S2 to S5 until the number of frames continuously estimated is greater than the number of frames P required for the received signal to be stable. Assume that the platform noise of P frames is estimated to obtain the baseband estimation result of P frames {f F,WLS }={f F1,WLS ,f F2,WLS ,...f Fp,WLS }; Take the maximum fundamental frequency f of the result F,WLSMax and the minimum value f F,WLSMin , and the interval [f F,WLSMax ,f F,WLSMin ] is divided into P-1 parts; Step S602: Take the center frequency of the interval with the highest probability of occurrence of the estimated fundamental frequency of the previous P frame as the robust estimated fundamental frequency result f of the P frame. F,Final ; When the current frame number is greater than P, f F,WLS The fundamental frequency result f of the previous frame data of the current estimation frame is robustly estimated F,Final , Based on the result of weighted least squares estimation of the current frame, a histogram-based sequential statistics is performed to perform robust estimation of the fundamental frequency; The estimated result of the previous P frame of the current frame falls within the previous robust f F,WLS is the probability of the interval with the center frequency, The robust f in the P frame falls in the front sequence F,WLS is the number of frames in the interval of the center frequency, The estimated base frequency result of the current estimated P frame falls within the current frame estimate is the probability of the interval with the center frequency, Estimated for the current frame The number of frames in the interval of the center frequency is calculated by executing the following formula to estimate the robust base frequency f of the current frame: F,Final renew: where p L and p H are the probabilities for evaluating the robustness of fundamental frequency estimation.

7. A computer device, characterized in that: The computer device includes a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, the steps of the platform multi-harmonic self-noise interference parameter estimation method as described in any one of claims 1 to 6 are implemented.

8. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores a computer program for executing the platform multi-harmonic self-noise interference parameter estimation method according to any one of claims 1 to 6.

Images