Correlation, Wavelets, and Gamma Ray Burst Timing

John Doty

MIT Center for Space Research
Room 37-541
77 Massachusetts Avenue
Cambridge, Massachusetts
USA 02139

jpd@space.mit.edu

Abstract

Cross correlation is useful as a method to estimate gamma ray burst arrival times. I show that selecting the most significant wavelets from a burst light curve is an effective way to obtain a template for cross correlation. One attractive feature of this method is that it requires the exchange of very little information to accurately compare burst arrival times across interplanetary distances.

On the Format

This report is a Mathematica notebook. If you are reading a text copy, the notebook itself may be obtained at <http://space.mit.edu/home/jpd/WakoPaper.html>. Although Mathematica notation is, in some cases, a bit harder to read than standard mathematical notation, it allows me to present the algorithmic content of the report in totality. Indeed, the electronic version, together with the electronic references, may be used as data reduction code as is.

One hazard of this format is that the thread of the exposition may be lost in programming detail. I have hidden Mathematica definitions that are not central to the exposition in "closed cells". Users of the electronic version may open and inspect these at will, and I can provide a fully opened reprint on request.

General Toolkit (hidden)

Off[General::spell1]
SetOptions[ListPlot,PlotRange->All];
<<Statistics`DescriptiveStatistics`
<<Statistics`DiscreteDistributions`
<<"Wavelet1Dv1.m"

The direct cross correlation function ccf will use the shift function to shift the data past the template, truncating it to the template length.

shift[data_,delay_,length_]:=Take[data,{delay+1,delay+length}]

The Fourier cross correlation method ftccf will use pad to extend the template with zeros to match the length of the data.

pad[template_,length_]:=
	Join[template,Table[0.,{length-Length[template]}]]

The function freqsq tabulates the square of the frequency in radians/sample period associated with the result of Fourier. The equivalent bandwidth calculation needs this.

freqsq[ len_ ] := freqsq[ len ] =
	Module[
		{ w1 = N[ 2 Pi / len ]},
	Join[ 
		Table[ (w1 k)^2, { k, 0, len/2 }], 
		Table[ (w1 (k-len))^2, 
			{ k, len/2 + 1, len - 1 }]
	]
]

The function peakdisplacement finds the displacement of the global peak of a smooth function (represented as a list) relative to the first position. Quadratic interpolation allows fractional displacements to be accurately determined. The function is assumed to be cyclic and the result corresponds to the branch of minimum absolute value. Thus, if the peak is at the last list element, the result will be -1.

findpeak[ y_List ] := 
	Position[y,Max[y],1,1][[1,1]]
zmod[ n_, m_ ] := Mod[ n + m/2 ,m] - m/2
interpeak[ y_List, ix_Integer ] :=
	zmod[ix-3+
		x/.FindMinimum[
			-Interpolation[snip3[y,ix],
				InterpolationOrder->2][x],{x,1}][[2]],
    Length[y]]
snip3[ y_List, x_Integer  ] := Take[ RotateLeft[y,x-2],3]
peakdisplacement[ y_List ] := interpeak[y,findpeak[y]]

I define wavf and wavr to be the forward and reverse wavelet transforms, normalized to preserve "power". In this normalization, the sum of the squares of the coefficients is preserved by the transform. If the input coefficients are independent random variables of equal variance, the output coefficients will be independent with the same variance as the input. The actual transform algorithm is defined in "Wavelet1Dv1.m".

wnormv[level_] := 
    wnormv[level] = 
      N[{Sqrt[ 2^level ], Table[Sqrt[2^n], 
      	{n, level, 1, -1}]}]
  wavf[data_, level_] := 
  	DPWT1D[data, basis, level] wnormv[level]
  wavr[data_] := 
    IDPWT1D[data / wnormv[Length[data[[2]]]], basis]

basis =Coif4={0.011587596739,-0.02932013798,-0.04763959031,0.273021046535,
      0.574682393857,0.294867193696,-0.054085607092,-0.042026480461,
      0.016744410163,0.003967883613,-0.001289203356,-0.000509505399};

The waveletfilter function returns a time series with all wavelets of significance less than sigma removed. It also removes the constant wavelet. The result is smoothed slightly with a convolutional filter. The smoothing step is necessary to make interpolation to continuous time possible. A better choice of wavelet basis might eliminate the need for this step.

smooth[ d_List ] :=
	0.3 RotateLeft[ d, 1 ] + 0.4 d + 0.3 RotateRight[ d, 1 ]

waveletfilter[ curve_List, sigma_Real ]:= 
  smooth[ wavr[ ReplacePart[ 
      Chop[ wavf[ curve, 
      	Floor[Log[2,Length[curve]]]],sigma],{0.},1]]]

The function poisson turns a list of rates into a randomly constructed list of counts for simulat.ion.

poisson[r_]:=Random[PoissonDistribution[r]]
SetAttributes[ poisson, Listable ]

simulate[delay_Real, burstcounts_Real, backcounts_Real]:=
  poisson[ burstcounts sampburst[128,delay,256]+
  	backcounts/256]

The function stat displays statistics of a list of measurements.

stat[list_] := {
		StringForm["Mean=``", NumberForm[Mean[list],4]], 
		StringForm["Sigma=``", NumberForm[StandardDeviation[list],4]], 
    StringForm["Range=``", NumberForm[{Min[list], Max[list]},4]]}

Introduction

Cross correlation has long been used for accurately measuring radar echo delays1. In this application, there is generally a clean template available for cross correlation: in the radar case, the pulse shape is determined by design.

The problem of determining gamma ray burst times is similar to the radar problem, so this method is potentially useful. However, determining a template is more difficult in this case because a clean light curve isn't usually available. Instead, the available light curve will contain large amounts of measurement noise. The presence of high frequency noise in the template light curve drastically increases the statistical uncertainty of the final result.

To use this method for gamma ray bursts, an effective smoothing method for the template is thus needed. It must preserve significant fast variability (since fast variability is most effective at indicating time delays) while rejecting high frequency noise. Wavelet analysis offers a means to this end.

A simulated burst

Definition (hidden)

sburst={0.,0.0204949,0.195477,0.593549,1.1791,1.77286,2.14865,2.18779,1.95901,
    1.65546,1.44947,1.38613,1.38654,1.3342,1.16573,0.905159,0.634113,0.434138,
    0.346066,0.364309,0.452186,0.556853,0.620056,0.594961,0.470532,0.287305,
    0.125439,0.0630876,0.127886,0.278703,0.440236,0.567214,0.676547,0.80953,
    0.964552,1.08334,1.11386,1.07443,1.03537,1.03673,1.04134,0.979782,
    0.833909,0.670944,0.603946,0.721705,1.03215,1.43981,1.77664,1.89383,
    1.76454,1.50584,1.28577,1.19986,1.23401,1.32896,1.44989,1.58296,1.69237,
    1.72039,1.64065,1.49272,1.34655,1.23607,1.13772,1.0124,0.858742,0.720738,
    0.645765,0.635596,0.636416,0.580201,0.444858,0.277844,0.155883,0.11692,
    0.134321,0.160692,0.191322,0.275641,0.464997,0.753505,1.06848,1.31081,
    1.40552,1.33154,1.12652,0.870301,0.649371,0.51569,0.467356,0.468911,
    0.493341,0.543776,0.632369,0.741288,0.814613,0.799422,0.698929,0.581218,
    0.528649,0.572443,0.67132,0.749705,0.757195,0.698916,0.619602,0.565419,
    0.559213,0.602726,0.689915,0.810863,0.949304,1.09516,1.27468,1.56263,
    2.03842,2.70179,3.419,3.95973,4.11016,3.78557,3.07033,2.17029,1.31489,
    0.664814,0.268542,0.0799396,0.0140864,0.};

simburst[ t_Real /; ( t >= 0. && t <= 1. ) ] =
  	Interpolation[ sburst ] [ 1 + t (Length[ sburst ]-1)];
simburst[ t_Real ] := 0.

sampburst[burstsamp_Integer, delay_Real, totsamp_Integer] := 
    Table[simburst[(n - delay)/burstsamp], {n, totsamp}] /
    	burstsamp

Here's what it looks like

Plot[simburst[t],{t,-0.2,1.2}]

Correlation

"Correlation" has several different meanings. Here, I'm using the "cross correlation function", which involves displacing one function of time past another, and integrating their product.

In the digital representation, measured functions are represented by samples, so integration is replaced by summation. The cross correlation function is  just the "dot product" of one set of samples with another:

ccf[template_,data_,delay_]:=
	template.shift[data,delay,Length[template]]

The function ccf computes an element of the non-cyclic discrete cross  correlation of a template with a longer data series.

If the whole function is needed, it is often more efficient to use FFT methods:

ftccf[template_,data_]:=
  Chop[InverseFourier[
      Fourier[data] Conjugate[Fourier[pad[template,Length[data]]]]]]

The Radar Problem

In radar range determination, the shape of the signal is known (although its  lag is not), and the measurement errors in samples of the received signal may often be reasonably assumed to be of constant variance and uncorrelated (at least over a wider bandwidth than the signal). In this case, it is known1 that cross-correlating the data stream against the known pulse shape and finding the peak in the cross correlation minimizes the effect of measurement noise to first order.

The lag estimate is just the location of the peak of the cross correlation (see the definition of peakdisplacement in the toolkit):

estlag[ template_List, data_List ]:= 
	peakdisplacement[ ftccf[template,data]]

To first order, the error in the estimated lag is inversely proportional to the product of the signal to noise ratio and the "effective bandwidth" of the  measurement. For a derivation in continuous time, see reference 1.

For the following discussion, I will assume that the signal is normalized so that the measurement variance of each sample is unity.  The received signal is presumed to be the product of an amplitude (amp) and a shape. If shape is normalized to represent unity signal to noise ratio (however defined), amp is then the (amplitude) signal to noise ratio (units are "sigmas"). The most theoretically reasonable normalization is probably shape.shape==1, but the lag error calculation is insensitive to this normalization.

The function ebw computes the effective bandwidth of shape. The output units are radians per sample period if shape is dimensionless: the result is proportional to the norm of shape. The expression under the square root represents the second moment of the Fourier power density (see the definition of freqsq in the toolkit).

ebw[shape_] := 
    Module[{fs = Fourier[shape]}, 
      Chop[Sqrt[(freqsq[Length[shape]] fs) . Conjugate[fs]]]]

The function deltat estimates the RMS error in the lag estimate. For this estimate to be accurate,the autocorrelation of shape should have a  maximum that is well approximated by a quadratic in the vicinity of zero lag, and the estimated error should be small compared with the width of the quadratic maximum.

deltat[template_, amp_] := 1/(ebw[template] amp)

A Radar-like Example (hidden)

In this case, the signal is proportional to the template, and the noise is simply added in. Here, I'm normalizing the template to make the amplitude equal to the signal to noise ratio, as determined by a matched filter1. This is reasonable because the pulse shape is known ahead of time, and a matched filter is therefore possible.

template = pad[sburst, 256]/Sqrt[sburst.sburst];
snr = 10.;

To save time, I'm using uniformly distributed noise. I can do this because the correlation method minimizes the variance of the result (to first order) for any uncorrelated noise of uniform variance. The special properties of the Gaussian distribution are not required. Also, here, I'm using a lag of zero.

noisy := snr template + 
	Sqrt[12.] Table[Random[] - 0.5, {256}];

Here's a typical pulse, at a signal to noise ratio of 10, with a noiseless profile superimposed:

Show[ListPlot[noisy,DisplayFunction->Identity],
  ListPlot[snr template,{PlotJoined->True,
      DisplayFunction->Identity}],{
    DisplayFunction->$DisplayFunction}]

Since the profile has a sharp peak, one might try estimating the lag by directly searching out this peak. Here are 100 simulation results:

lags = Table[ peakdisplacement[ noisy ],{100}]- 
	peakdisplacement[template] ;
stat[lags]

Just finding the peak works very poorly here, because even though there's a dominant peak in the profile, it isn't always the dominant peak in the data.

On the other hand, the theory suggests that the correlation method should be very effective:

deltat[ template,snr]

Simulation verifies this:

lags =Table[ estlag[template,noisy],{100}];
stat[lags]

The standard deviation is close to the theoretical prediction. 

One advantage of this method is that it's rather insensitive to distortions of the template, so long as the template has peaks and valleys in approximately the correct places:

distorted = Sqrt[template];
ListPlot[ distorted]

lags =Table[ estlag[distorted,noisy],{100}];
stat[lags]

The square root of the profile overemphasizes the downward peaks, and underemphasizes the upward peaks. Nevertheless, it is still a reasonably effective template, only slightly worse than the true profile.

The Gamma Ray Burst Problem

In the gamma ray burst problem, the pulse shape is not known beforehand, and the Poisson noise in the data is correlated with the pulse. For data in which the expected count rate exceeds one per bin, Poisson correlations may be effectively removed by a "square root transformation2":

sqtrans[x_] := N[Sqrt[x] + Sqrt[1 + x]]

The function sqtrans yields a list of data points with uncorrelated noise of unit variance.

To estimate the burst shape, I take a wavelet transform4 of the measured reference light curve, zero all wavelet coefficients with a significance of less than a given threshold, and transform back to the time domain. This yields a relatively smooth template, representing the "true" burst shape without noise. This representation is very crude when the signal to noise ratio is low, but the cross correlation method is so forgiving that it is nevertheless an effective template.

If the template actually represented the true noiseless light curve, the lag errors between the template and any measured light curve could be estimated as in the radar problem. For a quick estimate, I'll assume that the template is acuurate below. In actual fact, the template is usually better correlated with the reference light curve than with an independent measurement, and there are also systematic offsets due to the transform truncation, so the errors are unsymmetrical. This has little effect, however, on the lag difference between the reference and an independent measurement, and this is what counts in the actual analysis.

