Commit 01b428a0 authored by Yuanhao Zhang's avatar Yuanhao Zhang
Browse files

PulsarCrosscorr_v2: reduced redundant times of calling sinc function and...

PulsarCrosscorr_v2: reduced redundant times of calling sinc function and output sinc values into a list in DopplerShiftedInfo function
Original: 8b9f7daf4f01d1e9a95472bfb4b391cd6792455e
parent dc1d2cbe
......@@ -116,6 +116,7 @@ int main(int argc, char *argv[]){
UINT4Vector *lowestBins = NULL;
REAL8Vector *kappaValues = NULL;
REAL8Vector *signalPhases = NULL;
REAL8Vector *sincValueList = NULL;
PulsarDopplerParams XLAL_INIT_DECL(dopplerpos);
PulsarDopplerParams thisBinaryTemplate, binaryTemplateSpacings;
......@@ -142,7 +143,7 @@ int main(int argc, char *argv[]){
CrossCorrBinaryOutputEntry thisCandidate;
UINT4 checksum;
LogPrintf (LOG_CRITICAL, "starting time\n", 0);/*for debug convenience to record calculating time*/
LogPrintf (LOG_CRITICAL, "Starting time\n", 0);/*for debug convenience to record calculating time*/
/* initialize and register user variables */
if ( XLALInitUserVars( &uvar ) != XLAL_SUCCESS ) {
LogPrintf ( LOG_CRITICAL, "%s: XLALInitUserVars() failed with errno=%d\n", __func__, xlalErrno );
......@@ -448,6 +449,10 @@ int main(int argc, char *argv[]){
LogPrintf ( LOG_CRITICAL, "%s: XLALCreateREAL8Vector() failed with errno=%d\n", __func__, xlalErrno );
XLAL_ERROR( XLAL_EFUNC );
}
if ((sincValueList = XLALCreateREAL8Vector ( numSFTs ) ) == NULL){
LogPrintf ( LOG_CRITICAL, "%s: XLALCreateREAL8Vector() failed with errno=%d\n", __func__, xlalErrno );
XLAL_ERROR( XLAL_EFUNC );
}
/* args should be : spacings, min and max doppler params */
while ( (GetNextCrossCorrTemplate(&dopplerShiftFlag, &dopplerpos, &binaryTemplateSpacings, &minBinaryTemplate, &maxBinaryTemplate) == 0) )
......@@ -464,13 +469,12 @@ int main(int argc, char *argv[]){
}
}
if ( (XLALGetDopplerShiftedFrequencyInfo( shiftedFreqs, lowestBins, kappaValues, signalPhases, uvar.numBins, &dopplerpos, sftIndices, inputSFTs, multiBinaryTimes, Tsft ) != XLAL_SUCCESS ) ) {
if ( (XLALGetDopplerShiftedFrequencyInfo( shiftedFreqs, lowestBins, kappaValues, signalPhases, sincValueList, uvar.numBins, &dopplerpos, sftIndices, inputSFTs, multiBinaryTimes, Tsft ) != XLAL_SUCCESS ) ) {
LogPrintf ( LOG_CRITICAL, "%s: XLALGetDopplerShiftedFrequencyInfo() failed with errno=%d\n", __func__, xlalErrno );
XLAL_ERROR( XLAL_EFUNC );
}
if ( (XLALCalculatePulsarCrossCorrStatistic( &ccStat, &evSquared, curlyGUnshifted, signalPhases, lowestBins, kappaValues, uvar.numBins, sftPairs, sftIndices, inputSFTs, multiWeights) != XLAL_SUCCESS ) ) {
if ( (XLALCalculatePulsarCrossCorrStatistic( &ccStat, &evSquared, curlyGUnshifted, signalPhases, lowestBins, kappaValues, sincValueList, sftPairs, sftIndices, inputSFTs, multiWeights, uvar.numBins) != XLAL_SUCCESS ) ) {
LogPrintf ( LOG_CRITICAL, "%s: XLALCalculateAveCrossCorrStatistic() failed with errno=%d\n", __func__, xlalErrno );
XLAL_ERROR( XLAL_EFUNC );
}
......@@ -502,6 +506,7 @@ int main(int argc, char *argv[]){
XLALDestroyREAL8Vector ( kappaValues );
XLALDestroyUINT4Vector ( lowestBins );
XLALDestroyREAL8Vector ( shiftedFreqs );
XLALDestroyREAL8Vector ( sincValueList );
XLALDestroyMultiSSBtimes ( multiBinaryTimes );
XLALDestroyMultiSSBtimes ( multiSSBTimes );
XLALDestroyREAL8Vector ( curlyGUnshifted );
......@@ -526,9 +531,9 @@ int main(int argc, char *argv[]){
/* check memory leaks if we forgot to de-allocate anything */
LALCheckMemoryLeaks();
LogPrintf (LOG_CRITICAL, "The metric element g_ff=%.9f\n", diagff);
LogPrintf (LOG_CRITICAL, "The metric element g_aa=%.9f\n", diagaa);
LogPrintf (LOG_CRITICAL, "The metric element g_TT=%.9f\n", diagTT);
LogPrintf (LOG_CRITICAL, "The metric element g_ff=%.9f\n", diagff, 0);
LogPrintf (LOG_CRITICAL, "The metric element g_aa=%.9f\n", diagaa, 0);
LogPrintf (LOG_CRITICAL, "The metric element g_TT=%.9f\n", diagTT, 0);
LogPrintf (LOG_CRITICAL, "End time\n", 0);/*for debug convenience to record calculating time*/
return 0;
......
......@@ -31,6 +31,7 @@ int XLALGetDopplerShiftedFrequencyInfo
UINT4Vector *lowestBins, /**< Output list of bin indices */
REAL8Vector *kappaValues, /**< Output list of bin offsets */
REAL8Vector *signalPhases, /**< Output list of signal phases */
REAL8Vector *sincList, /**< Output list of sinc factors */
UINT4 numBins, /**< Number of frequency bins to use */
PulsarDopplerParams *dopp, /**< Doppler parameters for signal */
SFTIndexList *sftIndices, /**< List of indices for SFTs */
......@@ -47,7 +48,8 @@ int XLALGetDopplerShiftedFrequencyInfo
if ( signalPhases->length !=numSFTs
|| shiftedFreqs->length !=numSFTs
|| lowestBins->length !=numSFTs
|| kappaValues->length !=numSFTs ) {
|| kappaValues->length !=numSFTs
|| sincList->length !=numSFTs ) {
XLALPrintError("Lengths of SFT-indexed lists don't match!");
XLAL_ERROR(XLAL_EBADLEN );
}
......@@ -101,6 +103,10 @@ int XLALGetDopplerShiftedFrequencyInfo
lowestBins->data[sftNum]
= ceil( fminusf0 * Tsft - 0.5*numBins );
kappaValues->data[sftNum] = lowestBins->data[sftNum] - fminusf0 * Tsft;
sincList->data[sftNum]=1;
for (UINT8 l=0; l < numBins; l++) {
sincList->data[sftNum] *= gsl_sf_sinc(kappaValues->data[sftNum]+l);
}
/* printf("f=%.7f, f0=%.7f, Tsft=%g, numbins=%d, lowestbin=%d, kappa=%g\n",
shiftedFreqs->data[sftNum],
inputSFTs->data[detInd]->data[sftInd].f0,
......@@ -281,18 +287,20 @@ int XLALCalculatePulsarCrossCorrStatistic
REAL8Vector *signalPhases, /* Input: Phase of signal for each SFT */
UINT4Vector *lowestBins, /* Input: Bin index to start with for each SFT */
REAL8Vector *kappaValues, /* Input: Fractional offset of signal freq from best bin center */
UINT4 numBins, /* Input: Number of bins to include in calc */
REAL8Vector *sincList, /* Input: input the sinc factors*/
SFTPairIndexList *sftPairs, /* Input: flat list of SFT pairs */
SFTIndexList *sftIndices, /* Input: flat list of SFTs */
MultiSFTVector *inputSFTs, /* Input: SFT data */
MultiNoiseWeights *multiWeights /* Input: nomalizeation factor S^-1 & weights for each SFT */
MultiNoiseWeights *multiWeights, /* Input: nomalizeation factor S^-1 & weights for each SFT */
UINT4 numBins /**Input Number of frequency bins to be taken into calc */
)
{
UINT8 numSFTs = sftIndices->length;
if ( signalPhases->length !=numSFTs
|| lowestBins->length !=numSFTs
|| kappaValues->length !=numSFTs ) {
|| kappaValues->length !=numSFTs
|| sincList->length !=numSFTs ) {
XLALPrintError("Lengths of SFT-indexed lists don't match!");
XLAL_ERROR(XLAL_EBADLEN );
}
......@@ -355,7 +363,6 @@ int XLALCalculatePulsarCrossCorrStatistic
lowestBin1, numBins, lenDataArray1 );
for (UINT8 j=0; j < numBins; j++) {
COMPLEX16 data1 = dataArray1[lowestBin1+j];
REAL8 sincFactor = gsl_sf_sinc(kappaValues->data[sftNum1]+j);
/* Normalized sinc, i.e., sin(pi*x)/(pi*x) */
INT4 ccSign = baseCCSign;
UINT4 lowestBin2 = lowestBins->data[sftNum2];
......@@ -365,9 +372,9 @@ int XLALCalculatePulsarCrossCorrStatistic
lowestBin2, numBins, lenDataArray2 );
for (UINT8 k=0; k < numBins; k++) {
COMPLEX16 data2 = dataArray2[lowestBins->data[sftNum2]+k];
sincFactor *= gsl_sf_sinc(kappaValues->data[sftNum2]+k);
nume += /*multiWeights->data[detInd1]->data[sftNum1] * multiWeights->data[detInd2]->data[sftNum2] **/ creal ( GalphaCC * ccSign * sincFactor * conj(data1) * data2 );
REAL8 GalphaAmp = curlyGAmp->data[alpha] /** multiWeights->data[detInd1]->data[sftNum1] * multiWeights->data[detInd2]->data[sftNum2]*/ * sincFactor;
REAL8 sincFactor = sincList->data[sftNum1] * sincList->data[sftNum2];
nume += creal ( GalphaCC * ccSign * sincFactor * conj(data1) * data2 ); /*multiWeights->data[detInd1]->data[sftNum1] * multiWeights->data[detInd2]->data[sftNum2] **/
REAL8 GalphaAmp = curlyGAmp->data[alpha] * sincFactor ; /** multiWeights->data[detInd1]->data[sftNum1] * multiWeights->data[detInd2]->data[sftNum2]*/
curlyGSqr += SQUARE( GalphaAmp );
ccSign *= -1;
}
......
......@@ -109,6 +109,7 @@ int XLALGetDopplerShiftedFrequencyInfo
UINT4Vector *lowestBins,
REAL8Vector *kappaValues,
REAL8Vector *signalPhases,
REAL8Vector *sincList,
UINT4 numBins,
PulsarDopplerParams *dopp,
SFTIndexList *sfts,
......@@ -152,11 +153,12 @@ int XLALCalculatePulsarCrossCorrStatistic
REAL8Vector *signalPhases,
UINT4Vector *lowestBins,
REAL8Vector *kappaValues,
UINT4 numBins,
REAL8Vector *sincList,
SFTPairIndexList *sftPairs,
SFTIndexList *sftIndices,
MultiSFTVector *inputSFTs,
MultiNoiseWeights *multiWeights
MultiNoiseWeights *multiWeights,
UINT4 numBins
)
;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment