Instantaneous Spectrum Estimation of Event-based Densities

We present a method for obtaining a time-varying spectrum that is particularly suited when the data are in event-based form. This form arises in many areas of science and engineering, and especially in astronomy, where one has photon counting detectors. The method presented consists of three procedures. First, estimating the density using the kernel method; second, highpass filtering the manifestly positive density; finally, obtaining the time-frequency distribution with a modified Welch's method. For the sake of validation, event-based data are generated from a given distribution and the proposed method is used to construct the time-frequency spectrum and is compared to the original density. The results demonstrate the effectiveness of the method.


INTRODUCTION
An important problem in many fields is the study of densities and spectra where the data are event-based. We illustrate with a specific example of the study of a signal coming from the X-ray binary system XTE-J1550. The signal is made by a sequence t i of the photon arrival times, measured by an Xray detector orbiting the Earth [1,2]. The more concentrated are the photons in a given time interval, the higher will be the intensity of the radiation (in the X-ray region) for that particular time interval. Hence the density f (t) of the sequence t i is here intended as the normalized intensity of the radiation underlying the arriving photons.
Since a fundamental property of such a density is its instantaneous spectrum [3,4], because it is believed to carry important physical information on the system that generated the events, we present here a method that enables us to perform this estimation in an effective way. We want to build an estimation algorithm that performs the following operation where P (t, ω) is the estimated instantaneous spectrum of the (normalized) density f (t) underlying the sequence of events t i . The astronomy case mentioned earlier and described in detail in [1,2] can be considered as an excellent paradigm for instantaneous spectrum estimation of event-based densities.
Here the main physical quantity is time, but the identical method can be applied to any quantity x (with associated density f (x)).
The method consists of three steps. First, a density estimation is accomplished by using the kernel method, namely sliding a window that makes a weighted count of the events in the window support. Then highpass filtering allows one to eliminate the low frequency components that are always present in a density due to its positivity, but are often of no interest. Finally, the estimation of the time-frequency spectrum is performed by using a sliding Welch's estimator to recover instantaneous spectral components even in very noisy environments. The method is tested on simulated eventbased data, and the results show its robustness and effectiveness in estimating the time-varying spectral features of the initial density f (t).

METHOD
We now describe the basic steps of the method for the instantaneous spectrum estimation of event-based densities.

Density estimation
The original data represent a sequence of events t i that is typically recorded as a sequence of the time instants when the events occurred. This sequence of events has a density f (t). In the physical example that we mentioned in the introduction, the events are the photon arrival times in the X-ray detector mounted on the satellite, and the density is proportional to the strength of the incident radiation field. If no prior information is given on the stochastic process that underlies the events, the first step is to estimate the density f (t) underlying the data.
To accomplish such estimation we use the kernel method [5]. This method consists in sliding a function h(t) called kernel on the event data t i , and then operating a summation of the events in the kernel support using the kernel itself as weighting function. The kernel must be a positive function with integral equal to one. The estimated densityf (t) at a certain value t is the weighted summation obtained when the kernel is centered at that analysis value t. The kernel method generates a continuous estimated density. The shape and analytic properties of the kernel can control the properties of the estimated density. For example, differentiability off (t) holds if the kernel is differentiable.
In [1,2] we used a Hanning kernel, whereas in this paper results are obtained using a rectangular kernel. In general, the kernel must be chosen according to the nature of the data. 1 Also, in [7] a filter design approach has been discussed for kernel choice. Some general rules that we can follow in choosing an appropriate kernel are: • If we want the estimated densityf (t) to have certain analytical properties, these properties must be imposed on the kernel h(t). • As is the case in selecting a window for generating a spectrogram, tracking of fast variations in f (t) is better accomplished by more concentrated kernels, such as rectangular kernels. • The kernel method can be written as a convolution where the kernel h(t) plays the role of a filter [7]. In this sense one can choose the kernel by noting the effect on the spectrum of the estimated density (noise reduction, resolution, etc.).

Highpass filtering
The estimated densityf (t) is a positive function, and hence its spectrum always shows strong low frequency components and so will the time-frequency spectrum. Since the spectral components near zero are often not of interest, we highpass filterf (t) to remove them. The digital filter used in the operation can be a standard one, but care must be taken to not introduce nonlinear phase distortion in the filter process. This is particularly important in the time-frequency plane, where instantaneous spectral components are connected to the derivative of the instantaneous phase, and hence phase distortion directly changes the instantaneous spectrum. A forward-backward filtering is chosen, that consists in filtering, reversing, and filtering again the signal. The output densitŷ f h (t) can be proven to have zero phase distortion [8].

