Research on the Method of De-noising the Short-wave Time Signal at Low SNR

CHINESE ASTRONOMY AND ASTROPHYSICS Chinese Astronomy and Astrophysics 43 (2019) 579–589

Research on the Method of De-noising the Short-wave Time Signal at Low SNR †  XIE Liang National Time Service Centre, Chinese Academy of Sciences, Xi’an 710600

Abstract A voice enhancement algorithm based on the Empirical Mode Decomposition (EMD) and the improved spectral subtraction is proposed for the low-SNR (Signal Noise Ratio) shortwave time signal. This method is proposed to solve the problem that the shortwave time signal cannot be used for timing in complex noisy environments. The core idea of this method is to use the HilbertHuang Transform (HHT) algorithm to make the empirical mode decomposition on the noisy shortwave signal, and to select the intrinsic mode functions containing the shortwave signal information for the signal reconstruction by through the maximum correlation. Then, to make the spectral subtraction on the reconstructed signal to achieve the purpose of noise reduction. The experimental result shows that this method has a better noise reduction than the traditional methods. Key words short-wave—time service—Empirical Mode Decomposition (EMD)— spectral subtraction 1.

INTRODUCTION

The time-frequency ensuring system is the infrastructure of national important strategic resources. The all-direction whole-space three-dimensional time service network is the important measure to ensure the security and reliability of the national time-frequency system. In the 1960’s, considered the strategic necessity of this country, the first special radio time service station of our country was set up in Pucheng County, Shaanxi Province with the call sign BPM. It began its regular shortwave time service of standard time and standard †

Supported by the Project of Young Innovation Talents of National Time Service Centre, Chinese

Academy of Sciences (Y824SC1S11) Received 2019–01–21; revised version 2019–03–05  

A translation of Acta Astron. Sin. Vol. 60, No. 3, pp. 23.1–23.8, 2019 [email protected]

0275-1062/19/$-see front matter © 2019 B. V. AllScience rights reserved. c Elsevier 0275-1062/01/$-see front matter  2019 Elsevier B. V. All rights reserved. doi: 10.1016/j.chinastron.2019.11.008 PII:

580

XIE Liang / Chinese Astronomy and Astrophysics 43 (2019) 579–589

frequency of our country in 1981. The shortwave time service mainly relies on the sky-wave reflection, without the restrictions of the network hub and active repeater, the short-wave time service has the advantage of long-distance signal propagation, and its propagating channel possesses the wartime anti-destruction capability, which is unattainable by other means of time service. Especially in the area of weak satellite signals, the shortwave time service will give a full play of its unique priority. The shortwave time service signal adopts the conventional double-sideband amplitude modulation, of which the anti-interference ability is low. Meanwhile, because of the complexity and incessant variation of the ionosphere, the shortwave channel is a complicated, typical frequency-selective decay channel. Some factors, such as the multipath interference, the broadcasting stations of various forms, and the outbreak of atmospheric noise etc., may interfere significantly the BPM shortwave time signals. Hence, aiming at the characteristics of shortwave signals, to study the suitable filtering method will play a very important role for improving the reliability and accuracy of the shortwave time service. For the detection of shortwave time service signals the common-used methods are spectral subtraction, short-time energy detection, spectral entropy method and so on. For these signal processing methods it is usually assumed that the background noise is stationary, and these methods are generally very effective for stationary noises, but they are not suitable to the case of non-stationary noise. Because the distribution feature of the noise occurred in the shortwave communication is different from the Gaussian white noise, these methods often give a poor result in the shortwave noise treatment, especially, even no result may be obtained while the signal-to-noise ratio is less than 0 dB. Thus the timing accuracy of shortwave time service signals is severely affected. Aiming at the features of background noises of short-wave time service signals, this paper adopts the Hilbert-Huang Transform (HHT) algorithm to make decomposition with a basis function constructed from the noisy signal itself, and to select the intrinsic mode components in order to filter out a part of the noise not containing the shortwave signal, then the shortwave time service signal after pretreatment is obtained through a composition. The spectral subtraction is the most frequently-used one among the traditional methods to deal with the background noises of voice signals. The combination of both methods can restrain better the shortwave communication noise, so as to raise the output signal-to-noise ratio.

2.

ENCODING SCHEME OF THE BPM SHORTWAVE TIME SERVICE SIGNALS

The double-sideband amplitude modulation is adopted for the BPM shortwave time service signals of the National Time Service Centre. For the signals of the coordinated universal time (UTC), the carrier frequency signal is modulated by using 10 cycles of 1 kHz sine signal to generate the audio frequency signal of 10 ms duration, and for the universal time (UT1)

