Logarithmic scale adaptive multi-window power spectral density estimation method

Abstract

According to the logarithmic scale adaptive multi-window power spectral density estimation method provided by the invention, the conicity degree is adjusted through the equilibrium relationship between the deviation and the variance, the problems of low resolution and high algorithm complexity of the existing adaptive sine window algorithm are solved, and the improvement of the adaptive sine window algorithm is realized. The method comprises the following steps: firstly, initializing conicity degrees, setting the conicity degrees at all frequencies as fixed values, and calculating a power spectral density initial value; secondly, cone degree gradient processing is carried out on low-frequency signals, the cone degree is smoothed through filters with different intensities at different frequency bands, and then the optimal cone degree of the whole frequency domain is obtained; and finally, carrying out iterative operation on the conicity degree and the power spectrum density to obtain an optimal power spectrum density estimation value. Compared with an LPSD algorithm, an adaptive sine window algorithm and an LISA-LPSD algorithm, the method provided by the invention has the advantages of high spectrum resolution, low algorithm complexity and the like.

Description

Technical Field

[0001] The present invention belongs to the field of gravitational wave detection, and in particular relates to a logarithmic scale adaptive multi-window power spectral density estimation method. Background Art

[0002] In order to ensure the implementation of the low-frequency (0.1mHz-1Hz) detection project in space gravitational wave detection, spectral density estimation algorithms are widely used in detection system evaluation. The spectral density estimation algorithms mainly include the LPSD (logarithmic frequency axis power spectral density) algorithm, the LISA-LPSD (LISA improved logarithmic frequency axis power spectral density) algorithm included in the toolbox developed by the LISA team, and the adaptive sine window algorithm. The LPSD algorithm and the LISA-LPSD algorithm are widely used in the detection system evaluation. The high-frequency end close to 1Hz is processed with a wider resolution, and the low-frequency end is processed with a narrower resolution to adjust the spectrum details, but there is still a problem of low resolution of low-frequency signals and high computational complexity.

[0003] The existing adaptive sine window function adopts an optimal taper number acquisition method for the full frequency domain. , for the cost function ( The best taper number, is the true power spectral density estimate, is the second-order derivative, T represents the time length of the window function, For power spectral density estimation), the variance and deviation are equally weighted and balanced, the number of tapers is large, so the resolution is low and the running time is relatively long, so it has never been applied in this field. Summary of the invention

[0004] In order to overcome the shortcomings of existing algorithms and address the problems of low resolution, large taper number and relatively long running time, the present invention provides a logarithmic scale adaptive multi-window power spectral density estimation method, which is an adaptive sine window method with improved taper number. The taper number is adjusted by the balanced relationship between deviation and variance, and the taper number gradient processing is performed on the low-frequency signal. Finally, by comparing with the LPSD algorithm, the original adaptive sine window algorithm and the LISA-LPSD algorithm, it is shown that the method has the advantages of high spectral resolution and low algorithm complexity.

[0005] In order to achieve the above object, the present invention adopts the following technical scheme:

[0006] A logarithmic scale adaptive multi-window power spectral density estimation method comprises the following steps:

[0007] Step 1: Initialize the taper number. The taper number at all frequencies is 20. Zero-mean real-valued stationary random process ,in, Indicates the time The random process value when is a vector containing Quantity, It's in time hour independent random processes, the integration interval is , is the upper limit of the frequency range, defining the adaptive multi-window power spectral density estimate ;

[0008] Step 2: Estimation of the adaptive multi-window power spectral density defined in step 1 to process;

[0009] Step 3: When estimating the power spectral density of step 1, outliers are eliminated when the taper number suddenly changes to infinity. At the same time, linear interpolation is used to filter and smooth the taper number of step 2 for three iterations, while ensuring that when the frequency is greater than 0.1 Hz, the difference in the taper number of adjacent frequencies is less than or equal to 1.

[0010] Beneficial effects:

[0011] The present invention first performs the first two steps of initial power spectrum density estimation and optimizes the different frequency bands of the taper number solution. Adjustment The value is taken, a smaller taper number is used for lower frequencies, and taper number strong and weak filtering and smoothing are performed for different frequency ranges. The last third step is to iterate the taper number strong and weak filtering and smoothing, and iterative calculation of power spectrum density estimation. The single optimal taper number used in the full frequency domain of the original adaptive sine window algorithm is innovatively improved to achieve the solution of power spectrum density. The present invention solves the pain points of existing algorithms in this application field and the common difficulties of existing algorithms in the field. While improving the spectral resolution, it reduces the algorithm complexity and speeds up the calculation time. BRIEF DESCRIPTION OF THE DRAWINGS

