Pulsar signal denoising method based on EWT

Abstract

The application discloses a pulsar signal denoising method based on EWT, and specifically comprises the following steps: S1, acquiring a pulsar cumulative pulse profile x(t), wherein t represents a time variable; S2, performing EWT decomposition on the pulsar cumulative pulse profile x(t); S3, performing FFT on M EMFs {emf m (t), m=1, 2, …, M}; S4, calculating a frequency domain feature set {EMA m , m=1, 2, …, M}; S5, determining a starting layer number K of high-frequency EMF components; S6, directly removing the Mth layer EMF component emf M (t) as noise; S7, performing denoising processing on emf K (t)~emf M‑1 (t) by using a wavelet threshold method, and denoising results are recorded as S8, performing denoising processing on the low-frequency EMF components emf1(t)~emf K‑1 (t) by using a Savitzky-Golay smoothing filtering method, and denoising results are recorded as S9, and a denoised pulsar signal is obtained. The method can effectively improve the signal-to-noise ratio of the pulsar cumulative profile and improve the quality of the pulsar cumulative profile.

Description

TECHNICAL FIELD

[0001] The application belongs to the field of astronomical navigation signal processing, and particularly relates to a pulsar signal denoising method based on EWT. BACKGROUND

[0002] A pulsar is a kind of neutron star that rotates rapidly and can periodically emit pulse signals in the process of rotation. The pulse signals cover the electromagnetic wave frequency bands such as radio, infrared, visible light, X-ray and gamma ray. The rotation of the pulsar is not affected by external factors, and the rotation of the pulsar is very regular. The electromagnetic wave frequency swept out by the pulsar is very stable, and the stability can be comparable to that of a cesium atomic clock, and is called the most accurate astronomical clock in nature. Therefore, the pulse signals radiated by the pulsar can be compared to a lighthouse, which is used to realize autonomous navigation of a spacecraft and provide position, speed, attitude and other information for the spacecraft. In addition, the FAST telescope can observe very weak pulsar signals from very far away, which is used for astronomical and physical research, so as to reveal more mysteries of the universe. Therefore, the pulse signals radiated by the pulsar contain a large amount of useful information, which can provide data support for many researches. However, the pulsar is very far away from the earth, and in the process of propagation, it will be affected by many interferences. The very low energy flow density limits the signal-to-noise ratio of the pulsar signal, and is seriously completely submerged by noise. Therefore, the study of pulsar signal denoising can provide higher quality pulsar data for subsequent application research.

[0003] The denoising method based on FFT is the most commonly used traditional pulsar signal denoising method. The pulsar is a non-stationary signal, and the use of the FFT method is easy to be affected by harmonic interference. In order to overcome this problem, some scholars choose the wavelet-based denoising method. However, in the process of using the wavelet method, it is necessary to select a suitable wavelet basis and decomposition layer according to the characteristics of the signal. In view of the problems existing in the wavelet method, some scholars use EMD for denoising. The EMD method can adaptively decompose the signal, avoiding the shortcomings of wavelet decomposition. However, in the process of signal decomposition, envelope, under envelope, and different degrees of end effect and modal aliasing problems will occur, which will affect the quality of signal decomposition.

[0004] In recent years, some new time-frequency analysis methods have been gradually proposed, which provide a new adaptive time-frequency analysis idea for non-stationary signal processing. For example, the empirical wavelet transform (EWT) method, which inherits the advantages of adaptive decomposition of EMD and has a tight support framework of wavelet transform theory. In the process of decomposition, the EWT method can adaptively select the frequency band, overcoming the modal aliasing problem in the EMD method, and avoiding the problems of over envelope and under envelope; in addition, due to the completely reliable mathematical theoretical basis, the EWT has the advantage of low computational complexity. The EWT is widely used in many signal processing fields.

[0005] The observation signal can be decomposed into multiple empirical mode functions EMF by using EWT, and the specific process is as follows: firstly, the modes in the signal and the empirical wavelet base are determined by analyzing the Fourier spectrum of the observation signal; on this basis, the wavelet filter bank is constructed; finally, the observation signal is adaptively decomposed into multiple EMF by using the wavelet filter bank. The EMF is a time domain signal component at different frequencies. Among them, the noise is mainly concentrated in the high frequency region, so it is necessary to screen out the high frequency EMF component. On this basis, according to the different characteristics of noise and signal in the high frequency EMF component and the low frequency EMF component, a suitable method is selected for denoising processing, so as to realize the noise suppression of the signal. As described above, there are two problems when using EWT for signal denoising: ① distinguishing high frequency EMF components and low frequency EMF components; ② selecting a suitable denoising method for high frequency EMF components and low frequency EMF components. SUMMARY

[0006] The purpose of the present application is to provide an EWT-based pulsar signal denoising method, which can effectively improve the signal-to-noise ratio of the pulsar cumulative profile and improve the quality of the pulsar cumulative profile.

