Atomic clock signal denoising method

Abstract

An atomic clock signal denoising method comprises the following steps: loading an atomic clock original signal and carrying out EMD decomposition; determining a smoothing factor of the IMF component; carrying out improved wavelet denoising on each IMF component; selecting an optimal wavelet basis for the IMF components by adopting a self-adaptive strategy; performing wavelet decomposition on each IMF component based on the optimal wavelet basis to obtain a corresponding wavelet coefficient vector; performing threshold processing on the wavelet coefficient vector to obtain a de-noised IMF component; during threshold processing, firstly determining an optimal threshold of a wavelet coefficient vector, then constructing an improved threshold function to process a wavelet coefficient in the wavelet coefficient vector, and then performing wavelet reconstruction to obtain a de-noised IMF component; the denoised IMF component and the residual component are superposed and reconstructed into a denoised signal; and carrying out moving average filtering on the reconstructed de-noised signal to obtain a final de-noised signal. According to the method, noise can be effectively suppressed, and the detail features and smoothness of the signals can be reserved to the maximum extent.

Description

Technical Field

[0001] This invention belongs to the field of signal processing technology, and in particular relates to a method for denoising atomic clock signals. Background Technology

[0002] Wavelet thresholding denoising possesses excellent time-frequency localization characteristics and multi-resolution analysis capabilities, making it one of the mainstream methods widely used in the field of atomic clock signal denoising. Wavelet thresholding denoising works by performing wavelet transform on the original signal, decomposing it into wavelet coefficients of different scales. Based on the differences in the distribution of noise and the true signal in the wavelet domain—for example, noise coefficients have small amplitudes and are widely distributed, while the true signal is concentrated in a few coefficients with large amplitudes—an appropriate threshold is set to filter and process the wavelet coefficients. Finally, the denoised signal is reconstructed through inverse wavelet transform, achieving noise suppression.

[0003] In traditional wavelet thresholding denoising methods, hard thresholding and soft thresholding are the two most representative thresholding methods. Hard thresholding retains coefficients with absolute values ​​greater than the threshold and sets coefficients less than the threshold to zero. This method preserves signal details well, but discontinuities at the threshold points can easily lead to pseudo-Gibbs phenomena in the reconstructed signal. Soft thresholding, on the other hand, shrinks coefficients greater than the threshold by a fixed amount, maintaining a constant deviation between the processed and original coefficients. This method ensures the overall continuity of the function and avoids oscillations caused by discontinuities, but the introduction of a constant deviation can lead to over-smoothing of the signal and significant loss of detail features. Summary of the Invention

[0004] The purpose of this invention is to provide a method for denoising atomic clock signals that can effectively suppress noise while reducing signal distortion and better restoring the original form and detailed features of the signal.

[0005] To achieve the above objectives, the present invention adopts the following technical solution:

[0006] A method for denoising atomic clock signals includes the following steps:

[0007] S1, Load the original atomic clock signal;

[0008] S2. Perform EMD decomposition on the original signal, decomposing it into a series of IMF components and a residual component;

[0009] S3. Determine the smoothing factor for each IMF component;

[0010] S4. Perform improved wavelet denoising on each IMF component to obtain the denoised IMF components. The steps are as follows:

[0011] S401. An adaptive strategy is used to select the optimal wavelet basis for each IMF component.

[0012] S402. Based on the optimal wavelet basis, perform wavelet decomposition on each IMF component to obtain the wavelet coefficient vector and length vector under the optimal wavelet basis for each IMF component.

[0013] S403. Threshold the wavelet coefficient vectors of each IMF component under the optimal wavelet basis to obtain the denoised IMF components.

[0014] In step S403, the thresholding process for the wavelet coefficient vector under the optimal wavelet basis is as follows:

[0015] S403-1. Determine the optimal threshold for the wavelet coefficient vector under the optimal wavelet basis;

[0016] S403-2. Construct an improved threshold function. Apply the improved threshold function to the wavelet coefficient vector under the optimal wavelet basis. The improved threshold function is as follows: In the formula This represents the improved threshold function, and sgn(·) represents the sign function. T' is the wavelet coefficient vector under the optimal wavelet basis. best The optimal threshold is δ, which is the auxiliary smoothing parameter, and α is the value of α. i Let be the smoothing factor for the i-th IMF component;

[0017] S403-3. For the IMF component after thresholding, the wavelet coefficient vector and length vector of the component obtained by wavelet decomposition in step S402, and the optimal wavelet basis selected in step S401 are used to perform wavelet reconstruction to restore the time domain signal, which is the denoised IMF component.

[0018] S5. Superimpose all the denoised IMF components with the residual components to reconstruct the denoised signal;

[0019] S6. Perform moving average filtering on the reconstructed denoised signal to obtain the final denoised signal;