[0012] Figure 1 A flow chart of a logarithmic scale adaptive multi-window power spectral density estimation method of the present invention;

[0013] Figure 2 When the external sine signal is stimulated by 1mHz, the comparison diagram of the noise amplitude spectrum density algorithm of the mechanism angle measurement (the upper figure is the comparison diagram of the noise amplitude spectrum density estimation of the theoretical PSD, LPSD algorithm, the algorithm of the present invention, the LISA-LPSD algorithm, and the adaptive sine window algorithm, and the lower figure is the comparison diagram of the full frequency domain distribution of the optimal taper number of the algorithm of the present invention and the adaptive sine window algorithm);

[0014] Figure 3When the 1mHz external sinusoidal signal is used for stimulation, the algorithm fitting variance, deviation and mean square error diagram (variance, deviation and mean square error of amplitude spectrum density estimation of LPSD algorithm, algorithm of the present invention, LISA-LPSD algorithm and adaptive sinusoidal window algorithm);

[0015] Figure 4 When the 1mHz external sinusoidal signal is used for stimulation, the algorithm running time diagram for different data collection durations;

[0016] Figure 5 For the organization to swing to , the comparison diagram of the noise amplitude spectrum density algorithm of the mechanism angle measurement (the upper figure is the comparison diagram of the noise amplitude spectrum density estimation of the theoretical PSD, LPSD algorithm, the algorithm of the present invention, the LISA-LPSD algorithm, and the adaptive sine window algorithm, and the lower figure is the comparison diagram of the full frequency domain distribution of the optimal taper number of the algorithm of the present invention and the adaptive sine window algorithm);

[0017] Figure 6 For the organization to swing to When , the algorithm fitting variance, deviation and mean square error diagram (LPSD algorithm, the algorithm of the present invention, LISA-LPSD algorithm, and adaptive sine window algorithm amplitude spectrum density estimation variance, deviation and mean square error);

[0018] Figure 7 For the organization to swing to Algorithm running time diagram with different data collection durations. DETAILED DESCRIPTION

[0019] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention. In addition, the technical features involved in the various embodiments of the present invention described below can be combined with each other as long as they do not conflict with each other.

[0020] like Figure 1 As shown, a logarithmic scale adaptive multi-window power spectrum density estimation method of the present invention comprises the following steps:

