Iterative methods, media, and devices for correcting power spectrum sampling ambiguity.

Abstract

The present invention relates to a correction method, medium and device for power spectrum sampling confusion based on an iteration method, and belongs to the technical field of power spectrum calculation of ocean atmospheric process. The method obtains a time series of a physical process of interest, performs Fourier transform, obtains an initial power spectrum, removes spectrum peaks that significantly deviate from the spectrum slope, calculates the power spectrum slope by using the spectrum values ​​of the remaining frequency bands and a linear fit under exponential coordinates, and estimates the power spectrum slope at the sampling frequency and the Nyquist frequency ω according to the calculated power spectrum slope. N Calculate the spectrum value corresponding to the frequency ω outside the Nyquist frequency ω N The corresponding frequency after folding is finally calculated to obtain the corrected spectrum value. The present invention can correct the power spectrum feature estimation distortion caused by sampling confusion without using higher frequency sampling data.

Description

Technical Field

[0001] The present invention belongs to the technical field of ocean-atmosphere process power spectrum calculation, and in particular relates to a correction method, medium and equipment for power spectrum sampling confusion based on an iteration method. Background Art

[0002] In the field of digital signal processing, the frequency of half the sampling rate is called the Nyquist frequency. When the discrete sampling method is used to acquire a signal, the components with frequencies higher than the Nyquist frequency will be mistakenly identified as low-frequency components symmetrical to the Nyquist frequency. This misdiagnosed signal is called an aliased signal, which may distort the initial spectral distribution characteristics of the signal and obscure the physical laws behind it.

[0003] Geophysical processes usually follow a spectral slope from large-scale to small-scale energy levels. Affected by sampling aliasing energy, the spectral slope estimated from discrete sampling data often deviates from the actual spectral slope. As geophysical research shrinks in time and space, and the sampling frequency often cannot be further increased due to current technology and storage limitations, the impact of sampling aliasing on the physical process under study becomes more and more significant. For this reason, predecessors often regard the spectral values ​​close to one end of the Nyquist frequency as noise and only analyze the characteristics of the low-frequency process. However, there are obvious artificial factors in how to clearly define the range near the Nyquist frequency, which may cause deviations in the research conclusions. And as the scale of interest in the study approaches the Nyquist frequency scale, traditional practices are no longer feasible.

[0004] In atmospheric turbulence observations, instrument manufacturers use a pre-given correction function in the spectral domain to correct the spectral value deviation caused by sampling aliasing. However, this correction function often deviates from the true spectral value characteristics and may even further aggravate the spectral value aliasing error.

[0005] Wang et al. proposed in 2022 that sampling can be performed at a higher frequency, and then the results can be saved and analyzed at a sampling frequency that the storage conditions can withstand. The above approach can effectively reflect the power spectrum characteristics near the Nyquist end. However, in most cases, it is difficult to further increase the sampling frequency, so the effectiveness of this method is also limited. Summary of the invention

[0006] In view of the above problems, the present invention proposes a method, medium and device for correcting power spectrum sampling confusion based on an iterative method, which corrects the power spectrum feature estimation distortion caused by sampling confusion without using higher frequency sampling data.

[0007] The present invention achieves the above-mentioned object through the following schemes.

[0008] A method for correcting power spectrum sampling confusion based on an iterative method, the method comprising the following steps:

[0009] Step 1: Obtain the time series f(t) of the physical process of interest, perform Fourier transform F(ω), and obtain the initial power spectrum P raw (ω):

[0010]

[0011] Where N is the data sampling length, ω N is the Nyquist frequency;

[0012] Step 2: Eliminate the spectral peaks that significantly deviate from the spectral slope;

[0013] Step 3: Use the spectrum values ​​of the remaining frequency bands to perform a linear fit in the exponential coordinates to estimate the spectrum value P fit (f) = a + k n f calculates the power spectrum slope k n , where n represents the nth correction cycle, and f is the frequency of the exponential coordinate, that is, f = log 10 (ω); first-order linear fitting parameters a and k n The calculation formula is as follows:

[0014]

[0015]