[0007] The technical scheme adopted by the present application is an EWT-based pulsar signal denoising method, which is implemented according to the following steps:

[0008] Step 1, obtaining a pulsar cumulative pulse profile x(t), wherein t represents a time variable;

[0009] Step 2, EWT decomposition of the pulsar cumulative pulse profile x(t) to obtain {emf m (t), m = 1, 2, …, M};

[0010] Step 3, FFT of the {emf m (t), m = 1, 2, …, M} to obtain M frequency spectrum components {EMF m (f), m = 1, 2, …, M}, wherein f represents a frequency variable;

[0011] Step 4, calculating a frequency domain feature set {EMA m , m = 1, 2, …, M};

[0012] Step 5, using the {EMA m , m = 1, 2, …, M} to determine the starting layer number K of the high frequency EMF component;

[0013] Step 6, directly removing the Mth layer EMF component emf M (t) as noise;

[0014] Step 7, using a wavelet threshold method to process emf K (t) ~ emfM-1 (t) performing noise reduction processing on the low-frequency EMF components, and the noise reduction result is denoted as

[0015] Step 8, using the Savitzky-Golay smoothing filter method to process emf1(t)~emf K-1 (t) performing noise reduction processing on the low-frequency EMF components, and the noise reduction result is denoted as

[0016] Step 9, using and the approximate component emf0(t) to reconstruct, to obtain the denoised pulsar signal

[0017] The present application is also characterized in that,

[0018] The specific implementation of step 1 is: determining the type of pulsar, obtaining the pulsar data from the EPN database; setting the number of sampling points and the signal-to-noise ratio of the observed signal, and resampling and normalizing the obtained pulsar data in the database to obtain the pulsar standard profile s(t); according to the signal-to-noise ratio of the observed signal, noise n(t) is added to s(t) to obtain the pulsar cumulative pulse profile x(t), that is, x(t) = s(t) + n(t), where t represents the time variable, and n(t) is a Gaussian white noise with a mean of zero and a variance of σ 2 .

[0019] The specific implementation of step 2 is:

[0020] Step 2.1, performing Fourier transform on x(t) to obtain the signal spectrum X(ω), and the specific formula is as follows:

[0021]

[0022] In formula (1), j represents the imaginary unit; t represents time; ω represents frequency, and ω ∈ [0, π];

[0023] Step 2.2, normalizing the signal spectrum X(ω) to [0, π], and dividing [0, π] into M consecutive frequency bands, and the mth frequency band is Λ m = [ω m-1 ,ω m ], where ω m is the boundary value of the mth frequency band; M+1 endpoints are needed to divide M consecutive boundaries, and ω0=0, ω M = π;

[0024] Step 2.3, using the maximum value-based Fourier analysis method to determine M-1 endpoints, and the specific method is:

[0025] Firstly, find N local maximum values in the signal spectrum X(ω), i.e. N local peaks, where N≥M; the frequency corresponding to each local peak is called a local peak frequency, denoted as {Ω j , j = 1, 2, …, N}.

[0026] Secondly, arrange the N local peaks in descending order of frequency, retain the first M local peaks, and eliminate the redundant local peaks; then re-sort the local peak sequence in ascending order of frequency, and the frequency is denoted as {Ω j , j = 1, 2, …, M}.

[0027] Thirdly, take the median value of the frequencies corresponding to the adjacent two peak points as the boundary point, and formula (2) is the calculation formula of the mth frequency band for the boundary value ω m ; finally, divide [0, π] into M continuous intervals, and there are

[0028]

[0029] Finally, take ω m as the center frequency, define T m as the transition band, and T m = 2τ m , where τ m =γω m ; wherein, the parameter γ ∈ (0, 1), and satisfies the following conditions:

[0030]

[0031] Step 2.4, after determining Λ m , construct an empirical scale function and M empirical wavelet functions φ m (t); wherein, is a low-pass filter, and φ m (t) is a band-pass filter; their expressions are as follows:

[0032]

[0033]

[0034] The function β(x) in formula (4) and formula (5) satisfies the following conditions:

[0035]

[0036] There are multiple polynomials that satisfy the second condition in formula (6), and here we will use the following formula:

[0037] β(x) = x 4 (35-84x + 70x2 -20x 3 ) (7)

[0038] Step 2.5, inner product of x(t) and the empirical scaling function to get the approximation coefficient term W f (t); inner product of x(t) and the empirical wavelet function φ m (t) to get the mth detail coefficient term, the specific formula is as follows:

[0039]

[0040] W f (m,t) = <x(t), φ m (t) > (9)

[0041] In formula (8) and formula (9), <a, b> represents the inner product of a and b; m is 1 to M; t represents a time variable;

