Commit aa9a839d authored by Reinhard Prix's avatar Reinhard Prix

ResampFstat: fix truncated sinc-interpolation artifacts

- use Hamming-window on the sinc-interpolation, increasing
  sidebands (empty) to account for the transition band of ~ 4/(2Dterms+1)*Band
- use user input argument 'Dterms' to control length of sinc-interpolation window
- refs #5292
Original: 6138ebb3081913af227603f0600077ca19bfe713
parent ca67bd3d
......@@ -124,6 +124,7 @@ typedef struct tagResampWorkspace
typedef struct
{
UINT4 Dterms; // Number of terms to use (on either side) in Windowed-Sinc interpolation kernel
MultiCOMPLEX8TimeSeries *multiTimeSeries_DET; // input SFTs converted into a heterodyned timeseries
// ----- buffering -----
PulsarDopplerParams prev_doppler; // buffering: previous phase-evolution ("doppler") parameters
......@@ -252,11 +253,24 @@ XLALSetupFstatResamp ( void **method_data,
ResampMethodData *resamp = *method_data = XLALCalloc( 1, sizeof(*resamp) );
XLAL_CHECK( resamp != NULL, XLAL_ENOMEM );
resamp->Dterms = optArgs->Dterms;
// Set method function pointers
funcs->compute_func = XLALComputeFstatResamp;
funcs->method_data_destroy_func = XLALDestroyResampMethodData;
funcs->workspace_destroy_func = XLALDestroyResampWorkspace;
// Extra band needed for resampling: Hamming-windowed sinc used for interpolation has a transition bandwith of
// TB=(4/L)*fSamp, where L=2*Dterms+1 is the window-length, and here fSamp=Band (i.e. the full SFT frequency band)
// However, we're only interested in the physical band and we'll be throwing away all bins outside of this.
// This implies that we're only affected by *half* the transition band TB/2 on either side, as the other half of TB is outside of the band of interest
// (and will actually get aliased, i.e. the region [-fNy - TB/2, -fNy] overlaps with [fNy-TB/2,fNy] and vice-versa: [fNy,fNy+TB/2] overlaps with [-fNy,-fNy+TB/2])
// ==> therefore we only need to add an extra TB/2 on each side to be able to safely avoid the transition-band effects
REAL8 f0 = multiSFTs->data[0]->data[0].f0;
REAL8 dFreq = multiSFTs->data[0]->data[0].deltaF;
REAL8 Band = multiSFTs->data[0]->data[0].data->length * dFreq;
REAL8 extraBand = 2.0 / ( 2 * optArgs->Dterms + 1 ) * Band;
XLAL_CHECK ( XLALMultiSFTVectorResizeBand ( multiSFTs, f0 - extraBand, Band + 2 * extraBand ) == XLAL_SUCCESS, XLAL_EFUNC );
// Convert SFTs into heterodyned complex timeseries [in detector frame]
XLAL_CHECK ( (resamp->multiTimeSeries_DET = XLALMultiSFTVectorToCOMPLEX8TimeSeries ( multiSFTs )) != NULL, XLAL_EFUNC );
......@@ -1001,11 +1015,10 @@ XLALBarycentricResampleMultiCOMPLEX8TimeSeries ( ResampMethodData *resamp, // [
} // for alpha < numSFTsX
const UINT4 Dterms = 8;
XLAL_CHECK ( ti_DET->length >= TimeSeries_SRCX_a->data->length, XLAL_EINVAL );
UINT4 bak_length = ti_DET->length;
ti_DET->length = TimeSeries_SRCX_a->data->length;
XLAL_CHECK ( XLALSincInterpolateCOMPLEX8TimeSeries ( TimeSeries_SRCX_a->data, ti_DET, TimeSeries_DETX, Dterms ) == XLAL_SUCCESS, XLAL_EFUNC );
XLAL_CHECK ( XLALSincInterpolateCOMPLEX8TimeSeries ( TimeSeries_SRCX_a->data, ti_DET, TimeSeries_DETX, resamp->Dterms ) == XLAL_SUCCESS, XLAL_EFUNC );
ti_DET->length = bak_length;
// apply heterodyne correction and AM-functions a(t) and b(t) to interpolated timeseries
......
......@@ -32,6 +32,7 @@
#include <lal/ComputeFstat.h>
#include <lal/SinCosLUT.h>
#include <lal/Factorial.h>
#include <lal/Window.h>
/*---------- DEFINES ----------*/
#define MYMAX(x,y) ( (x) > (y) ? (x) : (y) )
......@@ -853,30 +854,39 @@ XLALCheckVectorComparisonTolerances ( const VectorComparison *result, ///< [in]
} // XLALCheckVectorComparisonTolerances()
/** Interpolate a given regularly-spaced COMPLEX8 timeseries 'ts_in = x_in(j * dt)' onto new samples
* 'y_out(t_out)' using Shannon sinc interpolation truncated to (2*Dterms+1) terms, namely
* 'y_out(t_out)' using *windowed* Shannon sinc interpolation, windowed to (2*Dterms+1) terms, namely
* \f[
* \newcommand{\Dterms}{\mathrm{Dterms}}
* \f]
*
* \f{equation}{
* x(t) = \sum_{j = j^* - \Delta j}^{j^* + \Delta j} x_j \,\, \frac{\sin(\pi\delta_j)}{\pi\delta_j}\,,\quad\text{with}\quad
* x(t) = \sum_{j = j^* - \Dterms}^{j^* + \Dterms} x_j \, w_j \, \frac{\sin(\pi\delta_j)}{\pi\delta_j}\,,\quad\text{with}\quad
* \delta_j \equiv \frac{t - t_j}{\Delta t}\,,
* \f}
* and where \f$j^* \equiv \mathrm{round}(t / \Delta t)\f$.
* where \f$j^* \equiv \mathrm{round}(t / \Delta t)\f$, and
* where \f$w_j\f$ is the window used (here: Hamming)
*
* In order to implement this more efficiently, we observe that \f$\sin(\pi\delta_j) = (-1)^{(j-j0)}\sin(\pi\delta_{j0})\f$ for integer \f$j\f$,
* and therefore
*
* \f{equation}{
* x(t) = \frac{\sin(\pi\,\delta_{j0})}{\pi} \, \sum_{j = j^* - \Delta j}^{j^* + \Delta j} (-1)^{(j-j0)}\frac{x_j}{\delta_j}\,,
* x(t) = \frac{\sin(\pi\,\delta_{j0})}{\pi} \, \sum_{j = j^* - \Dterms}^{j^* + \Dterms} (-1)^{(j-j0)}\frac{x_j \, w_j}{\delta_j}\,,
* \f}
*
* NOTE: Using Dterms=0 corresponds to closest-bin interpolation
*
* NOTE2: samples *outside* the original timespan are returned as 0
*
* NOTE3: we're using a Hamming window of length L = 2*Dterms + 1, which has a pass-band wiggle of delta_p ~ 0.0022,
* and a transition bandwidth of (4/L) * Bandwidth.
* You need to make sure to include sufficient effective sidebands to the input timeseries, so that the transition band can
* be safely ignored or 'cut out' at the end
*/
int
XLALSincInterpolateCOMPLEX8TimeSeries ( COMPLEX8Vector *y_out, ///< [out] output series of interpolated y-values [must be same size as t_out]
const REAL8Vector *t_out, ///< [in] output time-steps to interpolate input to
const COMPLEX8TimeSeries *ts_in,///< [in] regularly-spaced input timeseries
UINT4 Dterms ///< [in] truncate sinc kernel sum to +-Dterms around max
UINT4 Dterms ///< [in] window sinc kernel sum to +-Dterms around max
)
{
XLAL_CHECK ( y_out != NULL, XLAL_EINVAL );
......@@ -889,6 +899,10 @@ XLALSincInterpolateCOMPLEX8TimeSeries ( COMPLEX8Vector *y_out, ///< [out] outpu
REAL8 dt = ts_in->deltaT;
REAL8 tmin = XLALGPSGetREAL8 ( &(ts_in->epoch) ); // time of first bin in input timeseries
REAL8Window *win;
UINT4 winLen = 2 * Dterms + 1;
XLAL_CHECK ( (win = XLALCreateHammingREAL8Window ( winLen )) != NULL, XLAL_EFUNC );
const REAL8 oodt = 1.0 / dt;
for ( UINT4 l = 0; l < numSamplesOut; l ++ )
......@@ -911,8 +925,10 @@ XLALSincInterpolateCOMPLEX8TimeSeries ( COMPLEX8Vector *y_out, ///< [out] outpu
continue;
}
UINT8 jStart = MYMAX ( jstar - Dterms, 0 );
UINT8 jEnd = MYMIN ( jstar + Dterms, numSamplesIn - 1 );
INT4 jStart0 = jstar - Dterms;
UINT4 jEnd0 = jstar + Dterms;
UINT4 jStart = MYMAX ( jStart0, 0 );
UINT4 jEnd = MYMIN ( jEnd0, numSamplesIn - 1 );
REAL4 delta_jStart = (t_by_dt - jStart);
REAL4 sin0, cos0;
......@@ -923,7 +939,7 @@ XLALSincInterpolateCOMPLEX8TimeSeries ( COMPLEX8Vector *y_out, ///< [out] outpu
REAL8 delta_j = delta_jStart;
for ( UINT8 j = jStart; j <= jEnd; j ++ )
{
COMPLEX8 Cj = sin0oopi / delta_j;
COMPLEX8 Cj = win->data->data[j - jStart0] * sin0oopi / delta_j;
y_l += Cj * ts_in->data->data[j];
......@@ -935,6 +951,8 @@ XLALSincInterpolateCOMPLEX8TimeSeries ( COMPLEX8Vector *y_out, ///< [out] outpu
} // for l < numSamplesOut
XLALDestroyREAL8Window ( win );
return XLAL_SUCCESS;
} // XLALSincInterpolateCOMPLEX8TimeSeries()
......
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