56#include <lal/LALStdlib.h>
57#include <lal/LALHashTbl.h>
59#include <lal/SFTfileIO.h>
60#include <lal/SFTReferenceLibrary.h>
61#include <lal/LALPulsarVCSInfo.h>
71#define COMPARE_BY( x, y ) do { if ( (x) < (y) ) return -1; if ( (x) > (y) ) return +1; } while(0)
104 if ( seconds > 10 ) {
107 if ( seconds == 0 ) {
118 struct stat stat_buf;
120 if (
stat( fname, &stat_buf ) ) {
124 if ( ! S_ISDIR( stat_buf.st_mode ) ) {
165 fseek( fpin, move, SEEK_CUR );
170int main(
int argc,
char **argv )
178 char *cmdline = NULL;
182 char outdir0[] =
".";
183 char *
outdir = (
char * )outdir0;
185 int firstfile =
TRUE;
186 int allcomments =
FALSE;
188 unsigned int startBin = 0, endBin = 0;
189 unsigned int width = 0, overlap = 0;
190 double fMin = -1.0, fMax = -1.0;
191 double fWidth = -1.0, fOverlap = -1;
192 unsigned int nactivesamples;
198 int assumeSorted = 0;
200 LALHashTbl *nbsfts = NULL;
209 if ( ( argv[1] == NULL ) ||
210 ( strcmp( argv[1],
"-h" ) == 0 ) ||
211 ( strcmp( argv[1],
"--help" ) == 0 ) ) {
215 " Write this help message\n"
218 " [-a|--all-comments]\n"
219 " [-c|--add-comment 0|1|2]\n"
220 " [-v|--no-validation]\n"
221 " [-s|--start-bin <startbin>]\n"
222 " [-e|--end-bin <endbin (exclusively)>]\n"
223 " [-b|--width <sftbins>]\n"
224 " [-x|--overlap <overlap>]\n"
225 " [-fs|--start-frequency <startfrequency>]\n"
226 " [-fe|--end-frequency <endfrequency (exclusively)>]\n"
227 " [-fb|--frequency-bandwidth <frequencywidth>]\n"
228 " [-fx|--frequency-overlap <frequencyoverlap>]\n"
229 " [-ts|--minStartTime <minStartTime>]\n"
230 " [-te|--maxStartTime <maxStartTime (exclusively)>]\n"
231 " [-as|--assumeSorted 1|0]\n"
232 " [-m|--factor <factor>]\n"
233 " [-n|--output-directory <outputdirectory>]\n"
234 " [--] <inputfile> ...\n"
236 " This program reads in binary SFTs and writes out narrow-banded\n"
239 " The frequency bands of the ouput SFTs (first frequency bin of first output SFT,\n"
240 " last frequency bin of last output SFT, number of bins in each output SFT)\n"
241 " and a possible overlap of the output files can be specified\n"
242 " either in bins ('-s', '-e', '-b', -'x') or Hz ('-fs', '-fe', '-fb', '-fx')\n"
243 " (or mixed - if both are given, frequency values take precedence).\n"
245 " A 'mystery factor' can be specified with '-m' option.\n"
247 " '-v' skips validation of the CRC checksum of the input files\n"
249 " The name of the output SFT follows the standard SFT filename convention, with\n"
250 " information specific to narrow-band SFTs included in the description:\n"
251 " <outdir>/<obs>-<nSFT>_<det>_<timebase>SFT_NBF<firstbin>W<binwidth>[_<misc>]-<start>-<span>.sft\n"
253 " <outdir> = output directory given by the '-n' option, default is '.'\n"
254 " <nSFT> = number of SFTs in the file\n"
255 " <obs> = 1-character observatory prefix, e.g. 'H', 'L', 'V'\n"
256 " <det> = 2-character detector prefix, e.g. 'H1', 'L1', 'V1'\n"
257 " <timebase> = time base of the SFTs, rounded to the nearest second\n"
258 " <firstbin> = first bin of the SFTs, expressed in the form\n"
259 " <firstbin> = <freq>Hz<remainder>, where\n"
260 " <freq> = <firstbin> / <timebase>, rounded down to nearest Hz\n"
261 " <remainder> = <firstbin> - <freq> * <timebase>, number of remaining bins\n"
262 " <binwidth> = bin width of the SFTs, expressed in the form\n"
263 " <binwidth> = <freq>Hz<remainder>, where\n"
264 " <freq> = <binwidth> / <timebase>, rounded down to nearest Hz\n"
265 " <remainder> = <binwidth> - <freq> * <timebase>, number of remaining bins\n"
266 " <misc> = misc description field, if any, from the input files\n"
267 " <start> = start time of the first SFT as a GPS timestamp\n"
268 " <span> = time spanned by the SFT in seconds\n"
269 " If an output file already exists, the program will append the new SFTs to them, making\n"
270 " it possible to construct the final narrow-band SFTs by running the program multiple\n"
271 " times with different input SFTs. The GPS timestamps of the input SFTs need to be in\n"
272 " ascending order to get valid merged SFT files.\n"
274 " The '-c' options specifies how to deal with comments - 0 means no comment is written\n"
275 " at all, 1 means that the comment is taken unmodified from the input SFTs, 2 (default)\n"
276 " means that the program appends its RCS id and command-line to the comment.\n"
277 " By default a comment is written only to the first SFT of a merged SFT output 'block'\n"
278 " (i.e. call to this program). Adding the option '-a' to the command line specifies that\n"
279 " instead the comment is written into every SFT in the resulting file.\n"
281 " The options 'minStartTime' and 'maxStartTime' constrain input SFTs to be read\n"
282 " according to their time stamps. Following the XLALCWGPSinRange() convention,\n"
283 " these are interpreted as half-open intervals:\n"
284 " \tminStartTime <= time stamp < maxStartTime.\n"
285 " After encountering a time stamp outside of the given range, the program jumps to\n"
286 " the next input file as specified in the input argument; i.e. time stamp constraints\n"
287 " are imposed file-wise. For a multi-SFT file, headers are checked to ensure proper\n"
288 " sorting. Mind that all timestamps refer to SFT **start** times.\n"
290 " If 'assumeSorted' is set to 1, SFT input files will be assumed to be chronologically\n"
291 " sorted, which means the program will stop as soon as an SFT located after the\n"
292 " specified range is encountered.\n"
294 " After all options (and an optional '--' separator), the input files are given, as many\n"
295 " as you wish (or the OS supports - using xargs should be simple with this command-line\n"
298 " The program adds its own VCS ID and command-line to the comment of the written SFTs,\n"
299 " a mystery factor should show up as 'xxx' there.\n",
308 strcat( cmdline,
"\n" );
309 for ( arg = 0; arg < argc; arg++ ) {
310 if ( strcmp( argv[arg],
"-m" ) == 0 ) {
313 strcat( cmdline,
"-m xxx " );
317 strcat( cmdline, argv[arg] );
318 if ( arg == argc - 1 ) {
319 strcat( cmdline,
"\n" );
321 strcat( cmdline,
" " );
327 for ( arg = 1; arg < argc; arg++ ) {
328 if ( ( strcmp( argv[arg],
"-c" ) == 0 ) ||
329 ( strcmp( argv[arg],
"--add-comment" ) == 0 ) ) {
330 add_comment = atoi( argv[++arg] );
331 }
else if ( ( strcmp( argv[arg],
"-a" ) == 0 ) ||
332 ( strcmp( argv[arg],
"--all-comments" ) == 0 ) ) {
334 }
else if ( ( strcmp( argv[arg],
"-s" ) == 0 ) ||
335 ( strcmp( argv[arg],
"--start-bin" ) == 0 ) ) {
336 startBin = atoi( argv[++arg] );
337 }
else if ( ( strcmp( argv[arg],
"-e" ) == 0 ) ||
338 ( strcmp( argv[arg],
"--end-bin" ) == 0 ) ) {
339 endBin = atoi( argv[++arg] );
340 }
else if ( ( strcmp( argv[arg],
"-b" ) == 0 ) ||
341 ( strcmp( argv[arg],
"--width" ) == 0 ) ) {
342 width = atoi( argv[++arg] );
343 }
else if ( ( strcmp( argv[arg],
"-x" ) == 0 ) ||
344 ( strcmp( argv[arg],
"--overlap" ) == 0 ) ) {
345 overlap = atoi( argv[++arg] );
346 }
else if ( ( strcmp( argv[arg],
"-fs" ) == 0 ) ||
347 ( strcmp( argv[arg],
"--start-frequency" ) == 0 ) ) {
348 fMin = atof( argv[++arg] );
349 }
else if ( ( strcmp( argv[arg],
"-fe" ) == 0 ) ||
350 ( strcmp( argv[arg],
"--end-frequency" ) == 0 ) ) {
351 fMax = atof( argv[++arg] );
352 }
else if ( ( strcmp( argv[arg],
"-fb" ) == 0 ) ||
353 ( strcmp( argv[arg],
"--frequency-bandwidth" ) == 0 ) ) {
354 fWidth = atof( argv[++arg] );
355 }
else if ( ( strcmp( argv[arg],
"-fx" ) == 0 ) ||
356 ( strcmp( argv[arg],
"--frequency-overlap" ) == 0 ) ) {
357 fOverlap = atof( argv[++arg] );
358 }
else if ( ( strcmp( argv[arg],
"-ts" ) == 0 ) ||
359 ( strcmp( argv[arg],
"--minStartTime" ) == 0 ) ) {
361 minStartTime = &userMinStartTime;
362 }
else if ( ( strcmp( argv[arg],
"-te" ) == 0 ) ||
363 ( strcmp( argv[arg],
"--maxStartTime" ) == 0 ) ) {
365 maxStartTime = &userMaxStartTime;
366 }
else if ( ( strcmp( argv[arg],
"-as" ) == 0 ) ||
367 ( strcmp( argv[arg],
"--assumeSorted" ) == 0 ) ) {
368 assumeSorted = atoi( argv[++arg] );
369 }
else if ( ( strcmp( argv[arg],
"-m" ) == 0 ) ||
370 ( strcmp( argv[arg],
"--factor" ) == 0 ) ) {
371 factor = atof( argv[++arg] );
372 }
else if ( ( strcmp( argv[arg],
"-v" ) == 0 ) ||
373 ( strcmp( argv[arg],
"--no-validation" ) == 0 ) ) {
375 }
else if ( ( strcmp( argv[arg],
"-n" ) == 0 ) ||
376 ( strcmp( argv[arg],
"--output-directory" ) == 0 ) ) {
378 }
else if ( ( strcmp( argv[arg],
"-rb" ) == 0 ) ||
379 ( strcmp( argv[arg],
"--read-bandwidth" ) == 0 ) ) {
381 }
else if ( ( strcmp( argv[arg],
"-ror" ) == 0 ) ||
382 ( strcmp( argv[arg],
"--read-open-rate" ) == 0 ) ) {
384 }
else if ( ( strcmp( argv[arg],
"-wb" ) == 0 ) ||
385 ( strcmp( argv[arg],
"--write-bandwidth" ) == 0 ) ) {
387 }
else if ( ( strcmp( argv[arg],
"-wor" ) == 0 ) ||
388 ( strcmp( argv[arg],
"--write-open-rate" ) == 0 ) ) {
390 }
else if ( strcmp( argv[arg],
"--" ) == 0 ) {
393 }
else if ( strncmp( argv[arg],
"-", 1 ) == 0 ) {
394 fprintf( stderr,
"unknown option '%s', try '-h' for help\n", argv[arg] );
402 char *constraint_str;
403 const size_t constraint_str_len = 32;
405 XLAL_ENOMEM,
"out of memory allocating constraint_str" );
406 if ( minStartTime ) {
407 XLAL_CHECK_MAIN( snprintf( constraint_str, constraint_str_len,
"%d,", minStartTime->
gpsSeconds ) < (
int )constraint_str_len,
408 XLAL_ESYS,
"output of snprintf() was truncated" );
410 XLAL_CHECK_MAIN( snprintf( constraint_str, constraint_str_len,
"," ) < (
int )constraint_str_len,
411 XLAL_ESYS,
"output of snprintf() was truncated" );
413 if ( maxStartTime ) {
414 XLAL_CHECK_MAIN( snprintf( constraint_str + strlen( constraint_str ), constraint_str_len - strlen( constraint_str ),
"%d", maxStartTime->
gpsSeconds ) < (
int )( constraint_str_len - strlen( constraint_str ) ),
415 XLAL_ESYS,
"output of snprintf() was truncated" );
431 glob( pattern, GLOB_MARK | GLOB_NOSORT, NULL, &globbuf );
432 for (
size_t i = 0;
i < globbuf.gl_pathc; ++
i ) {
449 XLALPrintInfo(
"Existing SFT: det=%c%c timebase=%u firstbin=(%u,%u) binwidth=(%u,%u) filename=%s\n",
452 globbuf.gl_pathv[
i] );
456 globfree( &globbuf );
460 int stopInputCauseSorted = 0;
462 for ( ; ( arg < argc ) && !stopInputCauseSorted; arg++ ) {
470 XLAL_CHECK_MAIN( ( fpin = fopen( argv[arg],
"r" ) ) != NULL,
XLAL_EIO,
"could not open SFT file for reading" );
473 int firstread =
TRUE;
474 int nSFT_this_file = 0;
480 sfterrno =
ReadSFTHeader( fpin, &hd, &oldcomment, &swap, validate );
497 XLAL_CHECK_MAIN( !( firstread && ( inRange == 1 ) ),
XLAL_EIO,
"First timestamp %d in file was after user constraint [%s)!", hd.
gps_sec, constraint_str );
499 if ( inRange == -1 ) {
503 if ( inRange == 1 ) {
504 stopInputCauseSorted = assumeSorted;
519 if ( fWidth >= 0.0 ) {
522 if ( fOverlap >= 0.0 ) {
533 "ERROR: start bin (%d) is smaller than first bin in input SFT (%d)\n",
539 if ( startBin + width > endBin + 1 ) {
541 "ERROR: end bin (%d) is larger than last bin in input SFT (%d)\n",
547 if ( overlap >= width ) {
549 "ERROR: overlap (%d) is not smaller than the width (%d)\n",
568 }
else if ( add_comment ==
CMT_OLD ) {
580 for ( nactivesamples = 0; nactivesamples < endBin - startBin; nactivesamples += width - overlap );
582 nactivesamples += overlap + 1;
590 sfterrno =
ReadSFTData( fpin,
data, startBin, nactivesamples, NULL, NULL );
594 for ( bin = 0; bin < 2 * nactivesamples; bin++ ) {
599 for ( bin = startBin; bin < endBin; bin += width - overlap ) {
605 int last_output_bin = bin + width - 1;
606 int max_input_width = last_output_bin <= last_input_bin ? width : width - ( last_output_bin - last_input_bin );
607 int max_output_width = endBin - bin + 1;
608 int this_width = max_input_width < max_output_width ? max_input_width : max_output_width;
610 int timebase = (
int ) round( hd.
tbase );
611 int outfreq = (
int )floor( bin / timebase );
612 int outfreqbin = bin - outfreq * timebase;
613 int outwidth = (
int )floor( this_width / timebase );
614 int outwidthbin = this_width - outwidth * timebase;
619 .SFTtimebase = timebase,
620 .nbFirstBinFreq = outfreq,
621 .nbFirstBinRem = outfreqbin,
622 .nbBinWidthFreq = outwidth,
623 .nbBinWidthRem = outwidthbin,
625 XLALPrintInfo(
"Looking for SFT:det=%c%c timebase=%u firstbin=(%u,%u) binwidth=(%u,%u)\n",
626 key.detector[0],
key.detector[1],
key.SFTtimebase,
627 key.nbFirstBinFreq,
key.nbFirstBinRem,
key.nbBinWidthFreq,
key.nbBinWidthRem );
630 char *existingSFTfilename = NULL;
662 "New SFT timestamp %d is before startTime+span-timebase=%u+%u-%u=%u from existing file %s. Appending would yield an invalid merged SFT file.",
679 if ( existingSFTfilename != NULL && strcmp( existingSFTfilename, newSFTfilename ) != 0 ) {
680 XLALPrintInfo(
"Renaming SFT '%s' to '%s'\n", existingSFTfilename, newSFTfilename );
681 rename( existingSFTfilename, newSFTfilename );
685 XLALPrintInfo(
"New/updated SFT: det=%c%c timebase=%u firstbin=(%u,%u) binwidth=(%u,%u) filename=%s\n",
693 XLAL_CHECK_MAIN( ( fpout = fopen( newSFTfilename,
"a" ) ) != NULL,
XLAL_EIO,
"could not open SFT for writing" );
void LALCheckMemoryLeaks(void)
const LALVCSInfo lalPulsarVCSIdentInfo
Identable VCS and build information for LALPulsar.
#define XLAL_INIT_DECL(var,...)
int XLALPearsonHash(void *hval, const size_t hval_len, const void *data, const size_t data_len)
int XLALHashTblAdd(LALHashTbl *ht, void *x)
int XLALHashTblExtract(LALHashTbl *ht, const void *x, void **y)
LALHashTbl * XLALHashTblCreate(LALHashTblDtorFcn dtor, LALHashTblHashFcn hash, LALHashTblCmpFcn cmp)
void XLALHashTblDestroy(LALHashTbl *ht)
void * XLALMalloc(size_t n)
void * XLALCalloc(size_t m, size_t n)
void * XLALRealloc(void *p, size_t n)
char * XLALStringAppendFmt(char *s, const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(2
const char * SFTErrorMessage(int errorcode)
int CheckSFTHeaderConsistency(struct headertag2 *headerone, struct headertag2 *headertwo)
int ReadSFTHeader(FILE *fp, struct headertag2 *info, char **comment, int *swapendian, int validate)
int WriteSFT(FILE *fp, int gps_sec, int gps_nsec, double tbase, int firstfreqindex, int nsamples, const char *detector, unsigned short windowspec, const char *comment, float *data)
int ReadSFTData(FILE *fp, float *data, int firstbin, int nsamples, char **comment, struct headertag2 *info)
char * XLALBuildSFTFilenameFromSpec(const SFTFilenameSpec *spec)
Build an SFT file name from the given specification.
UINT4 XLALRoundFrequencyUpToSFTBin(const REAL8 freq, const REAL8 df)
Round a REAL8 frequency up to the nearest integer SFT bin number.
UINT4 XLALRoundFrequencyDownToSFTBin(const REAL8 freq, const REAL8 df)
Round a REAL8 frequency down to the nearest integer SFT bin number.
int XLALCWGPSinRange(const LIGOTimeGPS gps, const LIGOTimeGPS *minGPS, const LIGOTimeGPS *maxGPS)
Defines the official CW convention for whether a GPS time is 'within' a given range,...
int XLALParseSFTFilenameIntoSpec(SFTFilenameSpec *spec, const char *SFTpath)
Parse a SFT file path and return its specification.
int XLALFillSFTFilenameSpecStrings(SFTFilenameSpec *spec, const CHAR *path, const CHAR *extn, const CHAR *detector, const CHAR *window_type, const CHAR *privMisc, const CHAR *pubObsKind, const CHAR *pubChannel)
Convenience function for filling out the string fields in a SFTFilenameSpec.
int int int XLALPrintInfo(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_CHECK(assertion,...)
#define XLAL_CHECK_MAIN(assertion,...)
LIGOTimeGPS * XLALGPSSetREAL8(LIGOTimeGPS *epoch, REAL8 t)
static UINT8 hash_SFTFilenameSpec(const void *x)
int main(int argc, char **argv)
main program
static int is_directory(const char *fname)
UNIT_SOURCE read_bandwidth
throtteling settings
static void request_resource(UNIT_SOURCE *us, int units)
request a resurce.
UNIT_SOURCE write_open_rate
static int compare_SFTFilenameSpec(const void *x, const void *y)
UNIT_SOURCE write_bandwidth
static int move_to_next_SFT(FILE *fpin, int nsamples, int comment_length)
UNIT_SOURCE read_open_rate
const char *const vcsStatus
Structure specifying an SFT file name, following the convention in .
UINT4 SFTtimebase
Timebase in seconds of the SFT; set by XLALWriteSFT[Vector]2StandardFile()
UINT4 SFTspan
Total time interval in seconds covered by SFT file; set by XLALWriteSFT[Vector]2StandardFile()
UINT4 nbBinWidthRem
For narrow-band SFTs: remainder of division of SFT bandwidth by SFT time base.
UINT4 gpsStart
GPS time in seconds at the beginning of the first SFT in the file; set by XLALWriteSFT[Vector]2Standa...
UINT4 nbBinWidthFreq
For narrow-band SFTs: SFT bandwidth divided by SFT time base, rounded down.
CHAR detector[3]
2-character detector prefix (e.g.
UINT4 nbFirstBinRem
For narrow-band SFTs: remainder of division of SFT first bin frequency by SFT time base.
UINT4 numSFTs
Number of SFTs in the file; set by XLALWriteSFT[Vector]2StandardFile()
UINT4 nbFirstBinFreq
For narrow-band SFTs: SFT first bin frequency divided by SFT time base, rounded down.
int resource_rate
this many resource units become available in one second
int resource_units
number of resource units available for consumption
time_t last_checked
time we last checked