- Email: [email protected]
Wavelet Analysis Algorithm for Synthetic Pulsar Time
CHINESE ASTRONOMY AND ASTROPHYSICS
ELSEVIER
Chinese Astronomy and Astrophysics 31 (2007) 443–454
Wavelet Analysis Algorithm for Synthetic Pulsar Time† ZHONG Chong-xia 1
1,2
YANG Ting-gao1
National Time Service Center, Chinese Academy of Sciences, Xi’an 710600 2 Graduate School of the Chinese Academy of Sciences, Beijing 100049
Abstract The pulsar time defined by a single pulsar is affected by several noises and in order to weaken their effects and acquire a much more stable time scale we define a synthetic pulsar time from many single pulsar times. Synthesis of the two pulsars, PSR B1855 + 09 and PSR B1937 + 21 is implemented by two methods: the classical weighting algorithm and the wavelet decomposition algorithm. The results are compared. The classical weighting algorithm is unable to take into consideration the different degrees of stability at different frequencies while the wavelet algorithm can, and thereby get better results. Key words: astrometry—time—method: data analysis
1. INTRODUCTION Both universal time and ephemeris time are based on fundamental characteristics of the motion of celestial objects[1] . The invention of atomic clocks brought about an essential change in the traditional determination of astronomical time and people no longer depend on macroscopic celestial motions, but on microscopic physical regularities when measuring time. After the development over a few decades the technique of atomic clocks has achieved great progress. At present, the degree of stability of frequency of the atomic clock is progressing at a rate of one order of magnitude every 7 years and, at the same time, the discovery of the millisecond pulsar with a higher degree of stability in frequency has stimulated people to turn to the research on pulsar time scales. In 1982, Backer et al.[2] found the first millisecond pulsar and measured with a relative error of 10−14 of the period of rotation of the millisecond †
Supported by National Time Service Center of Chinese Academy of Sciences Received 2005–12–26; revised version 2006–02–20 A translation of Acta Astron. Sin. Vol. 48, No. 2, pp. 228–238, 2007
0275-1062/07/$-see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.chinastron.2007.10.009
444
ZHONG Chong-xia et al. / Chinese Astronomy and Astrophysics 31 (2007) 443–454
pulsar PSR B1937 + 21. This accuracy was not inferior to that of an atomic clock and the degree of stability was very good, thus providing a means of highly accurate timing by taking advantage of astronomical observation. A millisecond pulsar is a natural frequency source and the only one in the natural world that can rival the quantum frequency scale. It needs no maintenance and will not stop due to man-made causes. It is a natural clock and a frequency primary standard of which much is expected as comparison with the ground-based atomic time system. The pulsar time scale may become a hopeful time scale[1] . However, the pulsar time P T defined by a single pulsar is affected by several sources of noise, such as the atomic time error, uncertainty of planetary ephemeris, propagation through the interstellar medium, gravitational waves of the initial cosmic background and instability of the pulsar itself, etc. [3] . Except the noise of the atomic time itself, it can be considered that the other noise sources are independent of the individual pulsars. Thus, a synthetic pulsar time P Tens can be established through a weighted averaging of the pulsar times P Ti defined by different pulsars so as to weaken the influences of individual independent noise sources[4], in the same way as the synthetic atomic time is constructed to abate the effect of the noise from a single atomic clock, thereby increasing the degree of stability of the atomic time frequency standard. In this paper we propose a new algorithm based on wavelet analysis in addition to the classical weighting algorithm, in an effort to achieve a better synthetic pulsar time. 2. BASIC METHOD OF SYNTHETIC ATOMIC TIME CALCULATION Assume that there is a group of N atomic clocks and their readings at the epoch t are T1 (t), T2 (t), . . ., T (t), . . ., TN (t). The atomic time defined by this group of clocks should be the weighted average of these readings, with weights Pi , say, that is, T A(t) =
N i=1
Pi Ti (t)/
N
Pi .
(1)
i=1
Usual calculation of the time scale is made for a certain time interval I and in order to keep the long-term stability of T A(t) we should require that at the common epoch t0 of the time intervals (I − 1) and I, the phase T A(t0 ) and rate dT A(t)/dt|t0 remain unchanged, or T AI−1 (t0 ) = T AI (t0 ) , (2) dT AI−1 (t)/dt|t0 = dT AI (t)/dt|t0 where t0 is the initial epoch of the present time interval I and also the terminal epoch of the preceding time interval (I − 1). Since different weighting systems are adopted in the two time intervals (I − 1) and I, the T A(t0 ) calculated by means of Eq.(1) will not satisfy the two conditions of Eq.(2). So, Eq.(1) should be changed to T A(t) =
N i=1
Pi [Ti (t) + Ai + Di (t − t0 )]/
N
Pi ,
(3)
i=1
where Ai is the phase difference of T A(t0 )−Ti (t0 ), and Di is the rate difference dT A(t)/dt|t0 , both at the common epoch of the two time intervals t0 , caused by the difference in the
ZHONG Chong-xia et al. / Chinese Astronomy and Astrophysics 31 (2007) 443–454
445
weighting systems of the two time intervals. If Ai and Di are not taken into account, then we have the case (a) shown in Fig. 1. Here, the slope of l1 is the rate of T A(t) − Ti (t) at the end of the time interval (I − 1), and the slope of l2 is d[T A(t) − Ti (t)]/ dt| t0 , which is the rate of T A(t) − Ti (t) at the beginning of the time interval I, so we have l1 = l2 and Di = l2 − l1 . Case (b) shows what happens when Di = 0 and only Ai is taken into account. Case (c) is when both Ai and Di are considered. They are respectively a phase constant and a rate constant which ensure that T A(t) remains continuous both in phase and frequency at the epoch t0 , and they are to be applied to the whole time interval I.
Fig. 1
The representations of T A(t) − Ti (t) (a)Ai and Di both are 0, (b) Di = 0, (c) Ai and Di both = 0
3. CLASSICAL WEIGHTED ALGORITHM OF SYNTHETIC PULSAR TIME To obtain a high-quality synthetic pulsar time, attention should be paid to the following particulars: the participating pulsars should be as stable as possible, the observation of TOA should have a signal-noise-ratio as high as possible, there should be a sizeable number of pulsars to participate in the synthesis, the TOA data should be long enough, the reference atomic time standard should have a good long-term degree of stability, the time transmission between the observing station and reference clock should be as accurate as possible, and lastly, the most suitable model and theory should be adopted for the analytical model and correlative correction[5]. Discovery and observation of new millisecond pulsars or interruption of the observational programs for any reason will cause either increase or decrease in the number of pulsars participating in the synthesis and will result in discontinuity in the time series. This is similar to the discontinuities in the phase and rate of the synthetic atomic time caused by different values of the weights for the clocks in consecutive time intervals in the classical weighted algorithm of the synthetic atomic time. Therefore, in order to avoid such discontinuities in the synthetic pulsar time, a correction a should be added when the number of pulsars changes according to AT − P Tens = AT (t + δt) − [P Tens (t + δt) + a] ,
(4)
446
ZHONG Chong-xia et al. / Chinese Astronomy and Astrophysics 31 (2007) 443–454
where P Tens is the new synthetic pulsar time obtained from the calculation after the change in the number of pulsars, and a is the value of the difference between P Tens and its value calculated before the change. Unlike the algorithm of the atomic time scale, no frequency correction needs be made for the pulsar time because any systematic trend of variation in frequency has already been eliminated in P Ti . The starting point of the design of synthetic pulsar time algorithm is the acquisition of the best long-term degree of stability. Therefore, (a) the synthetic pulsar time P Tens should be the weighted average of the pulsar times P Ti ’s defined by the individual pulsars; (b) the selection of weights should be based on the long-term degree of stability of the pulsar times P Ti ’s defined by the individual pulsars; (c) Any unnecessary noise introduced by some systematic trend of the residuals should be avoided as far as possible. P Tens can be expressed by the following formula Wi · (AT − P Ti ) , (5) AT − P Tens = i
where the weight Wi is inversely proportional to the degree of instability of the pulsar. If the number of P Ti ’s participating in the synthesis is large enough, it is best to adopt the N -gonal cap method to calculate the degree of instability of P Ti [6] , and in this article the observed data of the two millisecond pulsars PSR B1855 + 09 and PSR B1937 + 21 are used to make the calculation and the degree of instability of P Ti by the method of σz (τ ). It is understood that all the weights of P Ti are normalized. 4. WAVELET DECOMPOSITION ALGORITHM OF SYNTHETIC PULSAR TIME Only one fixed weight is given to a given pulsar in the classical weighted algorithm for synthesizing the pulsar time, therefore consideration can not be given to both the longterm and short-term degrees of stability of the pulsar time. This is a defect in the classical algorithm. On the other hand, for a synthetic pulsar time algorithm based on wavelet analysis, the signal of a single pulsar clock is decomposed in the wavelet field, the components at different frequencies are extracted and then the weighted average of signals is taken by making use of the degrees of stability within different ranges of the signals characterized by the wavelet variance at the time of arrival of pulses of the pulsar. As far as time-frequency analysis is concerned, the Fourier transform is most often used. However, being universal, the Fourier transform can not provide information on the correlation of time and frequency. In order to overcome this shortcoming, Gabor introduced the “window” Fourier transform, which has the capability of localization in the time and frequency domains, with fixed sizes of the windows of the time and frequency domains. They lack adaptivity and are not suitable for the analysis of multiple scale signals and sudden changes. The wavelet transform inherits and develops the window Fourier transform. The action of a wavelet basis function is equivalent to that of a window function, and the wavelet translation is equal to the translation of window. The wavelet analysis is a rather ideal implementation of analysing time and frequency[7] .
ZHONG Chong-xia et al. / Chinese Astronomy and Astrophysics 31 (2007) 443–454
The wavelet transform of a signal is
447
[7]
∞ ¯ a,b (x)dx , f (x)Ψ
Wf (a,b) =
(6)
−∞
¯ a,b (x) = 1 Ψ x − b , and it is the result of the window function Ψ (x) under the where Ψ a |a| action of a time translation b with a scale expansion and contraction a. The signals separated by wavelet transform are divided according to time and scales, and this is the time-frequency domain analysis. The signals with different frequency components are successively separated out according to scales, there are high frequency signals in small scales and low frequency signals in large scales. The noises in pulsar time include error in the atomic time, slowing-down of the pulsar rotation itself, propagation in the interstellar medium, uncertainty of planetary ephemeris and background radiation of gravitational wave, and so forth[3] . These noises are different for different frequency components and the classical weighted algorithm can not deal with the different noises having different effects on the degrees of stability of the pulsar time at different frequencies. But this problem can be solved by the method of wavelet decomposition. Let the signal of the pulsar clock be expanded in wavelet series[8] : f (t) =
∞
< f, ϕj0 ,k > ϕj0 ,k (t) +
j0
∞
< f, φj,k > φj,k (t) .
(7)
j=−∞ k=−∞
k=−∞
From the point of view of wavelet scales, the signals are divided into two layers. The one above j0 is the extraction of fundamental characteristics and the one below j0 is the detail approximation. For the i-th signal measured, Eq.(7) can be written as i
f (x) =
∞
βji0 ,k ϕj0 ,k (x)
k=−∞
+
j0
∞
αij,k φj,k (x) .
The time-frequency energy distribution of pulsar time can be expressed as ∞ ∞ Ef = −∞ −∞
(8)
j=−∞ k=−∞
Wf (a, b)2 dadb . a2
[9]
(9)
For the binary wavelet transform and within a certain local frequency range, the energy of pulsar signals can be expressed as ⎧ 1 ⎪ ⎪ ⎪ n2 ⎪ ⎪ ⎪ (n2 − n1 )α2j,k ⎨ k=n1 . (10) Ej = 1 ⎪ ⎪ ⎪ n ⎪ 2 ⎪ 2 ⎪ (n2 − n1 )βj,k ⎩ k=n1
448
ZHONG Chong-xia et al. / Chinese Astronomy and Astrophysics 31 (2007) 443–454
The local energy distribution of signals and the variance have the same dimension, and without loss of generality, one may have the definition ⎧ 1 n2 ⎪ ⎪ ⎪ ⎪ (n −n1 )α2j,k 2 ⎨ k=n1 2 . (11) σj = 1 ⎪ n2 ⎪ ⎪ 2 ⎪ (n2 −n1 )βj,k ⎩ k=n1
Thus, σj is the multi-resolution weighting at the wavelet scale j. Generally, the weighted sum of l signals f i (x)(i = · · · , l) is W f i (x) i f (x) = , (12) Wi where Wi is the weight of signal. According to wavelet transform and its recomposition relation the above formula can also be written as ⎧ i i ⎫ ⎧ i i ⎫ σj βj0 ,k ϕj0 ,k (x) ⎬ σj αj,k φj,k (x) ⎬ j0 ∞ ⎨ ∞ ⎨ i i i i f (x) = + , (13) ⎩ ⎭ ⎩ ⎭ σj σj k=−∞
i
j=−∞ k=−∞
i
where σji indicates the multi-resolution weighting of signal f i (x) at the wavelet scale j. 5. CALCULATED RESULTS AND DISCUSSION ON THE ANALYSIS To the timing residuals obtained from the observation of the pulsars PSR B1855 + 09 and PSR B1937 + 21[10] , we first applied a sliding average over every five days. We then obtained 251 data points at 10-day intervals from MJD 46450 to MJD 48950, from the sliding averages by cubic spline interpolation. Figs. 2 and 3 show the original timing residuals and the timing residuals after the interpolation, respectively. Fig. 3 is much smoother than Fig. 2, having partially eliminated the noises by the sliding averaging. Figs. 4 and 5 are, respectively, the decomposition graph and coefficient recomposition graph obtained from a 3-scale decomposition of the interpolated timing residuals by means of the Daubechies wavelet[11] of the pulsar PSR B1855 + 09. Figs. 6 and 7 are the same for the pulsar PSRB1937+21. It can be seen from the figures that different frequencies in the signals are resolved by the wavelet decomposition and the discrepancy between the reconstructed and original signals is very small. Fig. 8 displays and compares the results of synthesizing the timing residuals of the pulsars PSR B1855 + 09 and PSR B 1937 + 21 using the classical weighting algorithm and the wavelet decomposition method. The weight used in the classical algorithm τ is approximately taken as the reciprocal of σz2 (τ ) at about the time of two years, which is nearly the square root of the total span of about 7 years, following the usual practice [12] . Fig. 8 shows that the classical weighting curve of residuals has a much greater amplitude of variation than the wavelet curve, where the residuals are generally less than 1 micro-second, with the influence of the noises much more effectively removed.
ZHONG Chong-xia et al. / Chinese Astronomy and Astrophysics 31 (2007) 443–454
Fig. 2 Timing residuals of the pulsars PSR B1855 + 09 and PSR B1937 + 21 obtained from the observation
Fig. 3 Timing residuals obtained after the observed data of the pulsars PSR B1855 + 09 and PSR B1937 + 21 are interpolated
449
450
ZHONG Chong-xia et al. / Chinese Astronomy and Astrophysics 31 (2007) 443–454
Fig. 4 Interpolated signal of the pulsar PSR B1855 + 09 and the decomposition coefficients of every layer after the Daubechies wavelet decomposition
Fig. 5 Recomposition graphs of decomposition coefficients of the layers (a3,d3, d2 and d1) of the pulsar PSR B1855 + 09 and the synthesized recomposition graph (a0)
ZHONG Chong-xia et al. / Chinese Astronomy and Astrophysics 31 (2007) 443–454
Fig. 6 Interpolated signal of the pulsar PSR B1937 + 21 and the decomposition coefficients of the layers after the Daubechies wavelet decomposition
Fig. 7 Recomposition graphs of the decomposition coefficients of the layers (a3, d3, d2 and d1) of the pulsar PSR B1937 + 21 and the synthesized recomposition graph (a0)
451
452
ZHONG Chong-xia et al. / Chinese Astronomy and Astrophysics 31 (2007) 443–454
Fig. 8 Timing residuals after the syntheses of the two pulsars with the wavelet decomposition algorithm and classical weighted algorithm, respectively
Fig. 9
A log-log plot between the degree of stability and the time for the two single pulsar times
ZHONG Chong-xia et al. / Chinese Astronomy and Astrophysics 31 (2007) 443–454
453
Fig. 10 Log-log plots between the degree of stability and time for the synthetic pulsar time obtained respectively by the classic weighting and wavelet decomposition algorithms
Fig. 9 shows, separately for the single pulsar times (after sliding averaging and interpolation) given by PSR B1855 + 09 and PSR B1937 + 21 as well as a log-log plot between the degree of stability σz (τ ) and the time τ . It can be seen that the short-term degree of stability in the PSR B1855 + 09 curve is not very good: it falls very sharply to a minimum of 10−14.15 , at τ = 1.08 × 108 s, then rises again. In comparison, the B1937 + 21 curve is much better: it falls slowly to its best value of 10−1403 , at τ = 5.4 × 107 s, and afterwards it begins to rise. Fig. 10 shows the same log-log plots (between stability and time) for the synthetic pulsar time obtained using the classical weighting and wavelet algorithms. Now, since the value of weight, τ , adopted in the synthesis by the classical weighted algorithm is the reciprocal of the degree of stability of of the two single pulsar times at about 3 year, so it is only when τ is close to this time that the result from the classical weighting method reaches its best value. However, in the classical weighting algorithm, consideration can not be given to both the short-term and long-term stability, while in the wavelet decomposition algorithm different weights in different time intervals are adopted, so giving consideration to both the short-term and long-term stability. From a comparison between Figs. 9 and 10, the conclusion is obvious that, no matter which method is used to synthesize the pulsar times, the synthetic pulsar time is always better than the pulsar time defined by a single pulsar in regard of stability. References 1
Pan L.D., Radio pulsars and time scales, Observation and research on pulsars, Urumchi: Xinjiang
2
Backer D.C., et al., Nature, 1982, 300, 615
People’s Publishing House, 2001, 168 3
Kaspi V.M., et al., Millisecond pulsar timing: recent advances, 1995
4
Yang T.G., Pan L.D., Ni G.R., et al., Research progress in metrological uses of high precision millisec-
454
ZHONG Chong-xia et al. / Chinese Astronomy and Astrophysics 31 (2007) 443–454
ond pulsar timing data, Observation and research on pulsars, Urumch: Xinjiang People’s Publishing House, 2001, 182 5
Appenzeller I., eds., Highlights of Astronomy, printed in the Netherlands, 1995, 10, 277
6
Petit G., Tavella P., et al., In: Proc. 6th European Frequency and Time Forum, 1992, 57
7
Stephane M., A Wavelet Tour of Signal Processing, translated by Yang L.H., Dai D.Q., Huang W.L.,
8
Cheng Z.X., Wavelet Analysis Algorithm and Applications, Xian: Xi’an Jiaotong University Press,
9
Ke X.Z., On Wavelet Variance, IEEE International Frequency Control Symposium, 1997, 515
et al., Beijing: Machinery Industry Press, 2002, 2 1998, 63 10
Kaspi V.M., Taylor J.H., Ryba M. F., ApJ, 1994, 428, 713
11
Hu C.H., Zhang J.B., Xia J., et al., System Analysis and Design Based on MATLAB-Wavelet Analysis,
12
Petit G., Tavella P., A&A, 1996, 308, 290
Xi’an: Xi’an Electronics University of Science and Technology Press, 2001, 112
ELSEVIER
Chinese Astronomy and Astrophysics 31 (2007) 443–454
Wavelet Analysis Algorithm for Synthetic Pulsar Time† ZHONG Chong-xia 1
1,2
YANG Ting-gao1
National Time Service Center, Chinese Academy of Sciences, Xi’an 710600 2 Graduate School of the Chinese Academy of Sciences, Beijing 100049
Abstract The pulsar time defined by a single pulsar is affected by several noises and in order to weaken their effects and acquire a much more stable time scale we define a synthetic pulsar time from many single pulsar times. Synthesis of the two pulsars, PSR B1855 + 09 and PSR B1937 + 21 is implemented by two methods: the classical weighting algorithm and the wavelet decomposition algorithm. The results are compared. The classical weighting algorithm is unable to take into consideration the different degrees of stability at different frequencies while the wavelet algorithm can, and thereby get better results. Key words: astrometry—time—method: data analysis
1. INTRODUCTION Both universal time and ephemeris time are based on fundamental characteristics of the motion of celestial objects[1] . The invention of atomic clocks brought about an essential change in the traditional determination of astronomical time and people no longer depend on macroscopic celestial motions, but on microscopic physical regularities when measuring time. After the development over a few decades the technique of atomic clocks has achieved great progress. At present, the degree of stability of frequency of the atomic clock is progressing at a rate of one order of magnitude every 7 years and, at the same time, the discovery of the millisecond pulsar with a higher degree of stability in frequency has stimulated people to turn to the research on pulsar time scales. In 1982, Backer et al.[2] found the first millisecond pulsar and measured with a relative error of 10−14 of the period of rotation of the millisecond †
Supported by National Time Service Center of Chinese Academy of Sciences Received 2005–12–26; revised version 2006–02–20 A translation of Acta Astron. Sin. Vol. 48, No. 2, pp. 228–238, 2007
0275-1062/07/$-see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.chinastron.2007.10.009
444
ZHONG Chong-xia et al. / Chinese Astronomy and Astrophysics 31 (2007) 443–454
pulsar PSR B1937 + 21. This accuracy was not inferior to that of an atomic clock and the degree of stability was very good, thus providing a means of highly accurate timing by taking advantage of astronomical observation. A millisecond pulsar is a natural frequency source and the only one in the natural world that can rival the quantum frequency scale. It needs no maintenance and will not stop due to man-made causes. It is a natural clock and a frequency primary standard of which much is expected as comparison with the ground-based atomic time system. The pulsar time scale may become a hopeful time scale[1] . However, the pulsar time P T defined by a single pulsar is affected by several sources of noise, such as the atomic time error, uncertainty of planetary ephemeris, propagation through the interstellar medium, gravitational waves of the initial cosmic background and instability of the pulsar itself, etc. [3] . Except the noise of the atomic time itself, it can be considered that the other noise sources are independent of the individual pulsars. Thus, a synthetic pulsar time P Tens can be established through a weighted averaging of the pulsar times P Ti defined by different pulsars so as to weaken the influences of individual independent noise sources[4], in the same way as the synthetic atomic time is constructed to abate the effect of the noise from a single atomic clock, thereby increasing the degree of stability of the atomic time frequency standard. In this paper we propose a new algorithm based on wavelet analysis in addition to the classical weighting algorithm, in an effort to achieve a better synthetic pulsar time. 2. BASIC METHOD OF SYNTHETIC ATOMIC TIME CALCULATION Assume that there is a group of N atomic clocks and their readings at the epoch t are T1 (t), T2 (t), . . ., T (t), . . ., TN (t). The atomic time defined by this group of clocks should be the weighted average of these readings, with weights Pi , say, that is, T A(t) =
N i=1
Pi Ti (t)/
N
Pi .
(1)
i=1
Usual calculation of the time scale is made for a certain time interval I and in order to keep the long-term stability of T A(t) we should require that at the common epoch t0 of the time intervals (I − 1) and I, the phase T A(t0 ) and rate dT A(t)/dt|t0 remain unchanged, or T AI−1 (t0 ) = T AI (t0 ) , (2) dT AI−1 (t)/dt|t0 = dT AI (t)/dt|t0 where t0 is the initial epoch of the present time interval I and also the terminal epoch of the preceding time interval (I − 1). Since different weighting systems are adopted in the two time intervals (I − 1) and I, the T A(t0 ) calculated by means of Eq.(1) will not satisfy the two conditions of Eq.(2). So, Eq.(1) should be changed to T A(t) =
N i=1
Pi [Ti (t) + Ai + Di (t − t0 )]/
N
Pi ,
(3)
i=1
where Ai is the phase difference of T A(t0 )−Ti (t0 ), and Di is the rate difference dT A(t)/dt|t0 , both at the common epoch of the two time intervals t0 , caused by the difference in the
ZHONG Chong-xia et al. / Chinese Astronomy and Astrophysics 31 (2007) 443–454
445
weighting systems of the two time intervals. If Ai and Di are not taken into account, then we have the case (a) shown in Fig. 1. Here, the slope of l1 is the rate of T A(t) − Ti (t) at the end of the time interval (I − 1), and the slope of l2 is d[T A(t) − Ti (t)]/ dt| t0 , which is the rate of T A(t) − Ti (t) at the beginning of the time interval I, so we have l1 = l2 and Di = l2 − l1 . Case (b) shows what happens when Di = 0 and only Ai is taken into account. Case (c) is when both Ai and Di are considered. They are respectively a phase constant and a rate constant which ensure that T A(t) remains continuous both in phase and frequency at the epoch t0 , and they are to be applied to the whole time interval I.
Fig. 1
The representations of T A(t) − Ti (t) (a)Ai and Di both are 0, (b) Di = 0, (c) Ai and Di both = 0
3. CLASSICAL WEIGHTED ALGORITHM OF SYNTHETIC PULSAR TIME To obtain a high-quality synthetic pulsar time, attention should be paid to the following particulars: the participating pulsars should be as stable as possible, the observation of TOA should have a signal-noise-ratio as high as possible, there should be a sizeable number of pulsars to participate in the synthesis, the TOA data should be long enough, the reference atomic time standard should have a good long-term degree of stability, the time transmission between the observing station and reference clock should be as accurate as possible, and lastly, the most suitable model and theory should be adopted for the analytical model and correlative correction[5]. Discovery and observation of new millisecond pulsars or interruption of the observational programs for any reason will cause either increase or decrease in the number of pulsars participating in the synthesis and will result in discontinuity in the time series. This is similar to the discontinuities in the phase and rate of the synthetic atomic time caused by different values of the weights for the clocks in consecutive time intervals in the classical weighted algorithm of the synthetic atomic time. Therefore, in order to avoid such discontinuities in the synthetic pulsar time, a correction a should be added when the number of pulsars changes according to AT − P Tens = AT (t + δt) − [P Tens (t + δt) + a] ,
(4)
446
ZHONG Chong-xia et al. / Chinese Astronomy and Astrophysics 31 (2007) 443–454
where P Tens is the new synthetic pulsar time obtained from the calculation after the change in the number of pulsars, and a is the value of the difference between P Tens and its value calculated before the change. Unlike the algorithm of the atomic time scale, no frequency correction needs be made for the pulsar time because any systematic trend of variation in frequency has already been eliminated in P Ti . The starting point of the design of synthetic pulsar time algorithm is the acquisition of the best long-term degree of stability. Therefore, (a) the synthetic pulsar time P Tens should be the weighted average of the pulsar times P Ti ’s defined by the individual pulsars; (b) the selection of weights should be based on the long-term degree of stability of the pulsar times P Ti ’s defined by the individual pulsars; (c) Any unnecessary noise introduced by some systematic trend of the residuals should be avoided as far as possible. P Tens can be expressed by the following formula Wi · (AT − P Ti ) , (5) AT − P Tens = i
where the weight Wi is inversely proportional to the degree of instability of the pulsar. If the number of P Ti ’s participating in the synthesis is large enough, it is best to adopt the N -gonal cap method to calculate the degree of instability of P Ti [6] , and in this article the observed data of the two millisecond pulsars PSR B1855 + 09 and PSR B1937 + 21 are used to make the calculation and the degree of instability of P Ti by the method of σz (τ ). It is understood that all the weights of P Ti are normalized. 4. WAVELET DECOMPOSITION ALGORITHM OF SYNTHETIC PULSAR TIME Only one fixed weight is given to a given pulsar in the classical weighted algorithm for synthesizing the pulsar time, therefore consideration can not be given to both the longterm and short-term degrees of stability of the pulsar time. This is a defect in the classical algorithm. On the other hand, for a synthetic pulsar time algorithm based on wavelet analysis, the signal of a single pulsar clock is decomposed in the wavelet field, the components at different frequencies are extracted and then the weighted average of signals is taken by making use of the degrees of stability within different ranges of the signals characterized by the wavelet variance at the time of arrival of pulses of the pulsar. As far as time-frequency analysis is concerned, the Fourier transform is most often used. However, being universal, the Fourier transform can not provide information on the correlation of time and frequency. In order to overcome this shortcoming, Gabor introduced the “window” Fourier transform, which has the capability of localization in the time and frequency domains, with fixed sizes of the windows of the time and frequency domains. They lack adaptivity and are not suitable for the analysis of multiple scale signals and sudden changes. The wavelet transform inherits and develops the window Fourier transform. The action of a wavelet basis function is equivalent to that of a window function, and the wavelet translation is equal to the translation of window. The wavelet analysis is a rather ideal implementation of analysing time and frequency[7] .
ZHONG Chong-xia et al. / Chinese Astronomy and Astrophysics 31 (2007) 443–454
The wavelet transform of a signal is
447
[7]
∞ ¯ a,b (x)dx , f (x)Ψ
Wf (a,b) =
(6)
−∞
¯ a,b (x) = 1 Ψ x − b , and it is the result of the window function Ψ (x) under the where Ψ a |a| action of a time translation b with a scale expansion and contraction a. The signals separated by wavelet transform are divided according to time and scales, and this is the time-frequency domain analysis. The signals with different frequency components are successively separated out according to scales, there are high frequency signals in small scales and low frequency signals in large scales. The noises in pulsar time include error in the atomic time, slowing-down of the pulsar rotation itself, propagation in the interstellar medium, uncertainty of planetary ephemeris and background radiation of gravitational wave, and so forth[3] . These noises are different for different frequency components and the classical weighted algorithm can not deal with the different noises having different effects on the degrees of stability of the pulsar time at different frequencies. But this problem can be solved by the method of wavelet decomposition. Let the signal of the pulsar clock be expanded in wavelet series[8] : f (t) =
∞
< f, ϕj0 ,k > ϕj0 ,k (t) +
j0
∞
< f, φj,k > φj,k (t) .
(7)
j=−∞ k=−∞
k=−∞
From the point of view of wavelet scales, the signals are divided into two layers. The one above j0 is the extraction of fundamental characteristics and the one below j0 is the detail approximation. For the i-th signal measured, Eq.(7) can be written as i
f (x) =
∞
βji0 ,k ϕj0 ,k (x)
k=−∞
+
j0
∞
αij,k φj,k (x) .
The time-frequency energy distribution of pulsar time can be expressed as ∞ ∞ Ef = −∞ −∞
(8)
j=−∞ k=−∞
Wf (a, b)2 dadb . a2
[9]
(9)
For the binary wavelet transform and within a certain local frequency range, the energy of pulsar signals can be expressed as ⎧ 1 ⎪ ⎪ ⎪ n2 ⎪ ⎪ ⎪ (n2 − n1 )α2j,k ⎨ k=n1 . (10) Ej = 1 ⎪ ⎪ ⎪ n ⎪ 2 ⎪ 2 ⎪ (n2 − n1 )βj,k ⎩ k=n1
448
ZHONG Chong-xia et al. / Chinese Astronomy and Astrophysics 31 (2007) 443–454
The local energy distribution of signals and the variance have the same dimension, and without loss of generality, one may have the definition ⎧ 1 n2 ⎪ ⎪ ⎪ ⎪ (n −n1 )α2j,k 2 ⎨ k=n1 2 . (11) σj = 1 ⎪ n2 ⎪ ⎪ 2 ⎪ (n2 −n1 )βj,k ⎩ k=n1
Thus, σj is the multi-resolution weighting at the wavelet scale j. Generally, the weighted sum of l signals f i (x)(i = · · · , l) is W f i (x) i f (x) = , (12) Wi where Wi is the weight of signal. According to wavelet transform and its recomposition relation the above formula can also be written as ⎧ i i ⎫ ⎧ i i ⎫ σj βj0 ,k ϕj0 ,k (x) ⎬ σj αj,k φj,k (x) ⎬ j0 ∞ ⎨ ∞ ⎨ i i i i f (x) = + , (13) ⎩ ⎭ ⎩ ⎭ σj σj k=−∞
i
j=−∞ k=−∞
i
where σji indicates the multi-resolution weighting of signal f i (x) at the wavelet scale j. 5. CALCULATED RESULTS AND DISCUSSION ON THE ANALYSIS To the timing residuals obtained from the observation of the pulsars PSR B1855 + 09 and PSR B1937 + 21[10] , we first applied a sliding average over every five days. We then obtained 251 data points at 10-day intervals from MJD 46450 to MJD 48950, from the sliding averages by cubic spline interpolation. Figs. 2 and 3 show the original timing residuals and the timing residuals after the interpolation, respectively. Fig. 3 is much smoother than Fig. 2, having partially eliminated the noises by the sliding averaging. Figs. 4 and 5 are, respectively, the decomposition graph and coefficient recomposition graph obtained from a 3-scale decomposition of the interpolated timing residuals by means of the Daubechies wavelet[11] of the pulsar PSR B1855 + 09. Figs. 6 and 7 are the same for the pulsar PSRB1937+21. It can be seen from the figures that different frequencies in the signals are resolved by the wavelet decomposition and the discrepancy between the reconstructed and original signals is very small. Fig. 8 displays and compares the results of synthesizing the timing residuals of the pulsars PSR B1855 + 09 and PSR B 1937 + 21 using the classical weighting algorithm and the wavelet decomposition method. The weight used in the classical algorithm τ is approximately taken as the reciprocal of σz2 (τ ) at about the time of two years, which is nearly the square root of the total span of about 7 years, following the usual practice [12] . Fig. 8 shows that the classical weighting curve of residuals has a much greater amplitude of variation than the wavelet curve, where the residuals are generally less than 1 micro-second, with the influence of the noises much more effectively removed.
ZHONG Chong-xia et al. / Chinese Astronomy and Astrophysics 31 (2007) 443–454
Fig. 2 Timing residuals of the pulsars PSR B1855 + 09 and PSR B1937 + 21 obtained from the observation
Fig. 3 Timing residuals obtained after the observed data of the pulsars PSR B1855 + 09 and PSR B1937 + 21 are interpolated
449
450
ZHONG Chong-xia et al. / Chinese Astronomy and Astrophysics 31 (2007) 443–454
Fig. 4 Interpolated signal of the pulsar PSR B1855 + 09 and the decomposition coefficients of every layer after the Daubechies wavelet decomposition
Fig. 5 Recomposition graphs of decomposition coefficients of the layers (a3,d3, d2 and d1) of the pulsar PSR B1855 + 09 and the synthesized recomposition graph (a0)
ZHONG Chong-xia et al. / Chinese Astronomy and Astrophysics 31 (2007) 443–454
Fig. 6 Interpolated signal of the pulsar PSR B1937 + 21 and the decomposition coefficients of the layers after the Daubechies wavelet decomposition
Fig. 7 Recomposition graphs of the decomposition coefficients of the layers (a3, d3, d2 and d1) of the pulsar PSR B1937 + 21 and the synthesized recomposition graph (a0)
451
452
ZHONG Chong-xia et al. / Chinese Astronomy and Astrophysics 31 (2007) 443–454
Fig. 8 Timing residuals after the syntheses of the two pulsars with the wavelet decomposition algorithm and classical weighted algorithm, respectively
Fig. 9
A log-log plot between the degree of stability and the time for the two single pulsar times
ZHONG Chong-xia et al. / Chinese Astronomy and Astrophysics 31 (2007) 443–454
453
Fig. 10 Log-log plots between the degree of stability and time for the synthetic pulsar time obtained respectively by the classic weighting and wavelet decomposition algorithms
Fig. 9 shows, separately for the single pulsar times (after sliding averaging and interpolation) given by PSR B1855 + 09 and PSR B1937 + 21 as well as a log-log plot between the degree of stability σz (τ ) and the time τ . It can be seen that the short-term degree of stability in the PSR B1855 + 09 curve is not very good: it falls very sharply to a minimum of 10−14.15 , at τ = 1.08 × 108 s, then rises again. In comparison, the B1937 + 21 curve is much better: it falls slowly to its best value of 10−1403 , at τ = 5.4 × 107 s, and afterwards it begins to rise. Fig. 10 shows the same log-log plots (between stability and time) for the synthetic pulsar time obtained using the classical weighting and wavelet algorithms. Now, since the value of weight, τ , adopted in the synthesis by the classical weighted algorithm is the reciprocal of the degree of stability of of the two single pulsar times at about 3 year, so it is only when τ is close to this time that the result from the classical weighting method reaches its best value. However, in the classical weighting algorithm, consideration can not be given to both the short-term and long-term stability, while in the wavelet decomposition algorithm different weights in different time intervals are adopted, so giving consideration to both the short-term and long-term stability. From a comparison between Figs. 9 and 10, the conclusion is obvious that, no matter which method is used to synthesize the pulsar times, the synthetic pulsar time is always better than the pulsar time defined by a single pulsar in regard of stability. References 1
Pan L.D., Radio pulsars and time scales, Observation and research on pulsars, Urumchi: Xinjiang
2
Backer D.C., et al., Nature, 1982, 300, 615
People’s Publishing House, 2001, 168 3
Kaspi V.M., et al., Millisecond pulsar timing: recent advances, 1995
4
Yang T.G., Pan L.D., Ni G.R., et al., Research progress in metrological uses of high precision millisec-
454
ZHONG Chong-xia et al. / Chinese Astronomy and Astrophysics 31 (2007) 443–454
ond pulsar timing data, Observation and research on pulsars, Urumch: Xinjiang People’s Publishing House, 2001, 182 5
Appenzeller I., eds., Highlights of Astronomy, printed in the Netherlands, 1995, 10, 277
6
Petit G., Tavella P., et al., In: Proc. 6th European Frequency and Time Forum, 1992, 57
7
Stephane M., A Wavelet Tour of Signal Processing, translated by Yang L.H., Dai D.Q., Huang W.L.,
8
Cheng Z.X., Wavelet Analysis Algorithm and Applications, Xian: Xi’an Jiaotong University Press,
9
Ke X.Z., On Wavelet Variance, IEEE International Frequency Control Symposium, 1997, 515
et al., Beijing: Machinery Industry Press, 2002, 2 1998, 63 10
Kaspi V.M., Taylor J.H., Ryba M. F., ApJ, 1994, 428, 713
11
Hu C.H., Zhang J.B., Xia J., et al., System Analysis and Design Based on MATLAB-Wavelet Analysis,
12
Petit G., Tavella P., A&A, 1996, 308, 290
Xi’an: Xi’an Electronics University of Science and Technology Press, 2001, 112