21 #include <lal/SFTClean.h>
22 #include <lal/SFTfileIO.h>
23 #include <lal/LogPrintf.h>
99 INT4 harmonicCount,
r, tempint;
100 REAL8 temp1, temp2, temp3, temp4;
117 fp = fopen( fname,
"r" );
124 r =
fscanf(
fp,
"%lf%lf%d%lf%lf%s\n", &temp1, &temp2,
125 &tempint, &temp3, &temp4, dump );
131 }
while (
r != EOF );
157 REAL8 *startFreq = NULL;
158 REAL8 *gapFreq = NULL;
159 INT4 *numHarmonics = NULL;
160 REAL8 *leftWing = NULL;
161 REAL8 *rightWing = NULL;
178 fp = fopen( fname,
"r" );
183 gapFreq = harmonicsInfo->
gapFreq;
217 INT4 count1, count2, nHarmonicSets, maxCount, position;
223 REAL8 f0, deltaf, leftDeltaf, rightDeltaf;
246 gapFreq = harmonicsInfo->
gapFreq;
252 for ( count1 = 0; count1 < nHarmonicSets; count1++ ) {
253 maxCount = *( numHarmonics + count1 );
254 f0 = *( startFreq + count1 );
255 deltaf = *( gapFreq + count1 );
256 leftDeltaf = *( leftWing + count1 );
257 rightDeltaf = *( rightWing + count1 );
258 for ( count2 = 0; count2 < maxCount; count2++ ) {
259 *( lineInfo->
lineFreq + count2 + position ) =
f0 + count2 * deltaf;
260 *( lineInfo->
leftWing + count2 + position ) = leftDeltaf;
261 *( lineInfo->
rightWing + count2 + position ) = rightDeltaf;
263 position += maxCount;
290 REAL8 temp1, temp2, temp3;
302 fp = fopen( fname,
"r" );
307 r =
fscanf(
fp,
"%lf%lf%lf%s\n", &temp1, &temp2, &temp3, dump );
313 }
while (
r != EOF );
315 lineInfo->
nLines = lineCount;
337 REAL8 *lineFreq = NULL;
338 REAL8 *leftWing = NULL;
339 REAL8 *rightWing = NULL;
354 fp = fopen( fname,
"r" );
357 nLines = lineInfo->
nLines;
394 INT4 nLinesIn, nLinesOut,
j;
395 REAL8 lineFreq, leftWing, rightWing;
409 nLinesIn = inLine->
nLines;
412 for (
j = 0;
j < nLinesIn;
j++ ) {
416 if ( ( lineFreq >=
fMin ) && ( lineFreq <=
fMax ) ) {
417 outLine->
lineFreq[nLinesOut] = lineFreq;
418 outLine->
leftWing[nLinesOut] = leftWing;
419 outLine->
rightWing[nLinesOut] = rightWing;
421 }
else if ( ( lineFreq <
fMin ) && ( lineFreq + rightWing >
fMin ) ) {
422 outLine->
lineFreq[nLinesOut] = lineFreq;
423 outLine->
leftWing[nLinesOut] = leftWing;
424 outLine->
rightWing[nLinesOut] = rightWing;
426 }
else if ( ( lineFreq >
fMax ) && ( lineFreq - leftWing <
fMax ) ) {
427 outLine->
lineFreq[nLinesOut] = lineFreq;
428 outLine->
leftWing[nLinesOut] = leftWing;
429 outLine->
rightWing[nLinesOut] = rightWing;
435 if ( nLinesOut == 0 ) {
436 outLine->
nLines = nLinesOut;
441 outLine->
nLines = nLinesOut;
471 REAL8 lineFreq, leftWing, rightWing, doppler;
484 for (
j = 0;
j < nLines;
j++ ) {
489 doppler =
VTOT * ( lineFreq + rightWing );
491 rightWing += doppler;
493 if ( (
freq <= lineFreq + rightWing ) && (
freq >= lineFreq - leftWing ) ) {
591 INT4 nLines,
count, leftCount, rightCount, lineBin, minBin, maxBin,
k, tempk;
592 INT4 leftWingBins, rightWingBins, length, sumBins;
594 REAL8 stdPow, medianPow;
595 REAL8 *tempDataPow = NULL;
596 REAL8 *lineFreq = NULL;
597 REAL8 *leftWing = NULL;
598 REAL8 *rightWing = NULL;
627 nLines = lineInfo->
nLines;
636 maxBin = minBin + length - 1;
655 tempSumBins = floor( tBase * leftWing[
count] ) + floor( tBase * rightWing[
count] );
656 sumBins += tempSumBins < 2 * width ? tempSumBins : 2 * width;
670 lineBin = lround( tBase * lineFreq[
count] );
671 leftWingBins = floor( tBase * leftWing[
count] );
672 rightWingBins = floor( tBase * rightWing[
count] );
675 if ( ( lineBin >= minBin ) && ( lineBin <= maxBin ) ) {
678 if ( ( leftWingBins > width ) || ( rightWingBins > width ) ) {
681 leftWingBins = leftWingBins < width ? leftWingBins : width;
682 rightWingBins = rightWingBins < width ? rightWingBins : width;
685 if ( 2 * window > ( maxBin - minBin ) ) {
686 LogPrintf(
LOG_NORMAL,
"%s: Window of +-%d bins around lineBin=%d does not fit within SFT range [%d,%d], noise floor estimation will be compromised.\n",
__func__, window, lineBin, minBin, maxBin );
688 if ( ( window < leftWingBins ) || ( window <= rightWingBins ) ) {
689 LogPrintf(
LOG_NORMAL,
"%s: Window of +-%d bins around lineBin=%d does not extend further than line wings [-%d,+%d], noise floor estimation will be compromised by the very same line that is supposed to be cleaned.\n",
__func__, window, lineBin, leftWingBins, rightWingBins );
691 for (
k = 0;
k < window ;
k++ ) {
692 if ( maxBin - lineBin - rightWingBins -
k > 0 ) {
693 inData = sft->
data->
data + lineBin - minBin + rightWingBins +
k + 1;
695 inData = sft->
data->
data + length - 1;
698 tempDataPow[
k] = ( crealf( *inData ) ) * ( crealf( *inData ) ) + ( cimagf( *inData ) ) * ( cimagf( *inData ) );
700 if ( lineBin - minBin - leftWingBins -
k > 0 ) {
701 inData = sft->
data->
data + lineBin - minBin - leftWingBins -
k - 1;
706 tempDataPow[
k + window] = ( crealf( *inData ) ) * ( crealf( *inData ) ) + ( cimagf( *inData ) ) * ( cimagf( *inData ) );
709 gsl_sort( tempDataPow, 1, 2 * window );
710 medianPow = gsl_stats_median_from_sorted_data( tempDataPow, 1, 2 * window );
711 stdPow = sqrt( medianPow / ( 2 * bias ) );
714 inData = sft->
data->
data + lineBin - minBin;
716 randVal = ranVector->
data + tempk;
717 *( inData ) =
crectf( stdPow * ( *randVal ), cimagf( *( inData ) ) );
720 randVal = ranVector->
data + tempk;
721 *( inData ) =
crectf( crealf( *( inData ) ), stdPow * ( *randVal ) );
727 for ( leftCount = 0; leftCount < leftWingBins; leftCount++ ) {
728 if ( ( lineBin - minBin - leftCount > 0 ) ) {
729 inData = sft->
data->
data + lineBin - minBin - leftCount - 1;
731 randVal = ranVector->
data + tempk;
732 *( inData ) =
crectf( stdPow * ( *randVal ), cimagf( *( inData ) ) );
735 randVal = ranVector->
data + tempk;
736 *( inData ) =
crectf( crealf( *( inData ) ), stdPow * ( *randVal ) );
742 for ( rightCount = 0; rightCount < rightWingBins; rightCount++ ) {
743 if ( ( maxBin - lineBin - rightCount > 0 ) ) {
744 inData = sft->
data->
data + lineBin - minBin + rightCount + 1;
746 randVal = ranVector->
data + tempk;
747 *( inData ) =
crectf( stdPow * ( *randVal ), cimagf( *( inData ) ) );
750 randVal = ranVector->
data + tempk;
751 *( inData ) =
crectf( crealf( *( inData ) ), stdPow * ( *randVal ) );
800 for (
k = 0;
k < sftVect->
length;
k++ ) {
839 for (
k = 0;
k < multVect->
length;
k++ ) {
864 INT4 nLines = 0, count1, nHarmonicSets;
888 if ( nHarmonicSets > 0 ) {
899 for ( count1 = 0; count1 < nHarmonicSets; count1++ ) {
900 nLines += *(
harmonics.numHarmonics + count1 );
974 numifos = MultiSFTVect->
length;
976 if ( linefiles != NULL ) {
982 for (
k = 0;
k < linefiles->
length;
k++ ) {
990 for (
j = 0;
j < numifos;
j++ ) {
994 width, window, linefiles->
data[
k], randPar ),
status );
#define __func__
log an I/O error, i.e.
#define ABORT(statusptr, code, mesg)
#define TRY(func, statusptr)
#define ATTATCHSTATUSPTR(statusptr)
#define ASSERT(assertion, statusptr, code, mesg)
#define DETATCHSTATUSPTR(statusptr)
#define INITSTATUS(statusptr)
#define RETURN(statusptr)
#define VTOT
Total detector velocity/c TO BE CHANGED DEPENDING ON DETECTOR.
void LogPrintf(LogLevel_t, const char *format,...) _LAL_GCC_PRINTF_FORMAT_(2
void LALNormalDeviates(LALStatus *status, REAL4Vector *deviates, RandomParams *params)
void LALRngMedBias(LALStatus *status, REAL8 *biasFactor, INT4 blkSize)
#define SFTCLEANH_EHEADER
#define SFTCLEANH_MSGENULL
void LALCleanMultiSFTVect(LALStatus *status, MultiSFTVector *multVect, INT4 width, INT4 window, LineNoiseInfo *lineInfo, RandomParams *randPar)
Function to clean a sft vector – calls LALCleanCOMPLEX8SFT repeatedly for each sft in vector.
#define SFTCLEANH_ELINENAME
void LALFindNumberLines(LALStatus *status, LineNoiseInfo *lineInfo, CHAR *fname)
Finds total number of spectral-lines contained in case the input file is a list of explicit spectral ...
#define SFTCLEANH_MSGEVAL
void LALRemoveKnownLinesInSFTVect(LALStatus *status, SFTVector *sftVect, INT4 width, INT4 window, CHAR *linefile, RandomParams *randPar)
function to remove lines from a sft vector given a file containing list of lines
void LALRemoveKnownLinesInMultiSFTVector(LALStatus *status, MultiSFTVector *MultiSFTVect, INT4 width, INT4 window, LALStringVector *linefiles, RandomParams *randPar)
top level function to remove lines from a multi sft vector given a list of files containing list of k...
void LALCleanSFTVector(LALStatus *status, SFTVector *sftVect, INT4 width, INT4 window, LineNoiseInfo *lineInfo, RandomParams *randPar)
Function to clean a sft vector – calls LALCleanCOMPLEX8SFT repeatedly for each sft in vector.
#define SFTCLEANH_MSGELINENAME
void LALChooseLines(LALStatus *status, LineNoiseInfo *outLine, LineNoiseInfo *inLine, REAL8 fMin, REAL8 fMax)
Takes a set of spectral lines and a frequency range as input and returns those lines which lie within...
void LALReadLineInfo(LALStatus *status, LineNoiseInfo *lineInfo, CHAR *fname)
Reads information from file containing list of explicit lines - obsolete.
void LALCheckLines(LALStatus *status, INT4 *countL, LineNoiseInfo *lines, REAL8 freq)
Function to count how many lines affect a given frequency.
void LALFindNumberHarmonics(LALStatus *status, LineHarmonicsInfo *harmonicInfo, CHAR *fname)
Looks into the input file containing list of lines, does some checks on the file format,...
void LALReadHarmonicsInfo(LALStatus *status, LineHarmonicsInfo *harmonicsInfo, CHAR *fname)
Reads in the contents of the input line-info file and fills up the LineHarmonicsInfo structure.
void LALHarmonics2Lines(LALStatus *status, LineNoiseInfo *lineInfo, LineHarmonicsInfo *harmonicsInfo)
Converts the list of harmonic sets into an explicit list of spectral lines.
#define SFTCLEANH_MSGEFILE
void LALCleanCOMPLEX8SFT(LALStatus *status, SFTtype *sft, INT4 width, INT4 window, LineNoiseInfo *lineInfo, RandomParams *randPar)
Function for cleaning a SFT given a set of known spectral disturbances.
CHAR * XLALGetChannelPrefix(const CHAR *name)
Find the valid CW detector prefix.
void LALDestroyVector(LALStatus *, REAL4Vector **)
void LALCreateVector(LALStatus *, REAL4Vector **, UINT4)
A vector of COMPLEX8FrequencySeries.
COMPLEX8FrequencySeries * data
Pointer to the data array.
UINT4 length
Number of elements in array.
structure for storing the contents of the input list of known spectral disturbances
REAL8 * gapFreq
frequency difference between adjacent harmonics in Hz
REAL8 * rightWing
width to the right in Hz
INT4 nHarmonicSets
number of sets of harmonics
REAL8 * leftWing
width to the left of each line in set in Hz
INT4 * numHarmonics
Number of harmonics.
REAL8 * startFreq
starting frequency of set in Hz
structure for storing list of spectral lines – constructed by expanding list of harmonics
REAL8 * lineFreq
central frequency of the line in Hz
REAL8 * leftWing
width to the left from central frequency in Hz
REAL8 * rightWing
width to the right in Hz
INT4 nLines
number of lines
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