[0020] As can be seen from the above technical solutions, the atomic clock signal denoising method of the present invention first applies empirical mode decomposition to atomic clock signal denoising, effectively separating the noise and real physical components aliased in the signal, laying the foundation for subsequent threshold denoising that is precisely matched with the characteristics of each frequency band. When performing wavelet denoising, an improved threshold function is constructed to introduce a smoothing factor. This smoothing factor controls the transition characteristics of the function near the threshold, thereby effectively suppressing noise while maximizing the preservation of signal details and smoothness. Applying this invention to atomic clock signal processing can effectively remove noise, providing a new technical approach to improving the data quality and timescale stability of atomic clocks.

[0021] Furthermore, in step S3, the smoothing factor of the IMF component is determined according to the following steps: for the i-th IMF component c i (t),

[0022] S301, Calculate IMF component c i The noise level estimate, kurtosis, energy concentration, and local variability of (t);

[0023] S302, Regarding IMF component c i The four features of noise level estimate, kurtosis, energy concentration and local fluctuation of (t) are normalized respectively to obtain the normalized values ​​of noise level feature, kurtosis feature, energy concentration feature and local fluctuation feature.

[0024] S303. Based on the four feature normalized values, a weighted fusion strategy is used to calculate the IMF component c. i The overall score S of (t) i ;

[0025] S304. Use a nonlinear mapping function to convert the IMF components c i The overall score S of (t) i Convert to IMF component c i Smoothing factor α of (t) i , α in the formula max α is the maximum smoothing factor. min It is the minimum smoothing factor.

[0026] In some embodiments, by constructing an adaptive smoothing factor selection strategy based on multi-feature fusion, multi-dimensional features such as noise level, kurtosis, energy concentration and local fluctuations are fused to dynamically adapt to the non-stationary and complex noise characteristics of atomic clock signals, thus solving the problem that traditional fixed parameter methods cannot balance noise suppression and detail preservation.

[0027] Furthermore, in step S401, the adaptive selection of the optimal wavelet basis is as follows:

[0028] S401-1. Initialize parameters and set the initial minimum entropy value;

[0029] S401-2. Determine the number of decomposition layers;

[0030] S401-3. Perform wavelet decomposition on each wavelet basis in the candidate wavelet basis set;

[0031] S401-4, Calculate the entropy index H. in, Indicates |c i The number of (t)|>ε, where ε is the reference value and a is the penalty coefficient. It is the number of decomposition levels;

[0032] S401-5. Traverse all candidate wavelet bases and select the wavelet base that minimizes the entropy index H as the optimal wavelet base.

[0033] Furthermore, in step S403-1, the step of determining the optimal threshold of the wavelet coefficient vector under the optimal wavelet basis is as follows:

[0034] a. Calculate the noise standard deviation of the wavelet coefficient vector under the optimal wavelet basis;

[0035] b. Generate a candidate threshold sequence Z = {T1, T2, ..., T...} M}(k=1,2,…,M), Z is the candidate threshold range, T1, T2,…,T k Let M be the first candidate threshold, the second candidate threshold, ..., the kth candidate threshold, and M be the number of candidate thresholds in the candidate threshold sequence.

[0036] c. Calculate the risk value SURE(T) for each candidate threshold. k ):

[0037] In the formula This represents the noise variance of the wavelet coefficient vector under the optimal wavelet basis. The wavelet coefficient vector under the optimal wavelet basis The number of wavelet coefficients, express The divergence;

[0038] d. Take the κ-th candidate threshold corresponding to the minimum risk value as the theoretically optimal threshold T. best ;

[0039] e. Shift the position of the theoretically optimal threshold to the right. The position after shifting to the right is κ' = min(κ + Δ, M), where Δ is the offset amount.

[0040] f. Select the κ'-th candidate threshold from the candidate threshold sequence after shifting to the right, and perform inverse normalization on the κ'-th candidate threshold to obtain the optimal threshold T'. best .

[0041] Furthermore, the noise standard deviation of the wavelet coefficient vector under the optimal wavelet basis. Med{} represents the median operator.

[0042] In some embodiments, an improved threshold selection strategy is employed, which enhances the robustness of denoising by introducing a right-biased conservative adjustment to adaptively increase the threshold. Attached Figure Description

[0043] To more clearly illustrate the embodiments of the present invention, the accompanying drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0044] Figure 1 This is a flowchart of the method of the present invention;

[0045] Figure 2 A flowchart for selecting a smoothing factor in an embodiment of the present invention;

[0046] Figure 3 This is a flowchart illustrating the thresholding process for wavelet coefficients according to an embodiment of the present invention;

[0047] Figure 4 This is a flowchart illustrating the process of determining the optimal threshold for the wavelet coefficient vector in an embodiment of the present invention.

[0048] Figure 5 Allen bias in simulation data from 5 atomic clocks;

[0049] Figure 6 The results of Cs1 EMD decomposition of cesium clock;

[0050] Figure 7 A comparison chart of data before and after noise reduction for five atomic clocks.

[0051] The specific embodiments of the present invention will be described in further detail below with reference to the accompanying drawings. Detailed Implementation