[0042] Step 2.6, calculate M+1 EMF components, including 1 approximation component and M detail components, the specific form is as follows:

[0043]

[0044] emf m (t) = W f (m,t) * φ m (t) (11)

[0045] In formula (10), emf0(t) is the approximation component; in formula (11), emf m (t) is the mth detail component, m = 1, 2, …, M;

[0046] Use the M+1 EMF components to reconstruct the signal, as follows:

[0047]

[0048] The pulsar cumulative pulse profile x(t) = s(t) + n(t), where s(t) is a standard pulse profile, n(t) is a Gaussian white noise with mean zero and variance σ 2 ; The noise here is additive Gaussian white noise, emf m (t) in formula (11) is expressed as:

[0049] emf m (t) = <s(t), φ m (t) > * φ m (t) + <n(t), φ m (t) > * φ m(t) = emf i s (t) + emf n i (t) (13)

[0050] In formula (13), emf i s (t) represents the mth EMF component of s(t); emf i n (t) represents the mth EMF component of n(t).

[0051] EMF m (f) is expressed as follows:

[0052] EMF m (f) = FFT[emf m (t)], m = 1, 2, …, M (14)

[0053] In formula (14), FFT[a] represents performing FFT operation on a; EMF m (f) is the FFT result corresponding to emf m (t); f represents a frequency variable.

[0054] In step 4, EMA m is specifically in the following form:

[0055] EMA m = |mean[EMF m (f)]|, m = 1, 2, …, M (15)

[0056] In formula (15), mean[a] represents taking the mean value of a, and |b| represents taking the modulus of b; EMA m is the eigenvalue corresponding to EMF m (f).

[0057] In step 5, the starting layer number K of the high-frequency EMF component is determined by using {EMA m , m = 1, 2, …, M}, and the expression of K is as follows:

[0058]

[0059] Formula (16) represents finding the maximum eigenvalue in the frequency domain feature set {EMA m , m = 1, 2, …, M}, and the layer number corresponding to the maximum eigenvalue is the value of K.

[0060] In step 7, the wavelet threshold method is used to process emf K (t) ~ emf M-1(t) The method of noise reduction processing is as follows: a Sym8 wavelet base is selected to perform 3-layer wavelet decomposition on the EMF component; after decomposition, a soft threshold function is used to remove noise from the detail component, and the threshold is determined by a Heursure method; finally, the high-frequency EMF component after noise removal is reconstructed, denoted as

[0061] In step 9, the pulsar signal after noise reduction is reconstructed The specific formula is as follows:

[0062]

[0063] The beneficial effects of the present application are:

[0064] (1) The method of the present application first decomposes the pulsar signal by using the EWT method to obtain a plurality of EMF components, clearly describing the relationship between the frequency of the pulsar signal and time; secondly, the high-frequency EMF component and the low-frequency EMF component are screened out by using the frequency domain characteristics of the EMF component; finally, the high-frequency EMF component and the low-frequency EMF component are respectively subjected to noise reduction processing by using the wavelet threshold denoising method and the Savitzky-Golay smoothing filter method, and a high-quality pulsar cumulative pulse profile is reconstructed.

[0065] (2) Compared with the commonly used denoising method, the method of the present application can accurately screen out the high-frequency EMF component and effectively suppress the noise component in each EMF component. The reason for the above advantages is that: first, because the FFT method is used to analyze the EMF component, it is beneficial to screen out the high-frequency EMF component; secondly, because two different denoising methods are used to denoise the high-frequency EMF component and the low-frequency EMF component respectively, the noise component in each EMF component is effectively suppressed. A high-quality pulsar cumulative pulse profile is reconstructed, overcoming the problem of low signal-to-noise ratio caused by low energy flow density. The method can effectively improve the signal-to-noise ratio of the pulsar cumulative profile, and provide effective data support for various applications of the pulsar signal. BRIEF DESCRIPTION OF DRAWINGS

[0066] Figure 1 is a flowchart of the method of the present application;

[0067] Figure 2 is a standard pulse profile of PSR B1508+55;

[0068] Figure 3 is a cumulative pulse profile of PSR B1508+55;

[0069] Figure 4 is an EWT decomposition diagram of the cumulative pulse profile of PSR B1508+55;

[0070] Figure 5 is the EWT decomposition of the cumulative pulse profile of PSR B1508+55 EWT;

[0071] Figure 6 is the frequency domain feature (EMA) curve;

[0072] Figure 7 is the denoised pulsar signal obtained by using the denoising method of the present application;

[0073] Figure 8 is the SNR comparison of the denoising results of PSR B1508+55 at different SNRs;

[0074] Figure 9 is the PrmsD comparison of the denoising results of PSR B1508+55 at different SNRs. DETAILED DESCRIPTION

[0075] The present application will be described in detail below in conjunction with the accompanying drawings and specific embodiments.