[0016] Among them, m represents the data length with a total of m data points, i represents the i-th data point, and f i represents the frequency value of the i-th point in the exponential coordinate system, P i Represents the power spectrum value of the i-th point in the exponential coordinate system;

[0017] Step 4: Estimate the sampling frequency and Nyquist frequency ω according to the calculated power spectrum slope N The fitted estimated spectrum value P corresponding to the frequency ω other than fit (ω);

[0018] Step 5: Calculate each frequency ω according to the Nyquist frequency ω N The corresponding frequency after folding is ω a , where ω a =ω-2nω N where n is the value that makes ω a Belong to [-ω N ,ω N ];

[0019] Step 6: Use the initial power spectrum P raw(ω) and the frequency domain correction coefficient Φ(ω) are used to calculate the corrected spectrum value P x (ω):

[0020]

[0021]

[0022] P am (ω) is the estimated spectrum value P fit (f) with ω N The folded spectrum value:

[0023]

[0024] Step 7: Use the convergence standard to determine whether the power spectrum before and after correction converges, that is, to determine whether the newly obtained corrected spectrum value P x The spectral slope k of (ω) n The corrected spectrum value P obtained in the previous step x The spectral slope k of (ω) n-1 Is the difference less than H? If not, use the new correction spectrum P x (ω), repeat steps 2 to 7; if converged, output the spectrum value, and use the output spectrum value to manage the physical process.

[0025] Furthermore, H in step seven is 0.05.

[0026] Furthermore, the physical process is characterized in that the power spectrum value is distributed according to a negative slope after taking the logarithm.

[0027] The present invention also provides a computer-readable storage medium, wherein the computer-readable storage medium stores a computer program, and the computer program is suitable for being loaded by a processor and executing the method for correcting power spectrum sampling confusion based on an iteration method.

[0028] The present invention also provides a computer device, which includes a memory and a processor, wherein the memory stores a computer program, and when the computer program is executed by the processor, the power spectrum sampling confusion correction method based on the iteration method is run.

[0029] The beneficial effects of the present invention compared with the prior art are as follows: the method only relies on the assumption that the observed physical process extends above the Nyquist frequency with a certain spectrum slope and the known power spectrum results, and can correct the error caused by sampling aliasing in the power spectrum, and give the true spectrum domain characteristics of the physical process near the Nyquist frequency band. The method does not require the prior knowledge of the high-frequency (above the Nyquist frequency) spectrum characteristics (such as slope, etc.) of the physical process, nor does it require an increase in the sampling frequency, and therefore has a wide range of practical value. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] Figure 1 The figure is the atmospheric turbulent velocity time series f(t), where the upper figure is the 5Hz sampling result and the lower figure is the 100Hz sampling result;

[0031] Figure 2 is the initial power spectrum of atmospheric turbulence P raw (ω) graph;

[0032] Figure 3 is the initial power spectrum P of atmospheric turbulence 5Hz sampling data raw (ω) and its fitting result diagram;

[0033] Figure 4 is the initial power spectrum P of atmospheric turbulence 5Hz sampling data raw (ω), power spectrum fitting results and high-frequency estimated spectrum value diagram;

[0034] Figure 5 is the initial power spectrum P of atmospheric turbulence 5Hz sampling data raw (ω), power spectrum fitting results and high-frequency estimated spectrum value diagram, where “●” and “+” represent two groups of frequencies and their corresponding folding frequencies;

[0035] Figure 6 is the initial power spectrum P of atmospheric turbulence 5Hz sampling data raw (ω) and the first corrected output power spectrum P x (ω) graph;

[0036] Figure 7 is the initial power spectrum P of atmospheric turbulence 5Hz sampling data raw (ω), and the power spectra P obtained from the first, third, and fifth correction cycles x (ω) graph;

[0037] Figure 8 The power spectrum of atmospheric turbulence 100Hz sampling data, 5Hz sampling power spectrum, and the final corrected output power spectrum P x (ω) graph;

[0038] Fig. 9 is the high-frequency observation power spectrum, low-frequency observation power spectrum and low-frequency observation corrected power spectrum P x (ω) graph;

[0039] Fig.10 Flow chart of the method of the present invention. DETAILED DESCRIPTION

[0040] The technical solution of the present invention is further explained below by way of embodiments, but the protection scope of the present invention is not limited in any form by the embodiments.

[0041] The present invention hereinafter only takes the physical processes of atmospheric turbulence and buoy velocity as examples to illustrate the technical solution of the present invention, but the method of the present invention is not limited to these two embodiments. Any power spectrum value that conforms to the physical process and is distributed according to a negative slope after taking the logarithm can be corrected by the method of the present invention.

[0042] Example 1: Atmospheric turbulence power spectrum correction.

[0043] A correction method for power spectrum sampling confusion based on an iterative method, such as Fig.10 As shown, the method comprises the following steps:

[0044] Step 1: Obtain the time series of atmospheric turbulent velocity f(t), such as Figure 1 As shown, Fourier transform F(f) is performed to obtain the initial power spectrum P raw (ω), the formula is as follows:

[0045]

[0046] Where N is the data sampling length, ω N is the Nyquist frequency;

[0047] like Figure 2 As shown, it can be seen that the power spectrum of the 5Hz sampling data deviates significantly from the power spectrum of the 100Hz sampling data, indicating that the power spectrum of the 5Hz sampling data is distorted due to sampling confusion, and its power spectrum needs to be corrected.

[0048] Step 2: Eliminate the spectral peaks that significantly deviate from the spectral slope. There are no spectral peaks that significantly deviate from the spectral slope in the data of this embodiment, so this step is ignored.

[0049] Step 3: Use the spectrum values ​​of the remaining frequency bands to perform a linear fit in the exponential coordinates to estimate the value P fit (f) = a + k n f calculates the power spectrum slope k n , where n represents the nth correction cycle, n is set to 0 in the initial cycle, and f is the frequency of the exponential coordinate, that is, f = log 10 (ω); first-order linear fitting parameters a and k n The calculation formula is as follows:

[0050]

[0051]

[0052] Among them, m represents the data length with a total of m data points, i represents the i-th data point, and f i represents the i-th frequency value in the exponential coordinate system, P irepresents the ith power spectrum value in the exponential coordinate system; the initial power spectrum P of the atmospheric turbulence 5Hz sampling data raw (ω) and its fitting results are shown in Figure 3 As shown, the power spectrum slope k0 = -1.16;

[0053] Step 4: Estimate the Nyquist frequency ω according to the power spectrum slope fitting result calculated in step 3 N The fitted estimated spectrum value P corresponding to the frequency ω other than fit (f); According to the fitting results, the spectrum values ​​at higher frequencies than the sampling frequency are estimated as follows Figure 4 shown.

[0054] Step 5: Calculate each frequency ω according to the Nyquist frequency ω N The corresponding frequency after folding is ω a , where ω a =ω-2nω N where n is the value that makes ω a Belong to [-ω N ,ω N ] integer; for each point of the atmospheric turbulence power spectrum, calculate its folding frequency and the power spectrum at the folding frequency, the result is as follows Figure 5 shown.

[0055] Step 6: Subtract the corrected spectrum value P obtained by folding each frequency from the initial power spectrum x (ω);

[0056]

[0057] Where Φ(ω) is the frequency domain correction coefficient, which is the ratio of each frequency minus the corrected spectrum value

[0058]

[0059] P am (ω) is the estimated spectrum value P fit (ω) N The folded spectrum value,

[0060]

[0061] The first corrected spectrum value P is obtained by subtracting the spectrum value at the corresponding folding frequency from the initial power spectrum of atmospheric turbulence. x The result of (ω) is as follows: Figure 6 As shown by the black line, k1 is equal to -0.57 at this time.

[0062] Step 7: Use the convergence standard to determine whether the power spectrum before and after correction converges, that is, to determine whether the newly obtained corrected spectrum value P x The spectral slope k of (ω) nThe corrected spectrum value P obtained in the previous step x The spectral slope k of (ω) n-1 Is the difference less than H? If not, use the new corrected spectrum value P x (ω), repeat steps 2 to 7; if converged, output the corrected spectrum value.