[0052] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0053] This invention's atomic clock denoising method first utilizes Empirical Mode Decomposition (EMD) to adaptively decompose the non-stationary atomic clock signal into a series of IMF (Intrinsic Mode Function) components and residual components, achieving multi-scale separation of the signal in the time-frequency domain. Then, for wavelet denoising of each IMF component, an improved threshold function based on the hyperbolic tangent function is constructed. This improved threshold function introduces a smoothing factor to achieve a continuous and smooth transition between traditional hard and soft thresholding characteristics. Finally, all denoised IMF components and the residual components obtained from the adaptive decomposition are reconstructed to obtain the final denoised signal, effectively suppressing noise while maximizing the preservation of signal details and smoothness. To dynamically adapt to signal characteristics, this invention employs a multi-feature fusion-based strategy to adaptively select the smoothing factor and utilizes an improved SURE criterion, combined with a right-shifting strategy, to determine the optimal threshold for each IMF component.

[0054] The following is combined Figure 1 The method for denoising atomic clock signals according to the present invention will be described. For example... Figure 1 As shown, the atomic clock signal denoising method in this embodiment includes the following steps:

[0055] S1. Load the original atomic clock signal x(t), read the atomic clock error observation sequence from the specified data file path, and identify the data dimension; construct a storage architecture for parameters and results to store performance indicators such as SNR and RMSE of the denoised and reconstructed signal sequence, as well as the value sequence of the adaptive smoothing factor.

[0056] S2. Perform EMD decomposition on the original signal x(t) to decompose it into a series of IMF components c. i (t) and a residual component r n (t), i = 1, ..., n, where n is the number of IMF components, i.e.

[0057] Each IMF component obtained from EMD decomposition must satisfy the following two conditions:

[0058] (1) In the entire data sequence, the number of local extrema and zero crossings are equal, or differ by at most 1;

