Time-Scale chirp detection AND estimation:
Application to gravity waves detection
Marcela Morvidone Bruno Torrésani
LATP, CMI, Université de Provence,LATP, CMI, Université de Provence
39 rue Joliot-Curie, 39 rue Joliot-Curie,
13453 Marseille cedex 13, 13453 Marseille cedex 13,
France. France.
Abstract:
We address the problem of detection of time-shifted and rescaled signals embedded into second order stationary random noise, whose covariance operator does not necessarily commute with rescalings. We show that in a certain number of situations (i.e for certain types of noise spectral densities), optimal detection may be achieved, and the corresponding parameter estimation problem may be solved. We formulate the problem in generic algebraic terms emphasizing the importance of symmetry groupe, and illustrate our results with the problem of gravity waves detection at interferometric detectors.
- Introduction:
The gravity waves emitted by coalescing binary stars are (at first order of approximation) a good example of parametric family of signals generated from a single one through the action of a symmetry group: in that case (first order Newtonian approximation), the signals are chirps, i.e. amplitude and frequency modulated signals, with power law amplitude and frequency modulation. The amplitude and frequency modulations are in fact characterized by a pair of physical parameters: the instantaneous frequency of the signal 1 second before the coalescence, which is in fact a (power of) scale parameter; and the coalescence date. It may be shown that the parametric family of signals under consideration may in fact be seen as the family of the images of a reference (chirp) signal by the action of a two parameters group, the group of affine transformations of the real line. Moreover, such an action is unitary in an appropriate Hilbert space.
Motivated by that problem, we develop an optimal filter theory adapted to the detection of such signals in random noise, in situations where the group action does not commute with the covariance of the noise (i.e. when the operators generating the parametric family of signals are not diagonal in the Karhunen-Loeve basis), and then ambiguity functions do not necessarily yield the optimal results.
While the general situation does not admit any simple solution, we show in the present paper that when the group action and the covariance of the noise satisfy some compatibility relations, it is possible to come out with an explicit optimal solution, which takes the form of a modified ambiguity function (a modified broad band
______ambiguity function in the case of scalings and time shifts).
This is in particular the case when the symmetry group is the group of affine transformations of the real line, and when the noise has a power law spectral density. In that case, even though the covariance operator of the noise (a convolution operator) and the group action do not commute, they are close enough to commutation to allow for the existence of an optimal strategy, which takes the form of a modified wide band ambiguity function.
We illustrate this approach with the example of the detection of chirp signals embedded in a second order stationary random noise with power law spectral density, and provide numerical illustrations. We also come back to the problem of gravitational waves detection, and propose several extensions of the proposed strategy to account for the specificities of the problem.
- Algebraic presentation of the problem
Before addressing the specific problem of chirp detection, let us start by describing the general setting of the problem, i.e. give general definitions for the space of signals and noise, and the corresponding group actions. The global setting is fairly general, and may be derived in an abstract framework. The latter applies to one-dimensional signals, as well as to higher dimensional generalizations.
a)Algebraic setting
The problem may be formulated in completely abstract and algebraic terms. We consider a Hilbert space H of signals (functions of one real variable), and suppose that a symmetry group G acts continuously and unitarily on H: for all element g of G, there is a linear operator (g) on H such that for all element x of H, ||(g)x|| = ||x||, and the mapping is continuous.
Consider a template signal x0 belonging to H, and suppose we are given an observation y, with the classical decision problem
H0:y(t)= b(t) ,
H1:y(t)=A (g)x0(t)+ b(t) ,
where A is an (unknown) attenuation factor (a real or complex number), g is an (unknown) element of the group G, and b is a random noise modelled as a (generalized) random process, characterized by its covariance operator (see below). The problem is to test the null hypothesis, and under the hypothesis H1, estimate the parameters A and g.
For that purpose, one seeks a family Tg of continuous linear forms on H, of the form
,
in such a way that the signal to noise ratio
be maximal, for all g in G, and the variance of the output noisebe independent of the group element g (here, E{X} denotes the expectation of the random variable X). At this point, notice that the template has to be such that Tgb is a second order random variable, for all group element g.
This problem is a classical one when G is the group of translations and b is a second order stationary random process, and yields the matched filter. The main reason is the fact that in such situations, the covariance function of the noise defines a covariance operator which is a convolution operator, and therefore commutes with the translations. In more general situations, the problem may or may not admit a simple solution, depending on the relationship of the group G with the noise.
More precisely, let K be the linear operator defined by its matrix elements as follows: for all f, h such that the following expressions make sense,
.
K is the covariance operator of the random process b. Then, by the classical argument, we may write, for all element of H ,
,
so that, by Cauchy-Schwarz inequality,
.
The inequality becomes an equality if and only if is such that there exists a scalar c satisfying
,
in other words, for a templateof the form
.
The constant c is to be fixed in such a way that the quantitybe independent of the group element g. Such a constraint may not be imposed in the most general situation. However, the constraint may be imposed if the covariance operator K satisfies a quasi-commutation relation with the group action: in other words, if there exists a character of G (a complex-valued mapping satisfying for all ) such that for all,
.
Indeed, in such a case, one may set
.
This ensures that
,
which is the desired property (the latter quantity being independent of g.) The corresponding value for the optimal signal to noise ratio is therefore given by
.
It is worth noticing that even though the quantity is independent of the group element g, the optimal signal to noise ratio depends on the value of g. In practice, this implies that detectability will depend on g.
b)Example: affine transformations and stationary random noise with power law spectral density:
Let us consider as a particular case the case where G is the group of affine transformations of the real line. The group G is isomorphic to the half plane (a generic group element consists in a pair (b,a) where b is a time variable and a -a positive number- is a scale variable), with group multiplication
Let Hr denote the space of distributions u whose Fourier transform vanishes for negative frequencies, and which satisfy
,
where U is the Fourier transform of u, defined as
It may be shown [7] that the above-defined norm || ||r and the corresponding scalar product turn Hr into a Hilbert space. The (unitary) action of G onto Hr is given as follows: if
,
whose Fourier W transform reads
.
Let us now turn to the description of the noise. The latter is modelled as a (possibly generalized) random process, stationary in the sense that its covariance operator K acts as a multiplier in the Fourier representation. More precisely, we use the following spectral representation:
,
where S is the (possibly unbounded) spectral density, and W(d) is a white noise measure of zero mean, and second order properties
.
We refer to [6,8] for more details on the spectral theory of random processes, and for more details on the significance of the latter equation (in the case of second order processes).
Let us assume that the spectral density has a (unbounded) power law behavior, in other words that there exists a real number such that
.
If we denote by the operator defined in the Fourier domain by
,
it follows that K is nothing but the covariance operator of the noise b in the considered Hilbert space. In addition, it is easy to verify that the covariance operator indeed satisfies the quasi-commutation relations:
,
with a character given by
,
(notice that the latter is a function of the scale parameter only.) Therefore, the discussion of the previous section applies, and we may introduce the family of filters defined by
.
Since the function is a function of the scale, the optimal signal to noise ratio depends on the scale. This implies that detectability of a signal (b,a)x0 will depend on its scale a.
Summarizing, we end up with the following expression for the optimal detector :
,
with X0 the Fourier transform of x0. The optimal detector may be interpreted as a matched filter bank, the main point being that the different matched filters (corresponding to the various values of the scale variable a) being normalized in such a way that none of them is privileged. Its output may also be interpreted as a generalized cross ambiguity function between the expected signal and the observation. Indeed, it coincides with the usual ambiguity function when the noise is a white noise and r=1. For that reason, we shall also call ambiguity function the output of the optimal filter.
Remark: The spaces of signals we have considered here are in fact spaces of analytic signals, i.e. signals whose Fourier transform vanish on the negative frequencies half axis. The modifications required for processing arbitrary signals are fairly simple, at least for real-valued signals (whose Fourier transform is completely characterized by its positive-frequencies part). In fact, it is enough to first project onto the space of interest, by cancelling the negative frequencies. Or alternatively, project the templates onto that space.
- Numerical examples:
To illustrate the proposed algorithm, we first study simple examples, namely modulated Gaussians or Gaussian chirps embedded into power law stationary noise.
Figure 1: the two reference templates: modulated gaussian (up), and gaussian chirp (down).
We display in Figure 1 the two reference signals we use in our tests, and in Figure 2 the same signals, after rescaling and shifting, embedded into Gaussian stationary noise, with power law spectral density. The parameter was set to = –1/3. We ran several simulations, with different values of the scaling parameter a0 and the amplitude parameter A.
Figure 2: the two reference templates embedded into Gaussian stationary noise, with power law spectral density.
The modulus of the corresponding (complex-valued) generalized ambiguity functions are displayed in Figures 3 and 4, for different values of the amplitude parameter A. As expected, the ambiguity functions present a peak at the correct values of the time and scale parameters. The output peak to RMS signal to noise ratio was equal to PSNR = 9.52dB and PSNR = 5.96dB respectively. Figure 5 shows another example with the modulated Gaussian with higher noise level, and output PSNR equal to 5.32dB. Further numerical tests (not presented here) confirm that the variance of the output noise indeed does not depend on the time and scale parameters (b,a).
Figure 3: Output of the detector in the case of a modulated Gaussian (figure 2 left). PSNR = 9.52dB.
Figure 4: Output of the detector in the case of a Gaussian chirp (figure 2 right). PSNR = 5.96dB.
Figure 5: Output of the detector in the case of a modulated Gaussian (figure 2 left). PSNR = 5.32dB.
We also display in Figure 6 the evolution of the (empirical) peak to RMS signal to noise ratio
,
as a function of the scale variable a0. Here, b0 and a0 denote the peak values of the time and scale variables (i.e. the values for which the output signal is maximal). denotes the empirical standard deviation of the output noise
averaged over all scales s and times . More precisely, we display numerical experiments for detection of a scaled modulated Gaussian embedded into power law noise, for three different values of the amplitude parameter A, and several values of the scale parameter a0. In each case, the result (full line in Figure 6) is an average over 40 realizations of the noise. As expected, the variations of the empirical signal to noise ratio, averaged over the 40 realizations of the noise, are found to follow closely the inverse square root of (a) (displayed with dotted lines in the figure.)
Figure 6: Empirical PSNR as a function of the scale a0
- Gravity waves at first Newtonian approximation
Let us now come back to the problem that motivated this work, namely the problem of detection of gravitational waves at interferometric detectors, and discuss the necessary modifications to the scheme proposed above.
Gravity waves are waves of deformation of the metric of space, generated by violent deformations of the matter density. It may be shown that at first order Post-Newtonian approximation, the gravity waves emitted by coalescing binary stars manifest themselves by signals of the form
,
where denotes the Heaviside step function, = 1/4 and = 3/8. Here, is the coalescence date, is a global phase, and A and F are related to physical parameters. In particular, F depends on the masses of the two stars, and A (besides being a function of the masses) is inversely proportional to the distance of the stars. More details on the specific physical problem may be found in Refs. [1,3,4].
The expected signal is then a chirp, to which local amplitude A( t) and frequencies (t) = F ( t) may be associated (notice that F is the value of the frequency one second before the coalescence). It is readily seen that the signal above is in fact of the form
,
where r 2 the scale parameter a is given by
,
and the reference signal x0 reads
.
Therefore, it makes sense to adapt to that situation the techniques we have developed in the previous sections.
The noise is generally assumed to be second order stationary and Gaussian. Its spectral measure has a continuous part, with spectral density assumed to be piecewise power law (within a bounded frequency range, say 5Hz–5 KHz), and a discrete part (corresponding to the resonances of the detector), which we shall disregard here for the sake of simplicity. An example of typical spectral density of detector noise is displayed in Figure 7. If resonances are not taken into account, the noise’s spectral density is globally well approximated by a piecewise power law function, with three distinct values of the exponent , in the low, intermediate and high frequencies domains. For such a reason it makes sense to adapt the detector described in the present paper to that situation.
Figure 7: Example of spectral density for detector noise in gravitational waves interferometric detection experiment.
Since there are at least three different values of in that case, it is necessary to adapt the strategy. We propose here two possible approaches.
- The first one is based on the following remark: because of the extremely low values of the input signal to noise ratio in the expected experimental situation, the only part of the spectrum which is likely to contribute to detection is the central part. This suggests to use the corresponding value of for defining the factor (a)in the detector.
- An alternative, already proposed by several authors, would amount to give up optimality and seek a more flexible and robust detection scheme. For example, by replacing the template signal x0 in the detector by a simpler waveform, with well-defined time and frequency locali-zation. In such a case, simple approximation techniques (see e.g. [2,4]) may yield approximate formulas for the (a) term. In fact, the corresponding detector would be in such a case closer to a modified wavelet transform (in the spirit of [5]) than to a cross-ambiguity function. The time-scale localization properties of such generalized wavelet transform are of course different from those of an ambiguity function: the wavelet transform is expected to be localized near a ridge in the time-scale plane, instead of a single point. Therefore, this other approach will have to be associated with other detection algorithms.
- Conclusions and perspectives:
We have developed and studied in this paper an algebraic approach to the problem of time-scale detection for signals embedded into non-white noise. We have shown that an optimal strategy may be derived in the case where the noise satisfies some scale covariance properties. As a result, our approach yields a time-scale transform, which may be viewed as an adapted version of a cross ambiguity function.
In the present paper, we have essentially focused on the algebraic aspects, and derived an abstract algebraic approach which covers the case of shifted and rescaled signals. Also, we have not paid attention to functional aspects. The latter will be studied elsewhere. We also plan to dig more deeply into the associated decision problems.
The application of such techniques to the problem of gravitational waves detection requires some adaptations, in particular because of the structure of the detector noise, which is in first approximation a superposition of several power law components. Another important problem is the fact that the approximation of gravity wave we have used here is only a first order approximation. Even though they are fairly close to the first order approximation, higher order post Newtonian approximations no longer take the form of shifted and rescaled copies of a reference signal. Thus, the robustness of the current approach to higher order extensions is also a serious issue.
It is worth noticing also the similarities of the proposed generalized ambiguity function with wavelet transforms, or the prewhitened wavelet transforms proposed in [5]. A natural question is also to see whether there is anything to be gained by improving the time-scale localization of the templates, as suggested for example in [2]. Such an improvement requires using non-linear time-scale transforms such as reassigned wavelet transforms. Those issues will be discussed elsewhere.