XIE Liang / Chinese Astronomy and Astrophysics 43 (2019) 579–589

581

the audio frequency signal of 100 ms duration is generated by the same way. For the integral minute signals of UTC and UT1, the carrier frequency signal is modulated using the 300 cycles of 1 KHz sine signal to generate the audio frequency signal of 300 ms duration. The starting points of the second signal and minute signal are respectively the starting points of the first cycle of their audio frequency signals. The signal format is shown in Figure 1.

UTC (10 ms)

UT1 (100 ms)

Minute (300 ms)

Fig. 1 The format of the BPM second signal and minute signal

In order to accommodate to the rapid development of time-frequency technology and to raise the level of shortwave time service, on the basis of the existing 5 MHz carrier frequency, the National Time Service Centre has additionally broadcasted the short-wave time codes, including year, month, day, hour, minute, second, UT1 correction and leap second warning. As the broadcasting program, the time information codes modulated by the 125 Hz subcarrier are inserted in the UTC time signal, and the initial edge of the time code pulse lags that of the UTC second by 20 ms[1] . The broadcasting program of the BPM time signal is shown in Figure 2. With 30 min as a period, in each period are broadcasted 20 min of UTC (in two separated intervals), 4 min of UT1, 5 min of carrier (standard frequency), 1 min of radio call sign alternatively at 2.5 MHz, 5 MHz, 10 MHz and 15 MHz, so that the lighting area of low carrier frequency can cover the quiescent region of high carrier frequency to realize a mutual complement, and to ensure the full coverage of shortwave time service signals on the whole land.

0

5

10

15

20

25 29 30

35

40

45

50

55 59 60 Time /min

UTC 00 min 00 s-09 min 59 s 15 min 00 s-24 min 59 s 30 min 00 s-39 min 59 s 45 min 00 s-54 min 59 s

Carrier 10 min 00 s14 min 59 s 40 min 00 s44 min 59 s

UT1 25 min 00 s28 min 59 s 55 min 00 s58 min 59 s

Radio call sign 29 min 00 s and 59 min 00 s start BPM morse code and female voice announcement

Fig. 2 The broadcasting program of the BPM time signal

582

XIE Liang / Chinese Astronomy and Astrophysics 43 (2019) 579–589

3.

ANALYSIS OF HILBERT-HUANG TRANSFORM (HHT) AND IMPROVED SPECTRAL SUBTRACTION ALGORITHMS

The Hilbert-Huang Transform (HHT) is an algorithm of integrating the empirical mode decomposition (EMD) with the Hilbert spectral analysis (HSA), and is a kind of adaptive method used especially to analyze the data of non-linear and non-stationary process[2] . The key part of HHT is EMD, any complex data can be decomposed into limited number of intrinsic mode functions (IMFs). This decomposition process is a smoothing treatment on the non-stationary signals, by through an iterative process, the signals of different scales are decomposed according to the sequence from high frequency to low frequency[3] . This method is adaptive. Viewing from the angle of signal processing, the EMD decomposition is analogous to an adaptive band-pass filter, while IMFs represent the result after making the adaptive filtering on the signal with the filters of different bandwidths. 3.1 The EMD (Empirical Mode Decomposition) Algorithm The process of EMD is a process of continuous selection, with the purpose to remove the overlapped waves so as to make the data waveform more symmetrical. The IMFs decomposed by the EMD must meet the two conditions as follows[4] : (1) for component signals the number of extremum points must be equal to that of zero-crossing points, or differ by one at most; (2) for any sampling points in the signal the upper and lower envelopes are defined by the maximum and minimum values, and the mean value is zero. As for the condition (1), it is similar to the definition on the narrow-band signal in the traditional stationary Gaussian process, while the condition (2) has the traditional global restriction converted into a local one, so as to reduce mainly the fluctuation of instantaneous frequency. The decomposing process is as follows. First, all the extremum points with the decomposition signal x(t) are derived, and adopting the 3rd-order spline fitting function to obtain the upper and lower envelope curves of the signal, which are averaged to obtain the mean value m1 (t) of x(t), then to have the original data set x(t) minus m1 (t) to obtain the a new data set h1 (t). Second, to judge whether h1 (t) meets the two conditions of the IMF, if yes, it is the first intrinsic mode component IMF; if no, h1 (t) is taken as a new initial signal to repeat the previous decomposition process until the IMF conditions are satisfied. Finally, the hk (t) which satisfies the conditions is taken as the first IMF component c1 (t), this is the component of x(t) with the highest frequency. To obtain the remainder signal by subtracting c1 (t) from x(t), and to repeat the previous decomposition process, the various IMFs will be obtained, as shown by Equation (1):

XIE Liang / Chinese Astronomy and Astrophysics 43 (2019) 579–589

x(t) − m1 = h1 ,

583

x(t) − c1 = r

h1 − m2 = h2 ,

r1 − c2 = r2

.. .

.. .

hk−1 − mk = hk ,

rn−1 − cn = rn

⇒ hk = c1

⇒ x(t) =

(1) n 

cj + rn ,

j=1

in which, x(t) is the signal to be treated; cj represents the IMF components obtained from the decomposition; rn indicates the remainders. The EMD algorithm decomposes a non-stationary signal into the sum of the different frequencies and tendency term. By through an iterative process, the component of highest frequency in the signal is obtained by the local feature selection of the signal, then in the same manner after the layer-by-layer selection, finally the frequency component whose signal period is longer than the signal duration, namely the minimal frequency, is obtained, this component is known as the tendency term, also the indicator of stopping selection. Because the decomposition is carried out according to the characteristic timescale of the signal in the EMD decomposition algorithm, along with the increase of the order number of the intrinsic mode component, the characteristic scale of the corresponding decomposed component increases as well. As shown by the spectrum, the signal is filtered in accordance with the sequence from high frequency to low frequency. But when the signal timescale has a stepped variation, this algorithm may lead to the problem of mode overlap, which is exhibited as that multiple IMFs contain the same frequency component. Figure 3 shows that by the EMD decomposition on the noisy shortwave BPM time signal recorded by an all-band receiver, the intrinsic mode components (IMFs) are obtained (taking the former 8 IMF components). We can find that after the EMD decomposition, for the intrinsic mode components (IMFs) of the decomposed noisy shortwave time signal, the fluctuation decreases gradually from the upper to the lower, namely, the frequency varies constantly from high frequency to low frequency, even though there exists mode overlap, the trend of frequency variation remains unchanged.

584

Amplitude Amplitude

2

0.1 0.05 0

Amplitude

1

0.04 0.02 0

Amplitude

0

0.04 0.02 0

0.02 0.01 0

Amplitude

IMF:number 1 0.5 0 −0.5

0.02 0.01 0

Amplitude

Amplitude

XIE Liang / Chinese Astronomy and Astrophysics 43 (2019) 579–589

0.02 0.01 0

IMF spectrum:number 1

0

0.5

1

Amplitude

Time /s 0.5 0 −0.5

IMF:number 2

0

1

2

Amplitude

IMF:number 3

0

1

2

0

0.5

1

Amplitude

IMF:number 4

0

1

2

0

0.5

1

Amplitude

IMF:number 5

0

2

1

0

0.5

1

Amplitude

IMF:number 6

0

1

2

0

0.5

1

1

2

Amplitude

0

0

0.5

1

0.01 0.005 0

0

0.5

1

Amplitude

Time /s IMF:number 8

0.5 0 −0.5 0

1

2

1.5 f /Hz

2

2.5 ×104

1.5 f /Hz

2

2.5 ×104

1.5 f /Hz

2

2.5 ×104

1.5 f /Hz

2

2.5 ×104

IMF spectrum:number 7 0.01 0.005 0

Amplitude

Amplitude

IMF:number 7

2.5 ×104

IMF spectrum:number 6

Time /s 0.5 0 −0.5

2

IMF spectrum:number 5

Time /s 0.5 0 −0.5

1.5 f /Hz

IMF spectrum:number 4

Time /s 0.5 0 −0.5

2.5 ×104

IMF spectrum:number 3

Time /s 0.5 0 −0.5

2

IMF spectrum:number 2

Time /s 0.5 0 −0.5

1.5 f /Hz

1.5 f /Hz

2

2.5 ×104

IMF spectrum:number 8

0

0.5

Time /s

1

1.5 f /Hz

2

2.5 ×104

Fig. 3 The EMD decomposition of the noisy shortwave time signal

3.2

The Improved Spectral Subtraction Algorithm

As a common-used method of voice signal treatment, the spectral subtraction method demonstrates itself mainly as the noise reduction in the frequency domain. Its basic idea is that assuming the noise be independent to the voice signal, to remove the power spectrum of noise from the power spectrum of the noisy signal, so as to obtain a relatively pure spectrum

XIE Liang / Chinese Astronomy and Astrophysics 43 (2019) 579–589

585

of the voice signal. This method is widely used because of its explicit physical meaning, less constraint conditions and less calculations. Assuming that the expression of the noisy voice signal in the time domain is

y(t) = d(t) + n(t) ,

(2)

in which, y(t) represents the noisy voice signal; d(t), the pure voice signal; n(t) is the additive noise. Then, the corresponding expression in the frequency domain is

Y (ω) = D(ω) + N (ω) ,

(3)

in which, Y (ω), D(ω) and N (ω) are the Fourier transforms of y(t), d(t) and n(t), respectively, and because d(t) and n(t) are not correlated with each other, hence 2

2

2

E[|Y (ω)| ] = E[|D(ω)| ] + E[|N (ω)| ] ,

(4)

2

in which E[∗] expresses the statistical mean value. E[|N (ω)| ] is obtained from the statistical mean value of the voice signal in the “quiescent duration” by means of end point detection on the voice signal, hence the estimate of voice amplitude spectrum of the original signal can be expressed as:

 |D(ω)| =

⎧ ⎨ ⎩0

2  (ω)|2 |Y (ω)| − |N

2  (ω)|2 |Y (ω)|  |N

,

(5)

other

  (ω)| indicate the estimated value of voice and that of noise, respecin which, |D(ω)| and |N tively. This traditional spectral substraction method has a drawback that in the estimation of noise, the estimated value of noise is obtained simply by the statistical average. For the noisy signal, the noise intensity fluctuates around the noise mean value, this leads to that after the spectral subtraction there is residual noise existed, this kind of residual noise is called the musical noise. In order to solve the problem of musical noise, an improvement is put forward on the traditional spectral subtraction method by adding in the over-subtraction factor α and the spectral lower-limit threshold parameter β (0 < β < 1). Let α > 1, then the filtering effect is better than the traditional spectral subtraction. By the adjustment of the β value, all the  (ω)| after the spectral subtraction are values with the amplitudes less than the value of β|N unified to be a fixed value, so as to reduce the influence of musical noise. This fixed value is also a broadband noise. The intensity of the broadband noise can be adjusted by tuning β,

586

XIE Liang / Chinese Astronomy and Astrophysics 43 (2019) 579–589

in general, the value of α is defined according to the signal-to-noise ratio of the audio frame, and the expression is:

 |X(ω)| =

4.

⎧ ⎨

 (ω)| |Y (ω)| − α|N ⎩β|N  (ω)| 2

2

 (ω)| |Y (ω)|  (α + β)|N 2

2

.

(6)

other

DENOISING ALGORITHM OF SHORTWAVE TIME SIGNALS BASED ON THE EMPIRICAL MODE DECOMPOSITION AND SPECTRAL SUBTRACTION

For the stationary sine signal its spectrum can be obtained from a Fourier transformation, while for the stationary non-sine signal its frequency can be expressed as a synthesis of several sine signals, and hereby the signal spectrum can be obtained through a Fourier transformation of the synthetic signal. The shortwave time signal is composed of the second signal and minute signal of the standard time expressed with the different durations of 1 kHz sine wave, and for the time code signal the 125 Hz sine wave is similarly adopted as subcarrier, so that the different pieces of time information can be expressed by the different durations. The pure shortwave time signal and time codes are stationary sine signals, but because of the feature of shortwave propagation, the various uncertain factors, such as the dispersion, absorption, decay, Doppler frequency drift, interference etc., caused by the ionospheric effect, will make the shortwave time signal vary randomly. These random variations result in that the initial points of time signal seconds can not be defined with the oscillographic timing method[5] . Hence, the noisy shortwave time signals with random variations are of non-stationarity. In the HHT analysis, the EMD decomposition algorithm constructs the basis function from the noisy voice signal itself, so that the inconvenience of prior selection on the basis function and the new noise caused by the mismatch of the fixed basis function with the signal can be avoided. Through the adaptive decomposition, the noisy signal is decomposed into a sum of intrinsic mode components, each intrinsic mode component makes the correlation detection with the sine wave containing the frequency information of pure shortwave time signal, and the IMFs those possess the frequency information of shortwave time signal are selected to make summation, so as to obtain the noisy shortwave time signal with the approximately Gaussian white noise. Finally, a voice enhancement is implemented by using the improved spectral subtraction method to achieve the filtering of the shortwave time signal. The principle is shown in Figure 4.

XIE Liang / Chinese Astronomy and Astrophysics 43 (2019) 579–589

587

BPM signal with noise

Sampling and framing

The decomposition of EMD

The generation of 1 kHz and 125 Hz signal

Correlation of each IMF to standard 1 kHz and 125 Hz sine signal

Highly correlated

NO

Discard

YES IMF summation and reconstruction signal

