An Efficient Algorithm for Improved Doppler Profile 1 Detection of MST Radar Signals 2

An efficient algorithm based on Empirical Mode Decomposition (EMD) de-noising using soft 10 threshold techniques for accurate doppler profile detection and Signal to Noise Ratio (SNR) improvement of 11 MST Radar Signals is discussed in this paper. Hilbert Huang Transform (HHT) is a time-frequency analysis 12 technique for processing radar echoes which constitutes EMD process that decomposes the non-stationary 13 signals into Intrinsic Mode Functions (IMFs). HHT process has been applied on the time series data of MST 14 (Mesosphere-Stratosphere–Troposphere) radar collected from NARL (National Atmospheric research 15 Laboratory), Gadanki, India. Further, spectral moments were estimated and signal parameters such as mean 16 doppler, signal power, noise power and SNR were calculated. Stacked doppler profile was plotted to observe the 17 improvement in doppler detection. It has been observed that there is a considerable improvement in recognition 18 of the doppler echo leading to improved Signal Power and SNR. The algorithm was tested for its efficacy on 19 various data sets for all the 6 beams and the results of two data sets are presented. 20


Introduction
Processing and analysis of radar echoes from MST region pose serious challenges to traditional signal processing techniques especially at higher altitudes above 15 Kms since the echo returns are very weak and buried in noise.The most common approach for analysis of atmospheric radar signals is the Fast Fourier Transform (FFT) and Wavelets.However, Fourier Transforms are not suitable for applications that involve nonlinear and non stationary signals as it has a serious drawback that in transforming to the frequency domain, time information is hidden and not explicitly visible.Wavelet Transform analysis, which is widely used method, overcomes some of the above limitations but there is a need for a-priori knowledge about the kind of scale elements present in the signals and the choice of a suitable mother wavelet for analysis since the accuracy of the results depends on the selection of the wavelet.Hence, an adaptive signal processing technique based on the Empirical Mode Decomposition (EMD) and de-noising using soft thresholding is proposed in this work.

Indian MST Radar System
The MST Radar facility at Gadanki (13.5°N, 79.2° E, 6.3° N mag.lat) is an excellent system used for atmospheric probing in the regions of Mesosphere, Stratosphere and Troposphere (MST) covering up to a height of about 100 km.It is used for coherent backscatter study of the ionospheric irregularities above 100 Km.MST radar is a state-of-the-art instrument capable of providing estimates of atmospheric parameters with very high

Hilbert Huang Transform
Hilbert Huang Transform (HHT) is one of the best time-frequency analysis technique for processing radar echoes.HHT basically comprises Empirical Mode Decomposition process based on numerical shifting to decompose the non-stationary signal into Intrinsic Mode Functions and obtain instantaneous frequency solution.
To obtain instantaneous frequency data, Hilbert Spectral Analysis (HSA) method is applied to the IMFs.Since the signal is decomposed in time domain and the IMFs length is the same as the original signal.HHT preserves the characteristics of the varying frequency.Also the, extraction of IMFs enables various de-noising techniques to be applied for accurate detection of doppler echo.This is the key benefit of HHT.It has been tested and validated comprehensively, but only empirically.In all cases studied, from any of the traditional analysis methods, HHT (Padmaja et al., 2011;Jing-tian et al., 2007) gave much sharper results in time-frequency-energy representations.HHT process was applied on time series MST radar data and was investigated for its efficacy.
The spectral moments were estimated and signal parameters such as mean doppler, signal power, noise power and SNR were calculated.

Empirical Mode Decomposition
The Empirical Mode Decomposition is an effective de-noising technique (Flandrin et al., 2004;Padmaja et al., 2017) that can be used to process non-linear and non-stationary data.This method is adaptive and intuitive.Each oscillatory mode is represented by an IMF that satisfies the conditions of an IMF (Huang et al., 2005).An IMF represents a simple oscillatory mode and it is a counterpart to the simple harmonic function, but it is much general instead of constant amplitude and frequency, as in a simple harmonic component.IMF can have a variable frequency and amplitude as function of time.

EMD Algorithm
Given a non-stationary signal x(t), the EMD algorithm (Rilling et al., 2003) can be summarized as follows: a) Find the local maxima and minima of the signal, then connect all the maxima and minima of signal X(t) and obtain the upper envelope X u (t) and the lower envelope X l (t) respctively.b) Compute the local mean value m 1 (t)=(X u (t)+X l (t))/2 of data X(t), subtract the mean value from signal X(t) and get the difference : h 1 (t)=X(t)-m 1 (t).c) Assume h 1 (t) as new data and repeat steps(1) and ( 2 is the mean value of h 1 (k-1)(t) and h 1 k(t).Step(c) is terminated untill the resulting data satisfies the two conditions of an IMF, defined as c 1 (t)=h 1 k.The residual data r 1 (t) is expressed as r 1 (t)= X(t)-c 1 (t).d) Assume r 1 (t) as new data and repeat steps (a-c) and extract all the IMFs.Terminate the sifting process untill n th residue r n (t) becomes less than a predetermined number or the residue becomes monotonic.
Geosci.Instrum.Method.Data Syst.Discuss., https://doi.org/10.5194/gi-2017-9Manuscript under review for journal Geosci.Instrum.Method.Data Syst.Discussion started: 19 June 2017 c Author(s) 2017.CC BY 3.0 License.e) Repeat steps (a-d) till the residual no longer contains any useful frequency information.The original signal is equal to the sum of its IMFs.If we have 'n' IMFs and a final residual r n (t), the original signal X(t) can be defined as follows

