next up previous contents
Next: Function: optimal_filter() Up: GRASP Routines: Stochastic background Previous: Function: extract_noise()   Contents


Function: extract_signal()

0

void extract_signal(int average, float *in1, float *in2, int n, float delta_t, double *whiten1, double *whiten2, double *signal12)
This function calculates the real-time cross-correlation spectrum $\tilde s_{12}(f)$ of the unwhitened data streams $s_1(t)$ and $s_2(t)$, using a Hann window and averaging the spectrum for two overlapped data sets, if desired.

The arguments of extract_signal() are:

average: Input. An integer variable that should be set equal to 1 if the values of the real-time cross-correlation spectra corresponding to two overlapped data sets are to be averaged.
in1: Input. in1[0..n-1] is an array of floating point variables containing the values of the assumed continuous-in-time whitened data stream $o_1(t)$ produced by the first detector. $o_1(t)$ is the convolution of detector whitening filter $W_1(t)$ with the data stream $s_1(t):=h_1(t)+n_1(t)$, where $h_1(t)$ is the gravitational strain and $n_1(t)$ is the noise intrinsic to the detector. The variables in1[] have units of rHz (or ${\rm sec}^{-1/2}$), which follows from the definition of $s_1(t)$ as a strain and $\tilde W_1(f)$ as the ``inverse'' of the square root of the noise power spectrum $P_1(f)$. in1[i] contains the value of $o_1(t)$ evaluated at the discrete time $t_i=i\Delta t$, where $i=0,1,\cdots,N-1$.
in2: Input. in2[0..n-1] is an array of floating point variables containing the values of the assumed continuous-in-time whitened data stream $o_2(t)$ produced by the second detector, in exactly the same format as the previous argument.
n: Input. The number $N$ of data points corresponding to an observation time $T:=N\ \Delta t$, where $\Delta t$ is the sampling period of the detectors, defined below. $N$ should equal an integer power of 2.
delta_t: Input. The sampling period $\Delta t$ (in sec) of the detectors.
whiten1: Input. whiten1[0..n-1] is an array of double precision variables containing the values of the real and imaginary parts of the spectrum $\tilde W_1(f)$ of the whitening filter of the first detector. These variables have units ${\rm rHz}/{\rm strain}$ (or ${\rm sec}^{-1/2}$), which are inverse to the units of the square root of the noise power spectrum $P_1(f)$. whiten1[2*i] and whiten1[2*i+1] contain, respectively, the values of the real and imaginary parts of $\tilde W_1(f)$ evaluated at the discrete frequency $f_i=i/(N\Delta t)$, where $i=0,1,\cdots,N/2-1$.
whiten2: Input. whiten2[0..n-1] is an array of double precision variables containing the values of the real and imaginary parts of the spectrum $\tilde W_2(f)$ of the whitening filter of the second detector, in exactly the same format as the previous argument.
signal12: Output. signal12[0..n/2-1] is an array of double precision variables containing the values of the real-time cross-correlation spectrum
\begin{displaymath}
\tilde s_{12}(f):=\left(\tilde s_1^*(f)\ \tilde s_2(f) + {\rm c.c.}\right)\ ,
\end{displaymath} (11.16.244)

where $\tilde s_1(f)$ and $\tilde s_2(f)$ are the Fourier transforms of the unwhitened data streams $s_1(t)$ and $s_2(t)$ produced by the two detectors. These variables have units of ${\rm strain}^2\cdot{\rm sec}^2$ (or simply ${\rm sec}^2$). signal12[i] contains the value of $\tilde s_{12}(f)$ evaluated at the discrete frequency $f_i=i/(N\Delta t)$, where $i=0,1,\cdots,N/2-1$.

extract_signal() calculates the real-time cross-correlation spectrum $\tilde s_{12}(f)$ as follows:

(i)
It first stores the input data streams $o_1(t)$ and $o_2(t)$ in the last two-thirds of internally-defined static buffers buf1[0..3*n/2-1] and buf2[0..3*n/2-1]. The first one-third of these buffers contains the input data left over from the previous call.
(ii)
It then multiplies the first two-thirds of these buffers by the Hann window function:
\begin{displaymath}
w(t):=\sqrt{{8\over 3}}\cdot{1\over 2}
\left[1-\cos\left({2\pi t\over T}\right)\right]\ .
\end{displaymath} (11.16.245)

The factor $\sqrt{8/3}$ is the ``window squared-and-summed'' factor described in Numerical Recipes in C, p.553. It is needed to offset the reduction in power that is introduced by the windowing.
(iii)
The windowed data is then Fourier transformed into the frequency domain, where it is unwhitened by dividing by the (complex) spectra $\tilde W_1(f)$ and $\tilde W_2(f)$, which represent the whitening filters of the two detectors. The resulting unwhitened frequency components are denoted by ${}^{(1)}\tilde s(f)$ and ${}^{(1)}\tilde s(f)$; the superscript $(1)$ indicates that we are analyzing the first of two overlapped data sets.
(iv)
The real-time cross-correlation spectrum is then calculated according to:
\begin{displaymath}
{}^{(1)}\tilde s_{12}(f):=\left[{}^{(1)}\tilde s_1^*(f)\ {}^{(1)}\tilde s_2(f) + {\rm c.c.}\ \right]\ .
\end{displaymath} (11.16.246)

(v)
The data contained in the last two-thirds of the buffers is then copied to the first two-thirds of the buffers, and steps (ii)-(iv) are repeated, yielding a second real-time cross-correlation spectrum ${}^{(2)}\tilde s_{12}(f)$.
(vi)
If average=1, $\tilde s_{12}(f)$ is given by:
\begin{displaymath}
\tilde s_{12}(f):={1\over 2}\left[\ {}^{(1)}\tilde s_{12}(f)+
\ {}^{(2)}\tilde s_{12}(f)\ \right]\ .
\end{displaymath} (11.16.247)

Otherwise, $\tilde s_{12}(f)={}^{(2)}\tilde s_{12}(f)$.
(vii)
Finally, the data contained in the last two-thirds of the buffers is again copied to the first two-thirds, in preparation for the next call to extract_sb(). The data saved in the first one-third of these buffers will match onto the next input data streams if the input data from one call of extract_sb() to the next is continuous.

Note: One should call extract_sb() with ${\tt average}\ne 1$, when one suspects that the current input data is not continuous with the data that was saved from the previous call. This is because a discontinuity between the ``old'' and ``new'' data sets has a tendency to introduce spurious large frequency components into the real-time cross-correlation spectrum, which should not be present. Since a single input data stream by itself is continuous, the cross-correlation spectrum ${}^{(2)}\tilde s_{12}(f)$ (which is calculated on the second pass through the data) will be free of these spurious large frequency components. This is why we set $\tilde s_{12}(f)$ equal to ${}^{(2)}\tilde s_{12}(f)$--and not equal to ${}^{(1)}\tilde s_{12}(f)$--when ${\tt average}\ne 1$.

Authors: Bruce Allen, ballen@dirac.phys.uwm.edu, and Joseph Romano, romano@csd.uwm.edu
Comments: Although it is possible and more efficient to write a single function to extract the real-time detector noise power and cross-correlation signal spectra simultaneously, we have chosen--for the sake of modularity--to write separate functions to perform these two tasks separately. (See also the comment at the end of Sec. [*].)


next up previous contents
Next: Function: optimal_filter() Up: GRASP Routines: Stochastic background Previous: Function: extract_noise()   Contents
Bruce Allen 2000-11-19