[0076] The present application provides a pulsar signal denoising method based on EWT, as shown in the following steps: Figure 1

[0077] Step 1, determine the type of pulsar, obtain the pulsar data from the EPN database; set the number of sampling points and the signal-to-noise ratio of the observed signal, resample and normalize the obtained pulsar data in the database to obtain the standard profile of the pulsar s(t); according to the signal-to-noise ratio of the observed signal, add noise n(t) to s(t) to obtain the cumulative pulse profile of the pulsar x(t), i.e. x(t) = s(t) + n(t), where t represents the time variable, n(t) is a Gaussian white noise with mean zero and variance σ 2 , the noise here is additive Gaussian white noise;

[0078] Step 2, EWT decomposition is performed on the cumulative pulse profile of the pulsar x(t), and the EWT algorithm will output M detail components and one approximation component according to the characteristics of the input signal, i.e. M+1 EMF components; each EMF component obtained is a time domain signal at different frequencies, denoted as {emf m (t), m = 0, 1, …, M}, where emf m (t) increases with m, the corresponding frequency becomes larger; in addition, emf0(t) is the approximation component, which does not need to be denoised;

[0079] The specific process of EWT decomposition of the cumulative pulse profile of the pulsar x(t) in step 2 is described as follows:

[0080] ​Step 2.1, Fourier transform of x(t) to get signal spectrum X(ω), the specific formula as follows:

[0081]

[0082] In formula (1), j represents the imaginary unit; t represents time; ω represents frequency, and ω ∈ [0, π];

[0083] Step 2.2, signal spectrum X(ω) is normalized to [0, π], and [0, π] is divided into M continuous frequency bands, and the mth frequency band is Λ m = [ω m-1 ,ω m ], where ω m is the boundary value of the mth frequency band; M+1 endpoints are needed to divide M continuous boundaries, and ω0=0, ω M = π; therefore, M-1 endpoints also need to be determined.

[0084] Step 2.3, M-1 endpoints are determined using the maximum value-based Fourier analysis method, and the specific method is as follows:

[0085] First, find N local maximum values in signal spectrum X(ω), that is, N local peaks, where N ≥ M. The frequency corresponding to each local peak is called a local peak frequency, denoted as {Ω j , j = 1, 2, …, N};

[0086] Second, arrange the N local peaks in descending order of frequency, retain the first M local peaks, and remove the redundant local peaks; then reorder the local peak sequence in ascending order of frequency, and the frequency is denoted as {Ω j , j = 1, 2, …, M};

[0087] Third, the median of the frequencies corresponding to the adjacent two peak points is taken as the boundary point, and formula (2) is the calculation formula of the mth frequency band for the boundary value ω m ; finally, [0, π] is divided into M continuous intervals, and

[0088]

[0089] Finally, take ω m as the center frequency, define T m as the transition band, T m = 2τ m , where τ m =γω m ; wherein, parameter γ ∈ (0, 1), and satisfies the following conditions:

[0090]

[0091] Step 2.4, after determining Λ m an empirical mode function and M empirical wavelet functions φ m (t); where, is a low-pass filter, φ m (t) is a band-pass filter; their expressions are as follows:

[0092]

[0093]

[0094] The function β(x) in formula (4) and formula (5) satisfies the following conditions:

[0095]

[0096] There are multiple polynomials that satisfy the second condition in formula (6), and here we will use the following formula:

[0097] β(x) = x 4 (35-84x+70x 2 -20x 3 ) (7)

[0098] Step 2.5, take the inner product of x(t) and the empirical mode function to get the approximation coefficient term W f (1,t); take the inner product of x(t) and the empirical wavelet function φ m (t) to get the mth detail coefficient term, the specific formula is as follows:

[0099]

[0100] W f (m,t) = <x(t), φ m (t)> (9)

[0101] In formula (8) and formula (9), <a,b> represents the inner product of a and b; m takes values from 1 to M; t represents the time variable;

[0102] Step 2.6, calculate M+1 EMF components, including 1 approximation component and M detail components, the specific form is as follows:

[0103]

[0104] emf m (t) = W f (m,t) * φ m (t) (11)

[0105] In formula (10), emf0(t) is an approximate component; in formula (11), emf m (t) is the mth detailed component, m = 1, 2, …, M.

[0106] Signal reconstruction is performed using the M+1 EMF components as follows:

[0107]

[0108] The pulsar accumulative pulse profile x(t) = s(t) + n(t), where s(t) is a standard pulse profile, and n(t) is a Gaussian white noise with a mean of zero and a variance of σ 2 The noise here is additive Gaussian white noise, and emf m (t) in formula (11) is expressed as:

[0109] emf m (t) = <s(t), φ m (t)> * φ m (t) + <n(t), φ m (t)> * φ m (t) = emf i s (t) + emf n i (t) (13)

