Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALPulsar 7.1.1.1-8a6b96f
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
LFTandTSutilsTest.c
Go to the documentation of this file.
1/*
2 * Copyright (C) 2014 Reinhard Prix
3 * Copyright (C) 2009 Reinhard Prix
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with with program; see the file COPYING. If not, write to the
17 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
18 * MA 02110-1301 USA
19 */
20
21/*********************************************************************************/
22/**
23 * \author R. Prix
24 * \file
25 * \brief Unit tests for LFTandTSutils module
26 *
27 */
28#include "config.h"
29
30/* System includes */
31#include <stdio.h>
32#include <time.h>
33
34/* GSL includes */
35#include <gsl/gsl_rng.h>
36#include <gsl/gsl_randist.h>
37
38
39/* LAL-includes */
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>
47#include <lal/Units.h>
48
49#include <lal/GeneratePulsarSignal.h>
50#include <lal/ComputeFstat.h>
51#include <lal/LFTandTSutils.h>
52
53/* local includes */
54
55/*---------- DEFINES ----------*/
56/* ----- Macros ----- */
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) )
60
61/* ---------- local types ---------- */
62
63/*---------- Global variables ----------*/
65
66/* ----- User-variables: can be set from config-file or command-line */
67/* ---------- local prototypes ---------- */
68int main( void );
69
70
71int test_XLALSFTVectorToLFT( void );
74
76int write_SFTdata( const char *fname, const SFTtype *sft );
77int XLALCompareSFTs( const SFTtype *sft1, const SFTtype *sft2, const VectorComparison *tol );
78static COMPLEX8 testSignal( REAL8 t, REAL8 sigmaN );
79
80/* These three functions were copied from the SFTutils module after being deprecated, to keep this test going. */
81int extractBandFromSFT( SFTtype **outSFT, const SFTtype *inSFT, REAL8 fMin, REAL8 Band );
84
85/*----------------------------------------------------------------------*/
86/* Main Function starts here */
87/*----------------------------------------------------------------------*/
88/**
89 * MAIN function
90 */
91int main( void )
92{
93
95
97
99
101
102 return XLAL_SUCCESS;
103
104} /* main() */
105
106/**
107 * Unit-Test for function XLALSFTVectorToLFT().
108 * Generates random data (timeseries + corresponding SFTs),
109 * then feeds the SFTs into XLALSFTVectorToLFT()
110 * and checks correctness of output Fourier transform by
111 * comparing to FT of original real-valued timeseries.
112 *
113 * \note This indirectly also checks XLALSFTVectorToCOMPLEX8TimeSeries()
114 * which is used by XLALSFTVectorToLFT().
115 *
116 * returns XLAL_SUCCESS on success, XLAL-error otherwise
117 */
118int
120{
121 // ----- generate real-valued random timeseries and corresponding SFTs
122 REAL4TimeSeries *tsR4 = NULL;
123 SFTVector *sfts0 = NULL;
125 UINT4 numSamplesR4 = tsR4->data->length;
126 REAL8 dt_R4 = tsR4->deltaT;
127 REAL8 TspanR4 = numSamplesR4 * dt_R4;
128 // ----- consider only the frequency band [3Hz, 4Hz]
129 REAL8 out_fMin = 3;
130 REAL8 out_Band = 1;
131 SFTVector *sftsBand;
132 XLAL_CHECK( ( sftsBand = extractBandFromSFTVector( sfts0, out_fMin, out_Band ) ) != NULL, XLAL_EFUNC );
133 XLALDestroySFTVector( sfts0 );
134
135 // ----- 1) compute FFT on original REAL4 timeseries -----
136 REAL4FFTPlan *planR4;
137 XLAL_CHECK( ( planR4 = XLALCreateForwardREAL4FFTPlan( numSamplesR4, 0 ) ) != NULL, XLAL_EFUNC );
138
139 SFTtype *lft0;
140 XLAL_CHECK( ( lft0 = XLALCreateSFT( numSamplesR4 / 2 + 1 ) ) != NULL, XLAL_EFUNC );
141 XLAL_CHECK( XLALREAL4ForwardFFT( lft0->data, tsR4->data, planR4 ) == XLAL_SUCCESS, XLAL_EFUNC );
142
143 strcpy( lft0->name, tsR4->name );
144 lft0->f0 = 0;
145 lft0->epoch = tsR4->epoch;
146 REAL8 dfLFT = 1.0 / TspanR4;
147 lft0->deltaF = dfLFT;
148
149 // ----- extract frequency band of interest
150 SFTtype *lftR4 = NULL;
151 XLAL_CHECK( extractBandFromSFT( &lftR4, lft0, out_fMin, out_Band ) == XLAL_SUCCESS, XLAL_EFUNC );
152 for ( UINT4 k = 0; k < lftR4->data->length; k ++ ) {
153 lftR4->data->data[k] *= dt_R4;
154 }
155
156 // ----- 2) compute LFT directly from SFTs ----------
157 SFTtype *lftSFTs0;
158 XLAL_CHECK( ( lftSFTs0 = XLALSFTVectorToLFT( sftsBand, 1 ) ) != NULL, XLAL_EFUNC );
159 XLALDestroySFTVector( sftsBand );
160
161 // ----- re-extract frequency band of interest
162 SFTtype *lftSFTs = NULL;
163 XLAL_CHECK( extractBandFromSFT( &lftSFTs, lftSFTs0, out_fMin, out_Band ) == XLAL_SUCCESS, XLAL_EFUNC );
164
165 if ( lalDebugLevel & LALINFO ) {
166 // ----- write debug output
167 XLAL_CHECK( XLALdumpREAL4TimeSeries( "TS_R4.dat", tsR4 ) == XLAL_SUCCESS, XLAL_EFUNC );
168
169 XLAL_CHECK( write_SFTdata( "LFT_R4T.dat", lftR4 ) == XLAL_SUCCESS, XLAL_EFUNC );
170 XLAL_CHECK( write_SFTdata( "LFT_SFTs.dat", lftSFTs ) == XLAL_SUCCESS, XLAL_EFUNC );
171
172 } // end: debug output
173
174 // ========== compare resulting LFTs ==========
176 XLALPrintInfo( "Comparing LFT with itself: should give 0 for all measures\n" );
177 XLAL_CHECK( XLALCompareSFTs( lftR4, lftR4, &tol0 ) == XLAL_SUCCESS, XLAL_EFUNC );
178 XLAL_CHECK( XLALCompareSFTs( lftSFTs, lftSFTs, &tol0 ) == XLAL_SUCCESS, XLAL_EFUNC );
179
181 tol.relErr_L1 = 4e-2;
182 tol.relErr_L2 = 5e-2;
183 tol.angleV = 5e-2; // rad
184 tol.relErr_atMaxAbsx = 1.2e-2;
185 tol.relErr_atMaxAbsy = 1.2e-2;
186
187 XLALPrintInfo( "Comparing LFT from REAL4-timeseries, to LFT from heterodyned COMPLEX8-timeseries:\n" );
188 XLAL_CHECK( XLALCompareSFTs( lftR4, lftSFTs, &tol ) == XLAL_SUCCESS, XLAL_EFUNC );
189
190 // ---------- free memory ----------
192 XLALDestroyREAL4FFTPlan( planR4 );
193
194 XLALDestroySFT( lft0 );
195 XLALDestroySFT( lftR4 );
196 XLALDestroySFT( lftSFTs );
197 XLALDestroySFT( lftSFTs0 );
198
199 return XLAL_SUCCESS;
200
201} // test_XLALSFTVectorToLFT()
202
203
204int
206{
207
208 COMPLEX8TimeSeries *tsIn;
209 REAL8 f0 = 100; // heterodyning frequency
210 REAL8 dt = 0.1; // sampling frequency = 10Hz
211 LIGOTimeGPS epoch = { 100, 0 };
212 REAL8 tStart = XLALGPSGetREAL8( &epoch );
213 UINT4 numSamples = 1000;
214 REAL8 Tspan = numSamples * dt;
215
216 XLAL_CHECK( ( tsIn = XLALCreateCOMPLEX8TimeSeries( "test TS_in", &epoch, f0, dt, &emptyLALUnit, numSamples ) ) != NULL, XLAL_EFUNC );
217 for ( UINT4 j = 0; j < numSamples; j ++ ) {
218 tsIn->data->data[j] = testSignal( tStart + j * dt, 0 );
219 } // for j < numSamples
220
221 // ---------- interpolate this onto new time-samples
222 UINT4 Dterms = 16;
223 REAL8 safety = ( Dterms + 1.0 ) * dt; // avoid truncated interpolation to minimize errors, set to 0 for seeing boundary-effects [they're not so bad...]
224 LIGOTimeGPS epochOut = epoch;
225 XLALGPSAdd( &epochOut, safety );
226 REAL8 TspanOut = Tspan - 2 * safety;
227
228 REAL8 dtOut = dt / 10;
229 UINT4 numSamplesOut = lround( TspanOut / dtOut );
230 COMPLEX8TimeSeries *tsOut;
231 XLAL_CHECK( ( tsOut = XLALCreateCOMPLEX8TimeSeries( "test TS_out", &epochOut, f0, dtOut, &emptyLALUnit, numSamplesOut ) ) != NULL, XLAL_EFUNC );
232
233
234 REAL8 tStartOut = XLALGPSGetREAL8( &epochOut );
235 REAL8Vector *times_out;
236 XLAL_CHECK( ( times_out = XLALCreateREAL8Vector( numSamplesOut ) ) != NULL, XLAL_EFUNC );
237 for ( UINT4 j = 0; j < numSamplesOut; j ++ ) {
238 REAL8 t_j = tStartOut + j * dtOut;
239 times_out->data[j] = t_j;
240 } // for j < numSamplesOut
241
242 XLAL_CHECK( XLALSincInterpolateCOMPLEX8TimeSeries( tsOut->data, times_out, tsIn, Dterms ) == XLAL_SUCCESS, XLAL_EFUNC );
243 XLALDestroyREAL8Vector( times_out );
244
245 // ---------- check accuracy of interpolation
246 COMPLEX8TimeSeries *tsFull;
247 XLAL_CHECK( ( tsFull = XLALCreateCOMPLEX8TimeSeries( "test TS_full", &epochOut, f0, dtOut, &emptyLALUnit, numSamplesOut ) ) != NULL, XLAL_EFUNC );
248 for ( UINT4 j = 0; j < numSamplesOut; j ++ ) {
249 tsFull->data->data[j] = testSignal( tStartOut + j * dtOut, 0 );
250 } // for j < numSamplesOut
251
252 // ----- out debug info
253 if ( lalDebugLevel & LALINFO ) {
255 XLAL_CHECK( XLALdumpCOMPLEX8TimeSeries( "TS_out.dat", tsOut ) == XLAL_SUCCESS, XLAL_EFUNC );
256 XLAL_CHECK( XLALdumpCOMPLEX8TimeSeries( "TS_full.dat", tsFull ) == XLAL_SUCCESS, XLAL_EFUNC );
257 } // if LALINFO
258
260 tol.relErr_L1 = 2e-2;
261 tol.relErr_L2 = 2e-2;
262 tol.angleV = 2e-2;
263 tol.relErr_atMaxAbsx = 2e-2;
264 tol.relErr_atMaxAbsy = 2e-2;
265
266 XLALPrintInfo( "Comparing sinc-interpolated timeseries to exact signal timeseries:\n" );
268 XLAL_CHECK( XLALCompareCOMPLEX8Vectors( &cmp, tsOut->data, tsFull->data, &tol ) == XLAL_SUCCESS, XLAL_EFUNC );
269
270 // ---------- free memory
274
275 return XLAL_SUCCESS;
276
277} // test_XLALSincInterpolateCOMPLEX8TimeSeries()
278
279int
281{
282 REAL8 f0 = 0; // heterodyning frequency
283 REAL8 sigmaN = 0.001;
284 REAL8 dt = 0.1; // sampling frequency = 10Hz
285 LIGOTimeGPS epoch = { 100, 0 };
286 REAL8 tStart = XLALGPSGetREAL8( &epoch );
287 UINT4 numSamples = 1000;
288 REAL8 Tspan = numSamples * dt;
289 REAL8 df = 1.0 / Tspan;
290
291 UINT4 numSamples0padded = 3 * numSamples;
292 REAL8 Tspan0padded = numSamples0padded * dt;
293 REAL8 df0padded = 1.0 / Tspan0padded;
294
295 UINT4 numBins = NhalfPosDC( numSamples );
296 UINT4 numBins0padded = NhalfPosDC( numSamples0padded );
297 // original timeseries
299 XLAL_CHECK( ( ts = XLALCreateREAL4TimeSeries( "test TS_in", &epoch, f0, dt, &emptyLALUnit, numSamples ) ) != NULL, XLAL_EFUNC );
300 for ( UINT4 j = 0; j < numSamples; j ++ ) {
301 ts->data->data[j] = crealf( testSignal( tStart + j * dt, sigmaN ) );
302 } // for j < numSamples
303
304 // zero-padded to double length
305 REAL4TimeSeries *ts0padded;
306 XLAL_CHECK( ( ts0padded = XLALCreateREAL4TimeSeries( "test TS_padded", &epoch, f0, dt, &emptyLALUnit, numSamples0padded ) ) != NULL, XLAL_EFUNC );
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] ) );
309
310 // compute FFT on ts and ts0padded
311 REAL4FFTPlan *plan, *plan0padded;
312 XLAL_CHECK( ( plan = XLALCreateForwardREAL4FFTPlan( numSamples, 0 ) ) != NULL, XLAL_EFUNC );
313 XLAL_CHECK( ( plan0padded = XLALCreateForwardREAL4FFTPlan( numSamples0padded, 0 ) ) != NULL, XLAL_EFUNC );
314 COMPLEX8Vector *fft, *fft0padded;
316 XLAL_CHECK( ( fft0padded = XLALCreateCOMPLEX8Vector( numBins0padded ) ) != NULL, XLAL_ENOMEM );
317
318 XLAL_CHECK( XLALREAL4ForwardFFT( fft, ts->data, plan ) == XLAL_SUCCESS, XLAL_EFUNC );
319 XLAL_CHECK( XLALREAL4ForwardFFT( fft0padded, ts0padded->data, plan0padded ) == XLAL_SUCCESS, XLAL_EFUNC );
321 XLALDestroyREAL4TimeSeries( ts0padded );
323 XLALDestroyREAL4FFTPlan( plan0padded );
324
325 SFTtype XLAL_INIT_DECL( tmp );
326 tmp.f0 = f0;
327 tmp.deltaF = df;
328 tmp.data = fft;
329
330 REAL8 Band = 0.5 / dt - f0;
331 SFTtype *sft = NULL;
334
335
336 tmp.f0 = f0;
337 tmp.deltaF = df0padded;
338 tmp.data = fft0padded;
339 SFTtype *sft0padded = NULL;
340 XLAL_CHECK( extractBandFromSFT( &sft0padded, &tmp, f0, Band ) == XLAL_SUCCESS, XLAL_EFUNC );
341 XLALDestroyCOMPLEX8Vector( fft0padded );
342
343 // ---------- interpolate input SFT onto frequency bins of padded-ts FFT for comparison
344 UINT4 Dterms = 16;
345 REAL8 safetyBins = ( Dterms + 1.0 ); // avoid truncated interpolation to minimize errors, set to 0 for seeing boundary-effects [they're not so bad...]
346
347 REAL8 fMin = f0 + safetyBins * df;
348 fMin = round( fMin / df0padded ) * df0padded;
349 UINT4 numBinsOut = numBins0padded - 3 * safetyBins;
350 REAL8 BandOut = ( numBinsOut - 1 ) * df0padded;
351
352 SFTtype *sftUpsampled = NULL;
353 XLAL_CHECK( extractBandFromSFT( &sftUpsampled, sft0padded, fMin + 0.5 * df0padded, BandOut - df0padded ) == XLAL_SUCCESS, XLAL_EFUNC );
354
355 SFTtype *sftInterpolated;
356 XLAL_CHECK( ( sftInterpolated = XLALSincInterpolateSFT( sft, fMin, df0padded, numBinsOut, Dterms ) ) != NULL, XLAL_EFUNC );
357
358 // ----- out debug info
359 if ( lalDebugLevel & LALINFO ) {
360 XLAL_CHECK( write_SFTdata( "SFT_in.dat", sft ) == XLAL_SUCCESS, XLAL_EFUNC );
361 XLAL_CHECK( write_SFTdata( "SFT_in0padded.dat", sft0padded ) == XLAL_SUCCESS, XLAL_EFUNC );
362 XLAL_CHECK( write_SFTdata( "SFT_upsampled.dat", sftUpsampled ) == XLAL_SUCCESS, XLAL_EFUNC );
363 XLAL_CHECK( write_SFTdata( "SFT_interpolated.dat", sftInterpolated ) == XLAL_SUCCESS, XLAL_EFUNC );
364 } // if LALINFO
365
366 // ---------- check accuracy of interpolation
368 tol.relErr_L1 = 8e-2;
369 tol.relErr_L2 = 8e-2;
370 tol.angleV = 8e-2;
371 tol.relErr_atMaxAbsx = 4e-3;
372 tol.relErr_atMaxAbsy = 4e-3;
373
374 XLALPrintInfo( "Comparing Dirichlet SFT interpolation with upsampled SFT:\n" );
375 XLAL_CHECK( XLALCompareSFTs( sftUpsampled, sftInterpolated, &tol ) == XLAL_SUCCESS, XLAL_EFUNC );
376
377 // ---------- free memory
378 XLALDestroySFT( sft );
379 XLALDestroySFT( sft0padded );
380 XLALDestroySFT( sftUpsampled );
381 XLALDestroySFT( sftInterpolated );
382
383 return XLAL_SUCCESS;
384
385} // test_XLALSincInterpolateSFT()
386
387
388/**
389 * function to generate random time-series with gaps, and corresponding SFTs
390 */
391int
393{
394 /* input sanity checks */
395 XLAL_CHECK( ( ts != NULL ) && ( ( *ts ) == NULL ), XLAL_EINVAL );
396 XLAL_CHECK( ( sfts != NULL ) && ( ( *sfts ) == NULL ), XLAL_EINVAL );
397
398 // test parameters
399 LIGOTimeGPS epoch0 = { 714180733, 0 };
400 UINT4 numSFTs = 20;
401 REAL8 Tsft = 1000.0;
402
403 // noise sigma
404 REAL8 sigmaN = 0.1;
405
406 /* prepare sampling constants */
407 REAL8 numR4SamplesPerSFT = 2 * 5000;
408 REAL8 dtR4 = ( Tsft / numR4SamplesPerSFT );
409
410 UINT4 numR4SamplesTS = numSFTs * numR4SamplesPerSFT;
411
412 /* ----- allocate timeseries ----- */
413 REAL4TimeSeries *outTS; // input timeseries
414 XLAL_CHECK( ( outTS = XLALCreateREAL4TimeSeries( "H1:test timeseries", &epoch0, 0, dtR4, &emptyLALUnit, numR4SamplesTS ) ) != NULL, XLAL_EFUNC );
415
416 REAL4 *TSdata = outTS->data->data;
417 /* initialize timeseries to zero (to allow for gaps) */
418 memset( TSdata, 0, outTS->data->length * sizeof( *TSdata ) );
419
420 /* also set up corresponding SFT timestamps vector */
421 LIGOTimeGPSVector *timestampsSFT;
422 XLAL_CHECK( ( timestampsSFT = XLALCalloc( 1, sizeof( *timestampsSFT ) ) ) != NULL, XLAL_ENOMEM );
423 XLAL_CHECK( ( timestampsSFT->data = XLALCalloc( numSFTs, sizeof( *timestampsSFT->data ) ) ) != NULL, XLAL_ENOMEM );
424 timestampsSFT->length = numSFTs;
425
426 /* ----- set up random-noise timeseries with gaps ---------- */
427 for ( UINT4 alpha = 0; alpha < numSFTs; alpha ++ ) {
428 /* record this SFT's timestamp */
429 timestampsSFT->data[alpha] = epoch0;
430 timestampsSFT->data[alpha].gpsSeconds += lround( alpha * Tsft );
431
432 /* generate all data-points of this SFT */
433 for ( UINT4 j = 0; j < numR4SamplesPerSFT; j++ ) {
434 UINT4 alpha_j = alpha * numR4SamplesPerSFT + j;
435 REAL8 ti = alpha * Tsft + j * dtR4;
436 /* unit-variance Gaussian noise + sinusoid */
437 TSdata[alpha_j] = crealf( testSignal( ti, sigmaN ) );
438 } // for js < numR4SamplesPerSFT
439
440 } /* for alpha < numSFTs */
441
442 /* ----- generate SFTs from this timeseries ---------- */
443 SFTParams XLAL_INIT_DECL( sftParams );
444 sftParams.Tsft = Tsft;
445 sftParams.timestamps = timestampsSFT;
446 sftParams.noiseSFTs = NULL;
447 sftParams.window = NULL;
448
449 SFTVector *outSFTs;
450 XLAL_CHECK( ( outSFTs = XLALSignalToSFTs( outTS, &sftParams ) ) != NULL, XLAL_EFUNC );
451
452 /* free memory */
453 XLALFree( timestampsSFT->data );
454 XLALFree( timestampsSFT );
455
456 /* return timeseries and SFTvector */
457 ( *ts ) = outTS;
458 ( *sfts ) = outSFTs;
459
460 return XLAL_SUCCESS;
461
462} // XLALgenerateRandomData()
463
464// compare two SFTs, return XLAL_SUCCESS if OK, error otherwise
465int
466XLALCompareSFTs( const SFTtype *sft1, const SFTtype *sft2, const VectorComparison *tol )
467{
468 // check input sanity
469 XLAL_CHECK( sft1 != NULL, XLAL_EINVAL );
470 XLAL_CHECK( sft2 != NULL, XLAL_EINVAL );
471 XLAL_CHECK( tol != NULL, XLAL_EINVAL );
472
473 XLAL_CHECK( XLALGPSCmp( &( sft1->epoch ), &( sft2->epoch ) ) == 0, XLAL_ETOL );
474 REAL8 tolFreq = 10 * LAL_REAL8_EPS;
475 REAL8 err_f0 = sft1->f0 - sft2->f0;
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 );
477 REAL8 relErr_df = RELERR( sft1->deltaF, sft2->deltaF );
478 XLAL_CHECK( relErr_df < tolFreq, XLAL_ETOL, "dFreq1 = %g, dFreq2 = %g, relErr_df = %g (tol = %g)", sft1->deltaF, sft2->deltaF, relErr_df, tolFreq );
479
480 XLAL_CHECK( XLALUnitCompare( &( sft1->sampleUnits ), &( sft2->sampleUnits ) ) == 0, XLAL_ETOL );
481
482 XLAL_CHECK( sft1->data->length == sft2->data->length, XLAL_ETOL, "sft1->length = %d, sft2->length = %d\n", sft1->data->length, sft2->data->length );
483
486
487 return XLAL_SUCCESS;
488
489} // XLALCompareSFTs()
490
491
492static COMPLEX8
494{
495 static gsl_rng *r = NULL;
496 static BOOLEAN firstCall = 1;
497 if ( firstCall ) {
498 XLAL_CHECK( ( r = gsl_rng_alloc( gsl_rng_taus2 ) ) != NULL, XLAL_EFAILED );
499 gsl_rng_set( r, 13 ); /* sets random-number seed */
500 firstCall = 0;
501 }
502
503 REAL8 amp1 = 0.001;
504 REAL8 amp2 = 1;
505 REAL8 freq1 = LAL_PI;
506 REAL8 freq2 = 2;
507
508 COMPLEX8 y = amp1 * sin( LAL_TWOPI * freq1 * t ) + I * amp2 * cos( LAL_TWOPI * freq2 * t );
509 y += gsl_ran_gaussian( r, sigmaN );
510 y += I * gsl_ran_gaussian( r, sigmaN );
511
512 return y;
513} // testSignal()
514
515int
516write_SFTdata( const char *fname, const SFTtype *sft )
517{
518 XLAL_CHECK( fname != NULL, XLAL_EINVAL );
519 XLAL_CHECK( sft != NULL, XLAL_EINVAL );
520
521 REAL8 f0 = sft->f0;
522 REAL8 df = sft->deltaF;
523
524 FILE *fp;
525 XLAL_CHECK( ( fp = fopen( fname, "wb" ) ) != NULL, XLAL_ESYS );
526 for ( UINT4 k = 0; k < sft->data->length; k ++ ) {
527 REAL8 fk = f0 + k * df;
528 fprintf( fp, "%.9f % 6e % 6e \n", fk, crealf( sft->data->data[k] ), cimagf( sft->data->data[k] ) );
529 } // for k < numFreqBins
530
531 fclose( fp );
532
533 return XLAL_SUCCESS;
534
535} // write_SFTdata()
536
537SFTVector *
538extractBandFromSFTVector( const SFTVector *inSFTs, ///< [in] input SFTs
539 REAL8 fMin, ///< [in] lower end of frequency interval to return
540 REAL8 Band ///< [in] band width of frequency interval to return
541 )
542{
543 XLAL_CHECK_NULL( inSFTs != NULL, XLAL_EINVAL, "Invalid NULL input SFT vector 'inSFTs'\n" );
544 XLAL_CHECK_NULL( inSFTs->length > 0, XLAL_EINVAL, "Invalid zero-length input SFT vector 'inSFTs'\n" );
545 XLAL_CHECK_NULL( fMin >= 0, XLAL_EDOM, "Invalid negative frequency fMin = %g\n", fMin );
546 XLAL_CHECK_NULL( Band > 0, XLAL_EDOM, "Invalid non-positive Band = %g\n", Band );
547
548 UINT4 numSFTs = inSFTs->length;
549
550 SFTVector *ret;
551 XLAL_CHECK_NULL( ( ret = XLALCreateSFTVector( numSFTs, 0 ) ) != NULL, XLAL_EFUNC );
552
553 for ( UINT4 i = 0; i < numSFTs; i ++ ) {
554 SFTtype *dest = &( ret->data[i] );
555 SFTtype *src = &( inSFTs->data[i] );
556
558
559 } /* for i < numSFTs */
560
561 /* return final SFT-vector */
562 return ret;
563
564} /* extractBandFromSFTVector() */
565
566int
567extractBandFromSFT( SFTtype **outSFT, ///< [out] output SFT (alloc'ed or re-alloced as required)
568 const SFTtype *inSFT, ///< [in] input SFT
569 REAL8 fMin, ///< [in] lower end of frequency interval to return
570 REAL8 Band ///< [in] band width of frequency interval to return
571 )
572{
573 XLAL_CHECK( outSFT != NULL, XLAL_EINVAL );
574 XLAL_CHECK( inSFT != NULL, XLAL_EINVAL );
575 XLAL_CHECK( inSFT->data != NULL, XLAL_EINVAL );
576 XLAL_CHECK( fMin >= 0, XLAL_EDOM, "Invalid negative frequency fMin = %g\n", fMin );
577 XLAL_CHECK( Band > 0, XLAL_EDOM, "Invalid non-positive Band = %g\n", Band );
578
579 REAL8 df = inSFT->deltaF;
581
582 REAL8 fMinSFT = inSFT->f0;
583 UINT4 numBinsSFT = inSFT->data->length;
584 UINT4 firstBinSFT = round( fMinSFT / df ); // round to closest bin
585 UINT4 lastBinSFT = firstBinSFT + ( numBinsSFT - 1 );
586
587 // find 'covering' SFT-band to extract
588 UINT4 firstBinExt, numBinsExt;
589 XLAL_CHECK( XLALFindCoveringSFTBins( &firstBinExt, &numBinsExt, fMin, Band, Tsft ) == XLAL_SUCCESS, XLAL_EFUNC );
590 UINT4 lastBinExt = firstBinExt + ( numBinsExt - 1 );
591
592 XLAL_CHECK( firstBinExt >= firstBinSFT && ( lastBinExt <= lastBinSFT ), XLAL_EINVAL,
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 );
595
596 INT4 firstBinOffset = firstBinExt - firstBinSFT;
597
598 if ( ( *outSFT ) == NULL ) {
599 XLAL_CHECK( ( ( *outSFT ) = XLALCalloc( 1, sizeof( *( *outSFT ) ) ) ) != NULL, XLAL_ENOMEM );
600 }
601 if ( ( *outSFT )->data == NULL ) {
602 XLAL_CHECK( ( ( *outSFT )->data = XLALCreateCOMPLEX8Vector( numBinsExt ) ) != NULL, XLAL_EFUNC );
603 }
604 if ( ( *outSFT )->data->length != numBinsExt ) {
605 XLAL_CHECK( ( ( *outSFT )->data->data = XLALRealloc( ( *outSFT )->data->data, numBinsExt * sizeof( ( *outSFT )->data->data[0] ) ) ) != NULL, XLAL_ENOMEM );
606 ( *outSFT )->data->length = numBinsExt;
607 }
608
609 COMPLEX8Vector *ptr = ( *outSFT )->data; // keep copy to data-pointer
610 ( *( *outSFT ) ) = ( *inSFT ); // copy complete header
611 ( *outSFT )->data = ptr; // restore data-pointer
612 ( *outSFT )->f0 = firstBinExt * df ; // set correct new fMin
613
614 /* copy the relevant part of the data */
615 memcpy( ( *outSFT )->data->data, inSFT->data->data + firstBinOffset, numBinsExt * sizeof( ( *outSFT )->data->data[0] ) );
616
617 return XLAL_SUCCESS;
618
619} // extractBandFromSFT()
int j
ProcessParamsTable * ptr
int k
void LALCheckMemoryLeaks(void)
int XLALCompareSFTs(const SFTtype *sft1, const SFTtype *sft2, const VectorComparison *tol)
int test_XLALSincInterpolateSFT(void)
static LALUnit emptyLALUnit
#define RELERR(x, y)
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)
REAL8 tol
#define fprintf
double e
SFTVector * XLALSignalToSFTs(const REAL4TimeSeries *signalvec, const SFTParams *params)
Turn the given time-series into SFTs and add noise if given.
#define LAL_REAL8_EPS
#define LAL_PI
#define LAL_TWOPI
unsigned char BOOLEAN
double REAL8
#define XLAL_INIT_DECL(var,...)
uint32_t UINT4
float complex COMPLEX8
int32_t INT4
float REAL4
#define lalDebugLevel
LALINFO
void * XLALCalloc(size_t m, size_t n)
void * XLALRealloc(void *p, size_t n)
void XLALFree(void *p)
SFTtype * XLALSFTVectorToLFT(SFTVector *sfts, REAL8 upsampling)
Turn the given multi-IFO SFTvectors into one long Fourier transform (LFT) over the total observation ...
Definition: LFTandTSutils.c:64
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.
#define NhalfPosDC(N)
Definition: LFTandTSutils.h:46
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)
static const INT4 r
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.
Definition: SFTtypes.c:152
void XLALDestroySFTVector(SFTVector *vect)
XLAL interface to destroy an SFTVector.
Definition: SFTtypes.c:300
REAL8 TSFTfromDFreq(REAL8 dFreq)
Definition: SFTtypes.c:1144
int XLALFindCoveringSFTBins(UINT4 *firstBin, UINT4 *numBins, REAL8 fMinIn, REAL8 BandIn, REAL8 Tsft)
Return the 'effective' frequency-band [fMinEff, fMaxEff] = [firstBin, lastBin] * 1/Tsft,...
Definition: SFTtypes.c:106
void XLALDestroySFT(SFTtype *sft)
Destructor for one SFT.
Definition: SFTtypes.c:176
SFTVector * XLALCreateSFTVector(UINT4 numSFTs, UINT4 numBins)
XLAL function to create an SFTVector of numSFT SFTs with SFTlen frequency-bins (which will be allocat...
Definition: SFTtypes.c:230
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,...)
XLAL_ENOMEM
XLAL_SUCCESS
XLAL_EFUNC
XLAL_ETOL
XLAL_EDOM
XLAL_ESYS
XLAL_EINVAL
XLAL_EFAILED
LIGOTimeGPS * XLALGPSAdd(LIGOTimeGPS *epoch, REAL8 dt)
int XLALGPSCmp(const LIGOTimeGPS *t0, const LIGOTimeGPS *t1)
REAL8 XLALGPSGetREAL8(const LIGOTimeGPS *epoch)
list y
src
ts
double alpha
CHAR name[LALNameLength]
COMPLEX8Sequence * data
A vector of COMPLEX8FrequencySeries.
COMPLEX8FrequencySeries * data
Pointer to the data array.
UINT4 length
Number of elements in array.
COMPLEX8Sequence * data
COMPLEX8 * data
A vector of 'timestamps' of type LIGOTimeGPS.
Definition: SFTfileIO.h:188
LIGOTimeGPS * data
array of timestamps
Definition: SFTfileIO.h:193
UINT4 length
number of timestamps
Definition: SFTfileIO.h:192
CHAR name[LALNameLength]
REAL4Sequence * data
LIGOTimeGPS epoch
REAL4 * data
REAL8 * data
Parameters defining the SFTs to be returned from LALSignalToSFTs().
Struct holding the results of comparing two floating-point vectors (real-valued or complex),...
Definition: LFTandTSutils.h:64
enum @1 epoch
double df