Intrinsic Mode Functions
After the application of EMD (George Tsolis et al., 2011;Flandrin et al., 2004;Norden et al., 1998;Wu et al., 2004;Dejie Yu et al., 2010), if the residue, r 1 still contains information of longer period components, then it is again treated as new data and subjected to the sifting process.The sifting process can be stopped by any of the predetermined criteria: either when the component value c n or the residue r n ,becomes less than the predetermined value or also when the residue, r n becomes a monotonic function from which no more IMFs can be extracted.The criterion for stopping the sifting process is based on limiting the size of the Standard Deviation (SD) ,which can be computed from the two consecutive sifting results as shown in the equation ( 2).

SD= ----------(2)
A typical value for SD is around 0.21 to 0.3.Hard and soft thresholding (similar to de-noising in Wavelets) is used to treat Intrinsic Mode Functions to achieve high SNR.

Denoising
Signal de-noising scheme based a multiresolution approach is referred to as empirical mode decomposition denoising (Padmaja et al., 2017;Donoho et al, 1995).A smooth version of the input signal can be obtained by thresholding the IMFs before signal reconstruction.Thresholding are of two types namely Hard and Soft Threshold.If Γ[τ j ] is a thresholding function, and τ j is the threshold parameter, the threshold can be determined in different ways.Donoho and Johnstone proposed a universal threshold, τ j for removing noise (Donoho et al, 1995).The method of soft threshold is applied to process the radar data.After extracting the Intrinsic Mode functions in each range bin, de-noising techniques are employed before reconstruction of the Doppler spectra by using threshold levels.

Hard Threshold
Hard threshold removes the corresponding IMFs depending on the frequencies if τ j is less than or equal to 1.The condition for hard threshold as shown in equation ( 3) and (4)

Soft Thresholding
Soft thresholding process tends to shrink noise towards zero.By taking the median values of IMFs, σ˜j and τ j were calculated using equations ( 5), (6),and (7).Where σ˜j is the estimation of the noise level of the j th IMF (scale level) and MAD j represents the absolute median deviation of the j th IMF.The soft thresholding shrinks the IMF samples by τ j towards zero as follows.
The signal can be reconstructed by adding all the IMFs which gives the de-noised signal.This procedure is applied for all the range bins.The processing steps are discussed in 4.21.

Processing steps for the Algorithm
(i) Read the time series data ( in ' .r'format).
(ii) Convert the .rfile into .matfile and read the data from .matfile.
(iii) Calculate IMFs by using EMD method for each range bin and apply Hilbert Transform on each IMF .
(iv) Apply Soft Threshold technique of de-noising as mentioned below.
 Calculate the noise level of the IMF viz.σ˜j (Donoho et al, 1995) by finding the median values of the IMFs.
 Calculate the value of universal Threshold (τ j ) by using the calculated value σ˜j .
 Taking universal Threshold as reference and the conditions proposed by Donoho and Johnstone, soft thresholding was done for each range bin.
 If the value of IMF is greater than or equal to τ j , the IMF will be subtracted by τ j .If the value of IMF is less than τ j then IMF value is made zero and if the value of IMF is less than or equal to -τ j , then the IMF will be incremented by τ j..
 This process is applied on all the range bins.
(v) Reconstruct the signal by adding all the IMFs for each range bin and apply three point running average method to each range bin.
(vi) Calculate the mean noise level for each range bin.(Hildebrand et al., 1974).
(vii) Subtract the mean noise level for each range bin and plot the stacked doppler spectrum.

Moments Calculations
Three lower order Spectral moments (zero, first and second) and SNR are calculated by using adaptive moments method (Anandan et al., 2004).These three spectral moment represents the signal strength (power), the weighted mean doppler shift and width of the spectrum (Woodman et al., 1985;Morse et al., 2002;Anandan et al., 2004).The moments were calculated for the data of 24 th July 2002 and 22 nd Jan 2007 data by using FFT and HHT.The expressions for the first three moments are as follows.
The 0 th moment representing the total signal power is -------------------------------------------( The 2 nd moment representing the variance, a measure of dispersion from the mean frequency is 148 Where m, n are the lower and upper limits of the Doppler bin of the spectral window.P i , f i are the powers and frequencies corresponding to the Doppler bins within the spectral window. Signal-to-noise ratio (SNR) in dB is calculated by equation ( 12).------------------------------(12) Where N and L are the total number of Doppler bins and mean noise level respectively which on multiplication gives the total noise over the whole bandwidth.Doppler width, which is taken to be the full width of the Doppler spectrum is calculated as:

6.Results And Discussion
The developed HHT algorithm was applied on various range bins of time series MST Radar data upto 25 Kms.

Conclusion
Thus an efficient algorithm based on Empirical Mode Decomposition de-noising using soft threshold technique for accurate doppler profile detection and improved SNR for MST Radar Signals was developed.Further, spectral moments were estimated and signal parameters such as mean doppler, signal power, noise power and SNR were calculated and the algorithm was tested on different radar data sets for its efficacy in comparison to FFT.It has been observed that there is a considerable improvement in recognition of the doppler echo and SNR.
Geosci.Instrum.Method.Data Syst.Discuss., https://doi.org/10.5194/gi-2017-9Manuscript under review for journal Geosci.Instrum.Method.Data Syst.Discussion started: 19 June 2017 c Author(s) 2017.CC BY 3.0 License.resolution on a continuous basis which contribute to study different dynamic process in the atmosphere.It is an important research tool in the investigation of prevailing winds, waves, turbulence and atmospheric stability and other phenomenon.The Indian MST radar is highly sensitive, pulse-coded, coherent VHF phased array radar operating at 53 MHz with a peak power-aperture product of 3x10 10 Wm 2 .

Figure 1 :
Figure 1:Doppler spectrum using FFT of height 10km Figure 2:Doppler spectrum using HHT of height 10km

Figure 5 :
Figure 5:Doppler spectrum using FFT of height 20.25km