[0110] In formula (13), emf i s (t) represents the mth EMF component of s(t); emf i n (t) represents the mth EMF component of n(t); the larger the value of m here, the higher the frequency corresponding to emf m (t). Therefore, as m increases, there is more and more noise in emf m (t). That is, there is more noise in high-frequency EMF components, and less noise in low-frequency EMF components. As can be seen from formula (13), the signal and noise in each EMF component are superimposed together, and it is difficult to filter out high-frequency EMF components directly using {emf m (t), m = 1, 2, …, M}.

[0111] Step 3, FFT analysis is performed on the M EMFs {emf m (t), m = 1, 2, …, M} to obtain M spectra, and the corresponding M spectral components are denoted as {EMF m (f), m = 1, 2, …, M}, where f represents the frequency variable;

[0112] The expression of EMF m (f) is as follows:

[0113] EMF m (f) = FFT[emf m (t)], m = 1, 2, …, M (14)

[0114] In formula (14), FFT[a] represents performing FFT operation on a; EMF m (f) is the corresponding FFT result of emf m (t); f represents a frequency variable.

[0115] Step 4, calculating the mean value of the mth layer EMF m (f) and taking the modulus of the mean value, denoted as EMA m , where EMA m is the frequency domain feature of the mth layer EMF component; the frequency domain features are calculated for emf1(t)~emf M (t) to obtain a frequency domain feature set {EMA m , m = 1, 2, …, M};

[0116] In step 4, the feature value EMA m of the mth layer EMF m (f) is calculated, and EMA m has the following form:

[0117] EMA m = |mean[EMF m (f)]|, m = 1, 2, …, M (15)

[0118] In formula (15), mean[a] represents calculating the mean value of a, and |b| represents taking the modulus of b; EMA m is the feature value corresponding to EMF m (f).

[0119] Step 5, using {EMA m , m = 1, 2, …, M} to determine the starting layer number K of the high-frequency EMF component, which is specifically: finding the maximum value in {EMA m , m = 1, 2, …, M}, and assuming that the layer number corresponding to this maximum value is K, then emf1(t)~emf K-1 (t) is a low-frequency component, and emf K (t)~emf M (t) is a high-frequency component.

[0120] In step 5, {EMA m , m = 1, 2, …, M} is used to determine the starting layer number K of the high-frequency EMF component, and the expression of K is as follows:

[0121]

[0122] Equation (16) represents finding the maximum eigenvalue in the frequency domain feature set {EMA m , m = 1, 2, …, M}, and the layer corresponding to the maximum eigenvalue is the K value.

[0123] Step 6, the high-frequency EMF contains a large amount of noise, and the Mth layer EMF contains the most noise; therefore, the Mth layer EMF component emf M (t) can be directly removed as noise;

[0124] Step 7, in emf K (t) ~ emf M-1 (t), there is a problem of signal and noise aliasing in these high-frequency EMF components. Here, the high-frequency EMF signal corresponds to the place where the pulsar signal changes dramatically, i.e., the detailed information of the pulsar signal. In order to effectively suppress the aliasing noise component in the high-frequency EMF signal and retain the high-frequency EMF signal component, the wavelet threshold method is used to denoise emf K (t) ~ emf M-1 (t), and the denoised high-frequency EMF is denoted as

[0125] In step 7, the wavelet threshold method is used to denoise emf K (t) ~ emf M-1 (t), and the method is as follows: a Sym8 wavelet basis is selected to perform 3-layer wavelet decomposition on the EMF component; after decomposition, a soft threshold function is used to denoise the detail component, and the threshold is determined by the Heursure method; finally, the denoised high-frequency EMF component is reconstructed, and the denoising result is denoted as

[0126] Step 8, in emf1(t) ~ emf K-1 (t), there is a problem of signal and noise aliasing in these high-frequency EMF components. Here, the high-frequency EMF signal corresponds to the place where the pulsar signal changes dramatically, i.e., the detailed information of the pulsar signal. In order to effectively suppress the aliasing noise component in the high-frequency EMF signal and retain the high-frequency EMF signal component, the wavelet threshold method is used to denoise emf1(t) ~ emf K-1 (t), and the denoised high-frequency EMF is denoted as K-1 (t).

[0127] In step 8, when filtering, a 3rd order polynomial is used to fit the data in the window, and the filter window width is 41.

[0128] Step 9, using and the approximate component emf0(t) to reconstruct the de-noised pulsar signal

[0129] In step 9, the de-noised pulsar signal is reconstructed The specific formula is as follows:

[0130]

[0131] Embodiment

[0132] The computer simulation is performed by taking the pulsar data in the EPN database as an example. The pulsar PSR B1508+55 is selected for simulation analysis. The observation frequency of PSR B1508+55 is 0.61 GHz, and the sampling point number is 1024. Figure 2 is the standard pulse profile of PSR B1508+55; Figure 3 is the cumulative pulse profile of PSR B1508+55, and the signal-to-noise ratio of the observation signal is 20 dB.

[0133] The EWT decomposition is performed on the noisy signal in Figure 3 The EMF components obtained after decomposition are shown in Figure 4 It can be seen from the figure that 7 layers of EMF, i.e., emf0(t)~emf6(t), are obtained after decomposition, wherein emf0(t) is the approximate component, which shows the profile information of the observation signal and does not need to be de-noised. The 6 layers of EMF components emf1(t)~emf6(t) are detail components, Figure 4 The observation frequency of the EMF components from top to bottom in

[0134] Next, the high-frequency EMF components are screened out by using the frequency domain features. First, the frequency spectrum of each EMF component in Figure 4 is calculated, and the results are shown in Figure 5 By comparing the abscissa of each frequency spectrum in Figure 5 it can be seen that the frequency bands corresponding to the frequency spectra of the EMF components from top to bottom are increasingly larger. In addition, the frequency band widths of EMF4(f)~EMF6(f) are larger.

[0135] The frequency domain features EMA1~EMA6 of EMF1(f)~EMF6(f) in Figure 5 are calculated. Figure 6 The frequency domain feature curves are drawn, wherein the abscissa is the layer number 1~6, and the ordinate is the frequency domain features EMA1~EMA6 corresponding to the 1~6 layers of EMF components. From Figure 6It can be seen intuitively that the maximum value of the six frequency domain features is EMA4, i.e. K=4. Therefore, emf1(t)~emf3(t) are low frequency components, and emf4(t)~emf6(t) are high frequency components. This result is consistent with the analysis in Figure 4 and Figure 5 , proving that the high frequency EMF components can be effectively screened out by using the frequency domain features.

[0136] From Figure 4 , it can be seen that, among the three high frequency EMF components emf4(t)~emf6(t), emf6(t) contains a large amount of noise components and can be directly removed as noise components; and emf4(t) and emf5(t) have noise and signal mixed. Since the signal components in the high frequency EMF components emf4(t) and emf5(t) contain the detailed information of the pulsar signal, when suppressing the mixed noise components in emf4(t) and emf5(t), the signal components in them also need to be preserved. Here, the wavelet threshold method is selected to denoise emf4(t) and emf5(t). Among them, the wavelet basis is selected as Sym8, the decomposition level is 3, the threshold function is selected as the soft threshold function, and the threshold is determined by the heuristic threshold method. The high frequency EMF obtained after denoising is denoted as

[0137] From Figure 4 , it can be seen that, among the three low frequency EMF components emf1(t)~emf3(t), only a small amount of noise components are contained. Since the signal components in the low frequency EMF components correspond to the outline information of the pulsar signal, here the Savitzky-Golay smoothing filter method is selected to denoise emf1(t)~emf3(t). When denoising, a 3rd order polynomial is used to fit the data in the window, and the filter window width is 41. The high frequency EMF obtained after denoising is denoted as

[0138] Finally, the denoised pulsar signal is reconstructed by using the approximate component emf0(t), and the time domain waveform of the denoised pulsar signal is shown in Figure 7 . The signal-to-noise ratio of the denoised pulsar signal is 30.49 dB, which is increased by 10.49 dB.

[0139] Comparing the observation signal of PSR B1508+55 in Figure 3 with the denoised waveform in Figure 7 , it can be seen that the noise in the pulse part and the non-pulse part is effectively removed; comparing the standard pulse outline of PSR B1508+55 in Figure 2 with the denoised waveform in Figure 7The comparison of the denoised waveforms in FIG. 8 shows that the pulse information in the PSR B1508+55 signal is well preserved in the denoised signal, Figure 7 The denoising results in FIG. 8 are similar to the standard pulse profile of the PSR B1508+55 in FIG. 7. In summary, the method can effectively suppress noise while preserving signal detail information, effectively improving the quality and signal-to-noise ratio of the observed signal. Figure 2

[0140] The performance of the method is verified below using Monte Carlo simulation. The standard pulse profile of the PSR B1508+55 in FIG. 7 is added with different noise, and the signal-to-noise ratio of the noisy signal is set to be -10dB, -5dB, 0dB, 5dB, 10dB, 15dB, 20dB, 25dB and 30dB. At each specified signal-to-noise ratio, 100 noisy samples are taken, and each sample is denoised using the method. The signal-to-noise ratio SNR and the percentage root mean square difference PrmsD of each denoised sample signal are counted, and the SNR mean and PrmsD mean at each specified signal-to-noise ratio are calculated. The results are shown in FIG. 9 and FIG. 10. Figure 2 Figure 8 Figure 9