Simulations

Simulate a burst of 300 expected counts against a background of 500 expected counts:

SeedRandom[0];
ref = simulate[0.,300.,500.];
ListPlot[ ref ]

Now, make a template from the wavelet transform, ignoring the constant wavelet and all wavelets of less than 4.5 sigma significance:

template = waveletfilter[sqtrans[ref],4.5];
ListPlot[ template, PlotRange->All]

Only four wavelet coefficients met the significance criterion in this case, so this template could be easily encoded in a few tens of bits. This is not a terribly high fidelity representation of the burst, but it does show the dominant peaks at the ends of a less variable plateau. Many of the details are lost, but this template is still effective. The template is displaced a bit from the reference:

reflag = estlag[ template,sqtrans[ref]]

Next, simulate an independent measurement at a true lag of 55.55, and reduce it:

lagged = simulate[55.55,300.,500.];
ListPlot[ lagged ]

lag = estlag[ template,sqtrans[lagged]];
lag-reflag-55.55

The error relative to the true lag is less than 1% of the burst duration, even at this rather poor signal to noise. This was with an arbitrarily chosen random seed  (0), and further simulation confirms that this result is typical.

I  can also estimate the error using the radar error estimate. I estimate the signal to noise ratio by taking the dot product of a normalized template with the reference, and then use the radar method. A factor of Sqrt[2.] is needed because  noise in both the reference and the lagged curve contributes.

normtemp = template/Sqrt[template.template];
estsnr = normtemp.sqtrans[ref]
Sqrt[2.] deltat[ normtemp, estsnr]

By this estimate, the burst is thus ~11 sigma, and the RMS (1 sigma) error estimate is a bit less than observed in the example above.

More Simulations (hidden in proceedings)

runsim[ bcount_ ] := Module[ {	ref,template,reflag,lagged,lag,
		lagerror,normtemp,estsnr,estrms},
	ref = simulate[0.,bcount,500.];
	template = waveletfilter[sqtrans[ref],4.5];
	reflag = estlag[ template,sqtrans[ref]];
	lagged = simulate[55.55,bcount,500.];
	lag = estlag[ template,sqtrans[lagged]];
	lagerror = lag-reflag-55.55;
	normtemp = template/Sqrt[template.template];
	estsnr = normtemp.sqtrans[ref];
	estrms = Sqrt[2.] deltat[ normtemp, estsnr];
	{lagerror,estrms}
]

Make a matrix containing the results of 25 simulations:

sims = Transpose[ Table[ runsim[300.], {25}]];

The first row is the actual errors in the lag calculation:

stat[sims[[1]]]

The second row is the corresponding error estimates:

stat[sims[[2]]]

The error estimates cover a wide range (not surprising at this low signal to noise ratio). Note, however, that the mean estimated error is not far from the observed value of sigma for the actual errors.

Increasing the burst amplitude improves the precision, as expected:

sims = Transpose[ Table[ runsim[600.], {25}]];

stat[sims[[1]]]
stat[sims[[2]]]

sims = Transpose[ Table[ runsim[1200.], {25}]];

stat[sims[[1]]]
stat[sims[[2]]]

I've also compared the wavelet method with direct correlation of the observed light curves3:

rundirect[ bcount_ ] := Module[ {template,lagged,lag,},
	template = simulate[0.,bcount,500.];
	lagged = simulate[55.55,bcount,500.];
	lag = estlag[ N[template],N[lagged]];
	lag-55.55
	]

sims = Table[ rundirect[1200.], {25}];

stat[sims]

The accuracy of the two methods is comparable, but the direct method requires the exchange of about 10 times as much information in this case.

Conclusions

Wavelet transform filtering is an effective technique for generating correlation templates for burst time lag analysis. The templates are characterized by a small number of coefficients, and may therefore be exchanged across interplanetary distances at low cost.

This method is not fully developed. I do not know what the best wavelet basis for this method might be. In the simulations I have run, this method shows very similar accuracy to direct correlation of the full profiles3, but I find it difficult to believe that this is a general feature.

Acknowledgements

I acknowledge stimulating conversations with Kevin Hurley and George Ricker. This work was supported by NASA grant NAGW4763.

Bibliography

1. Skolnik, Merrill I., Introduction to Radar Systems, 462-474 (McGraw-Hill 1962).

2. Pollard J. H., Numerical and Statistical Techniques, 243 (Cambridge University Press 1977).

3. Hurley,  K., AIP Conference Proceedings 307, 687 (1994).

4. Scholl,  James F., <http://www.mathsource.com/cgi-bin/MathSource/Applications/Engineering/Electrical/Signals/0207-144>.