[0021] Step 1: Initialize the taper number. The taper number at all frequencies is 20. Zero-mean real-valued stationary random process ( Indicates the time The random process value at time , where is a vector containing Quantity, It's in time hour independent random processes), the integration interval is ( is the upper limit of the frequency range), and the adaptive multi-window power spectral density estimation is defined as:

[0022] (1)

[0023] in, Indicates the frequency The power spectral density estimate at is the weighting factor corresponding to the window being processed by window i, is the periodogram obtained with the selected data window, is the optimal taper number.

[0024] The calculation formula is:

[0025] (2)

[0026] Here, K is the number of cone periodograms after averaging, expressed as:

[0027] (3)

[0028] in, Represents the overall window, is a sinusoidal window function, r represents the sample sequence number from 1 to L, L is the number of observations, represents the sampling interval, j is an imaginary unit, and is defined as follows:

[0029] (4)

[0030] In order to ensure the normalization of frequency domain energy, the overall window function The following conditions are met:

[0031] (5)

[0032] in, , , T represents the time length of the window function.

[0033] Step 2: Adaptive multi-window power spectral density estimation defined in step 1 Perform the following processing:

[0034] When the frequency is less than 1mHz, the taper number is set to 1. A strong filter is used to smooth the taper number between 1mHz and 0.01Hz. When the frequency is greater than 0.01Hz, a parabola weighting is used at each frequency to obtain the optimal taper number. , and use a weak filter to smooth the taper number from 0.01Hz to 1Hz.

[0035] Analyze the deviation between the power spectral density estimate and the true spectral density estimate:

[0036] (6)

[0037] Among them, E is the expected value, f is the frequency, and the variable replacement , also satisfies the normalization condition: , then the deviation can be expressed as:

[0038] (7)

[0039] In the analysis When , the window function In the time domain there are The window function frequency response main lobe width is proportional to the ratio of the window length T, following the time domain uncertainty principle, so (W is the main lobe width) To simplify the analysis, since the frequency domain capability of the sine window is concentrated in a narrow range, it is assumed that the sine window The corresponding frequency domain window function Small area The inner distribution is close to uniform:

[0040] (8)

[0041] in, is the true power spectral density estimate It changes slowly and can be expressed approximately as:

[0042] (9)

[0043] Therefore, the deviation Further simplified to:

[0044] ;

[0045] because :

[0046] (10)

[0047] Power spectral density estimation variance calculation requires covariance calculate:

[0048] ;

[0049] in, are the power spectral density estimates of the same signal x under m and n windows respectively. Assuming that the correlation between adaptive sine window spectrum estimates is small, the variance of each periodogram is approximately , averaging the periodograms, the variance is reduced to , then further simplify

[0050] (11)

[0051] Based on the mean square error criterion of the adaptive sine window algorithm, the optimal algorithm evaluation criterion is proposed. Adjusting the relationship between deviation and variance is essentially a matter of adjusting the number of tapers. Then the cost function

[0052] (12)

[0053] in, is a weight parameter used to balance the influence of bias and variance on the cost function. The mean square error MSE is:

[0054] (13)

[0055] To calculate the expected value, let MSE = 0. It can be seen that:

[0056] (14)

[0057] Parameters for adjusting the taper number of sub-bass frequencies The value of is calculated according to the dichotomy method, preferably, .

[0058] Step 3: When estimating the power spectrum density of step 1, outliers should be removed when the taper number suddenly changes to infinity, and linear interpolation should be used. The taper number filtering and smoothing of step 2 should be iterated three times to solve, and it should be ensured that when the frequency is greater than 0.1 Hz, the difference between the taper numbers of adjacent frequencies is less than or equal to 1. The power spectrum density estimation of step 1 is iterated 15 times, and logarithmic operations are performed to obtain the logarithmic scale optimal power spectrum density estimate and logarithmic scale amplitude spectrum density estimate.

[0059] In addition, the algorithm complexity formula of the LPSD algorithm itself is O(N J), when J increases, the operation time increases greatly. The complexity of the LPSD algorithm improved by the LISA project is O(N log2(N / K)+J)+O(M), the algorithm complexity is reduced and the running time is shortened. However, calling the toolbox O(M) increases the algorithm complexity and takes a long time. (The signal length is , the number of segments is , the frequency points are , M is the number of all .jar files loaded and output). The complexity of the adaptive sine window algorithm is O(5N log2(N)), which is obviously much smaller than the complexity of the LISA-LPSD algorithm and runs faster.

[0060] Example:

[0061] The embodiment of the present invention includes the following steps:

[0062] like Figure 1 As shown, a logarithmic scale adaptive multi-window power spectrum density estimation method of the present invention comprises the following steps:

[0063] Step 1: Use the data as the adaptive multi-window power spectral density estimation input. First, according to the initial taper number 20, the adaptive multi-window power spectral density estimation is obtained according to formula (1): ,Right now Figure 1 Medium PSD.

[0064] (1)

[0065] Step 2: Adaptive multi-window power spectral density estimation defined in step 1 Perform the following processing:

[0066] When the frequency is less than 1mHz, the taper number is set to 1. When the frequency is between 1mHz and 0.01Hz, a strong filter is used to smooth the taper number. When the frequency is greater than 0.01Hz, a parabola weighting is used at each frequency to obtain the optimal taper number according to formula (14). , and use a weak filter to smooth the taper number from 0.01Hz to 1Hz.

[0067] (14)

[0068] Step 3: When estimating the power spectrum density of step 1, outliers should be removed when the taper number suddenly changes to infinity. At the same time, linear interpolation is used. The taper number filtering and smoothing of step 2 are iterated 3 times to solve. At the same time, it is necessary to ensure that the difference between the taper numbers of adjacent frequencies greater than 0.1Hz is less than or equal to 1. The power spectrum density estimation of step 1 is iterated 15 times to obtain the optimal power spectrum density estimation value and amplitude spectrum density estimation value. The amplitude spectrum density estimation effect of the angle value is compared with the LPSD algorithm, LISA-LPSD algorithm and the original adaptive sine window algorithm through the improved adaptive sine window algorithm, and the algorithm running time is recorded.

[0069] Specifically, the present invention is applied to the evaluation of the pointing noise of a gravitational wave laser pointing mechanism. When a sinusoidal stimulus is applied, the angle measurement noise characteristics of the mechanism are as follows: Figure 2 , Figure 3As shown in the figure, near 1mHz, the mean square error of the algorithm of the present invention is reduced compared with the adaptive sine window algorithm, and the mean square error is greatly reduced compared with the LPSD algorithm. Taking 1mHz as the measurement standard, the mean square error of the algorithm of the present invention is reduced by 82% compared with the adaptive sine window algorithm. Compared with the LPSD algorithm, the mean square error is reduced by 90%, which is an order of magnitude, and the performance is significantly improved. The mean square error of the LISA-LPSD algorithm is reduced by 55% compared with the LPSD algorithm, and the mean square error of the algorithm of the present invention is reduced by 79% compared with the LISA-LPSD algorithm, and the performance is greatly improved. Figure 2 The upper figure is a comparison diagram of noise amplitude spectrum density estimation obtained by performing square root and logarithm operations on the theoretical PSD, LPSD algorithm, the algorithm of the present invention, the LISA-LPSD algorithm, and the adaptive sine window algorithm. The lower figure is a comparison diagram of the full frequency domain distribution of the optimal taper number of the algorithm of the present invention and the adaptive sine window algorithm. Figure 3 In the figure, square root and logarithm operations are performed for the LPSD algorithm, the algorithm of the present invention, the LISA-LPSD algorithm, and the adaptive sine window algorithm to obtain the variance, deviation, and mean square error of the amplitude spectral density estimation.

[0070] like Figure 4 As shown in the figure, when the data sampling time is less than 1 hour, the running time of the adaptive sine window algorithm is not much different from that of the algorithm of the present invention, but as the data sampling time increases, the advantage of the LISA-LPSD algorithm becomes more prominent. Considering the influence of temperature and other factors on the experiment, the longest data sampling time is 5 hours. Based on this, the algorithm operation rate is evaluated. Compared with the adaptive sine window algorithm, the operation time of the algorithm of the present invention is reduced to 47%, compared with the LPSD algorithm, the operation time is reduced to 11%, and compared with the LISA-LPSD algorithm, the operation time is reduced to 65%.

[0071] Similarly, the organization is biased towards When the mechanism angle measurement noise characteristics are as follows Figure 5 , Figure 6 As shown, the mean square error of the algorithm of the present invention is reduced compared with the adaptive sine window algorithm, and the mean square error is also significantly smaller than that of the LPSD algorithm. Taking 1mHz as the measurement standard, the mean square error of the algorithm of the present invention is reduced by about 45% compared with the adaptive sine window algorithm. In addition, compared with the LPSD algorithm, the mean square error is reduced by about 54%, and the performance is significantly improved. The mean square error of the LASA-LPSD algorithm is reduced by about 24% compared with the LPSD algorithm. Furthermore, the mean square error of the algorithm of the present invention is reduced by about 39% compared with the LASA-LPSD algorithm, and the performance is significantly improved. Figure 5 The upper figure is a comparison diagram of noise amplitude spectrum density estimation of theoretical PSD, LPSD algorithm, algorithm of the present invention, LISA-LPSD algorithm, and adaptive sine window algorithm, and the lower figure is a comparison diagram of full frequency domain distribution of optimal taper number of the algorithm of the present invention and the adaptive sine window algorithm. Figure 6 The variance, deviation and mean square error of the amplitude spectrum density estimation of the LPSD algorithm, the algorithm of the present invention, the LISA-LPSD algorithm and the adaptive sine window algorithm are shown in the figure.

[0072] At the same time, it can be obtained from FIG. 7 that the operation time of the algorithm of the present invention is reduced to 33% compared with the operation time of the adaptive sine window algorithm, reduced to 14% compared with the LPSD algorithm, and reduced to 70% compared with the LISA-LPSD algorithm.

[0073] Therefore, based on the two experimental simulations, the algorithm of the present invention has a reduced mean square error, reduced algorithm complexity, and faster operation speed compared with the adaptive sine window algorithm, LPSD algorithm, and LISA-LPSD algorithm. It should be emphasized that compared with the LISA-LPSD algorithm, its mean square error is reduced by about 40%-80%, which is about 20%-60% of its original mean square error. When generating steep signals such as sinusoidal noise, the improved adaptive sine window algorithm has a better performance, and the mean square error is reduced by even 20%. In contrast, the mean square error is reduced by nearly 40% in a relatively flat state. The running time is reduced to 65%-70%, and even the running time is reduced by an order of magnitude when processing less data. Overall, there is no large fluctuation and the performance is stable.

[0074] The adaptive sine window algorithm with improved taper number at low frequency can achieve better analysis results and higher resolution when applied in this field. Compared with the LISA-LPSD algorithm, the algorithm complexity is reduced and the running time is faster.

Claims

1. A logarithmic scale adaptive multi-window power spectral density estimation method, characterized in that: The steps include: Step 1: Initialize the taper number. The taper number at all frequencies is 20. Zero-mean real-valued stationary random process ,in, Indicates the time The random process value when is a vector containing Quantity, It's in time hour independent random processes, the integration interval is , is the upper limit of the frequency range, defining the adaptive multi-window power spectral density estimate ; Step 2: Estimation of the adaptive multi-window power spectral density defined in step 1 to process; Step 3: When estimating the power spectral density of step 1, outliers are eliminated when the taper number suddenly changes to infinity. At the same time, linear interpolation is used to filter and smooth the taper number of step 2 for three iterations, while ensuring that when the frequency is greater than 0.1 Hz, the difference in the taper number of adjacent frequencies is less than or equal to 1.

2. A logarithmic scale adaptive multi-window power spectral density estimation method according to claim 1, characterized in that: In step 1, adaptive multi-window power spectral density estimation for: （1） in, represents the power spectral density estimate at frequency f, is the weighting factor corresponding to the window being processed by window i, is the periodogram obtained with the selected data window, is the optimal taper number.

3. A logarithmic scale adaptive multi-window power spectral density estimation method according to claim 2, characterized in that: Window weighting factor The calculation formula is: （2） Where K is the number of cone periodograms after averaging, Represents the overall window function, which is expressed as: （3） Among them, r represents the sample sequence number from 1 to L, L is the number of observations, represents the sampling interval, j is the imaginary unit, The sine window function is defined as follows: （4） Where T represents the time length of the window function.

4. A logarithmic scale adaptive multi-window power spectral density estimation method according to claim 3, characterized in that: In order to ensure the normalization of frequency domain energy, the overall window function The following conditions are met: （5） Where f is the frequency, , , is the sampling frequency.

5. A logarithmic scale adaptive multi-window power spectral density estimation method according to claim 4, characterized in that: The step 2 comprises: When the frequency is less than 1mHz, the taper number is set to 1. A strong filter is used to smooth the taper number between 1mHz and 0.01Hz. When the frequency is greater than 0.01Hz, a parabola weighting is used at each frequency to obtain the optimal taper number. , and use a weak filter to smooth the taper number from 0.01Hz to 1Hz.

6. A logarithmic scale adaptive multi-window power spectral density estimation method according to claim 5, characterized in that: Analyze the estimated deviation between the power spectral density estimate and the true spectral density estimate for: （6） Among them, E is the expected value, f is the frequency, and the variable replacement , which also satisfies the normalization condition: ; The estimated bias It is expressed as: （7） ； Where W is the main lobe width.

7. A logarithmic scale adaptive multi-window power spectral density estimation method according to claim 6, characterized in that: Assuming a sinusoidal window function The corresponding frequency domain window function In a small area The distribution is close to uniform, that is: （8） in, is the true power spectral density estimate It changes slowly and can be expressed approximately as: （9） Therefore, the deviation is further simplified to: ； because , we get: （10）。 8. A logarithmic scale adaptive multi-window power spectral density estimation method according to claim 7, characterized in that: Perform covariance Calculation: ； in, They are the power spectral density estimates of the same signal x under m and n windows respectively; Assuming that the correlation between adaptive sine window spectrum estimates is small, the variance of each periodogram is approximately , averaging the periodograms, the variance is reduced to , then further simplifying it to get: （11） Based on the mean square error criterion of the adaptive sine window algorithm, the optimal algorithm evaluation criterion is proposed. Adjusting the relationship between deviation and variance is essentially a matter of adjusting the number of tapers. Then the cost function for: （12） in, is the weight parameter.

9. A logarithmic scale adaptive multi-window power spectral density estimation method according to claim 8, characterized in that: The mean square error MSE is: （13） To obtain the expected operation, let MSE = 0. get: （14） Parameters for adjusting the taper number of sub-bass frequencies The value of .

10. The logarithmic scale adaptive multi-window power spectral density estimation method according to claim 1, characterized in that: The power spectral density estimation in step 1 is iterated 15 times and solved, and logarithmic operations are performed to obtain the optimal logarithmic scale power spectral density estimation value and logarithmic scale amplitude spectral density estimation value.

Images