Synopsis

#include <lal/LALDemod.h>

The following is a brief synopsis of the demodulation, or 'Coherent Transform', procedure.

In order to remove frequency and amplitude modulation of a time series $x_a$, we need two basic components:

Frequency modulation information

This is given through a phase model $\Phi$.

Amplitude modulation information

This is given through two functions $\hat{a}$ and $\hat{b}$, which are derived from the beam-pattern functions $F_{+}$ and $F_{\times}$.

Given these, the F statistic in the $b^{th}$ frequency bin is

$$\mathcal{F}_{b} = \frac{4}{S_{h}(f_{0})T_{0}} \frac{B\vert\h... ...}}\vert^{2}+A\vert F_{b}\vert^{2} - 2C \Re(F_{a}F_{b}^{*})}{D}$$ (10.1)

where

$$\hat{F_{\hat{a}}}=\sum_{a=0}^{\it {NM-1}} x_{a} \hat{a} e^{-2{\pi}i{\Phi}_{ab}(\vec{\lambda})}$$ (10.2)

$$\hat{F_{\hat{b}}}=\sum_{a=0}^{\it {NM-1}} x_{a} \hat{b} e^{-2{\pi}i{\Phi}_{ab}(\vec{\lambda})}$$ (10.3)

$T_{0}$ is the observation time, $S_{h}$ is the noise power spectral density, and $A$, $B$, $C$, and $D$ are constants.

In writing the previous equation we have assumed that there is a total of $M\cdot N$ data samples and $0\leq a<MN$. $\Phi_{ab}$ is the expected phase at time $a$ for an intrinsic emission frequency $b\over T_{DeFT}$ (where the denominator is the DeFT time baseline). $\Phi$ depends on $\vec\lambda$, a vector of parameters that defines the phase model. Typically these are the source location and the spin-down parameter values of the template source for which one is demodulating. For simplicity, we will focus only on $F_{a}$; the analysis for $F_{b}$ is identical. Let us now suppose that the time series $x_a$ is composed of $M$ chunks, each of $N$ samples. If we introduce a short-time index $0\leq j<N-1$ and a short time-series index $0\leq \alpha <M-1$, so that $a=N\alpha+j$, we can rewrite the above sum as

$$\hat{F_{\hat{a}}}({\vec{\lambda}})=\sum_{\alpha=0}^{M-1}\sum... ...x_{\alpha j}a_{\alpha j}e^{-2{\pi}i{\Phi}_{ab}(\vec{\lambda})}$$ (10.4)

Note that $\hat{a}(t)$ is a periodic function with period equal to one sidereal day. Since the sum over $N$ is on a timescale much shorter than that (say, 1 hour), then $\hat{a}(t)$ won't change significantly, and thus can be taken outside of that summation, and then is evaluated at the midpoint of each SFT time. Now, If $\tilde{x}_{\alpha k}$ is the matrix of FTs formed along the short time index $j$

$$x_{\alpha j}=\frac{1}{N}\sum_{k=0}^{N-1}\tilde{x}_{\alpha k}e^{2\pi{i}\frac{jk}{N}},$$ (10.5)

making the appropriate substitutions, Eq. becomes

$$\hat{F_{\hat{a}}}({\vec{\lambda}})=\sum_{\alpha=0}^{M-1}\hat... ...{-2\pi i(\Phi_{\alpha jb}(\vec{\lambda})-\frac{jk}{N})}\right]$$ (10.6)

We assume that the phase evolution can be described as linear in $t$ during the time duration $T_{SFT}$; thus we can Taylor-expand $\Phi$ around the temporal midpoint of every SFT time data chunk. For large values of $N$, the summation over $j$ in eq. () can be expressed in closed form, thus saving computations, and eq. () can be rewritten as

$$\hat{F_{\hat{a}}}=\sum_{\alpha=0}^{M-1}\hat{a}_{\alpha}e^{i ... ...0}^{N-1}\tilde{x}_{\alpha\beta} P_{\alpha k}(b,\vec{\lambda}),$$ (10.7)

with

$$P_{\alpha k}(b,\vec{\lambda})={\sin{x'}\over x'}-i{1-\cos{x'}\over x'}$$ (10.8)

$$x'=\sum_{s} f_s B_{s\alpha} - k$$ (10.9)

$$y_\alpha=\sum_{s} f_s A_{s\alpha}.$$ (10.10)

In the previous expressions $f_s$ indicate the spin-down parameters of different orders (labeled by the index $s$^10.1), and $A_{s\alpha}$ and $B_{s\alpha}$ are functions that depend on the phase evolution, whose values depend on $\alpha$ and on $\vec\lambda$. The values of these functions are calculated by the ComputeSky() routine, also in this package. Incidentally, in the code, these are the values contained in the variable skyConst. Note that the function $P_{\alpha k}$ is peaked around $x'=0$. Thus in the summation over $k$ in eq. () one only needs to consider a few values (NTERMS) of $k$ around $k^*$ such that $x'(k^*)\approx 0$. This approximation again saves computations. Eq. () can then be rewritten as

$$\hat{F_{\hat{a}}}=\sum_{\alpha=0}^{M-1}\hat{a}_{\alpha}e^{i ... ...NTERMS} \tilde x_{\alpha\beta} P_{\alpha k}(b,\vec{\lambda}).$$ (10.11)

If $NTERMS$ is 8 the power loss due to this approximation is less than $\sim 5\%$.

Now, computing $\hat{F_{\hat{a}}}$ and $\hat{F_{\hat{b}}}$ can be done in parallel; given the approximations we have made, for each iteration of the $\alpha$ loop, one computes first $P_{\alpha k}$ (through the k-loop), multiplies by $\tilde{x}_{\alpha k}$, and then forms the statistics of at the same time. After all the iterations of the $\alpha$ loop are complete, that is, when all SFTs have been exhausted, the final statistic is computed.