Satellite navigation system’s main function is to realize the positioning, velocity and timing, Global Position System is an important part of GNSS, GPS L5 signal spectrum and B2 signal common Beidou Satellite Navigation System in our country, because there have been many radio navigation systems in the work in this band, so the GPS L5 receiver is facing the complex electromagnetic environment, DME equipment is one of the largest interference sources to the interference. In severe cases, GPS receiver cannot be accurately positioned, which poses a potential threat to the implementation of PBN of civil aviation. Currently, there are three main methods to deal with interference suppression of DME equipment, namely, time domain processing method, frequency domain processing method, and time-frequency domain combined processing method .
The time domain processing method is the most common pulse interference suppression method, which is to suppress the interference by setting a threshold in the time domain consistent with the noise level and zeroing the received signal samples that exceed the threshold. In the case of high pulse interference density, this method will lose a large number of useful satellite signals.
Compared with the time-domain processing method, the frequency-domain filtering method can filter out the interference more effectively, but at the same time, there will also be cases of filtering out useful signals . Due to the variability of aircraft flight environment, the frequency-domain processing method based on notch filter requires the construction of different notch filters, which is not feasible in the changeable flight environment.
Compared with the former two methods, the time-frequency mixing method can save more useful signals to suppress the interference. How to detect the center position of the pulse is the key of this method, but if there is a very bad interference environment, it is difficult to achieve this method.
To sum up, the traditional time-frequency analysis is unable to filter out the interference signals while preserving the useful signals as much as possible. Therefore, based on the detection and suppression of DME interference in the domain of wavelet packet transform, this paper further systematically studies and analyzes the function type, decomposition order and the selection of threshold value.
2. Signal Model
2.1. How the DME System Works
DME is the ranging part of TACAN adopted by the International Civil Organization. As a new international standard ranging system, DME has important applications in navigation, such as positioning service, speed and flight time, airway spacing and regional navigation.
DME consists of two parts, an airborne interrogator and a ground transponder, and calculates the distance by measuring the time of electromagnetic wave transmission in space. The block diagram of its basic working principle is shown in Figure 1.
After the pulse output by the oscillator passes the frequency division, a chatter rectangular pulse with random intervals is obtained; after the coding by the encoder, a Gauss pulse pair with fixed time intervals is generated; after the RF modulation, the RF query signal is transmitted through a non-directional antenna. As shown in Figure 2, after the signal is DME ground receiving platform,
Figure 1. Block diagram of basic working principle of DME system.
Figure 2. DME measurement of oblique distance between aircraft and DME platform.
the ground of the transponder decoder will according to the rf frequency and pulse interval between rf interrogation signal of acceptable, if the effective inquiry, decoder can produce a decoding pulse, the radio frequency modulated by omnidirectional antenna radiation also get the corresponding plane to ask signals of question and answer. The oblique distance R from the aircraft to the ground DME platform is
where, C is the speed of light, T0 refers to the time between the ground transponder receiving the question signal and sending the question signal, and T refers to the time difference recorded by the interrogator from sending the question signal to receiving the answer signal.
2.2. DME Pulse Signal Model
The coding format of the interrogation signal and the response signal of DME is Gaussian pulse pair. The mathematical model of a baseband DME pulse pair signal is as follows 
Among them, the half-amplitude pulse width of both Gauss pulses is 3.5 μs, and the modulated DME signal model can be expressed as
Figure 3 shows the DME/TACAN standard pulse pair generated according to the mathematical model. The transmitted DME/TACAN signal is composed of such pulse pair continuously , and the arrival time of the pulse pair is random.
2.3. DME Signal Quality Index
The main performance indexes of DME beacon transmitter are as follows.
1) Output power: the output power of the terminal DME beacon is more than 100 W, and the output power of the airway DME beacon is more than 1 KW.
2) Signal form: it has a Gaussian envelope. Since the shape of the envelope will affect the signal spectrum, the interference size of the adjacent channel will be affected.
3) Transmission rate: The frontal response pulse repetition frequency range of DME beacon transmitting base is 2700 - 700 pairs/s, the main performance index of DME beacon receiver is sensitivity, which refers to the receiver’s ability to accept weak signals, and the minimum detectable signal power is its measurement data. ICAO’s requirement for DME transponder receiver sensitivity is the need for query pulse pairs with correct intervals and nominal frequencies. When its peak power density at the transponder antenna is −103 dB W/m2, the transponder responds with at least 70% response efficiency. Under the condition of the highest sensitivity and good station environment, if the peak pulse power emitted by the DME beacon station is 1 kW, its operating distance can reach 200 nm (370 km).
Figure 3. Baseband DME waveform.
2.4. Effects of DME Pulse Interference on GNSS
The transmitting frequency of the response pulse pair in the DME platform working in 64× - 126× channel overlaps with the frequency band range of the signal received by the airborne GNSS receiver . These interference pulses will enter the GNSS receiver and reduce the signal-to-noise ratio of the signal, thus reducing the sensitivity of the user when capturing GNSS signals . This will further increase the GNSS signal tracking failure rate. When it is serious, the availability of observation data will be reduced, and even cycle slip will be caused, resulting in greater errors in positioning results.
The second effect of interference on the GNSS receiver is that it will increase the observed code noise and phase noise, thereby reducing the accuracy of the GNSS user terminal. In addition, the interference will also cause the acquisition time of the GNSS receiver to prolong , resulting in GNSS the user’s first positioning time is extended. The mathematical formula of the received signal of a GNSS signal receiver affected by a DME station can be shown as 
Among them, is the received signal, is the DME baseband signal , is the total logarithm of the impulse transmitted from the m-th DME station received during the observation time, is the received power of the m-th DME station signal , and is the arrival time of the U pulse pair generated by the m-th DME platform.
3. DME Pulse Interference Suppression Method Based on Wavelet Packet Decomposition
3.1. The Difference between Wavelet and Wavelet Packet
Wavelet analysis and wavelet packet analysis are both suitable for the nonstationary signal analysis. Compared with smaller wave analysis , it can not only decompose the low-frequency part of the signal, but also the high-frequency part . The wavelet packet analysis can be used to analyze the signal containing a large amount of medium and high frequency information. Perform better time frequency localization analysis, and adaptively select the best wavelet basis function according to the characteristics of the signal, thereby improving the time frequency resolution.
In the multi-resolution analysis based on wavelet transform, only the low frequency sub-band of the signal is decomposed at each level , and the high frequency sub-band is not decomposed. Therefore, the high frequency band resolution is poor. If the space corresponding to the low frequency band in the wavelet analysis is regarded as , the space corresponding to the high frequency part is regarded as , the high and low frequency bands here are relative to the frequency band of the upper level space; the wavelet transform only divides each space into two. The wavelet packet transform can divide the low-frequency space and the high-frequency space at the same time . Relative to wavelet decomposition, the signal can be decomposed more accurately and the time-frequency resolution can be improved.
Set as a multi-resolution analysis of , the corresponding scaling function is , and the wavelet is , satisfying the scaling equation 
, which called wavelet decomposition of space, , , , , etc, which are called wavelet decomposition of function.
Wavelet decomposition refers to decomposing on the basis of decomposing in a certain way.
Figure 4 shows the three-level wavelet packet decomposition tree diagram of the original signal, where S represents the low-frequency coefficient obtained by the signal through the low-pass filter , and D represents the high frequency coefficient obtained by the signal through the high-pass filter. Each layer of wavelet packet can divide the original frequency band into two, and the three-layer wavelet packet can divide the original frequency band into 2K sub-bands, thereby realizing frequency band subdivision and improving frequency domain resolution. For the entire wavelet packet, it is a binary-organized band-pass filter bank from wide to narrow, and various applications can find the best combination that meets the needs.
3.2. Principle of Interference Suppression Algorithm Based on Wavelet Packet Transform
The essence of wavelet transform is the similarity between the original signal and
Figure 4. Schematic diagram of three-level wavelet packet decomposition.
the wavelet basis function. The wavelet coefficient is the coefficient of the wavelet basis function similar to the original signal. In the problem of interference signal suppression, since the GNSS signal is submerged in noise, the similarity between noise and wavelet is less than that between wavelet and DME signal, so the received signal interfered by DME impulse can be passed through wavelet. The packet change is transformed into the wavelet packet coefficient domain. The coefficient with larger absolute value is expressed as DME interference signal , and the coefficient with smaller absolute value is expressed as noise and satellite signal. At the same time, setting an appropriate threshold in the wavelet packet domain will make the absolute value larger the coefficient is set to zero, so as to achieve the purpose of suppressing the pulse interference signal. Finally, the processed wavelet packet coefficient is inversely transformed to obtain the signal in the time domain. The basic steps are as follows :
1) Wavelet packet decomposition
Choose the appropriate wavelet and the appropriate decomposition level (denoted as K), and then perform K-level decomposition on the received signal interfered by DME.
2) Threshold quantization of high-frequency coefficients after decomposition
Set an appropriate threshold in the wavelet packet coefficient domain, and zero the coefficients whose absolute value is greater than this threshold to suppress impulse interference signals .
3) Inverse wavelet packet transform.
The K-th wavelet packet coefficient is inverse transformed to get the time domain waveform.
Because the satellite signal contains a lot of noise, the satellite signal and the noise can be viewed together as Gaussian noise with zero mean value ε.
Suppose the variance of noise ε is , and its power spectral density is . The variance of the wavelet coefficients after the first-level wavelet packet transformation of the filter is . The frequency response of is 
Among them, is the unit impulse response sequence, and the power spectral density of is:
The autocorrelation function of is:
Suppose the probability density function of white noise with zero mean variance is , then there is:
Suppose the probability density function of white noise with zero mean variance Th is , then there is:
According to the above formula, the threshold value is:
3.4. The Choice of Wavelet Packet Decomposition Series
Each layer of wavelet packet divides the original frequency band into two, and the k-layer wavelet packet can divide the original frequency band into 2K sub-bands, thereby realizing frequency band subdivision and improving frequency domain resolution. That is, the larger the decomposition level, the finer the division of the signal, and the better the interference suppression effect. When the decomposition gets to a certain level, continuing the decomposition will not make the performance significantly improved. In addition, Due to the wavelet packet filter is not ideal, the increase in the number of decomposition stages will cause aliasing of the subband spectrum. Therefore, the choice of series of the wavelet packet should take into account the bandwidth of the interference signal and the amount of calculation. That is to say, wavelet packet resolution after decomposition is equivalent to the frequency band of the suppressed interference signal.
3.5. Selection of Wavelet Function Type
In numerical analysis, we knowed that the essence of wavelet transform is the similarity between the original signal and the wavelet basis function. The wavelet coefficient is the coefficient of the wavelet basis function similar to the original signal. Therefore, the more similar the generator function is to the analyzed signal, the more local the coefficients after the wavelet transform can be, and the better the interference suppression effect. In the actual process, there is no display expression in a few wavelet functions, and it is difficult to give a certain degree of similarity between the wavelet function and the analyzed function. This experiment adopts the general selection criteria of wavelet function, in order to be able to perform multi-resolution analysis and apply the fast algorithm of wavelet transform.
In the paper, we consider the Daubechies wavelet function family, the Symlet wavelet function family, the Coifle wavelet function family, the bioorthogonal wavelet function family , the rbio inverse biorthogonal wavelet function family and the discrete Meyer wavelet function “dmey”. According to the capture performance of the signal after interference suppression, the effects of various wavelet functions in suppressing interference are compared.
4. Simulation Results and Analysis
This experiment is based on wavelet packet decomposition algorithm to suppress DME pulse interference. Under certain false alarm rate condition, uniform thresholds are set for comparison of the impact of different wavelet function types and levels on the interference suppression performance. The experimental equipment of this experiment is shown in Figure 5. In the experiment, a more real DME interference environment is simulated. The repetition frequency of each DME pulse pair is 2000 pairs/sec. The interval between pulses obeys the Poisson distribution, and the interference duty cycle is 37%. And add the generated DME pulse interference signal to the real GPS samples collected by the GNSS receiver.
Figure 8 shows the GNSS time domain waveform after impulse interference suppression. Obviously, The algorithm proposed in this paper can effectively suppress pulse jamming.
This paper uses the alpha_mean index (the ratio of the maximum value to the mean value after removing the maximum value) to measure the suppression performance of the pulse interference signal. Figure 9 is a graph showing the
Figure 5. Experimental setup of the measurement test in the lab.
Figure 6. GNSS signal waveform without impulse interference.
Figure 7. GNSS signal waveform with DME impulse interference.
Figure 8. GNSS signal after pulse interference suppression.
Figure 9. Comparison of interference suppression performance of wavelet packet transform.
results of interference suppression performance using different decomposition levels and wavelet functions.
This paper mainly studies the wavelet packet decomposition and reconstruction and the actual algorithm applied to suppress DME pulse interference. Through real experimental results, the correctness of the impulse interference suppression algorithm based on wavelet packet transform is verified. In addition, we compare the interference suppression performance under different decomposition orders and wavelet packet function, and the optimal wavelet function and corresponding to suppress DME impulse interference are found. The number of decomposition stages provides a reference for the parameter selection of wavelet packet transform for DME pulse interference suppression.
 Cheng, W. and Zhou, C.D. (2012) The Feature Extraction of Rolling Bearing Fault Based on Wavelet Packet-EMD Energy Distribution. Applied Mechanics and Materials, 233, 234-238.
 Liu, J.B., Wu, D.P., Wang, Z.M., Jin, X.Y., Dong, F., Jiang, L.R. and Cai, C.Y. (2020) Automatic Sleep Staging Algorithm Based on Random Forest and Hidden Markov Model. Computer Modeling in Engineering & Sciences, 123, 401-426.
 Zhao, J. (2010) Fault Diagnosis of Inter-Turn Short-Circuit in Rotor Windings Based on Artificial Intelligence. 2010 2nd IEEE International Conference on Information and Financial Engineering, Chongqing, 17-19 September 2010, 617-620.
 Zuo, L.-Q., Sun, H.-M., Mao, Q.-C., Liu, X.-Y. and Jia, R.-S. (2019) Noise Suppression Method of Microseismic Signal Based on Complementary Ensemble Empirical Mode Decomposition and Wavelet Packet Threshold. IEEE Access, 7, 176504-176513.
 Wu, R.B., Wang, W.Y., Li, L.L., Lu, D. and Wang, L. (2016) Distance Measuring Equipment Interference Suppression Based on Parametric Estimation and Wavelet-Packet Transformation for Global Navigation Satellite Systems. IEEE Transactions on Aerospace and Electronic Systems, 52, 1607-1617.
 Liu, X.F., Zhang, J., Liu, R., Shu, X.Y., Chen, M. and Tian, T.T. (2018) Identification of Epileptic Electroencephalograph in Frequency Bands. 2018 14th International Conference on Natural Computation, Fuzzy Systems and Knowledge Discovery (ICNC-FSKD), Huangshan, 368-372.
 Yu, B., Lou, L.F., Li, S., Zhang, Y.S., Qiu, W.Y., Wu, X., Wang, M.H. and Tian, B.G. (2017) Prediction of Protein Structural Class for Low-Similarity Sequences Using Chou’s Pseudo Amino Acid Composition and Wavelet Denoising. Journal of Molecular Graphics and Modelling, 76, 260-273.
 Liu, Y.L., Jiang, Y.Z., Jiang, W. and Feng, Y.Q. (2011) Estimation Algorithm of Doppler Shift and Time Delay in HF Channel Sounding. IEEE 2011 10th International Conference on Electronic Measurement & Instruments, Chengdu, 16-18 August 2011, 350-354.
 Gassama, S., Sonnendrücker, é., Schneider, K., Farge, M. and Domingues, M.O. (2007) Wavelet Denoising for Postprocessing of a 2D Particle-in-Cell Code. ESAIM: Proceedings, 16, 195-210.
 Hegarty, C., Kim, T., Ericson, S., et al. (1999) Methodology for Determining Compatibility of GPS L5 with Existing Systems and Preliminary Results. Proceedings of the Institute of Navigation Annual Meeting, Cambridge, June 1999, 635-644.
 He, H., Cheng, S., Zhang, Y. and Nguimbis, J. (2005) Home Network Power-Line Communication Signal Processing Based on Wavelet Packet Analysis. IEEE Transactions on Power Delivery, 20, 1879-1885.
 Sun, X.D. (2010) A Blind Digital Watermarking for Color Medical Images Based on PCA. 2010 IEEE International Conference on Wireless Communications Networking and Information Security, Beijing, 25-27 June 2010, 421-427.
 Zhang, J.P., Wang, X.F., Chen, T. and Zhang, Y. (2005) Change Detection for the Urban Area Based on Multiple Sensor Information Fusion. Proceedings 2005 IEEE International Geoscience and Remote Sensing Symposium, Seoul, 25-29 July 2005, 4.
 Yan, H.M., Mu, H.N., Yi, X.J., Yang, Y.Y. and Chen, G.L. (2019) Fault Diagnosis of Rolling Bearing with Small Samples Based on Wavelet Packet Theory and Random Forest. 2019 International Conference on Sensing, Diagnostics, Prognostics, and Control (SDPC), Beijing, 15-17 August 2019, 305-310.