[0059] (2) At any time t, the upper envelope e defined by the local maximum point max (t) and the lower envelope e defined by the local minimum point. min The average value of (t) is zero, i.e., [e max (t)+e min [(t)] / 2=0;

[0060] EMD decomposition is a conventional method, and this invention does not improve this part. For specific decomposition steps, please refer to the description of empirical mode decomposition of clock difference data in patent application No. 2015103920575, which will not be repeated here.

[0061] S3 represents each IMF component obtained after EMD decomposition, c. i (t) Select the smoothing factor α i ;

[0062] S4. After selecting the smoothing factor, for each IMF component c i (t) Perform improved wavelet denoising separately to obtain the denoised IMF components.

[0063] S5. All denoised IMF components The residual component r obtained from EMD decomposition n The signals (t) are superimposed to reconstruct a denoised signal.

[0064] S6. The reconstructed denoised signal... A moving average filter is applied to further improve the smoothness and continuity of the denoised signal; the moving average filter employs an adaptive window size strategy. The signal is the final denoised signal after moving average filtering, and K is the size of the moving window. λ is the denoised signal Length, This is the floor function.

[0065] To dynamically adapt to the non-stationary and complex noise characteristics of atomic clock signals, this invention determines a smoothing factor for each IMF component. This embodiment constructs an adaptive smoothing factor selection strategy based on multi-feature fusion to select the smoothing factor (step S3). The process for determining the smoothing factor is as follows: Figure 2 As shown, for the i-th IMF component c i (t), select the smoothing factor according to the following steps:

[0066] S301, Calculate IMF component c i Noise level estimate of (t) Kudo Energy Concentration and local volatility

[0067] S302, Regarding IMF component c iThe four characteristics of noise level estimate, kurtosis, energy concentration, and local variability of (t) are normalized to obtain the normalized values ​​of the noise level characteristics. Kurtosis characteristic normalized value Normalized value of energy concentration characteristics and the normalized value of local volatility characteristics

[0068] S303. Based on the four feature normalized values, a weighted fusion strategy is used to calculate the IMF component c. i The overall score S of (t) i , The weighting reflects the concept of prioritizing smoothness. This invention sets the weight of noise level and energy dispersion to 60% to ensure that smoothing effect is given priority in complex noise environments.

[0069] S304. Use a nonlinear mapping function to transform the IMF component c i The overall score S of (t) i Convert to this IMF component c i Smoothing factor α of (t) i , This invention employs a nonlinear mapping function to convert the composite score of IMF components into a smoothing factor, adapting to the sensitivity differences across different score intervals, α max α min These are the maximum and minimum smoothing factors, respectively. The values ​​of the maximum and minimum smoothing factors are empirical values. In this embodiment, α... max 15, α min The value is 3.

[0070] Noise level estimation is the direct basis for determining denoising intensity. Kurtosis is a dimensionless parameter that reflects the numerical statistical characteristics of the random variable distribution of a signal, and also represents the normalized fourth central moment of the signal. Local variability is a time-domain analysis index used to quantify the degree of drastic change in a signal between adjacent sampling points, and can directly reflect the short-term instability and high-frequency noise level of the signal.

[0071] This embodiment uses an estimation method based on the median absolute deviation to estimate the noise level, and the IMF component c i Noise level estimate of (t) Med{} represents the median operator.

[0072] The IMF component c in this embodiment i kurtosis of (t) In the formula, E(·) represents the expectation. These are IMF components c i The mean and variance of (t).

[0073] This embodiment calculates the IMF component c. i Energy concentration of (t) At that time, first analyze the IMF component c i (t) Perform discrete wavelet decomposition to separate the IMF components c i (t) is decomposed into approximate coefficients and detail coefficients at different scales, yielding the wavelet coefficient vector. and length vector IMF component c i Energy concentration of (t) In the formula, Γ(·) represents The number of Represents the wavelet coefficient vector The number of wavelet coefficients, For IMF component c i Noise level estimate of (t). Wavelet coefficient vector It contains wavelet coefficients at all levels, arranged from low to high frequency, with a length vector. It includes the length of coefficients at each level, as well as the currently processed IMF component c. i The signal length (t) can be used. In practical applications, the built-in wavedec function in MATLAB can be used to perform discrete wavelet decomposition on the IMF components. In the formula This indicates the name of the wavelet basis function; in this embodiment, the db4 wavelet basis is selected. The number of decomposition layers is adaptively determined based on the signal strength. N i It is an IMF component c i The signal length is (t), floor(·) represents the floor function, and min(·) represents the minimum of the two.

[0074] The IMF component c in this embodiment i Local fluctuations of (t) The IMF component c is calculated from the ratio of the standard deviation of the first-order difference sequence of the signal to that of the original sequence. i Local fluctuations of (t) In the formula The standard deviation of the first-order difference sequence; For IMF component c i The standard deviation of (t).

[0075] Obtain IMF component c i After estimating the noise level, kurtosis, energy dispersion, and local variability of (t), the normalized values ​​of the noise level characteristics, kurtosis characteristics, energy dispersion characteristics, and local variability characteristics are calculated using the following formulas:

[0076] Noise level characteristic normalization value In the formula, max(·) represents the maximum value of the two, and min(·) represents the minimum value of the two. R = max(c i (t))-min(c i (t)), max(c i (t) represents the IMF component c obtained from EMD decomposition. i The maximum value in (t), min(c) i (t) represents the IMF component c obtained from EMD decomposition. i The minimum value in (t);

[0077] Kurtosis characteristic normalized value

[0078] Normalized value of energy concentration characteristics

[0079] Local volatility characteristic normalized value

[0080] In step S4, the IMF component c i The steps for improving wavelet denoising are as follows:

[0081] S401. Adaptive strategy is adopted for all IMF components to adaptively select the optimal wavelet basis for each component.

[0082] S402. After selecting the optimal wavelet basis through an adaptive strategy, apply the optimal wavelet basis to each IMF component c. i (t) Perform discrete wavelet decomposition to separate the IMF components c i (t) is decomposed into approximation coefficients and detail coefficients at different scales, yielding the wavelet coefficient vector under the optimal wavelet basis corresponding to each IMF component. wavedec represents wavelet decomposition. Represents the optimal wavelet basis. This indicates that wavelet decomposition is performed using the optimal wavelet basis. For optimal wavelet basis The wavelet coefficient vector below, For optimal wavelet basis The length vector below, i.e. The i-th IMF component is based on the optimal wavelet basis. The wavelet coefficient vector and length vector obtained by wavelet decomposition;

[0083] S403. Thresholding is applied to the wavelet coefficient vector obtained under the optimal wavelet basis to suppress noise, ultimately yielding the denoised IMF component signal. The specific steps for thresholding the wavelet coefficients in this embodiment are as follows:

[0084] S403-1. Determine the wavelet coefficient vector under the optimal wavelet basis obtained from step S402. The optimal threshold T' best ;

[0085] S403-2, Obtaining the optimal threshold T' best Then, an improved threshold function based on hyperbolic tangent is constructed, and the improved threshold function based on hyperbolic tangent is applied to each layer of wavelet coefficients in the wavelet coefficient vector for processing. During the thresholding process, a smoothing factor is introduced for each IMF component through the improved threshold function. The improved threshold function is as follows: T' represents the improved threshold function. best For the optimal threshold, sgn(·) represents the sign function. For optimal wavelet basis The wavelet coefficient vector below, Represents the wavelet coefficient vector The modulus, α i Let δ be the smoothing factor for the i-th IMF component, and δ be the auxiliary smoothing parameter. In this embodiment, δ = T' best ;

[0086] S403-3. For the IMF component after thresholding, use the wavelet coefficient vector under the optimal wavelet basis obtained by wavelet decomposition in step S402. and length vector and the optimal wavelet basis selected in step S401 Perform wavelet reconstruction to restore the time-domain signal. That is, the denoised IMF component signal; in practical applications, wavelet reconstruction can be performed using the discrete wavelet inverse transform function waverec() in MATLAB. This is the denoised IMF component signal.

[0087] Traditional hard and soft thresholding functions suffer from the Gibbs phenomenon caused by discontinuous threshold points and excessive signal smoothing due to constant bias, respectively. To overcome these shortcomings, this embodiment employs an improved thresholding function based on the hyperbolic tangent function during wavelet denoising of the IMF components. By introducing a smoothing factor, a continuous and smooth transition within the threshold neighborhood is achieved. Combined with an improved threshold selection strategy, this effectively suppresses noise while preserving signal details.

[0088] Traditional wavelet thresholding denoising methods all use a single wavelet basis for wavelet denoising, which cannot adapt to the differences in time-frequency characteristics of different IMF components. To overcome these problems, this embodiment improves upon traditional wavelet thresholding denoising when performing wavelet denoising on IMF components. During wavelet denoising, an adaptive wavelet basis selection strategy based on the minimum entropy criterion with a penalty term is constructed to perform optimal adaptive wavelet basis selection. Then, the wavelet basis that can produce the sparsest representation of the current IMF component is selected from the candidate wavelet basis set, achieving adaptive wavelet basis selection and improving the effect of wavelet thresholding denoising. Figure 4 As shown, in step S401, the process of adaptively selecting the optimal wavelet basis includes the following steps:

[0089] S401-1. Initialize parameters; set the initial minimum entropy value H. min =∞, the number of coefficients in the coefficient vector whose absolute value is greater than the reference value ε, this number is the entropy value.

[0090] S401-2. Determine the number of decomposition levels; the number of wavelet decomposition levels determines the level of detail in signal analysis, directly affecting denoising performance and computational efficiency. The maximum number of layers is set to 4 to prevent the Gibbs phenomenon from occurring due to too many layers;

[0091] S401-3. Perform Discrete Wavelet Decomposition: To evaluate the suitability of different wavelet bases for the current signal, decompose each wavelet base in the candidate wavelet base set. The candidate wavelet base set is wave_list = {'db4','db5','db6','db8','sym6','sym8','coif3','coif4','coif5'}. For each wavelet base in the set, perform discrete wavelet decomposition. In specific applications, the wavedec function built into MATLAB can also be used to perform discrete wavelet decomposition on the IMF components. It is the number of decomposition levels;

[0092] S401-4. Calculate the entropy index; To quantify the sparsity of coefficients generated by each wavelet basis and comprehensively consider model complexity, this embodiment defines an entropy index with a penalty term as a selection criterion. in, Indicates |c i The number of (t)|>ε, where ε is the reference value. The decomposition layer is 'a', which is a penalty coefficient used to adjust the weight of the influence of the layer. In this invention, a = 0.1.

[0093] S401-5. Selecting the optimal wavelet basis: After traversing all candidate wavelet bases, select the wavelet base that minimizes the entropy index H as the optimal wavelet basis.

[0094] In some embodiments, the wavelet coefficient vector under the optimal wavelet basis is adaptively selected by combining the Stein unbiased risk estimation criterion with a right-shifting strategy. The optimal threshold T' best ,like Figure 4 As shown, in step S403-1, the optimal threshold T' of the wavelet coefficient vector under the optimal wavelet basis is determined. best The steps are as follows:

[0095] a. Calculate the wavelet coefficient vector under the optimal wavelet basis using the median absolute deviation estimation method. noise standard deviation The wavelet coefficients are normalized based on the noise standard deviation to unify the calculation scale. The i-th IMF component c i Normalized wavelet coefficients of (t) The method of estimating the absolute deviation of the median to calculate the noise standard deviation is robust to outliers and can accurately reflect the noise level.

[0096] b. Generate a candidate threshold sequence Z = {T1, T2, ..., T...} M}(k=1,2,…,M), Z is the candidate threshold range, T1, T2,…,T k These are the first candidate threshold, the second candidate threshold, ..., the kth candidate threshold, respectively, and M is the number of candidate thresholds. In this embodiment, the threshold range is [0.5, 4], and M = 100.

[0097] c. Calculate the risk value SURE(T) for each candidate threshold based on the SURE criterion. k ), In the formula This represents the noise variance of the wavelet coefficient vector under the optimal wavelet basis. The wavelet coefficient vector under the optimal wavelet basis The number of wavelet coefficients, Represents the improved threshold function The divergence;

