27 #include <lal/LFTandTSutils.h>
28 #include <lal/SFTfileIO.h>
29 #include <lal/LogPrintf.h>
30 #include <lal/TimeSeries.h>
31 #include <lal/ComplexFFT.h>
32 #include <lal/ComputeFstat.h>
33 #include <lal/SinCosLUT.h>
34 #include <lal/Factorial.h>
35 #include <lal/Window.h>
38 #define MYMAX(x,y) ( (x) > (y) ? (x) : (y) )
39 #define MYMIN(x,y) ( (x) < (y) ? (x) : (y) )
40 #define SQ(x) ((x)*(x))
42 #define FINITE_OR_ONE(x) (((x) != 0) ? (x) : 1.0)
43 #define cRELERR(x,y) ( cabsf( (x) - (y) ) / FINITE_OR_ONE( 0.5 * (cabsf(x) + cabsf(y)) ) )
44 #define fRELERR(x,y) ( fabsf( (x) - (y) ) / FINITE_OR_ONE( 0.5 * (fabsf(x) + fabsf(y)) ) )
45 #define OOTWOPI (1.0 / LAL_TWOPI)
46 #define OOPI (1.0 / LAL_PI)
47 #define LD_SMALL4 (2.0e-4)
83 REAL8 Tspan0 = numSamples0 *
dt;
92 UINT4 numSFTsFit = lround( ( Tspan0 * upsampling ) /
Tsft );
97 if ( numSamples > numSamples0 ) {
100 memset( lTS->
data->
data + numSamples0, 0, ( numSamples - numSamples0 ) *
sizeof( lTS->
data->
data[0] ) );
106 REAL8 fHet = f0SFT + 1.0 * NnegSFT * dfSFT;
116 strcpy( outputLFT->
name, firstSFT->
name );
117 strncat( outputLFT->
name,
":long Fourier transform",
sizeof( outputLFT->
name ) - 1 - strlen( outputLFT->
name ) );
119 outputLFT->
f0 = f0LFT;
124 COMPLEX8FFTPlan *LFTplan;
180 REAL8 fHet = f0SFT + 1.0 * NnegSFT * dfSFT;
183 COMPLEX8FFTPlan *SFTplan;
205 nudge_n = 1
e-9 * round( nudge_n * 1e9 );
214 REAL8 hetCycles = fmod( fHet * offsetEff, 1 );
216 if ( nudge_n != 0 ) {
217 XLALPrintInfo(
"n = %d, offset_n = %g, nudge_n = %g, offset = %g, offsetEff = %g, hetCycles = %g\n",
218 n, offset_n, nudge_n, offset, offsetEff, hetCycles );
221 REAL4 hetCorrection_re, hetCorrection_im;
223 COMPLEX8 hetCorrection =
crectf( hetCorrection_re, hetCorrection_im );
234 hetCorrection *= dfSFT;
244 UINT4 binsLeft = numSamples - bin0_n;
245 UINT4 copyLen =
MYMIN( numFreqBinsSFT, binsLeft );
307 memcpy( tmp, X->
data,
N *
sizeof( *tmp ) );
310 memcpy( X->
data, tmp + Npos_and_DC, Nneg *
sizeof( *tmp ) );
313 memcpy( X->
data + Nneg, tmp, Npos_and_DC *
sizeof( *tmp ) );
339 memcpy( tmp, X->
data,
N *
sizeof( *tmp ) );
342 memcpy( X->
data + Npos_and_DC, tmp, Nneg *
sizeof( *tmp ) );
345 memcpy( X->
data, tmp + Nneg, Npos_and_DC *
sizeof( *tmp ) );
374 REAL8 shiftCyles = shift * fk;
375 REAL4 fact_re, fact_im;
439 for (
UINT4 k = 0;
k <
x->data->length;
k++ ) {
441 REAL8 shiftCycles = shift * tk;
442 REAL4 fact_re, fact_im;
448 x->data->data[
k] *= fact;
484 while ( ( nspins > 0 ) && ( doppler->
fkdot[nspins] == 0 ) ) {
490 for (
UINT4 k = 0;
k < numSamples;
k++ ) {
491 REAL8 tk_pow_jp1 = tk;
493 for (
UINT4 j = 1;
j <= nspins;
j++ ) {
499 REAL4 cosphase, sinphase;
528 if ( ! multiTimes ) {
531 if ( multiTimes->
data != NULL ) {
587 memcpy(
out, ttimes,
sizeof( *ttimes ) );
639 ( *ts_out ) = ( *ts_in );
667 REAL8 x_L1 = 0, x_L2 = 0;
669 REAL8 y_L1 = 0, y_L2 = 0;
671 REAL8 diff_L1 = 0, diff_L2 = 0;
674 REAL8 maxAbsx = 0, maxAbsy = 0;
675 COMPLEX8 x_atMaxAbsx = 0, y_atMaxAbsx = 0;
676 COMPLEX8 x_atMaxAbsy = 0, y_atMaxAbsy = 0;
679 UINT4 numSamples =
x->length;
680 for (
UINT4 i = 0;
i < numSamples;
i ++ ) {
683 REAL8 xAbs2_i = x_i * conj( x_i );
684 REAL8 yAbs2_i = y_i * conj( y_i );
685 REAL8 xAbs_i = sqrt( xAbs2_i );
686 REAL8 yAbs_i = sqrt( yAbs2_i );
687 XLAL_CHECK( isfinite( xAbs_i ),
XLAL_EFPINVAL,
"non-finite element: x(%d) = %g + I %g\n",
i, creal( x_i ), cimag( x_i ) );
688 XLAL_CHECK( isfinite( yAbs_i ),
XLAL_EFPINVAL,
"non-finite element: y(%d) = %g + I %g\n",
i, creal( y_i ), cimag( y_i ) );
690 REAL8 absdiff = cabs( x_i - y_i );
692 diff_L2 +=
SQ( absdiff );
699 scalar += x_i * conj( y_i );
701 if ( xAbs_i > maxAbsx ) {
706 if ( yAbs_i > maxAbsy ) {
715 x_L2 = sqrt( x_L2sq );
716 y_L2 = sqrt( y_L2sq );
717 diff_L2 = sqrt( diff_L2 );
718 REAL8 sinAngle2 = ( x_L2sq * y_L2sq - scalar * conj( scalar ) ) /
FINITE_OR_ONE( x_L2sq * y_L2sq );
719 REAL8 angle = asin( sqrt( fabs( sinAngle2 ) ) );
752 REAL8 x_L1 = 0, x_L2 = 0;
754 REAL8 y_L1 = 0, y_L2 = 0;
756 REAL8 diff_L1 = 0, diff_L2 = 0;
759 REAL4 maxAbsx = 0, maxAbsy = 0;
760 REAL4 x_atMaxAbsx = 0, y_atMaxAbsx = 0;
761 REAL4 x_atMaxAbsy = 0, y_atMaxAbsy = 0;
763 UINT4 numSamples =
x->length;
764 for (
UINT4 i = 0;
i < numSamples;
i ++ ) {
769 REAL8 xAbs_i = fabs( x_i );
770 REAL8 yAbs_i = fabs( y_i );
774 REAL8 absdiff = fabs( x_i - y_i );
776 diff_L2 +=
SQ( absdiff );
785 if ( xAbs_i > maxAbsx ) {
790 if ( yAbs_i > maxAbsy ) {
799 x_L2 = sqrt( x_L2sq );
800 y_L2 = sqrt( y_L2sq );
801 diff_L2 = sqrt( diff_L2 );
807 result->
angleV = asin( sqrt( fabs( sinAngle2 ) ) );
832 const char *
names[5] = {
"relErr_L1",
"relErr_L2",
"angleV",
"relErr_atMaxAbsx",
"relErr_atMaxAbsy" };
834 const REAL4 tols[5] = { localTol.relErr_L1, localTol.relErr_L2, localTol.angleV, localTol.relErr_atMaxAbsx, localTol.relErr_atMaxAbsy };
840 if ( results[
i] > tols[
i] ) {
846 ( results[
i] > tols[
i] ) ?
"FAILED. Exceeded tolerance." :
"OK." );
908 UINT4 winLen = 2 * Dterms + 1;
913 for (
UINT4 l = 0;
l < numSamplesOut;
l ++ ) {
917 if ( (
t < 0 ) || (
t > ( numSamplesIn - 1 )*
dt ) ) {
923 INT8 jstar = lround( t_by_dt );
925 if ( fabs( t_by_dt - jstar ) <
LD_SMALL4 ) {
930 INT4 jStart0 = jstar - Dterms;
931 UINT4 jEnd0 = jstar + Dterms;
933 UINT4 jEnd =
MYMIN( jEnd0, numSamplesIn - 1 );
935 REAL4 delta_jStart = ( t_by_dt - jStart );
941 REAL8 delta_j = delta_jStart;
942 for (
UINT8 j = jStart;
j <= jEnd;
j ++ ) {
943 COMPLEX8 Cj =
win->data->data[
j - jStart0] * sin0oopi / delta_j;
947 sin0oopi = -sin0oopi;
951 y_out->
data[
l] = y_l;
998 for (
UINT4 l = 0;
l < numBinsOut;
l ++ ) {
1002 if ( (
f < 0 ) || (
f > ( numBinsIn - 1 )*
df ) ) {
1007 INT8 kstar = lround(
f * oodf );
1009 if ( fabs(
f - kstar *
df ) < 1
e-6 ) {
1017 REAL4 delta_kstar = (
f - kstar *
df ) * oodf;
1018 REAL4 sin2pidelta, cos2pidelta;
1021 COMPLEX8 prefact =
OOTWOPI * ( sin2pidelta + I * ( cos2pidelta - 1.0f ) );
1024 REAL4 delta_k = delta_kstar + Dterms;
1030 y_out->
data[
l] = prefact * y_l;
1058 for (
UINT4 k = 0;
k < numBinsOut;
k ++ ) {
1059 f_out->
data[
k] = f0Out +
k * dfOut;
1064 ( *out ) = ( *sft_in );
1066 out->deltaF = dfOut;
1096 newLen = oldLen * refineby;
1103 for (
l = 0;
l < newLen;
l++ ) {
1106 REAL8 remain, kstarREAL;
1107 UINT4 kstar, kmin, kmax,
k;
1109 REAL8 Yk_re, Yk_im, Xd_re, Xd_im;
1111 kstarREAL = 1.0 *
l / refineby;
1112 kstar = lround( kstarREAL );
1113 kstar =
MYMIN( kstar, oldLen - 1 );
1114 remain = kstarREAL - kstar;
1118 kmax =
MYMIN( oldLen, kstar + Dterms );
1121 if ( fabs( remain ) > 1
e-5 ) {
1124 coskm1 = cos(
LAL_TWOPI * remain ) - 1.0;
1127 for (
k = kmin;
k < kmax;
k++ ) {
1128 REAL8 Plk_re, Plk_im;
1130 Xd_re = crealf( in->
data[
k] );
1131 Xd_im = cimagf( in->
data[
k] );
1133 kappa_l_k = kstarREAL -
k;
1135 Plk_re = sink / kappa_l_k;
1136 Plk_im = coskm1 / kappa_l_k;
1138 Yk_re += Plk_re * Xd_re - Plk_im * Xd_im;
1139 Yk_im += Plk_re * Xd_im + Plk_im * Xd_re;
static LALUnit emptyLALUnit
COMPLEX8FFTPlan * XLALCreateReverseCOMPLEX8FFTPlan(UINT4 size, int measurelvl)
int XLALCOMPLEX8VectorFFT(COMPLEX8Vector *_LAL_RESTRICT_ output, const COMPLEX8Vector *_LAL_RESTRICT_ input, const COMPLEX8FFTPlan *plan)
COMPLEX8FFTPlan * XLALCreateForwardCOMPLEX8FFTPlan(UINT4 size, int measurelvl)
void XLALDestroyCOMPLEX8FFTPlan(COMPLEX8FFTPlan *plan)
static const REAL8 LAL_FACT_INV[]
#define XLAL_INIT_DECL(var,...)
void * XLALCalloc(size_t m, size_t n)
void * XLALMalloc(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 XLALReorderFFTWtoSFT(COMPLEX8Vector *X)
Change frequency-bin order from fftw-convention to a 'SFT' ie.
MultiCOMPLEX8TimeSeries * XLALMultiSFTVectorToCOMPLEX8TimeSeries(const MultiSFTVector *multisfts)
Turn the given multiSFTvector into multiple long COMPLEX8TimeSeries, properly dealing with gaps.
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 XLALFrequencyShiftMultiCOMPLEX8TimeSeries(MultiCOMPLEX8TimeSeries *x, const REAL8 shift)
Multi-detector wrapper for XLALFrequencyShiftCOMPLEX8TimeSeries NOTE: this modifies the MultiCOMPLEX8...
int XLALTimeShiftSFT(SFTtype *sft, REAL8 shift)
Time-shift the given SFT by an amount of 'shift' seconds, using the frequency-domain expression ,...
int XLALFrequencyShiftCOMPLEX8TimeSeries(COMPLEX8TimeSeries *x, const REAL8 shift)
Freq-shift the given COMPLEX8Timeseries by an amount of 'shift' Hz, using the time-domain expression ...
COMPLEX8TimeSeries * XLALSFTVectorToCOMPLEX8TimeSeries(const SFTVector *sftsIn)
Turn the given SFTvector into one long time-series, properly dealing with gaps.
int XLALCopyMultiCOMPLEX8TimeSeries(MultiCOMPLEX8TimeSeries *multiTimesOut, MultiCOMPLEX8TimeSeries *multiTimesIn)
Copies a MultiCOMPLEX8TimeSeries structure, output must be allocated of identical size as input!
int XLALCheckVectorComparisonTolerances(const VectorComparison *result, const VectorComparison *tol)
Check VectorComparison result against specified tolerances, to allow for standardized comparison and ...
int XLALReorderSFTtoFFTW(COMPLEX8Vector *X)
Change frequency-bin order from 'SFT' to fftw-convention ie.
int XLALCompareREAL4Vectors(VectorComparison *result, const REAL4Vector *x, const REAL4Vector *y, const VectorComparison *tol)
Compare two REAL4 vectors using various different comparison metrics.
int XLALCompareCOMPLEX8Vectors(VectorComparison *result, const COMPLEX8Vector *x, const COMPLEX8Vector *y, const VectorComparison *tol)
Compare two COMPLEX8 vectors using various different comparison metrics.
COMPLEX8TimeSeries * XLALDuplicateCOMPLEX8TimeSeries(COMPLEX8TimeSeries *ttimes)
Duplicates a COMPLEX8TimeSeries structure.
int XLALSpinDownCorrectionMultiTS(MultiCOMPLEX8TimeSeries *multiTimeSeries, const PulsarDopplerParams *doppler)
Apply a spin-down correction to the complex8 timeseries using the time-domain expression y(t) = x(t) ...
MultiCOMPLEX8TimeSeries * XLALDuplicateMultiCOMPLEX8TimeSeries(MultiCOMPLEX8TimeSeries *multiTimes)
Duplicates a MultiCOMPLEX8TimeSeries structure.
int XLALCopyCOMPLEX8TimeSeries(COMPLEX8TimeSeries *ts_out, COMPLEX8TimeSeries *ts_in)
Copies a COMPLEX8TimeSeries structure, output must be allocated of identical size as input!
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 XLALSincInterpolateCOMPLEX8FrequencySeries(COMPLEX8Vector *y_out, const REAL8Vector *f_out, const COMPLEX8FrequencySeries *fs_in, UINT4 Dterms)
Interpolate a given regularly-spaced COMPLEX8 frequency-series 'fs_in = x_in( k * df)' onto new sampl...
COMPLEX8Vector * XLALrefineCOMPLEX8Vector(const COMPLEX8Vector *in, UINT4 refineby, UINT4 Dterms)
Interpolate frequency-series to newLen frequency-bins.
void XLALDestroyMultiCOMPLEX8TimeSeries(MultiCOMPLEX8TimeSeries *multiTimes)
Destroy a MultiCOMPLEX8TimeSeries structure.
#define PULSAR_MAX_SPINS
maximal number of spin-parameters (Freq + spindowns) we can handle
SFTtype * XLALCreateSFT(UINT4 numBins)
XLAL function to create one SFT-struct.
SFTVector * XLALDuplicateSFTVector(const SFTVector *sftsIn)
Create a complete copy of an SFT vector.
void XLALDestroySFTVector(SFTVector *vect)
XLAL interface to destroy an SFTVector.
int XLALSinCosLUT(REAL4 *sinx, REAL4 *cosx, REAL8 x)
Calculate sin(x) and cos(x) to roughly 1e-7 precision using a lookup-table and Taylor-expansion.
int XLALSinCos2PiLUT(REAL4 *sin2pix, REAL4 *cos2pix, REAL8 x)
Calculate sin(2*pi*x) and cos(2*pi*x) to roughly 1e-7 precision using a lookup-table and Taylor-expan...
int XLALSinCos2PiLUTtrimmed(REAL4 *s, REAL4 *c, REAL8 x)
A function that uses the lookup tables to evaluate sin and cos values of 2*Pi*x, but relies on x bein...
COMPLEX8TimeSeries * XLALCreateCOMPLEX8TimeSeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaT, const LALUnit *sampleUnits, size_t length)
void XLALDestroyCOMPLEX8TimeSeries(COMPLEX8TimeSeries *series)
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
void XLALDestroyREAL8Vector(REAL8Vector *vector)
COMPLEX8Vector * XLALCreateCOMPLEX8Vector(UINT4 length)
void XLALDestroyREAL8Window(REAL8Window *window)
REAL8Window * XLALCreateHammingREAL8Window(UINT4 length)
int int int XLALPrintInfo(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_CHECK(assertion,...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_CHECK_NULL(assertion,...)
LIGOTimeGPS * XLALGPSAdd(LIGOTimeGPS *epoch, REAL8 dt)
REAL8 XLALGPSGetREAL8(const LIGOTimeGPS *epoch)
REAL8 XLALGPSDiff(const LIGOTimeGPS *t1, const LIGOTimeGPS *t0)
A vector of COMPLEX8FrequencySeries.
COMPLEX8FrequencySeries * data
Pointer to the data array.
UINT4 length
Number of elements in array.
Multi-IFO container for COMPLEX8 resampled timeseries.
COMPLEX8TimeSeries ** data
array of COMPLEX8TimeSeries (pointers)
UINT4 length
number of IFOs
A collection of SFT vectors – one for each IFO in a multi-IFO search.
UINT4 length
number of ifos
SFTVector ** data
sftvector for each ifo
Type containing the 'Doppler-parameters' affecting the time-evolution of the phase.
PulsarSpins fkdot
Intrinsic spins: [Freq, f1dot, f2dot, ...
LIGOTimeGPS refTime
Reference time of pulsar parameters (in SSB!)
Struct holding the results of comparing two floating-point vectors (real-valued or complex),...
REAL4 relErr_L2
relative error between vectors using L2 norm
REAL4 angleV
angle between the two vectors, according to
REAL4 relErr_L1
relative error between vectors using L1 norm
REAL4 relErr_atMaxAbsy
single-sample relative error at maximum |sample-value| of second vector 'x'
REAL4 relErr_atMaxAbsx
single-sample relative error at maximum |sample-value| of first vector 'x'