[0141] Figure 8 The horizontal axis of FIG. 9 is the SNR of the noisy sample signal, and the vertical axis is the SNR mean of the different denoised sample signals. The SNR mean is obtained by projecting each blue solid circle to the vertical axis. In FIG. 9, the SNR mean obtained under the condition of -10dB is 2.12dB, and the signal-to-noise ratio is improved by 12.12dB; under the condition of -5dB, the SNR mean obtained is 7.13dB, and the signal-to-noise ratio is improved by 12.13dB; under the condition of 0dB, the SNR mean obtained is 12.11dB, and the signal-to-noise ratio is improved by 12.11dB; under the condition of 5dB, the SNR mean obtained is 17.08dB, and the signal-to-noise ratio is improved by 12.08dB; under the condition of 10dB, the SNR mean obtained is 21.97dB, and the signal-to-noise ratio is improved by 11.97dB; under the condition of 15dB, the SNR mean obtained is 26.58dB, and the signal-to-noise ratio is improved by 11.58dB; under the condition of 20dB, the SNR mean obtained is 30.5dB, and the signal-to-noise ratio is improved by 10.5dB; under the condition of 25dB, the SNR mean obtained is 33.14dB, and the signal-to-noise ratio is improved by 8.14dB; under the condition of 30dB, the SNR mean obtained is 34.47dB, and the signal-to-noise ratio is improved by 4.47dB. It can be seen that under the condition of low signal-to-noise ratio, the method can effectively improve the signal-to-noise ratio of the observed signal. Figure 8

[0142] Figure 9 ​​​​​The horizontal axis of the figure is the SNR of the noisy sample signal, and the vertical axis is the PrmsD mean value of the different denoised sample signals. The PrmsD mean value is obtained by projecting each blue solid circle to the vertical axis. Figure 9 In the-10 dB condition, the obtained PrmsD mean value is 0.61; in the-5 dB condition, the obtained PrmsD mean value is 0.19 dB; in the 0 dB condition, the obtained PrmsD mean value is 0.061; in the 5 dB condition, the obtained PrmsD mean value is 0.019; in the 10 dB condition, the obtained PrmsD mean value is 0.006; in the 15 dB condition, the obtained PrmsD mean value is 0.002; in the 20 dB condition, the obtained PrmsD mean value is 0.0008; in the 25 dB condition, the obtained PrmsD mean value is 0.0004; and in the 30 dB condition, the obtained PrmsD mean value is 0.0003. It can be seen that when the SNR is-5 dB, the noise is effectively suppressed by the method. When the SNR is 0 dB, the corresponding PrmsD mean value is very small, and the noise content in the denoised signal is very small. It can be seen that under the condition of low SNR, the method can effectively suppress the noise component in the observation signal while retaining the signal component.

[0143] In summary, when the flux density is low and the SNR of the pulsar observation signal is low, the method can effectively improve the SNR of the pulsar cumulative profile and improve the quality of the pulsar cumulative profile, thereby providing effective data support for various applications of the pulsar signal.

Claims