[0063] In this application case, H in step 7 is selected as 0.05. k1 is equal to -0.57, which does not meet the convergence condition, so it needs to be corrected again. The spectral slope k of the next 2-5 cycles is n They are k2=-0.81, k3=-0.11, k4=-1.20, k5=-1.24 respectively. According to the convergence standard, the power spectrum output after the fifth cycle is the power spectrum result after the correction is completed. The convergence process is as follows: Figure 7 As shown. Figure 8 As shown in the figure, the power spectrum characteristics of the final corrected output are basically consistent with the high-frequency sampling power spectrum characteristics, indicating that this method corrects the power spectrum distortion caused by sampling confusion and increases the credibility of the results near the Nyquist frequency end of the power spectrum.

[0064] Example 2

[0065] The method described in Example 1 is also used for the correction of the power spectrum of the buoy velocity. Example 2 is the correction of the power spectrum of the buoy velocity. The correction result is as follows: Fig. 9 shown. Fig. 9 The middle gray solid line is the power spectrum of the horizontal flow velocity data observed at 0.5 second intervals, the gray dotted line is the power spectrum of the horizontal flow velocity data observed at 1 hour intervals, and the black solid line is the corrected power spectrum of the horizontal flow velocity data observed at 1 hour intervals. It can be seen that the spectral slope is similar to the power spectrum of the horizontal flow velocity data observed at 0.5 second intervals.

Claims

1. A method for correcting power spectrum sampling confusion based on an iterative method, characterized in that: The method comprises the following steps: Step 1: Obtain the time series f(t) of the physical process of interest, perform Fourier transform F(ω), and obtain the initial power spectrum P raw (ω): Where N is the data sampling length, ω N is the Nyquist frequency; Step 2: Eliminate the spectral peaks that significantly deviate from the spectral slope; Step 3: Use the spectrum values ​​of the remaining frequency bands to perform a linear fit in the exponential coordinates, P fit (f) = a + k n f calculates the power spectrum slope k n , where n represents the nth correction cycle, and f is the frequency of the exponential coordinate, that is, f = log 10 (ω); first-order linear fitting parameters a and k n The calculation formula is as follows: Among them, m represents the data length with a total of m data points, i represents the i-th data point, and f i represents the frequency value of the i-th point in the exponential coordinate system, P i Represents the power spectrum value of the i-th point in the exponential coordinate system; Step 4: Estimate the sampling frequency and Nyquist frequency ω according to the calculated power spectrum slope N The fitted estimated spectrum value P corresponding to the frequency ω other than fit (ω); Step 5: Calculate each frequency ω according to the Nyquist frequency ω N The corresponding frequency after folding is ω a , where ω a =ω-2nω N where n is the value that makes ω a Belong to [-ω N ,ω N ]; Step 6: Use the initial power spectrum P raw (ω) and the frequency domain correction coefficient Φ(ω) are used to calculate the corrected spectrum value P x (ω): P am (ω) is the estimated spectrum value P fit (f) with ω N The folded spectrum value: Step 7: Use the convergence standard to determine whether the power spectrum before and after correction converges, that is, to determine whether the newly obtained corrected spectrum value P x The power spectrum slope k of (ω) n The corrected spectrum value P obtained in the previous step x The power spectrum slope k of (ω) n-1 Is the difference less than H? If not, use the new correction spectrum P x (ω), repeat steps 2 to 7; if converged, output the corrected spectrum value, and use the output spectrum value to manage the physical process.

2. The correction method according to claim 1, characterized in that: H in step 7 is 0.

05.

3. The correction method according to claim 1, characterized in that: The physical process described in step 1 is that the power spectrum value is distributed according to a negative slope after taking the logarithm.

4. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores a computer program, and the computer program is suitable for being loaded by a processor and executing the method for correcting power spectrum sampling confusion based on an iteration method as described in any one of claims 1-3.

5. A computer device, characterized in that: The device includes a memory and a processor, wherein the memory stores a computer program, and when the computer program is executed by the processor, the method for correcting power spectrum sampling confusion based on the iteration method as described in any one of claims 1-3 is executed.

Images