29 #include <gsl/gsl_math.h>
30 #include <gsl/gsl_sort_double.h>
32 #include <lal/Units.h>
33 #include <lal/FrequencySeries.h>
34 #include <lal/NormalizeSFTRngMed.h>
36 #include <lal/PSDutils.h>
83 if (
psd->data->data ) {
104 if ( multvect == NULL ) {
130 multiWeights->
length = length;
157 for (
UINT4 X = 0; X < multiWeights->
length; X++ ) {
165 return copiedWeights;
202 if ( !multiPSDVect ) {
206 if ( !multiSFTVect ) {
210 if ( lastBin < firstBin ) {
218 if ( numIFOs != multiSFTVect->
length ) {
234 if ( ( firstBin == 0 ) && ( lastBin ==
numBins - 1 ) ) {
240 for ( X = 0; X < numIFOs; X ++ ) {
246 for ( iTS = 0; iTS < numTS; iTS ++ ) {
251 XLALPrintError(
"%s: inconsistent number of frequency-bins across multiPSDVector: X=%d, iTS=%d: numBins = %d != %d\n",
257 XLALPrintError(
"%s: inconsistent number of frequency-bins across multiSFTVector: X=%d, iTS=%d: numBins = %d != %d\n",
262 UINT4 numNewBins = lastBin - firstBin + 1;
287 UINT4 excludePercentile )
300 UINT4 numSFTsTot = 0;
301 REAL8 sumWeights = 0;
303 for (
UINT4 X = 0; X < numIFOs; X++ ) {
316 UINT4 halfBlock = blocksRngMed / 2;
324 UINT4 length = lengthsft - blocksRngMed + 1;
325 UINT4 halfLength = length / 2;
328 UINT4 excludeIndex = excludePercentile * halfLength ;
331 REAL8 Tsft_avgS2 = 0.0;
332 for (
UINT4 k = halfBlock + excludeIndex;
k < lengthsft - halfBlock - excludeIndex;
k++ ) {
335 Tsft_avgS2 /= lengthsft - 2 * halfBlock - 2 * excludeIndex;
337 REAL8 wXa = 1.0 / Tsft_avgS2;
349 REAL8 TsftS2_inv = sumWeights / numSFTsTot;
352 for (
UINT4 X = 0; X < numIFOs; X ++ ) {
380 if ( multiPSDVect == NULL ) {
384 if ( multiPSDVect->
length == 0 || multiPSDVect->
data == 0 ) {
410 memset(
Q->data->data, 0,
Q->data->length *
sizeof(
Q->data->data[0] ) );
415 for ( X = 0; X < numIFOs; X ++ ) {
420 UINT4 numSFTsInSeg = 0;
426 for ( iTS = 0; iTS < numTS; iTS++ ) {
430 if ( (
f0 != thisPSD->
f0 ) || ( dFreq != thisPSD->
deltaF ) || ( numFreqs != thisPSD->
data->
length ) ) {
431 XLALPrintError(
"%s: %d-th timestamp %f: inconsistent PSDVector: f0 = %g : %g, dFreq = %g : %g, numFreqs = %d : %d \n",
449 for ( iFreq = 0; iFreq < numFreqs; iFreq++ ) {
458 REAL8 duty_X = numSFTsInSeg *
Tsft / Tseg;
460 if ( ( duty_X < 0 ) && ( duty_X > 1 ) ) {
467 for ( iFreq = 0; iFreq < numFreqs; iFreq ++ ) {
468 Q->data->data[iFreq] += duty_X * SXinv->
data->
data[iFreq] / numSFTsInSeg;
504 for (
i = 0;
i < length; ++
i ) {
512 for (
i = 0;
i < length; ++
i ) {
515 res /= (
REAL8 )length;
521 for (
i = 0;
i < length; ++
i ) {
522 res += 1.0 / *(
data++ );
530 for (
i = 0;
i < length; ++
i ) {
531 res += 1.0 / *(
data++ );
533 res = (
REAL8 )length / res;
539 for (
i = 0;
i < length; ++
i ) {
542 res = 1.0 / sqrt( res );
548 for (
i = 0;
i < length; ++
i ) {
551 res = 1.0 / sqrt( res / (
REAL8 )length );
564 if ( ( sortidx =
XLALMalloc( length *
sizeof(
size_t ) ) ) == NULL ) {
568 gsl_sort_index( sortidx,
data, 1, length );
574 if ( length / 2 == ( length + 1 ) / 2 ) {
575 res = (
data[sortidx[length / 2 - 1]] +
data[sortidx[length / 2]] ) / 2;
577 res =
data[sortidx[length / 2]];
584 res =
data[sortidx[0]];
589 res =
data[sortidx[length - 1]];
628 const UINT4 totalNumSFTs,
635 switch ( optypeSFTs ) {
648 factor = totalNumSFTs;
652 factor = sqrt( totalNumSFTs );
682 const BOOLEAN returnMultiPSDVector,
684 const UINT4 blocksRngMed,
689 const BOOLEAN normalizeByTotalNumSFTs,
691 const REAL8 FreqBand,
692 const BOOLEAN normalizeSFTsInPlace
701 XLAL_CHECK( ( FreqMin >= 0 && FreqBand >= 0 ) || ( FreqMin == -1 && FreqBand == 0 ),
XLAL_EINVAL,
"Need either both Freqmin>=0 && FreqBand>=0 (truncate PSD band to this range) or FreqMin=-1 && FreqBand=0 (use full band as loaded from SFTs, including rngmed-sidebands." );
712 if ( normalizeSFTsInPlace ) {
716 for (
UINT4 X = 0; X < numIFOs; ++X ) {
720 for (
UINT4 X = 0; X < numIFOs; ++X ) {
728 REAL8 dFreq = ( *multiPSDVector )->data[0]->data[0].deltaF;
732 UINT4 firstBin, lastBin;
735 REAL8 fmaxSFTs = fminSFTs + numBinsSFTs * dFreq;
736 XLAL_CHECK( ( FreqMin >= fminSFTs ) && ( FreqMin + FreqBand <= fmaxSFTs ),
XLAL_EDOM,
"Requested frequency range [%f,%f]Hz not covered by available SFT band [%f,%f]Hz", FreqMin, FreqMin + FreqBand, fminSFTs, fmaxSFTs );
738 UINT4 rngmedSideBandBins = blocksRngMed / 2 + 1;
739 REAL8 rngmedSideBand = rngmedSideBandBins * dFreq;
745 if ( ( minBinPSD < minBinSFTs ) || ( maxBinPSD > maxBinSFTs ) ) {
746 XLAL_ERROR(
XLAL_ERANGE,
"Requested frequency range [%f,%f)Hz means rngmedSideBand=%fHz (%d bins) would spill over available SFT band [%f,%f)Hz",
747 FreqMin, FreqMin + FreqBand, rngmedSideBand, rngmedSideBandBins, fminSFTs, fmaxSFTs );
749 UINT4 binsBand = maxBinPSD - minBinPSD + 1;
750 firstBin = minBinPSD - minBinSFTs;
751 lastBin = firstBin + binsBand - 1;
754 lastBin = numBinsSFTs - 1;
756 XLAL_PRINT_INFO(
"loaded SFTs have %d bins, effective PSD output band is [%d, %d]", numBinsSFTs, firstBin, lastBin );
758 XLAL_CHECK(
XLALCropMultiPSDandSFTVectors( *multiPSDVector, SFTs, firstBin, lastBin ) ==
XLAL_SUCCESS,
XLAL_EFUNC,
"Failed call to XLALCropMultiPSDandSFTVectors (multiPSDVector, SFTs, firstBin=%d, lastBin=%d)", firstBin, lastBin );
761 UINT4 numBins = ( *multiPSDVector )->data[0]->data[0].data->length;
771 UINT4 maxNumSFTs = 0;
772 for (
UINT4 X = 0; X < numIFOs; ++X ) {
773 maxNumSFTs = GSL_MAX( maxNumSFTs, ( *multiPSDVector )->data[X]->length );
780 REAL8 normPSD = 2.0 * dFreq;
788 REAL8 totalNumSFTsNormalizingFactor = 1.;
789 if ( normalizeByTotalNumSFTs ) {
790 UINT4 totalNumSFTs = 0;
791 for (
UINT4 X = 0; X < numIFOs; ++X ) {
792 totalNumSFTs += ( *multiPSDVector )->data[X]->length;
802 for (
UINT4 X = 0; X < numIFOs; ++X ) {
809 ( *multiPSDVector )->data[X]->data[
alpha].data->data[
k] *= normPSD;
810 overSFTs->
data[
alpha] = ( *multiPSDVector )->data[X]->data[
alpha].data->data[
k];
823 if ( normalizeByTotalNumSFTs ) {
824 ( *finalPSD )->data[
k] *= totalNumSFTsNormalizingFactor;
831 if ( returnNormSFT ) {
839 for (
UINT4 X = 0; X < numIFOs; ++X ) {
847 overSFTs->
data[
alpha] = crealf( bin ) * crealf( bin ) + cimagf( bin ) * cimagf( bin );
867 if ( !returnMultiPSDVector ) {
869 *multiPSDVector = NULL;
872 if ( !normalizeSFTsInPlace ) {
897 const UINT4 blocksRngMed,
900 const BOOLEAN normalizeByTotalNumSFTs,
919 normalizeByTotalNumSFTs,
941 if ( outbname == NULL ) {
945 if ( multiPSDVect == NULL ) {
949 if ( multiPSDVect->
length == 0 || multiPSDVect->
data == 0 ) {
957 UINT4 len = strlen( outbname ) + 4;
965 for ( X = 0; X < numIFOs; X ++ ) {
969 sprintf( fname,
"%s-%c%c", outbname, thisPSDVect->
data[0].
name[0], thisPSDVect->
data[0].
name[1] );
971 if ( (
fp = fopen( fname,
"wb" ) ) == NULL ) {
984 fprintf(
fp,
"%%%% first line holds frequencies [Hz] of PSD-Columns\n" );
987 for ( iFreq = 0; iFreq < numFreqs; iFreq ++ ) {
988 sprintf( buf,
"f%d [Hz]", iFreq + 1 );
995 for ( iFreq = 0; iFreq < numFreqs; iFreq++ ) {
1002 for ( iFreq = 0; iFreq < numFreqs; iFreq ++ ) {
1003 sprintf( buf,
"PSD(f%d)", iFreq + 1 );
1011 for ( iTS = 0; iTS < numTS; iTS++ ) {
1019 if ( (
f0 != thisPSD->
f0 ) || ( dFreq != thisPSD->
deltaF ) || ( numFreqs != thisPSD->
data->
length ) ) {
1020 XLALPrintError(
"%s: %d-th timestamp %f: inconsistent PSDVector: f0 = %g : %g, dFreq = %g : %g, numFreqs = %d : %d \n",
1028 for ( iFreq = 0; iFreq < numFreqs; iFreq ++ ) {
1072 XLAL_CHECK( ( PSDVect != NULL ) && ( PSDVect->
data != NULL ) && ( PSDVect->
length > 0 ),
XLAL_EINVAL,
"PSDVect must be allocated and not zero length." );
1073 if ( outputNormSFT ) {
1074 XLAL_CHECK( ( normSFTVect != NULL ) && ( normSFTVect->
data != NULL ),
XLAL_EINVAL,
"If requesting outputNormSFT, normSFTVect must be allocated" );
1086 fprintf( fpOut,
"%%%% columns:\n%%%% FreqBinStart" );
1087 if ( outFreqBinEnd ) {
1088 fprintf( fpOut,
" FreqBinEnd" );
1091 if ( outputNormSFT ) {
1092 fprintf( fpOut,
" normSFTpower" );
1102 if ( outFreqBinEnd ) {
1111 if ( outputNormSFT ) {
1115 fprintf( fpOut,
" %f", nsft );
#define __func__
log an I/O error, i.e.
Internal SFT types and functions.
REAL8FrequencySeries * XLALCreateREAL8FrequencySeries(const CHAR *name, const LIGOTimeGPS *epoch, REAL8 f0, REAL8 deltaF, const LALUnit *sampleUnits, size_t length)
REAL8FrequencySeries * XLALResizeREAL8FrequencySeries(REAL8FrequencySeries *series, int first, size_t length)
COMPLEX8FrequencySeries * XLALResizeCOMPLEX8FrequencySeries(COMPLEX8FrequencySeries *series, int first, size_t length)
void XLALDestroyREAL8FrequencySeries(REAL8FrequencySeries *series)
void * XLALCalloc(size_t m, size_t n)
void * XLALMalloc(size_t n)
MultiPSDVector * XLALNormalizeMultiSFTVect(MultiSFTVector *multsft, UINT4 blockSize, const MultiNoiseFloor *assumeSqrtSX)
Function for normalizing a multi vector of SFTs in a multi IFO search and returns the running-median ...
int XLALDumpMultiPSDVector(const CHAR *outbname, const MultiPSDVector *multiPSDVect)
Dump complete multi-PSDVector over IFOs, timestamps and frequency-bins into per-IFO ASCII output-file...
int XLALWritePSDtoFilePointer(FILE *fpOut, REAL8Vector *PSDVect, REAL8Vector *normSFTVect, BOOLEAN outputNormSFT, BOOLEAN outFreqBinEnd, INT4 PSDmthopBins, INT4 nSFTmthopBins, INT4 binSize, INT4 binStep, REAL8 Freq0, REAL8 dFreq)
Write a PSD vector as an ASCII table to an open output file pointer.
int XLALCropMultiPSDandSFTVectors(MultiPSDVector *multiPSDVect, MultiSFTVector *multiSFTVect, UINT4 firstBin, UINT4 lastBin)
Function that truncates the PSD in place to the requested frequency-bin interval [firstBin,...
MultiNoiseWeights * XLALCopyMultiNoiseWeights(const MultiNoiseWeights *multiWeights)
Create a full copy of a MultiNoiseWeights object.
const UserChoices MathOpTypeChoices
int XLALComputePSDfromSFTs(REAL8Vector **finalPSD, MultiSFTVector *inputSFTs, const UINT4 blocksRngMed, const MathOpType PSDmthopSFTs, const MathOpType PSDmthopIFOs, const BOOLEAN normalizeByTotalNumSFTs, const REAL8 FreqMin, const REAL8 FreqBand)
Compute the PSD (power spectral density) over a MultiSFTVector.
void XLALDestroyPSDVector(PSDVector *vect)
Destroy a PSD-vector.
void XLALDestroyMultiPSDVector(MultiPSDVector *multvect)
Destroy a multi PSD-vector.
int XLALComputePSDandNormSFTPower(REAL8Vector **finalPSD, MultiPSDVector **multiPSDVector, REAL8Vector **normSFT, MultiSFTVector *inputSFTs, const BOOLEAN returnMultiPSDVector, const BOOLEAN returnNormSFT, const UINT4 blocksRngMed, const MathOpType PSDmthopSFTs, const MathOpType PSDmthopIFOs, const MathOpType nSFTmthopSFTs, const MathOpType nSFTmthopIFOs, const BOOLEAN normalizeByTotalNumSFTs, const REAL8 FreqMin, const REAL8 FreqBand, const BOOLEAN normalizeSFTsInPlace)
Compute the PSD (power spectral density) and the "normalized power" over a MultiSFTVector.
MultiNoiseWeights * XLALComputeMultiNoiseWeights(const MultiPSDVector *rngmed, UINT4 blocksRngMed, UINT4 excludePercentile)
Computes weight factors arising from MultiSFTs with different noise floors.
void XLALDestroyMultiNoiseWeights(MultiNoiseWeights *weights)
Destroy a MultiNoiseWeights object.
MultiNoiseWeights * XLALCreateMultiNoiseWeights(const UINT4 length)
Create a MultiNoiseWeights of the given length.
MathOpType
common types of mathematical operations over an array
REAL8 XLALMathOpOverArray(const REAL8 *data, const size_t length, const MathOpType optype)
Compute various types of "math operations" over the entries of an array.
REAL8FrequencySeries * XLALComputeSegmentDataQ(const MultiPSDVector *multiPSDVect, LALSeg segment)
Compute the "data-quality factor" over the given SFTs.
REAL8 XLALGetMathOpNormalizationFactorFromTotalNumberOfSFTs(const UINT4 totalNumSFTs, const MathOpType optypeSFTs)
Compute an additional normalization factor from the total number of SFTs to be applied after a loop o...
@ MATH_OP_POWERMINUS2_SUM
@ MATH_OP_ARITHMETIC_MEDIAN
@ MATH_OP_POWERMINUS2_MEAN
@ MATH_OP_ARITHMETIC_MEAN
void XLALDestroyMultiSFTVector(MultiSFTVector *multvect)
Destroy a multi SFT-vector.
REAL8 TSFTfromDFreq(REAL8 dFreq)
int XLALCopySFT(SFTtype *dest, const SFTtype *src)
Copy an entire SFT-type into another.
UINT4 XLALRoundFrequencyUpToSFTBin(const REAL8 freq, const REAL8 df)
Round a REAL8 frequency up to the nearest integer SFT bin number.
MultiSFTVector * XLALCreateEmptyMultiSFTVector(UINT4Vector *numsft)
Create an empty multi-IFO SFT vector with a given number of SFTs per IFO (which are not allocated).
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,...
const LALUnit lalHertzUnit
void XLALDestroyUINT4Vector(UINT4Vector *vector)
UINT4Vector * XLALCreateUINT4Vector(UINT4 length)
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
void XLALDestroyREAL8Vector(REAL8Vector *vector)
#define XLAL_ERROR_REAL8(...)
#define XLAL_PRINT_INFO(...)
#define XLAL_ERROR_NULL(...)
#define XLAL_CHECK(assertion,...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
#define XLAL_CHECK_NULL(assertion,...)
#define XLAL_IS_REAL8_FAIL_NAN(val)
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.
One noise-weight (number) per SFT (therefore indexed over IFOs and SFTs.
REAL8 Sinv_Tsft
normalization factor used: (using single-sided PSD!)
UINT4 length
number of detectors
REAL8Vector ** data
weights-vector for each detector
BOOLEAN isNotNormalized
if true: weights are saved unnormalized (divide by Sinv_Tsft to get normalized version).
A collection of PSD vectors – one for each IFO in a multi-IFO search.
PSDVector ** data
sftvector for each ifo
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
A vector of REAL8FrequencySeries.
REAL8FrequencySeries * data
Pointer to the data array.
UINT4 length
Number of elements in array.