Header EPSearch.h

Detailed Description

A set of functions to implement the excess power search technique which was suggested in Ref. [4] and later independently invented in Ref. [3] .

The implementation here is described in detail in Ref. [2] .

This package provides a suite for functions for computing time-frequency representations of data using stacked Fourier transforms.

The first few functions simply combine functionality from the packages fft and Window in a convenient way. They are designed to streamline the task of setting up structures to prepare for taking many discrete Fourier transforms (DFTs), including windows and plans for FFTW.

A general description of the time-frequency (TF) transform provided by TFTransform is as follows. Suppose one starts with some data \(h_j\), \(0 \le j < n\) in the time domain, with sampling time \(\Delta t\), so that the data point \(h_j\) corresponds to a time \(t_j = t_\textrm{start} + j \Delta t\). Taking the standard DFT yields complex data

\begin{equation} \label{standarddft} {\tilde h}_\gamma = \sum_{j=0}^{n-1} \, e^{-2 \pi i j \gamma / n} \, h_j \end{equation}

in the Fourier domain, for \(0 \le \gamma \le [n/2]+1\). Here the data point \({\tilde h}_\gamma\) corresponds to a frequency \(f_\gamma = \gamma \Delta f\), where \(\Delta f= 1/(n \Delta t)\) is the frequency resolution.

Now suppose that we can factorize the number \(n\) of data points as

\begin{equation} n = 2 N_T N_F. \end{equation}

Then, by a time-frequency plane we shall mean a set of \(N_T N_F\) complex numbers \(H_{I\Gamma}\) with \(0 \le I < N_T\) and \(0 \le \Gamma < N_F\), obtained by an invertible linear transformation from the original data, such that the data point \(H_{I\Gamma}\) corresponds approximately to a time \(t_I = t_\textrm{start} + I {\overline{\Delta t}}\) and to a frequency \(f_\Gamma = \Gamma {\overline{\Delta f}}\). Here \(N_F\) is the number of frequency bins in the TF plane, and \(N_T\) is the number of time bins. The time resolution \({\overline {\Delta t}}\) and frequency resolution \({\overline {\Delta f}}\) are related by \({\overline {\Delta t}} \ {\overline {\Delta f}} =1\), and are given by \({\overline {\Delta t}} = 2 N_F \Delta t\) and \({\overline {\Delta f}} = N_T \Delta f\). Note that there are many other time-frequency representations of data that are not of this type; see [1] .

There are many possible choices of linear transformations from the data \(h_j\) to data \(H_{J\Gamma}\) satisfying the above properties. Here we have implemented two simple choices. The first choice consists of dividing the time-domain data \(h_j\) into \(N_T\) equal-sized chunks, each of length \(n/N_T\), and then taking the forward DFT of each chunk. Then, \(H_{J\Gamma}\) is just the \(\Gamma\)th element of the \(J\)th chunk. In terms of formulae this corresponds to

\begin{equation} \label{verticalTFP} H_{J\Sigma} = \sum_{k=0}^{2 N_F-1} \, \exp \left[ 2 \pi i k \Sigma / (2 N_F) \right] \, h_{2 N_F J + k}, \end{equation}

for \(0 \le J < N_T\) and \(0 \le \Sigma < N_F\). We call this first type of TF plane a vertical TF plane, since it corresponds to a series of vertical lines if the time axis is horizontal and the frequency axis vertical.

The second type of TF plane is obtained by first taking a DFT of all the time domain data to obtain frequency domain data, then dividing the frequency domain data into \(N_F\) equal-sized chunks, then taking the inverse DFT of each chunk. We call the resulting TF plane a horizontal TF plane. In terms of formulae the TF plane elements are

\begin{equation} \label{horizontalTFP} H_{J\Sigma} = \sum_{\gamma=0}^{N_T-1} \, \exp \left[ -2 \pi i J \gamma / N_T \right] \, {\tilde h}_{N_T \Sigma + \gamma}, \end{equation}

for \(0 \le J < N_T\) and \(0 \le \Sigma < N_F\), where \({\tilde h}_\gamma\) is given by Eq. \eqref{standarddft}.

Prototypes
double	XLALExcessPowerFilterInnerProduct (const COMPLEX16FrequencySeries *filter1, const COMPLEX16FrequencySeries *filter2, const REAL8Sequence *correlation, const REAL8FrequencySeries *psd)
	Compute the magnitude of the inner product of two arbitrary channel filters. More...
COMPLEX16FrequencySeries *	XLALCreateExcessPowerFilter (double channel_flow, double channel_width, const REAL8FrequencySeries *psd, const REAL8Sequence *correlation)
	Generate the frequency domain channel filter function. More...
int	XLALOverlappedSegmentsCommensurate (int target_length, int segment_length, int segment_shift)
	Round target_length down so that an integer number of intervals of length segment_length, each shifted by segment_shift with respect to the interval preceding it, fits into the result. More...
int	XLALEPGetTimingParameters (int window_length, int max_tile_length, double fractional_tile_stride, int *psd_length, int *psd_shift, int *window_shift, int *window_pad, int *tiling_length)
	Compute and return the timing parameters for an excess power analysis. More...
SnglBurst *	XLALEPSearch (LIGOLwXMLStream *diagnostics, const REAL8TimeSeries *tseries, REAL8Window *window, double flow, double bandwidth, double confidence_threshold, double fractional_stride, double maxTileBandwidth, double maxTileDuration)
	Generate a linked list of burst events from a time series. More...

Function Documentation

XLALExcessPowerFilterInnerProduct()

double XLALExcessPowerFilterInnerProduct	(	const COMPLEX16FrequencySeries *	filter1,
		const COMPLEX16FrequencySeries *	filter2,
		const REAL8Sequence *	correlation,
		const REAL8FrequencySeries *	psd
	)

Compute the magnitude of the inner product of two arbitrary channel filters.

Note that the sums are done over only the positive frequency components, so this function multiplies by the required factor of 2. The result is the full inner product, not the half inner product. It is safe to pass the same filter as both arguments. If the PSD is set to NULL then no PSD weighting is applied. PSD weighting is only used in reconstructing h_rss.

The return value is NaN if the input frequency series have incompatible parameters. Note that the two-point spectral correlation function does not carry enough metadata to determine if it is compatible with the filters or PSD, for example it does not carry a deltaF parameter. It is left as an excercise for the calling code to ensure the two-point spectral correlation is appropriate.

Parameters

filter1	frequency-domain filter
filter2	frequency-domain filter
correlation	two-point spectral correlation function. see XLALREAL8WindowTwoPointSpectralCorrelation().
psd	power spectral density function. see XLALREAL8AverageSpectrumWelch() and friends.

Definition at line 73 of file EPFilters.c.

XLALCreateExcessPowerFilter()

COMPLEX16FrequencySeries* XLALCreateExcessPowerFilter	(	double	channel_flow,
		double	channel_width,
		const REAL8FrequencySeries *	psd,
		const REAL8Sequence *	correlation
	)

Generate the frequency domain channel filter function.

The filter corresponds to a frequency band [channel_flow, channel_flow + channel_width]. The filter is nominally a Hann window twice the channel's width, centred on the channel's centre frequency. This makes a sum across channels equivalent to constructing a Tukey window spanning the same frequency band. This trick is one of the ingredients that allows us to accomplish a multi-resolution tiling using a single frequency channel projection (*).

The filter is normalized so that its "magnitude" as defined by the inner product function XLALExcessPowerFilterInnerProduct() is N. Then the filter is divided by the square root of the PSD frequency series prior to normalilization. This has the effect of de-emphasizing frequency bins with high noise content, and is called "over whitening".

Note: the number of samples in the window is odd, being one more than the number of frequency bins in twice the channel width. This gets the Hann windows to super-impose to form a Tukey window. (you'll have to draw yourself a picture).

(*) Really, there's no need for the "effective window" resulting from summing across channels to be something that has a name, any channel filter at all would do, but this way the code's behaviour is more easily understood — it's easy to say "the channel filter is a Tukey window of variable central width".

Parameters

channel_flow	Hz
channel_width	Hz
psd	power spectral density function. see XLALREAL8AverageSpectrumWelch() and friends.
correlation	two-point spectral correlation function. see XLALREAL8WindowTwoPointSpectralCorrelation().

Definition at line 150 of file EPFilters.c.

XLALOverlappedSegmentsCommensurate()

int XLALOverlappedSegmentsCommensurate	(	int	target_length,
		int	segment_length,
		int	segment_shift
	)

Round target_length down so that an integer number of intervals of length segment_length, each shifted by segment_shift with respect to the interval preceding it, fits into the result.

Definition at line 191 of file EPSearch.c.

XLALEPGetTimingParameters()

int XLALEPGetTimingParameters	(	int	window_length,
		int	max_tile_length,
		double	fractional_tile_shift,
		int *	psd_length,
		int *	psd_shift,
		int *	window_shift,
		int *	window_pad,
		int *	tiling_length
	)

Compute and return the timing parameters for an excess power analysis.

Pass NULL for any optional pointer to not compute and return that parameter. The return is 0 on success, negative on failure.

Parameters

window_length	Number of samples in a window used for the time-frequency plane
max_tile_length	Number of samples in the tile of longest duration
fractional_tile_shift	Number of samples by which the start of the longest tile is shifted from the start of the tile preceding it, as a fraction of its length
psd_length	(optional) User's desired number of samples to use in computing a PSD estimate. Will be replaced with actual number of samples to use in computing a PSD estimate (rounded down to be comensurate with the windowing).
psd_shift	(optional) Number of samples by which the start of a PSD is to be shifted from the start of the PSD that preceded it in order that the tiling pattern continue smoothly across the boundary.
window_shift	Number of samples by which the start of a time-frequency plane window is shifted from the window preceding it in order that the tiling pattern continue smoothly across the boundary.
window_pad	How many samples at the start and end of each window are treated as padding, and will not be covered by the tiling.
tiling_length	How many samples will be covered by the tiling.

Definition at line 234 of file EPSearch.c.

XLALEPSearch()

SnglBurst* XLALEPSearch	(	LIGOLwXMLStream *	diagnostics,
		const REAL8TimeSeries *	tseries,
		REAL8Window *	window,
		double	flow,
		double	bandwidth,
		double	confidence_threshold,
		double	fractional_stride,
		double	maxTileBandwidth,
		double	maxTileDuration
	)

Generate a linked list of burst events from a time series.

Definition at line 969 of file EPSearch.c.