Time-frequency estimation
The last step is the estimation of the instantaneous spectrum P (t, ω) from the estimated densityf h (t) obtained from the previous step. Since the density is always corrupted with noise, we have to adopt noise reduction techniques. This is accomplished by using the sliding estimator, a time-frequency estimation technique [9,10,11,12] (see also [13]). The sliding estimator is simply obtained by sliding a window on the signal centered at every analysis time, and then taking the Welch's periodogram of the windowed signal at that time. The result is the extension of the well-known periodogram to the nondeterministic case of the time-frequency spectrum estimation of a signal embedded in noise.

VALIDATION OF THE METHOD
To prove the validity of the method, we test it on event-based densities generated by numerical algorithms. Since the generation of events with a fixed density is not a common problem in signal analysis, we provide a description of the method adopted for such purpose. After that we present and discuss the results of the time-frequency estimation method applied to the generated densities. We point out that our method has been motivated primarily by timing problems arising in astronomy as discussed in the introduction [1,2,14]. The only current method (to the best of our knowledge) that allows time-frequency spectrum estimation of event-based densities is the so-called "Dynamic Power Spectrum" [14,15]. We have already discussed this method and shown comparisons with our algorithm in [1,2] and hence we will not repeat them here. Our method shows an improvement in the details of the time-varying spectrum.
Before describing how to generate events, we first choose a continuous density f (t) by fixing a signal with a given timefrequency law. In particular, we generate a signal that has a time-frequency behavior similar to the QPO (Quasi-Periodic Oscillation, see [16]) that we investigated in [1,2]. In Figure 1, we show the time-frequency spectrum of the chosen density f (t) obtained with the sliding estimator technique described in Section 2.3. The density in this example is still regularly sampled as a standard digital signal. A better time-frequency distribution could be used to represent the signal, for example a Choi-Williams distribution [17] that exhibits greater resolution, but here we use the sliding estimator because we will compare against Figure 1 all the obtained results, and so we need to compute them in a similar way.

Density generation
We now describe a method to generate a series of events, t i , that fit the chosen density f (t). These, for example, may be time of arrival events. Among the several methods available, we choose the mixing method [18], for its simplicity and reliability. The mixing method is based on the possibility of writing the desired distribution f (t) as a linear combination of distributions, that is, where the constants p 1 , p 2 , . . . , p n play the role of weights and, due to the density property of all the functions f k involved, also probabilities. We generate the sequence t i by mixing m sequences t k i that are selected in the following way: where u i is a sequence with uniform distribution between 0 and 1. This decomposition is useful if the basic densities f k (t) are easier to generate by numerical algorithms. The easiest case is, of course, when we choose the f k to be uniform distributions. By applying the mixing method one can approximate any given distribution to any desired degree of accuracy.
In particular, we choose the f k densities to be linear densities. In this case we can approximate f (t) with greater convergence speed than by using uniform distributions. Hence

Estimation results
We now present the time-frequency spectrum obtained applying the proposed method to the chosen density f (t). We consider the noise free case and then we add to the original density a white Gaussian noise with decreasing SNR, to simulate measurement noise and other stochastic processes embedded in the data. We point out that often in astronomy the only source of noise are the statistical fluctuations associated with counting a finite number of discrete events. Nevertheless, in the general framework of signal analysis, it is a basic and important approach to consider the case where additive white noise is present.

Noise-free case
In Figure 2 we represent the estimated instantaneous spectrum of f (t) obtained from a sequence t i of events generated by the mixing method. The number of events used is 2N, N being the length of the original sampled density f (t). The quality of the result is surprising, if one considers that only two events are generated for each of the N samples of the noise free density. To give an idea on how rough is the approximation of f (t) obtained using 2N, we plot in Figure 3 a comparison between the original density f (t) and a reconstructed versionf (t) obtained from the kernel method. The quality of the reconstruction increases if one generates a sequence t i made of 4N events, as Figure 4 shows. Here a direct comparison with Figure 1 shows the quality of the estimation.

SNR = 2.5 dB
In Figure 5 we represent the estimation done when SNR = 2.5 dB. The reconstruction is very good, despite only 4N events have been used to generate the sequence t i .

SNR = −7 dB
In Figure 6 we plot the estimation obtained when SNR = −7 dB and 40N events are used. The noise has now increased consistently, and 4N points are not enough to retrieve the time-frequency signature of the original signal. The increased number of events allows us to clearly understand the structure of the time-varying frequency component. Also a "bumpy" behavior is noticeable in the signal component. This is due mainly to the interaction between the high noise level and the signal itself, more than to the event-based nature of the data. The latter can be in fact excluded considering the results obtained in the noise free case, where the event-based timefrequency spectrum represents an excellent approximation of the original density, with no visible "bump."

CONCLUSION
We have presented a new method that allows the estimation of the instantaneous spectrum of a density f (t) when only a sequence of events t i (with density f (t)) is given. We considered the case of a "noisy density," that means we developed the estimation algorithm taking care of noise reduction by using a Welch's algorithm. The results presented show that the estimation algorithm is robust and effective, and that has excellent performance even when the number of events t i available is small and a high noise level is present.