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


Function: extract_noise()

0

void extract_noise(int average, int which, float *in, int n, float delta_t, double *whiten_out, double *power)
This function calculates the real-time noise power spectrum $P(f)$ of a detector, using a Hann window and averaging the spectrum for two overlapped data sets, if desired.

The arguments of extract_noise() are:

average: Input. An integer variable that should be set equal to 1 if the values of the real-time noise power spectra corresponding to two overlapped data sets are to be averaged.
which: Input. An integer variable specifying which internally-defined static buffer should be used when overlapping the new input data set with data saved from a previous call. The allowed values are $1\le{\tt which}\le 16$.
in: Input. in[0..n-1] is an array of floating point variables containing the values of the assumed continuous-in-time whitened data stream $o(t)$ produced by the detector. $o(t)$ is the convolution of detector whitening filter $W(t)$ with the data stream $s(t):=h(t)+n(t)$, where $h(t)$ is the gravitational strain and $n(t)$ is the noise intrinsic to the detector. The variables in[] have units of rHz (or ${\rm sec}^{-1/2}$), which follows from the definition of $s(t)$ as a strain and $\tilde W(f)$ as the ``inverse'' of the square root of the noise power spectrum $P(f)$. in[i] contains the value of $o(t)$ evaluated at the discrete time $t_i=i\Delta t$, where $i=0,1,\cdots,N-1$.
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 detector, defined below. $N$ should equal an integer power of 2.
delta_t: Input. The sampling period $\Delta t$ (in sec) of the detector.
whiten_out: Input. whiten_out[0..n-1] is an array of double precision variables containing the values of the real and imaginary parts of the spectrum $\tilde W(f)$ of the whitening filter of the 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(f)$. whiten_out[2*i] and whiten_out[2*i+1] contain, respectively, the values of the real and imaginary parts of $\tilde W(f)$ evaluated at the discrete frequency $f_i=i/(N\Delta t)$, where $i=0,1,\cdots,N/2-1$.
power: Output. power[0..n/2-1] is an array of double precision variables containing the values of the real-time noise power spectrum $P(f)$ of the detector. Explicitly,
\begin{displaymath}
P(f):={2\over T}\ \tilde s^*(f)\tilde s(f)\ ,
\end{displaymath} (11.15.240)

where $\tilde s(f)$ is the Fourier transform of the unwhitened data stream $s(t)$ produced by the detector. These variables have units of ${\rm strain}^2/{\rm Hz}$ (or seconds). power[i] contains the value of $P(f)$ evaluated at the discrete frequency $f_i=i/(N\Delta t)$, where $i=0,1,\cdots,N/2-1$.

extract_noise() calculates the real-time noise power spectrum $P(f)$ as follows:

(i)
It first stores the input data stream $o(t)$ in the last two-thirds of an appropriately chosen static buffer buf[0..3*n/2-1]. The first one-third of this buffer contains the input data left over from the previous call.
(ii)
It then multiplies the first two-thirds of this buffer 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.15.241)

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) spectrum $\tilde W(f)$ of the whitening filter of the detector. The resulting unwhitened frequency components are denoted by ${}^{(1)}\tilde s(f)$; the superscript $(1)$ indicates that we are analyzing the first of two overlapped data sets.
(iv)
The real-time noise power spectrum is then calculated according to:
\begin{displaymath}
{}^{(1)}P(f):={2\over T}\ {}^{(1)}\tilde s^*(f)\ {}^{(1)}\tilde s(f)\ .
\end{displaymath} (11.15.242)

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

Otherwise, $P(f)={}^{(2)}P(f)$.
(vii)
Finally, the data contained in the last two-thirds of the buffer is again copied to the first two-thirds, in preparation for the next call to extract_noise(). The data saved in the first one-third of this buffer will match onto the next input data stream if the input data from one call of extract_noise() to the next is continuous.

Note: One should call extract_noise() 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 noise power spectrum, which should not be present. Since a single input data stream by itself is continuous, the noise power spectrum ${}^{(2)}P(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 $P(f)$ equal to ${}^{(2)}P(f)$--and not equal to ${}^{(1)}P(f)$--when ${\tt average}\ne 1$.

Authors: Bruce Allen, ballen@dirac.phys.uwm.edu, and Joseph Romano, romano@csd.uwm.edu
Comments: In the context of stochastic background simulations, it would be more efficient to extract the real-time noise power spectra at two detectors simultaneously. However, for modularity of design, and to allow this function to be used possibly for ``single-detector'' gravity-wave searches, we decided to write the above routine instead.


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