Header TimeFreqFFT_h

Detailed Description

Performs real-to-complex, complex-to-real FFTs and average power spectrum estimation.

Synopsis

#include <lal/TimeFreqFFT.h>

Perform time-to-frequency and frequency-to-time fast Fourier transforms. Also provides a function to compute mean and median power spectra with user specified windowning.

Description

The routines LALTimeFreqRealFFT() and LALTimeFreqComplexFFT() transform time series \(h_j\), \(0\le j<n\), into a frequency series \(\tilde{h}_k\). For LALTimeFreqRealFFT(),

\[ \tilde{h}_k = \Delta t \times H_k \; \mbox{for } 0\le k\le\lfloor n/2\rfloor \]

The packing covers the range from dc (inclusive) to Nyquist (inclusive if \(n\) is even). For LALTimeFreqComplexFFT(),

\[ \tilde{h}_k = \Delta t \left\{ \begin{array}{ll} H_{k+\lfloor(n+1)/2\rfloor} & \mbox{for } 0\le k<\lfloor n/2\rfloor, \\ H_{k-\lfloor n/2\rfloor} & \mbox{for } \lfloor n/2\rfloor\le k<n. \\ \end{array} \right. \]

The packing covers the range from negative Nyquist (inclusive if \(n\) is even) up to (but not including) positive Nyquist. Here \(H_k\) is the DFT of \(h_j\):

\[ H_k = \sum_{j=0}^{n-1} h_j e^{-2\pi ijk/n}. \]

The units of \(\tilde{h}_k\) are equal to the units of \(h_j\) times seconds.

The routines LALFreqTimeRealFFT() and LALFreqTimeComplexFFT() perform the inverse transforms from \(\tilde{h}_k\) back to \(h_j\). This is done by shuffling the data, performing the reverse DFT, and multiplying by \(\Delta f\).

The routine LALREAL4AverageSpectrum() uses Welch's method to compute the average power spectrum of the time series stored in the input structure tSeries and return it in the output structure fSeries. A Welch PSD estimate is defined by an FFT length, overlap length, choice of window function and averaging method. These are specified in the parameter structure; the FFT length is obtained from the length of the REAL4Window in the parameters.

On entry the parameter structure params must contain a valid REAL4Window, an integer that determines the overlap as described below and a forward FFT plan for transforming data of the specified window length into the time domain. The method used to compute the average must also be set.

If the length of the window is \(N\), then the FFT length is defined to be \(N/2-1\). The input data of length \(M\) is divided into \(i\) segments which overlap by \(o\), where

\begin{equation} i = \frac{M-o}{N-o}. \end{equation}

The PSD of each segment is obtained. The Welch PSD estimate is the average of these \(i\) sub-estimates. The average is computed using the mean or median method, as specified in the parameter structure.

Note: the return PSD estimate is a one-sided power spectral density normalized as defined in the conventions document. When the averaging method is choosen to be mean and the window type Hann, the result is the same as returned by the LDAS datacondAPI psd() action for a real sequence without detrending.

Operating Instructions

const UINT4 n  = 65536;
const REAL4 dt = 1.0 / 16384.0;
static LALStatus status; compute average power spectrum
static REAL4TimeSeries         x;
static COMPLEX8FrequencySeries X;
static COMPLEX8TimeSeries      z;
static COMPLEX8FrequencySeries Z;
RealFFTPlan    *fwdRealPlan    = NULL;
RealFFTPlan    *revRealPlan    = NULL;
ComplexFFTPlan *fwdComplexPlan = NULL;
ComplexFFTPlan *revComplexPlan = NULL;
 
LALSCreateVector( &status, &x.data, n );
LALCCreateVector( &status, &X.data, n / 2 + 1 );
LALCCreateVector( &status, &z.data, n );
LALCCreateVector( &status, &Z.data, n );
LALCreateForwardRealFFTPlan( &status, &fwdRealPlan, n, 0 );
LALCreateReverseRealFFTPlan( &status, &revRealPlan, n, 0 );
LALCreateForwardComplexFFTPlan( &status, &fwdComplexPlan, n, 0 );
LALCreateReverseComplexFFTPlan( &status, &revComplexPlan, n, 0 );
 
x.f0 = 0;
x.deltaT = dt;
x.sampleUnits = lalMeterUnit;
strncpy( x.name, "x", sizeof( x.name ) );
 
z.f0 = 0;
z.deltaT = dt;
z.sampleUnits = lalVoltUnit;
strncpy( z.name, "z", sizeof( z.name ) );
 
<assign data>
 
LALTimeFreqRealFFT( &status, &X, &x, fwdRealPlan );
LALFreqTimeRealFFT( &status, &x, &X, revRealPlan );
LALTimeFreqComplexFFT( &status, &Z, &z, fwdComplexPlan );
LALFreqTimeComplexFFT( &status, &z, &Z, revComplexPlan );
 
LALDestroyRealFFTPlan( &status, &fwdRealPlan );
LALDestroyRealFFTPlan( &status, &revRealPlan );
LALDestroyComplexFFTPlan( &status, &fwdComplexPlan );
LALDestroyComplexFFTPlan( &status, &revComplexPlan );
LALCDestroyVector( &status, &Z.data );
LALCDestroyVector( &status, &z.data );
LALCDestroyVector( &status, &X.data );
LALSDestroyVector( &status, &x.data );

Notes

1. The routines do not presently work properly with heterodyned data, i.e., the original time series data should have f0 equal to zero.

Prototypes
int	XLALREAL4TimeFreqFFT (COMPLEX8FrequencySeries *freq, const REAL4TimeSeries *tser, const REAL4FFTPlan *plan)
int	XLALREAL4FreqTimeFFT (REAL4TimeSeries *tser, const COMPLEX8FrequencySeries *freq, const REAL4FFTPlan *plan)
int	XLALREAL8TimeFreqFFT (COMPLEX16FrequencySeries *freq, const REAL8TimeSeries *tser, const REAL8FFTPlan *plan)
int	XLALREAL8FreqTimeFFT (REAL8TimeSeries *tser, const COMPLEX16FrequencySeries *freq, const REAL8FFTPlan *plan)
int	XLALCOMPLEX8TimeFreqFFT (COMPLEX8FrequencySeries *freq, const COMPLEX8TimeSeries *tser, const COMPLEX8FFTPlan *plan)
int	XLALCOMPLEX8FreqTimeFFT (COMPLEX8TimeSeries *tser, const COMPLEX8FrequencySeries *freq, const COMPLEX8FFTPlan *plan)
int	XLALCOMPLEX16TimeFreqFFT (COMPLEX16FrequencySeries *freq, const COMPLEX16TimeSeries *tser, const COMPLEX16FFTPlan *plan)
int	XLALCOMPLEX16FreqTimeFFT (COMPLEX16TimeSeries *tser, const COMPLEX16FrequencySeries *freq, const COMPLEX16FFTPlan *plan)
int	XLALREAL4ModifiedPeriodogram (REAL4FrequencySeries *periodogram, const REAL4TimeSeries *tseries, const REAL4Window *window, const REAL4FFTPlan *plan)
	Compute a "modified periodogram," i.e., the power spectrum of a windowed time series. More...
int	XLALREAL8ModifiedPeriodogram (REAL8FrequencySeries *periodogram, const REAL8TimeSeries *tseries, const REAL8Window *window, const REAL8FFTPlan *plan)
	Compute a "modified periodogram," i.e., the power spectrum of a windowed time series. More...
int	XLALREAL4AverageSpectrumWelch (REAL4FrequencySeries *spectrum, const REAL4TimeSeries *tseries, UINT4 seglen, UINT4 stride, const REAL4Window *window, const REAL4FFTPlan *plan)
	Use Welch's method to compute the average power spectrum of a time series. More...
int	XLALREAL8AverageSpectrumWelch (REAL8FrequencySeries *spectrum, const REAL8TimeSeries *tseries, UINT4 seglen, UINT4 stride, const REAL8Window *window, const REAL8FFTPlan *plan)
	Use Welch's method to compute the average power spectrum of a time series. More...
REAL8	XLALMedianBias (UINT4 nn)
	compute the median bias * See arXiv: gr-qc/0509116 appendix B for details More...
REAL8	XLALLogMedianBiasGeometric (UINT4 nn)
	Geometric mean median bias factor. More...
int	XLALREAL4AverageSpectrumMedian (REAL4FrequencySeries *spectrum, const REAL4TimeSeries *tseries, UINT4 seglen, UINT4 stride, const REAL4Window *window, const REAL4FFTPlan *plan)
	Median Method: use median average rather than mean. More...
int	XLALREAL8AverageSpectrumMedian (REAL8FrequencySeries *spectrum, const REAL8TimeSeries *tseries, UINT4 seglen, UINT4 stride, const REAL8Window *window, const REAL8FFTPlan *plan)
	Median Method: use median average rather than mean. More...
int	XLALREAL4AverageSpectrumMedianMean (REAL4FrequencySeries *spectrum, const REAL4TimeSeries *tseries, UINT4 seglen, UINT4 stride, const REAL4Window *window, const REAL4FFTPlan *plan)
	Median-Mean Method: divide overlapping segments into "even" and "odd" segments; compute the bin-by-bin median of the "even" segments and the "odd" segments, and then take the bin-by-bin average of these two median averages. More...
int	XLALREAL8AverageSpectrumMedianMean (REAL8FrequencySeries *spectrum, const REAL8TimeSeries *tseries, UINT4 seglen, UINT4 stride, const REAL8Window *window, const REAL8FFTPlan *plan)
	Median-Mean Method: divide overlapping segments into "even" and "odd" segments; compute the bin-by-bin median of the "even" segments and the "odd" segments, and then take the bin-by-bin average of these two median averages. More...
int	XLALREAL4SpectrumInvertTruncate (REAL4FrequencySeries *spectrum, REAL4 lowfreq, UINT4 seglen, UINT4 trunclen, REAL4FFTPlan *fwdplan, REAL4FFTPlan *revplan)
	UNDOCUMENTED. More...
int	XLALREAL8SpectrumInvertTruncate (REAL8FrequencySeries *spectrum, REAL8 lowfreq, UINT4 seglen, UINT4 trunclen, REAL8FFTPlan *fwdplan, REAL8FFTPlan *revplan)
	UNDOCUMENTED. More...
REAL4TimeSeries *	XLALRespFilt (REAL4TimeSeries *strain, COMPLEX8FrequencySeries *transfer)
	Apply transfer function to time series. More...
REAL4TimeSeries *	XLALREAL4Convolution (REAL4TimeSeries *strain, REAL4TimeSeries *transfer)
	SHOULD Convolve two time series, but doesn't. More...
COMPLEX8FrequencySeries *	XLALWhitenCOMPLEX8FrequencySeries (COMPLEX8FrequencySeries *fseries, const REAL4FrequencySeries *psd)
	Normalize a COMPLEX8 frequency series to a REAL4 average PSD. More...
COMPLEX16FrequencySeries *	XLALWhitenCOMPLEX16FrequencySeries (COMPLEX16FrequencySeries *fseries, const REAL8FrequencySeries *psd)
	Double-precision version of XLALWhitenCOMPLEX8FrequencySeries(). More...
REAL8Sequence *	XLALREAL8WindowTwoPointSpectralCorrelation (const REAL8Window *window, const REAL8FFTPlan *plan)
	Compute the two-point spectral correlation function for a whitened frequency series from the window applied to the original time series. More...
LALPSDRegressor *	XLALPSDRegressorNew (unsigned average_samples, unsigned median_samples)
	Allocate and initialize a LALPSDRegressor object. More...
void	XLALPSDRegressorFree (LALPSDRegressor *r)
	Free all memory associated with a LALPSDRegressor object. More...
void	XLALPSDRegressorReset (LALPSDRegressor *r)
	Reset a LALPSDRegressor object to the newly-allocated state. More...
int	XLALPSDRegressorSetAverageSamples (LALPSDRegressor *r, unsigned average_samples)
	Change the average_samples of a LALPSDRegressor object. More...
unsigned	XLALPSDRegressorGetAverageSamples (const LALPSDRegressor *r)
	Return the average_samples of a LALPSDRegressor object. More...
unsigned	XLALPSDRegressorGetNSamples (const LALPSDRegressor *r)
	Return the number of frequency series samples over which the running average has been computed. More...
int	XLALPSDRegressorSetMedianSamples (LALPSDRegressor *r, unsigned median_samples)
	Change the median_samples of a LALPSDRegressor object. More...
unsigned	XLALPSDRegressorGetMedianSamples (const LALPSDRegressor *r)
	Return the median_samples of a LALPSDRegressor object. More...
int	XLALPSDRegressorAdd (LALPSDRegressor *r, const COMPLEX16FrequencySeries *sample)
	Update a LALPSDRegressor object from a frequency series object. More...
REAL8FrequencySeries *	XLALPSDRegressorGetPSD (const LALPSDRegressor *r)
	Retrieve a copy of the current PSD estimate. More...
int	XLALPSDRegressorSetPSD (LALPSDRegressor *r, const REAL8FrequencySeries *psd, unsigned weight)
	Initialize the running average of a LALPSDRegressor to the given PSD, and record it has having counted as weight samples of the running history. More...

Data Structures
struct	LALPSDRegressor
	UNDOCUMENTED. More...

Files
file	TimeFreqFFTTest.c
	Tests the routines in TimeFreqFFT.h.

Function Documentation

XLALREAL4TimeFreqFFT()

int XLALREAL4TimeFreqFFT	(	COMPLEX8FrequencySeries *	freq,
		const REAL4TimeSeries *	tser,
		const REAL4FFTPlan *	plan
	)

Definition at line 38 of file TimeFreqFFT.c.

XLALREAL4FreqTimeFFT()

int XLALREAL4FreqTimeFFT	(	REAL4TimeSeries *	tser,
		const COMPLEX8FrequencySeries *	freq,
		const REAL4FFTPlan *	plan
	)

Definition at line 74 of file TimeFreqFFT.c.

XLALREAL8TimeFreqFFT()

int XLALREAL8TimeFreqFFT	(	COMPLEX16FrequencySeries *	freq,
		const REAL8TimeSeries *	tser,
		const REAL8FFTPlan *	plan
	)

Definition at line 117 of file TimeFreqFFT.c.

XLALREAL8FreqTimeFFT()

int XLALREAL8FreqTimeFFT	(	REAL8TimeSeries *	tser,
		const COMPLEX16FrequencySeries *	freq,
		const REAL8FFTPlan *	plan
	)

Definition at line 153 of file TimeFreqFFT.c.

XLALCOMPLEX8TimeFreqFFT()

int XLALCOMPLEX8TimeFreqFFT	(	COMPLEX8FrequencySeries *	freq,
		const COMPLEX8TimeSeries *	tser,
		const COMPLEX8FFTPlan *	plan
	)

Definition at line 206 of file TimeFreqFFT.c.

XLALCOMPLEX8FreqTimeFFT()

int XLALCOMPLEX8FreqTimeFFT	(	COMPLEX8TimeSeries *	tser,
		const COMPLEX8FrequencySeries *	freq,
		const COMPLEX8FFTPlan *	plan
	)

Definition at line 261 of file TimeFreqFFT.c.

XLALCOMPLEX16TimeFreqFFT()

int XLALCOMPLEX16TimeFreqFFT	(	COMPLEX16FrequencySeries *	freq,
		const COMPLEX16TimeSeries *	tser,
		const COMPLEX16FFTPlan *	plan
	)

Definition at line 333 of file TimeFreqFFT.c.

XLALCOMPLEX16FreqTimeFFT()

int XLALCOMPLEX16FreqTimeFFT	(	COMPLEX16TimeSeries *	tser,
		const COMPLEX16FrequencySeries *	freq,
		const COMPLEX16FFTPlan *	plan
	)

Definition at line 388 of file TimeFreqFFT.c.

XLALREAL4ModifiedPeriodogram()

int XLALREAL4ModifiedPeriodogram	(	REAL4FrequencySeries *	periodogram,
		const REAL4TimeSeries *	tseries,
		const REAL4Window *	window,
		const REAL4FFTPlan *	plan
	)

Compute a "modified periodogram," i.e., the power spectrum of a windowed time series.

Definition at line 46 of file AverageSpectrum.c.

XLALREAL8ModifiedPeriodogram()

int XLALREAL8ModifiedPeriodogram	(	REAL8FrequencySeries *	periodogram,
		const REAL8TimeSeries *	tseries,
		const REAL8Window *	window,
		const REAL8FFTPlan *	plan
	)

Compute a "modified periodogram," i.e., the power spectrum of a windowed time series.

Definition at line 122 of file AverageSpectrum.c.

XLALREAL4AverageSpectrumWelch()

int XLALREAL4AverageSpectrumWelch	(	REAL4FrequencySeries *	spectrum,
		const REAL4TimeSeries *	tseries,
		UINT4	seglen,
		UINT4	stride,
		const REAL4Window *	window,
		const REAL4FFTPlan *	plan
	)

Use Welch's method to compute the average power spectrum of a time series.

See: Peter D. Welch "The Use of Fast Fourier Transform for the Estimation of Power Spectra: A Method Based on Time Averaging Over Short, Modified Periodograms" IEEE Transactions on Audio and Electroacoustics, Vol. AU-15, No. 2, June 1967.

Definition at line 203 of file AverageSpectrum.c.

XLALREAL8AverageSpectrumWelch()

int XLALREAL8AverageSpectrumWelch	(	REAL8FrequencySeries *	spectrum,
		const REAL8TimeSeries *	tseries,
		UINT4	seglen,
		UINT4	stride,
		const REAL8Window *	window,
		const REAL8FFTPlan *	plan
	)

Use Welch's method to compute the average power spectrum of a time series.

See: Peter D. Welch "The Use of Fast Fourier Transform for the Estimation of Power Spectra: A Method Based on Time Averaging Over Short, Modified Periodograms" IEEE Transactions on Audio and Electroacoustics, Vol. AU-15, No. 2, June 1967.

Definition at line 289 of file AverageSpectrum.c.

XLALMedianBias()

REAL8 XLALMedianBias

(

UINT4

nn

)

compute the median bias * See arXiv: gr-qc/0509116 appendix B for details

Definition at line 378 of file AverageSpectrum.c.

XLALLogMedianBiasGeometric()

REAL8 XLALLogMedianBiasGeometric

(

UINT4

nn

)

Geometric mean median bias factor.

For a sample of nn 2-degree of freedom \(\chi^{2}\)-distributed random variables, compute and return the natural logarithm of the ratio of the sample median to the geometric mean. See Section III.A. of LIGO-T0900462 for the derivation.

Definition at line 1423 of file AverageSpectrum.c.

XLALREAL4AverageSpectrumMedian()

int XLALREAL4AverageSpectrumMedian	(	REAL4FrequencySeries *	spectrum,
		const REAL4TimeSeries *	tseries,
		UINT4	seglen,
		UINT4	stride,
		const REAL4Window *	window,
		const REAL4FFTPlan *	plan
	)

Median Method: use median average rather than mean.

Note: this will cause a bias if the segments overlap, i.e., if the stride is less than the segment length – even though the median bias for Gaussian noise is accounted for – because the segments are not independent and their correlation is non-zero.

Definition at line 444 of file AverageSpectrum.c.

XLALREAL8AverageSpectrumMedian()

int XLALREAL8AverageSpectrumMedian	(	REAL8FrequencySeries *	spectrum,
		const REAL8TimeSeries *	tseries,
		UINT4	seglen,
		UINT4	stride,
		const REAL8Window *	window,
		const REAL8FFTPlan *	plan
	)

Median Method: use median average rather than mean.

Note: this will cause a bias if the segments overlap, i.e., if the stride is less than the segment length – even though the median bias for Gaussian noise is accounted for – because the segments are not independent and their correlation is non-zero.

Definition at line 572 of file AverageSpectrum.c.

XLALREAL4AverageSpectrumMedianMean()

int XLALREAL4AverageSpectrumMedianMean	(	REAL4FrequencySeries *	spectrum,
		const REAL4TimeSeries *	tseries,
		UINT4	seglen,
		UINT4	stride,
		const REAL4Window *	window,
		const REAL4FFTPlan *	plan
	)

Median-Mean Method: divide overlapping segments into "even" and "odd" segments; compute the bin-by-bin median of the "even" segments and the "odd" segments, and then take the bin-by-bin average of these two median averages.

Definition at line 741 of file AverageSpectrum.c.

XLALREAL8AverageSpectrumMedianMean()

int XLALREAL8AverageSpectrumMedianMean	(	REAL8FrequencySeries *	spectrum,
		const REAL8TimeSeries *	tseries,
		UINT4	seglen,
		UINT4	stride,
		const REAL8Window *	window,
		const REAL8FFTPlan *	plan
	)

Median-Mean Method: divide overlapping segments into "even" and "odd" segments; compute the bin-by-bin median of the "even" segments and the "odd" segments, and then take the bin-by-bin average of these two median averages.

Definition at line 910 of file AverageSpectrum.c.

XLALREAL4SpectrumInvertTruncate()

int XLALREAL4SpectrumInvertTruncate	(	REAL4FrequencySeries *	spectrum,
		REAL4	lowfreq,
		UINT4	seglen,
		UINT4	trunclen,
		REAL4FFTPlan *	fwdplan,
		REAL4FFTPlan *	revplan
	)

UNDOCUMENTED.

Definition at line 1073 of file AverageSpectrum.c.

XLALREAL8SpectrumInvertTruncate()

int XLALREAL8SpectrumInvertTruncate	(	REAL8FrequencySeries *	spectrum,
		REAL8	lowfreq,
		UINT4	seglen,
		UINT4	trunclen,
		REAL8FFTPlan *	fwdplan,
		REAL8FFTPlan *	revplan
	)

UNDOCUMENTED.

Definition at line 1177 of file AverageSpectrum.c.

XLALRespFilt()

REAL4TimeSeries* XLALRespFilt	(	REAL4TimeSeries *	strain,
		COMPLEX8FrequencySeries *	transfer
	)

Apply transfer function to time series.

This function returns the convolution of the time series with the frequency domain transfer function that has been supplied by the user. It zero pads the input data by a factor of two to alleviate wraparound from the FFT. This means that the transfer function must have

\[ \texttt{deltaF} = 1.0 / ( 2.0 \times \texttt{strain->data->length} \times \texttt{strain->data->deltaT} ) \]

Definition at line 55 of file Convolution.c.

XLALREAL4Convolution()

REAL4TimeSeries* XLALREAL4Convolution	(	REAL4TimeSeries *	strain,
		REAL4TimeSeries *	transfer
	)

SHOULD Convolve two time series, but doesn't.

This function does nothing yet

Definition at line 146 of file Convolution.c.

XLALWhitenCOMPLEX8FrequencySeries()

COMPLEX8FrequencySeries* XLALWhitenCOMPLEX8FrequencySeries	(	COMPLEX8FrequencySeries *	fseries,
		const REAL4FrequencySeries *	psd
	)

Normalize a COMPLEX8 frequency series to a REAL4 average PSD.

If the frequency series is the Fourier transform of (coloured) Gaussian random noise, and the PSD is of the same noise, and both have been computed according to the LAL technical specifications (LIGO-T010095-00-Z), then the output frequency series' bins will be complex Gaussian random variables with mean squares of 1 (the real and imaginary components, individually, have variances of 1/2). Fourier transforms computed by XLALREAL4ForwardFFT(), and PSDs computed by XLALREAL4AverageSpectrumMedian() and friends conform to the LAL technical specifications.

PSDs computed from high- or low-passed data can include 0s at one or the other end of the spectrum, and these would normally result in divide-by-zero errors. This routine avoids PSD divide-by-zero errors by ignoring those frequency bins, setting them to zero in the frequency series. The justification is that if the PSD is 0 then there must be nothing in that frequency bin to normalize anyway.

The input PSD is allowed to span a larger frequency band than the input frequency series, but the frequency resolutions of the two must be the same. The frequency series is modified in place. The return value is the input frequency series' address on success, or NULL on error. When an error occurs, the contents of the input frequency series are undefined.

Reviewed: f16747afde9f92ede1c25b0b537d21160b0389a3 on 2014-08-10 by J. Creighton, B. Sathyaprakash, K. Cannon.

Definition at line 1309 of file AverageSpectrum.c.

XLALWhitenCOMPLEX16FrequencySeries()

COMPLEX16FrequencySeries* XLALWhitenCOMPLEX16FrequencySeries	(	COMPLEX16FrequencySeries *	fseries,
		const REAL8FrequencySeries *	psd
	)

Double-precision version of XLALWhitenCOMPLEX8FrequencySeries().

Reviewed: f16747afde9f92ede1c25b0b537d21160b0389a3 on 2014-08-10 by J. Creighton, B. Sathyaprakash, K. Cannon.

Definition at line 1362 of file AverageSpectrum.c.

XLALREAL8WindowTwoPointSpectralCorrelation()

REAL8Sequence* XLALREAL8WindowTwoPointSpectralCorrelation	(	const REAL8Window *	window,
		const REAL8FFTPlan *	plan
	)

Compute the two-point spectral correlation function for a whitened frequency series from the window applied to the original time series.

If \(x_{j}\) is a stationary process then the components of its Fourier transform, \(X_{k}\), are independent random variables, and let their mean square be \(\langle|X_{k}|^{2}\rangle = 1\). If \(x_{j}\) is multiplied by the window function \(w_{j}\) then it is no longer stationary and the components of its Fourier transform are no longer independent. Their correlations are

\[\langle X_{k} X*_{k'}\rangle\]

and depend only on \(|k - k'|\).

Given the window function \(w_{j}\), this function computes and returns a sequence containing \(\langle X_{k} X*_{k'}\rangle\). The sequence's indices are \(|k - k'|\). A straight-forward normalization factor can be applied to convert this for use with a sequence \(x_{j}\) whose Fourier transform does not have bins with equal mean square.

The FFT plan argument must be a forward plan (time to frequency) whose length equals that of the window. If the window has length N the return value is the address of a newly-allocated sequence of length floor(N/2 + 1), or NULL on error. This function assumes the window is symmetric about its midpoint (is an even function of the sample index if the midpoint is index 0), so that its Fourier transform is real-valued.

Parameters

window	window function used to prevent leakage when measuring PSD. see XLALCreateHannREAL8Window() and friends.
plan	forward FFT plan. see XLALCreateREAL8FFTPlan().

Definition at line 1930 of file AverageSpectrum.c.

XLALPSDRegressorNew()

LALPSDRegressor* XLALPSDRegressorNew	(	unsigned	average_samples,
		unsigned	median_samples
	)

Allocate and initialize a LALPSDRegressor object.

The LALPSDRegressor object implements an algorithm for iteratively estimating a moving mean power spectral density from a sequence of frequency series objects containing Fourier-transformed snippets of a time series. The algorithm is designed to converge quickly even with high dynamic range data and to be stable against periods of non-Gaussianity ("glitches"). These things are achieved, respectively, by using a running geometric mean instead of an arithmetic mean and updating the running mean from a median of recent frequency series objects. Obtaining a correctly normalized PSD from these steps requires the application of certain correction factors that have been derived from the assumption that the underlying statistics are Gaussian. If the noise is not Gaussian then the resulting PSD estimate will be biased — it will differ from the true arithmetic mean of the spectrum by some unknown factor.

Two parameters control the behaviour of the LALPSDRegressor object: the average_samples and median_samples. When new data is supplied, the median of the most recent median_samples frequency series objects is used to update the moving average. For each update, the new PSD estimate is computed from 1 part the new data and (average_samples - 1) parts the old data. This makes the regressor implement a moving exponentially-distributed weighted average, with average_samples setting the scale of the distribution.

Changes to the PSD must persist for approximately 1/2 the number of median samples before they begin to be reflected in the running average, and so the median_samples parameter controls the robustness of the PSD estimate to glitches. Increasing the median_samples also increases the time required for the PSD to converge.

Increasing the average_samples parameter decreases the sample noise in the resulting PSD by averaging over more frequency series samples but increases the time required to converge to a new PSD should the spectrum changes.

Definition at line 1467 of file AverageSpectrum.c.

XLALPSDRegressorFree()

void XLALPSDRegressorFree

(

LALPSDRegressor *

r

)

Free all memory associated with a LALPSDRegressor object.

The object must not be used again after calling this function.

Definition at line 1521 of file AverageSpectrum.c.

XLALPSDRegressorReset()

void XLALPSDRegressorReset

(

LALPSDRegressor *

r

)

Reset a LALPSDRegressor object to the newly-allocated state.

This resets the internal frequency series parameters.

Definition at line 1500 of file AverageSpectrum.c.

XLALPSDRegressorSetAverageSamples()

int XLALPSDRegressorSetAverageSamples	(	LALPSDRegressor *	r,
		unsigned	average_samples
	)

Change the average_samples of a LALPSDRegressor object.

Definition at line 1535 of file AverageSpectrum.c.

XLALPSDRegressorGetAverageSamples()

unsigned XLALPSDRegressorGetAverageSamples

(

const LALPSDRegressor *

r

)

Return the average_samples of a LALPSDRegressor object.

Definition at line 1547 of file AverageSpectrum.c.

XLALPSDRegressorGetNSamples()

unsigned XLALPSDRegressorGetNSamples

(

const LALPSDRegressor *

r

)

Return the number of frequency series samples over which the running average has been computed.

This initially counts the number of frequency series updates the LALPSDRegressor has received, until the count reaches average_samples and then it stops increasing.

Definition at line 1558 of file AverageSpectrum.c.

XLALPSDRegressorSetMedianSamples()

int XLALPSDRegressorSetMedianSamples	(	LALPSDRegressor *	r,
		unsigned	median_samples
	)

Change the median_samples of a LALPSDRegressor object.

Any existing median history is preserved to the extent possible: if median_samples is being reduced then the most recent samples are preserved.

Definition at line 1568 of file AverageSpectrum.c.

XLALPSDRegressorGetMedianSamples()

unsigned XLALPSDRegressorGetMedianSamples

(

const LALPSDRegressor *

r

)

Return the median_samples of a LALPSDRegressor object.

Definition at line 1612 of file AverageSpectrum.c.

XLALPSDRegressorAdd()

int XLALPSDRegressorAdd	(	LALPSDRegressor *	r,
		const COMPLEX16FrequencySeries *	sample
	)

Update a LALPSDRegressor object from a frequency series object.

The frequency series is assumed to have been computed from real-valued data and to contain only non-negative frequency components. The contents of the frequency series object are digested, and this code does not take ownership of if; the calling code is responsible for freeing the frequency series object when it is done with it.

The properties of the first frequency series object added to the regressor set certain internal parameters like the frequency resolution and Nyquist frequency of the PSD. Once those parameters have been set, the only mechanism by which they can be changed is to call XLALPSDRegressorReset() and reset the regressor to the newly-allocated state.

Definition at line 1632 of file AverageSpectrum.c.

XLALPSDRegressorGetPSD()

REAL8FrequencySeries* XLALPSDRegressorGetPSD

(

const LALPSDRegressor *

r

)

Retrieve a copy of the current PSD estimate.

The return value is a newly-allocated frequency series object. The calling code is responsible for freeing it when it no longer needs it.

Definition at line 1773 of file AverageSpectrum.c.

XLALPSDRegressorSetPSD()

int XLALPSDRegressorSetPSD	(	LALPSDRegressor *	r,
		const REAL8FrequencySeries *	psd,
		unsigned	weight
	)

Initialize the running average of a LALPSDRegressor to the given PSD, and record it has having counted as weight samples of the running history.

The median history is erased, and repopulated with values derived from the PSD. The PSD data is copied, this function does not take ownership of the PSD object, the calling code is responsible for freeing it when it no longer needs it.

Definition at line 1831 of file AverageSpectrum.c.