Noise reduction by spectral subtraction

End

BPM time signal after noise removed

Fig. 4 The schematic diagram for the shortwave time signal noise reduction based on the EMD and spectral subtraction

5.

EXPERIMENTAL RESULT

Because the spectral distribution feature of the shortwave communication noise is different from that of the Gaussian white noise, in order to demonstrate more truly the filtering effect, in this paper the noisy shortwave time signal is collected by using the all-band Tecsun radio receiver, and the signal-to-noise ratio of the original noisy shortwave time signal is about −8 dB. In order to show the advantage of the method of this paper, the experimental result is compared with those of other filtering and denoising algorithms. In Figure 5, from top to bottom are respectively the original noisy shortwave signal, the shortwave signal after the denoising of spectral subtraction, the result of the short-time spectral entropy detection on the shortwave time signal endpoints, the result of the shorttime energy and zero-crossing rate detection on the shortwave time signal endpoints, and the shortwave time signal after the EMD and spectral subtraction denoising. Because the duration of the shortwave time signal is 10 ms, while that of the experimental result is 3 s, the complete waveform of the shortwave time signal can not be demonstrated. In order to show more clearly the experimental result, the final panel in Figure 5 illustrates one amplified shortwave time signal after the voice enhancement by the EMD plus spectral subtraction to display intuitively the modulated waveform of the shortwave time signal after filtering. It can be seen that the amplitude of the shortwave time signal is some time small and some

588

XIE Liang / Chinese Astronomy and Astrophysics 43 (2019) 579–589

time large, this is caused by the double refraction and random fluctuation of the ionosphere and its inhomogeneous structure, leading to the fading of the amplitude of shortwave time signal.

Amplitude

Detection value

Detection value

Amplitude

Amplitude

The original waveform of speech 2 0 −2

0

0.5

1

1.5 Time /s

2

2.5

3

The voice waveform filtered by spectral subtraction

1 0 −1

0

0.5

1

1.5 2 2.5 Time /s The short time spectral entropy detection results

0

0.5

1

3

3 2 −1

1.5 Time /s

2

2.5

3

The result of energy and zero rate detection

0.4 0.2 0

0

0.5

1

1.5 2 2.5 3 Time /s The voice waveform filtered by EMD reconstruction and spectral subtraction 1

0

0.5

1

0 −1

1.5 Time /s

2

2.5

3

Fig. 5 The comparison of filtering effects for multiple filtering algorithms

It may be seen from Figure 5 that the spectral subtraction, and the short-time energy and zero-crossing rate detection on the shortwave time signal endpoints are almost invalid, when the signal-to-noise ratio of the input signal is −8 dB. The spectral entropy method has a certain effect, the signal-to-noise ratio after filtering is about 0 dB. The algorithm adopted in this paper achieves the signal-to-noise ratio of 2.3297 dB after filtering. The experimental result demonstrates that the denoising algorithm of shortwave time signals

XIE Liang / Chinese Astronomy and Astrophysics 43 (2019) 579–589

589

based on the EMD decomposition and spectral subtraction has some advantages in the shortwave time signal extraction from the shortwave communication noise in comparison with the conventional signal processing methods, such as the spectral subtraction, spectral entropy method, the short-time energy and zero crossing rate method, and so on. 6.

CONCLUSION

According to the characteristics of shortwave time signals this paper adopts a voice treatment algorithm in combination of the non-stationary signal treatment method with the stationary signal treatment method, first to make the adaptive decomposition on the noisy shortwave time signal with the EMD empirical mode decomposition algorithm, then according to the maximum similarity to select the necessary intrinsic mode components (IMFs) for the voice signal reconstruction, and finally to make filtering and to obtain the enhanced shortwave time signal by using the improved spectral subtraction algorithm. In this algorithm the input signal is collected with an ordinary shortwave radio receiver so as to bring the simulation more close to the reality. The experimental result indicates that this algorithm can substantially improve the signal-to-noise ratio of the input noisy shortwave time signal, and therefore the timing accuracy of the shortwave time signals in use. References 1

Meng Z. M., Journal of Time and Frequency, 2014, 37, 145

2

Huang N. E., Zhen S., Steven R. L., et al., RSPSA, 1998, 454, 903

3

Huang N. E., Chen X. Y., Lo M. T., et al., Advances in Adaptive Data Analysis, 2011, 3, 63

4

Zou X. J., Li X. Y., Lo M. T., et al., International Multi-Symposiums on Computer and Computational

5

Qi G. R., Foundation of Time Science, Beijing: Publishers of High Education, 2006, 75

Sciences, 2006, 1, 208