[0098] d. Determine the theoretically optimal threshold T by minimizing the risk value. best Theoretical optimal threshold T best T is the candidate threshold corresponding to the minimum risk value. best =minSURE(T kThe theoretically optimal threshold is located at position κ in the candidate threshold sequence, which is the κth candidate threshold. The κth candidate threshold with the smallest risk value is taken as the theoretically optimal threshold.

[0099] e. Considering the denoising requirements in complex noise environments, a right-shifting strategy is adopted to ensure sufficient noise suppression. The position of the theoretically optimal threshold is shifted to the right, and the position after the right shift is κ' = min(κ + Δ, M), where Δ is the shift amount. In this embodiment, Δ = 5.

[0100] f. Select the κ'-th candidate threshold from the candidate threshold sequence after shifting to the right, and perform inverse normalization on the κ'-th candidate threshold to obtain the optimal threshold.

[0101] To verify the performance of the method of the present invention, simulations were performed on the method of the present invention, the traditional wavelet thresholding denoising method based on hard threshold function, and the wavelet thresholding denoising method based on soft threshold function, and the denoising results of the three methods on atomic clock signals were compared.

[0102] During the simulation, two different models of cesium atomic clocks (cesium clock Cs1 and cesium clock Cs2) were first selected as frequency standards. Cesium atomic clocks have good long-term frequency stability, and the linear frequency drift due to aging is extremely low and negligible in practical applications, thus ensuring a relatively small frequency offset. The deviation between the output signal of the cesium atomic clock and the ideal reference signal can be systematically decomposed into three core parameters: phase deviation, frequency deviation, and frequency drift. These three parameters constitute a complete dynamic system capable of capturing the main noise characteristics and deterministic trends of the cesium clock.

[0103] Based on the aforementioned physical mechanism, the following discrete-time state-space model is established to characterize the properties of the cesium atomic clock. Three key state variables are defined in the model: p t f represents the clock difference at time t, i.e., the cumulative phase deviation; t d represents the frequency deviation at time t; t Let t represent the frequency drift at time t. Let the time interval be τ, then the state evolution equation of the model is:

[0104] In the formula, ε, μ, and σ are the phase deviation, frequency deviation, and random noise in the frequency drift process, respectively.

[0105] The output deviation of an atomic clock is a complex stochastic process resulting from the superposition of various unrelated noise types. To accurately simulate atomic clock data, it is essential to first determine the main noise types and their parameters. These parameters can be obtained by analyzing the Allan variance, a typical indicator of atomic clock frequency stability. The variation of the Allan variance over different averaging times τ exhibits a clear functional relationship with specific noise types; therefore, by fitting the measured Allan variance curve, the intensity coefficients of various noise types can be separated and solved.

[0106] The overall frequency stability of an atomic clock is composed of four main noise components: phase white noise (WPM), frequency white noise (WFM), frequency flicker noise (FFM), and frequency random walk noise (RWFM). These noises are uncorrelated, and the Allan variance can be expressed as the sum of the variances of each noise term:

[0107] The functional relationship between each noise term and the Allan standard deviation is as follows: A in the formula WP For phase white noise at τ = 1s, f h The Allan standard deviation at 1 Hz, A WF A FF A RW These are the Allan standard deviations of frequency white noise, frequency flicker noise, and frequency random walk noise at τ = 1s, respectively.

[0108] To solve for these four noise parameters A WP A WF A FF A RW It needs to be done over multiple average times τ i The Allan standard deviation σ is measured on (i = 1, 2, ..., N). y,m (τ i And establish the following system of equations:

[0109]

[0110] This system of equations can be written in matrix form: M·A=Σ(2), where,

[0111] By solving this system of equations using the least squares method, the optimal noise parameter vector A can be obtained:

[0112] A = (M T M) -1 M T Σ.

[0113] A WP A WF AFF A RW The noise parameters required to simulate atomic clock noise are used for subsequent clock difference data synthesis.

[0114] In addition, two hydrogen atomic clocks (H1 and H2) and one rubidium atomic clock (Rb1) were selected for comparison. Their frequency stability and drift rate were used to simulate the clock error data of the atomic clocks. The performance indicators of the five atomic clocks are shown in Table 1.

[0115] Table 1 Performance Indicators of Atomic Clocks

[0116]

[0117] The performance indicators of the five atomic clocks are used as inputs and substituted into equation (1) or equation (2). The least squares method is used to solve for the four noise figures: WPM, WFM, FFM and RWFM. The specific noise figures obtained are shown in Table 2.

[0118] Table 2 Noise parameters of atomic clocks

[0119]

[0120] The clock difference data of five atomic clocks were simulated using Stable32 software as the raw signal of the atomic clocks. 5000 data points were collected from each atomic clock at a time interval of 1 second.

[0121] To verify the validity of the simulation data, the Allan Deviation (ADEV) of the data from the five simulated atomic clocks was calculated. Figure 5 Allen deviation curves for five analog atomic clocks (cesium clocks Cs1 and Cs2, hydrogen clocks H1 and H2, and rubidium clock Rb1) are displayed. Figure 5 Overall, the Allen bias curves of the two hydrogen clocks (H1, H2) almost overlap and are located at the bottom of the graph, with values ​​approximately two orders of magnitude lower than the cesium clock. This aligns with the theory that the hydrogen maser exhibits excellent short-term stability due to atomic interactions. The stability of the cesium clock Cs1 is superior to that of the cesium clock Cs2 across all average times, consistent with the superior frequency stability index of Cs1 shown in Table 1. The Allen bias of all simulated clock bias data gradually decreases with increasing average time, consistent with the characteristics of typical atomic clock noise processes. The stability of the simulated atomic clock data is essentially consistent with that of the NTSC atomic clock; therefore, the clock bias data from this simulation are considered valid.

[0122] For atomic clock signals, EMD can naturally separate noise mixed in different frequency bands from the real signal components based on its own time-scale characteristics, laying an ideal foundation for subsequent implementation of precise denoising strategies that match the frequency band characteristics. To intuitively verify the effectiveness of EMD in atomic clock signal processing, simulated clock bias data of a cesium clock Cs1 is used as an example, and its decomposition results are as follows: Figure 6 As shown.

[0123] After EMD decomposition of the original atomic clock signal, an improved wavelet threshold denoising method is applied to each IMF component obtained from the decomposition. The denoising performance of the algorithm is evaluated using two metrics: SNR and RMSE. SNR measures the power ratio of the effective component to the residual noise in the denoised clock error signal; a higher SNR indicates a better separation and removal of noise energy from the observation data. RMSE reflects the average deviation between the denoised clock error and the original observed clock error; a lower RMSE indicates a closer similarity in waveform and amplitude between the denoised signal and the original signal. The expressions for these two metrics are as follows:

[0124]

[0125] Where, x(n) and These are the original noisy signal and the denoised reconstructed signal, respectively. and These represent the total power of the original signal and the total power of the noise, respectively. To test the denoising performance of each denoising method, a rubidium atomic clock, model PRS10, from the inventor's laboratory was selected for comparison and denoted as rubidium clock Rb2. 5000 data points were collected at a time interval of 1 second. Figure 7 Figures (a)-(f) show the comparison of denoising data for cesium clocks Cs1 and Cs2, hydrogen clocks H1 and H2, and rubidium clocks Rb1 and Rb2 before and after denoising. The gray curves represent the original clock difference data, and the black curves represent the clock difference data after denoising. Tables 3 and 4 respectively show the denoising performance of traditional wavelet thresholding denoising based on hard thresholding function, wavelet thresholding denoising based on soft thresholding function, and the method of this invention on the clock difference data of six different types of atomic clocks.

[0126] Table 3. SNR (dB) of various denoising methods under different atomic clocks.

[0127]

[0128] Table 4. RMSE (×10⁻¹¹ s) of each threshold method under different atomic clocks.

[0129]

[0130] As shown in Tables 3 and 4, the method of this invention achieved the highest signal-to-noise ratio (SNR) on all six atomic clocks in terms of noise suppression. Specifically, compared to the hard thresholding method, the method of this invention improved noise suppression by approximately 25% and 35% on cesium and hydrogen clocks, respectively; compared to the soft thresholding method, the method of this invention improved noise suppression by approximately 14% on cesium clocks, and even for the superior-performing hydrogen clock, the method of this invention achieved a further improvement of approximately 5%. For rubidium clocks, compared to the hard thresholding method, the soft thresholding method and the method of this invention improved noise suppression by approximately 6% and 18% on simulated data, respectively; however, in measured data, the soft thresholding method only improved noise suppression by 0.5%, indicating limited improvement when processing measured signals; in contrast, the method of this invention improved noise suppression by approximately 26%, demonstrating stronger adaptability and advantages to measured data. This result verifies that the improved threshold function of the method of this invention, by introducing a smooth transition mechanism, can more accurately distinguish between noise and signal energy, retain more large-amplitude wavelet coefficients that reflect the real physical process, thereby retaining a higher proportion of useful power in the output signal.

[0131] Regarding signal fidelity, the method of this invention achieved the lowest RMSE values ​​in all test scenarios, indicating that its denoising results are closest to the original observation data in the time domain. Specifically, for cesium clocks, the method of this invention reduces RMSE by approximately 41% and 28% compared to hard and soft thresholding methods, respectively. For hydrogen clocks, the RMSE is reduced to 8.7 × 10⁻⁶. -14 Compared to the soft thresholding method, the error is reduced by approximately 38%, achieving an almost order-of-magnitude improvement in accuracy. For rubidium clocks, the method of this invention reduces RMSE by approximately 25% in simulated data. In measured data, compared to the hard thresholding method, the soft thresholding method and the method of this invention achieve reductions of approximately 1.5% and 26%, respectively, further demonstrating that the soft thresholding method offers limited improvement in real-world testing environments, while the method of this invention offers significant advantages. The above simulation results show that the method of this invention effectively suppresses noise while minimizing signal distortion, better preserving the original form and detailed characteristics of atomic clock error data.

[0132] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention in any way. Although the present invention has been disclosed above with reference to preferred embodiments, it is not intended to limit the present invention. Any person skilled in the art can make some modifications or alterations to the above-disclosed technical content to create equivalent embodiments without departing from the scope of the present invention. Any simple modifications, equivalent changes, and alterations made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the scope of the present invention.

Claims

1. A method for denoising atomic clock signals, characterized in that, Includes the following steps: S1, Load the original atomic clock signal; S2. Perform EMD decomposition on the original signal, decomposing it into a series of IMF components and a residual component; S3. Determine the smoothing factor for each IMF component; S4. Perform improved wavelet denoising on each IMF component to obtain the denoised IMF components. The steps are as follows: S401. An adaptive strategy is used to select the optimal wavelet basis for each IMF component. S402. Based on the optimal wavelet basis, perform wavelet decomposition on each IMF component to obtain the wavelet coefficient vector and length vector under the optimal wavelet basis for each IMF component. S403. Threshold the wavelet coefficient vectors of each IMF component under the optimal wavelet basis to obtain the denoised IMF components. In step S403, the thresholding process for the wavelet coefficient vector under the optimal wavelet basis is as follows: S403-1. Determine the optimal threshold for the wavelet coefficient vector under the optimal wavelet basis; S403-2. Construct an improved threshold function. Apply the improved threshold function to the wavelet coefficient vector under the optimal wavelet basis. The improved threshold function is as follows: In the formula This represents the improved threshold function, and sgn(·) represents the sign function. T′ is the wavelet coefficient vector under the optimal wavelet basis. best The optimal threshold is δ, which is the auxiliary smoothing parameter, and α is the value of α. i Let be the smoothing factor for the i-th IMF component; S403-3. For the IMF component after thresholding, the wavelet coefficient vector and length vector of the component obtained by wavelet decomposition in step S402, and the optimal wavelet basis selected in step S401 are used to perform wavelet reconstruction to restore the time domain signal, which is the denoised IMF component. S5. Superimpose all the denoised IMF components with the residual components to reconstruct the denoised signal; S6. Perform moving average filtering on the reconstructed denoised signal to obtain the final denoised signal.

2. The atomic clock signal denoising method as described in claim 1, characterized in that: In step S3, the smoothing factor of the IMF component is determined according to the following steps: for the i-th IMF component c i (t), S301, Calculate IMF component c i The noise level estimate, kurtosis, energy concentration, and local variability of (t); S302, Regarding IMF component c i The four features of noise level estimate, kurtosis, energy concentration and local fluctuation of (t) are normalized respectively to obtain the normalized values ​​of noise level feature, kurtosis feature, energy concentration feature and local fluctuation feature. S303. Based on the four feature normalized values, a weighted fusion strategy is used to calculate the IMF component c. i The overall score S of (t) i ; S304. Use a nonlinear mapping function to convert the IMF components c i The overall score S of (t) i Convert to IMF component c i Smoothing factor α of (t) i , α in the formula max α is the maximum smoothing factor. min It is the minimum smoothing factor.

3. The atomic clock signal denoising method as described in claim 1, characterized in that: In step S401, the steps for adaptively selecting the optimal wavelet basis are as follows: S401-1. Initialize parameters and set the initial minimum entropy value; S401-2. Determine the number of decomposition layers; S401-3. Perform wavelet decomposition on each wavelet basis in the candidate wavelet basis set; S401-4, Calculate the entropy index H. in, Indicates |c i The number of (t)|>ε, where ε is the reference value and a is the penalty coefficient. It is the number of decomposition levels; S401-5. Traverse all candidate wavelet bases and select the wavelet base that minimizes the entropy index H as the optimal wavelet base.

4. The atomic clock signal denoising method as described in claim 1, characterized in that: In step S403-1, the steps for determining the optimal threshold of the wavelet coefficient vector under the optimal wavelet basis are as follows: a. Calculate the noise standard deviation of the wavelet coefficient vector under the optimal wavelet basis; b. Generate a candidate threshold sequence Z = {T1, T2, ..., T...} M }(k=1,2,…,M), Z is the candidate threshold range, T1, T2,…,T k Let M be the first candidate threshold, the second candidate threshold, ..., the kth candidate threshold, and M be the number of candidate thresholds in the candidate threshold sequence. c. Calculate the risk value SURE(T) for each candidate threshold. k ): In the formula This represents the noise variance of the wavelet coefficient vector under the optimal wavelet basis. The wavelet coefficient vector under the optimal wavelet basis The number of wavelet coefficients, express The divergence; d. Take the κ-th candidate threshold corresponding to the minimum risk value as the theoretically optimal threshold T. best ; e. Shift the position of the theoretically optimal threshold to the right. The position after shifting to the right is κ' = min(κ + Δ, M), where Δ is the offset amount. f. Select the κ'-th candidate threshold from the candidate threshold sequence after shifting to the right, and perform inverse normalization on the κ'-th candidate threshold to obtain the optimal threshold T′. best .

5. The atomic clock signal denoising method as described in claim 4, characterized in that: Noise standard deviation of wavelet coefficient vector under optimal wavelet basis Med{} represents the median operator.

Images