1. A method for pulsar signal denoising based on EWT, characterized in that, The method is implemented according to the following steps: Step 1, obtaining a cumulative pulse profile x(t) of a pulsar, where t represents a time variable; The specific implementation of step 1 is: determining the type of pulsar, obtaining the pulsar data from the EPN database; setting the number of sampling points and the signal-to-noise ratio of the observation signal, resampling and normalizing the obtained pulsar data in the database to obtain a standard profile s(t) of the pulsar; adding noise n(t) to s(t) according to the signal-to-noise ratio of the observation signal to obtain a cumulative pulse profile x(t) of the pulsar, that is, x(t) = s(t) + n(t), wherein t represents a time variable, n(t) is Gaussian white noise with a mean of zero and a variance of σ 2 . Step 2, EWT decomposition of the pulsar accumulated pulse profile x(t) to get {emf m (t), m = 1, 2,..., M} ; The specific implementation of step 2 is: Step 2.1, performing Fourier transform on x(t) to obtain a signal spectrum X(ω), and the specific formula is as follows: In formula (1), j represents a virtual unit; t represents time; ω represents frequency, and ω ∈ [0, π]; Step 2.2, the signal spectrum X(ω) is normalized to [0, π], and [0, π] is divided into M continuous frequency bands, the mth frequency band is Λ m = [ω m-1 ,ω m ], where ω m is the boundary value of the mth frequency band; M+1 endpoints are needed to divide M continuous boundaries, and ω0=0, ω M =π; Step 2.3, determining M-1 endpoints using a maximum value-based Fourier analysis method, and the specific method is as follows: First, N local maximum values in the signal spectrum X(ω) are found, i.e. N local peaks, where N≥M; the frequency corresponding to each local peak is called a local peak frequency, denoted as {Ω j , j = 1, 2, …, N}; Secondly, arrange the N local peaks in descending order of frequency, retain the first M local peaks, and eliminate the redundant local peaks; then re-sort the local peak sequence in ascending order of frequency, and the frequency is denoted as {Ω j j = 1, 2, …, M}. Again, the median of the frequencies corresponding to the adjacent 2 peak points is taken as the boundary point, and formula (2) is the calculation formula of the mth frequency band for the boundary value ω m ; finally, [0, π] is divided into M continuous intervals, and there are Finally, with ω m as the center frequency, define T m as the transition band, T m = 2τ m , where τ m =γω m ; where the parameter γ ∈ (0, 1) and satisfies the following conditions: Step 2.4, after determining Λ m an empirical scaling function and M empirical wavelet functions φ m (t); where, is a low-pass filter, φ m (t) is a band-pass filter; their expressions are as follows: The function β(x) in formula (4) and formula (5) satisfies the following conditions: There are multiple polynomials that satisfy the second condition in formula (6), and the following formula will be used here: β(x) = x 4 (35 - 84x + 70x 2 - 20x 3 ) (7) Step 2.5, take the inner product of x(t) with the empirical scaling function to obtain the approximation coefficient term W f (1,t); take the inner product of x(t) with the empirical wavelet function φ m (t) to obtain the mth detail coefficient term, which is given by W f (m,t) = <x(t), φ m (t)> (9) In formula (8) and formula (9), <a, b> represents the inner product of a and b; m takes a value of 1 to M; t represents a time variable; Step 2.6, calculating M+1 EMF components, including 1 approximate component and M detail components, and the specific form is as follows: emf m (t) = W f (m, t) * φ m (t) (11) In Equation (10), emf0(t) is an approximate component; in Equation (11), emf m (t) is the mth detailed component, m = 1, 2,..., M; The signal is reconstructed using the M+1 EMF components, as shown below: The pulsar accumulates the pulse profile x(t) = s(t) + n(t), where s(t) is the standard pulse profile and n(t) is the Gaussian white noise with zero mean and variance σ 2 Here, the noise is additive Gaussian white noise, and emf(t) in equation (11) is expressed as: m (t) = emf(t) + n(t) emf m (t) = <s(t), φ m (t)> * φ m (t) = <s(t), φ m (t)> * φ m (t) = emf i s (t) = emf n i (t) (13) In equation (13), emf i s (t) represents the mth EMF component of s(t); emf i n (t) represents the mth EMF component of n(t); Step 3, on {emf m (f), m = 1, 2, …, M} with M EMFs, to obtain M spectral components {EMF m (f), m = 1, 2, …, M}, where f represents the frequency variable; Step 4, calculate the frequency domain feature set {EMA m m = 1, 2,..., M} Step 5, determine the starting layer number K of the high-frequency EMF component using {EMA m m = 1, 2,..., M} Step 6, divide the Mth layer EMF component emf M (t) as noise direct rejection; Step 7, denoising emf (t) using wavelet thresholding method K (t) ~ emf M-1 (t) using wavelet thresholding method Step 8, denoising the emf1(t)~emf K-1 (t) these low frequency EMF components, and the denoising result is denoted as Step 9, reconstructing using and the approximate component emfo(t) to obtain the denoised pulsar signal 2. The EWT-based pulsar signal denoising method of claim 1, wherein, In step 3, EMF m The expression of (f) is as follows: EMF m (f) = FFT[emf m (t)], m = 1, 2,..., M (14) In Equation (14), FFT[a] represents an FFT operation on a; EMF m (f) emf m (t) corresponding FFT result; f represents a frequency variable.

3. The EWT-based pulsar signal denoising method of claim 2, wherein, In step 4, EMA m is in the following form: EMA m = |mean[EMF m (f)]|, m = 1, 2,..., M (15) In formula (15), mean[a] represents averaging a, and |b| represents taking the modulus of b; EMA m for EMF m (f) the corresponding eigenvalue.

4. The EWT-based pulsar signal denoising method of claim 3, characterized in that, In Step 5, the starting layer number K of the high-frequency EMF component is determined using {EMA m , m = 1, 2, …, M} and the expression of K is as follows: Equation (16) represents finding the largest eigenvalue in the frequency domain feature set {EMA m , m = 1, 2, …, M}, and the layer number corresponding to the largest eigenvalue is the K value.

5. The EWT-based pulsar signal denoising method of claim 4, characterized in that, In step 7, the wavelet thresholding method is used to evaluate the EMF. K (t)~emf M-1 The noise reduction method for (t) is as follows: The Sym8 wavelet basis is selected to perform a 3-level wavelet decomposition on the EMF components; after decomposition, a soft thresholding function is used to denoise the detail components, with the threshold determined heuristically; finally, the noise-removed high-frequency EMF components are reconstructed, denoted as...

6. The EWT-based pulsar signal denoising method according to claim 5, characterized in that, In step 9, the de-noised pulsar signal is reconstructed The specific formula is as follows:

Images