35#include <gsl/gsl_rng.h>
36#include <gsl/gsl_randist.h>
40#include <lal/AVFactories.h>
41#include <lal/LALStdlib.h>
42#include <lal/LogPrintf.h>
43#include <lal/TimeSeries.h>
44#include <lal/RealFFT.h>
45#include <lal/ComplexFFT.h>
46#include <lal/SFTfileIO.h>
49#include <lal/GeneratePulsarSignal.h>
50#include <lal/ComputeFstat.h>
51#include <lal/LFTandTSutils.h>
57#define MYMAX(x,y) ( (x) > (y) ? (x) : (y) )
58#define MYMIN(x,y) ( (x) < (y) ? (x) : (y) )
59#define RELERR(x,y) ( ( (x) - (y) ) / ( 0.5 * (fabs(x) + fabs(y)) + 2) )
127 REAL8 TspanR4 = numSamplesR4 * dt_R4;
136 REAL4FFTPlan *planR4;
146 REAL8 dfLFT = 1.0 / TspanR4;
176 XLALPrintInfo(
"Comparing LFT with itself: should give 0 for all measures\n" );
181 tol.relErr_L1 = 4
e-2;
182 tol.relErr_L2 = 5
e-2;
184 tol.relErr_atMaxAbsx = 1.2e-2;
185 tol.relErr_atMaxAbsy = 1.2e-2;
187 XLALPrintInfo(
"Comparing LFT from REAL4-timeseries, to LFT from heterodyned COMPLEX8-timeseries:\n" );
213 UINT4 numSamples = 1000;
217 for (
UINT4 j = 0;
j < numSamples;
j ++ ) {
223 REAL8 safety = ( Dterms + 1.0 ) *
dt;
229 UINT4 numSamplesOut = lround( TspanOut / dtOut );
237 for (
UINT4 j = 0;
j < numSamplesOut;
j ++ ) {
238 REAL8 t_j = tStartOut +
j * dtOut;
239 times_out->
data[
j] = t_j;
248 for (
UINT4 j = 0;
j < numSamplesOut;
j ++ ) {
260 tol.relErr_L1 = 2
e-2;
261 tol.relErr_L2 = 2
e-2;
263 tol.relErr_atMaxAbsx = 2
e-2;
264 tol.relErr_atMaxAbsy = 2
e-2;
266 XLALPrintInfo(
"Comparing sinc-interpolated timeseries to exact signal timeseries:\n" );
283 REAL8 sigmaN = 0.001;
287 UINT4 numSamples = 1000;
291 UINT4 numSamples0padded = 3 * numSamples;
292 REAL8 Tspan0padded = numSamples0padded *
dt;
293 REAL8 df0padded = 1.0 / Tspan0padded;
300 for (
UINT4 j = 0;
j < numSamples;
j ++ ) {
307 memcpy( ts0padded->
data->
data,
ts->data->data, numSamples *
sizeof( ts0padded->
data->
data[0] ) );
308 memset( ts0padded->
data->
data + numSamples, 0, ( numSamples0padded - numSamples ) *
sizeof( ts0padded->
data->
data[0] ) );
311 REAL4FFTPlan *plan, *plan0padded;
337 tmp.deltaF = df0padded;
338 tmp.data = fft0padded;
345 REAL8 safetyBins = ( Dterms + 1.0 );
348 fMin = round(
fMin / df0padded ) * df0padded;
349 UINT4 numBinsOut = numBins0padded - 3 * safetyBins;
350 REAL8 BandOut = ( numBinsOut - 1 ) * df0padded;
368 tol.relErr_L1 = 8
e-2;
369 tol.relErr_L2 = 8
e-2;
371 tol.relErr_atMaxAbsx = 4
e-3;
372 tol.relErr_atMaxAbsy = 4
e-3;
374 XLALPrintInfo(
"Comparing Dirichlet SFT interpolation with upsampled SFT:\n" );
407 REAL8 numR4SamplesPerSFT = 2 * 5000;
408 REAL8 dtR4 = (
Tsft / numR4SamplesPerSFT );
418 memset( TSdata, 0, outTS->
data->
length *
sizeof( *TSdata ) );
433 for (
UINT4 j = 0;
j < numR4SamplesPerSFT;
j++ ) {
437 TSdata[alpha_j] = crealf(
testSignal( ti, sigmaN ) );
444 sftParams.Tsft =
Tsft;
445 sftParams.timestamps = timestampsSFT;
446 sftParams.noiseSFTs = NULL;
447 sftParams.window = NULL;
476 XLAL_CHECK( err_f0 < tolFreq,
XLAL_ETOL,
"f0_1 = %.16g, f0_2 = %.16g, relErr_f0 = %g (tol = %g)\n", sft1->
f0, sft2->
f0, err_f0, tolFreq );
495 static gsl_rng *
r = NULL;
499 gsl_rng_set(
r, 13 );
509 y += gsl_ran_gaussian(
r, sigmaN );
510 y += I * gsl_ran_gaussian(
r, sigmaN );
584 UINT4 firstBinSFT = round( fMinSFT /
df );
585 UINT4 lastBinSFT = firstBinSFT + ( numBinsSFT - 1 );
588 UINT4 firstBinExt, numBinsExt;
590 UINT4 lastBinExt = firstBinExt + ( numBinsExt - 1 );
593 "Requested frequency-bins [%f,%f]Hz = [%d, %d] not contained within SFT's [%f, %f]Hz = [%d,%d].\n",
594 fMin,
fMin +
Band, firstBinExt, lastBinExt, fMinSFT, fMinSFT + ( numBinsSFT - 1 ) *
df, firstBinSFT, lastBinSFT );
596 INT4 firstBinOffset = firstBinExt - firstBinSFT;
598 if ( ( *outSFT ) == NULL ) {
601 if ( ( *outSFT )->data == NULL ) {
604 if ( ( *outSFT )->data->length != numBinsExt ) {
606 ( *outSFT )->data->length = numBinsExt;
610 ( *( *outSFT ) ) = ( *inSFT );
611 ( *outSFT )->data =
ptr;
612 ( *outSFT )->f0 = firstBinExt *
df ;
615 memcpy( ( *outSFT )->data->data, inSFT->
data->
data + firstBinOffset, numBinsExt *
sizeof( ( *outSFT )->data->data[0] ) );
void LALCheckMemoryLeaks(void)
int XLALCompareSFTs(const SFTtype *sft1, const SFTtype *sft2, const VectorComparison *tol)
int test_XLALSincInterpolateSFT(void)
static LALUnit emptyLALUnit
static COMPLEX8 testSignal(REAL8 t, REAL8 sigmaN)
SFTVector * extractBandFromSFTVector(const SFTVector *inSFTs, REAL8 fMin, REAL8 Band)
int write_SFTdata(const char *fname, const SFTtype *sft)
int test_XLALSFTVectorToLFT(void)
Unit-Test for function XLALSFTVectorToLFT().
int main(void)
MAIN function.
int test_XLALSincInterpolateCOMPLEX8TimeSeries(void)
int XLALgenerateRandomData(REAL4TimeSeries **ts, SFTVector **sfts)
function to generate random time-series with gaps, and corresponding SFTs
int extractBandFromSFT(SFTtype **outSFT, const SFTtype *inSFT, REAL8 fMin, REAL8 Band)
SFTVector * XLALSignalToSFTs(const REAL4TimeSeries *signalvec, const SFTParams *params)
Turn the given time-series into SFTs and add noise if given.
#define XLAL_INIT_DECL(var,...)
void * XLALCalloc(size_t m, size_t n)
void * XLALRealloc(void *p, size_t n)
SFTtype * XLALSFTVectorToLFT(SFTVector *sfts, REAL8 upsampling)
Turn the given multi-IFO SFTvectors into one long Fourier transform (LFT) over the total observation ...
int XLALSincInterpolateCOMPLEX8TimeSeries(COMPLEX8Vector *y_out, const REAL8Vector *t_out, const COMPLEX8TimeSeries *ts_in, UINT4 Dterms)
Interpolate a given regularly-spaced COMPLEX8 timeseries 'ts_in = x_in(j * dt)' onto new samples 'y_o...
int XLALCompareCOMPLEX8Vectors(VectorComparison *result, const COMPLEX8Vector *x, const COMPLEX8Vector *y, const VectorComparison *tol)
Compare two COMPLEX8 vectors using various different comparison metrics.
SFTtype * XLALSincInterpolateSFT(const SFTtype *sft_in, REAL8 f0Out, REAL8 dfOut, UINT4 numBinsOut, UINT4 Dterms)
(Complex)Sinc-interpolate an input SFT to an output SFT.
int XLALdumpREAL4TimeSeries(const char *fname, const REAL4TimeSeries *series)
int XLALdumpCOMPLEX8TimeSeries(const char *fname, const COMPLEX8TimeSeries *series)
REAL4FFTPlan * XLALCreateForwardREAL4FFTPlan(UINT4 size, int measurelvl)
void XLALDestroyREAL4FFTPlan(REAL4FFTPlan *plan)
int XLALREAL4ForwardFFT(COMPLEX8Vector *output, const REAL4Vector *input, const REAL4FFTPlan *plan)
SFTtype * XLALCreateSFT(UINT4 numBins)
XLAL function to create one SFT-struct.
void XLALDestroySFTVector(SFTVector *vect)
XLAL interface to destroy an SFTVector.
REAL8 TSFTfromDFreq(REAL8 dFreq)
int XLALFindCoveringSFTBins(UINT4 *firstBin, UINT4 *numBins, REAL8 fMinIn, REAL8 BandIn, REAL8 Tsft)
Return the 'effective' frequency-band [fMinEff, fMaxEff] = [firstBin, lastBin] * 1/Tsft,...
void XLALDestroySFT(SFTtype *sft)
Destructor for one SFT.
SFTVector * XLALCreateSFTVector(UINT4 numSFTs, UINT4 numBins)
XLAL function to create an SFTVector of numSFT SFTs with SFTlen frequency-bins (which will be allocat...
REAL4TimeSeries * XLALCreateREAL4TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
COMPLEX8TimeSeries * XLALCreateCOMPLEX8TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
void XLALDestroyCOMPLEX8TimeSeries(COMPLEX8TimeSeries *series)
void XLALDestroyREAL4TimeSeries(REAL4TimeSeries *series)
int XLALUnitCompare(const LALUnit *unit1, const LALUnit *unit2)
COMPLEX8Vector * XLALCreateCOMPLEX8Vector(UINT4 length)
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
void XLALDestroyREAL8Vector(REAL8Vector *vector)
void XLALDestroyCOMPLEX8Vector(COMPLEX8Vector *vector)
int int int XLALPrintInfo(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_CHECK(assertion,...)
#define XLAL_CHECK_NULL(assertion,...)
LIGOTimeGPS * XLALGPSAdd(LIGOTimeGPS *epoch, REAL8 dt)
int XLALGPSCmp(const LIGOTimeGPS *t0, const LIGOTimeGPS *t1)
REAL8 XLALGPSGetREAL8(const LIGOTimeGPS *epoch)
A vector of COMPLEX8FrequencySeries.
COMPLEX8FrequencySeries * data
Pointer to the data array.
UINT4 length
Number of elements in array.
A vector of 'timestamps' of type LIGOTimeGPS.
LIGOTimeGPS * data
array of timestamps
UINT4 length
number of timestamps
Parameters defining the SFTs to be returned from LALSignalToSFTs().
Struct holding the results of comparing two floating-point vectors (real-valued or complex),...