next up previous contents
Next: Structure: struct removed_lines Up: GRASP Routines: General purpose Previous: Function: multitaper_spectrum()   Contents

Function: multitaper_cross_spectrum()

0 void multitaper_cross_spectrum(float *o1, float *o2, int npoints, int padded_length, float delta_t, int nwin, float nwdt, double *ReImSpec12)

This function calculates the high resolution multitaper estimate of the (complex-valued) cross-correlation spectrum $\tilde o_1^*(f)\tilde o_2(f)$ of the two input time-series $o_1(t)$ and $o_2(t)$.

The arguments are:

o1: Input. o1[0..npoints-1] is an array of floating point variables containing the values of the input time-series $o_1(t)$. o1[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$.
o2: Input. o2[0..npoints-1] is an array of floating point variables containing the values of the input time-series $o_2(t)$, in exactly the same format as the previous argument.
npoints: Input. The total number $N$ of data points contained in the two input time-series.
padded_length: Input. The padded length $N_p$ is an integer power of 2, greater than (or equal to) the total number of data points. The tapered data sets are zero-padded out to this length. The total number of frequency bins (including DC and Nyquist) in the output cross-correlation spectrum is $N_f=1+N_p/2$.
delta_t: Input. The sampling period $\Delta t$ (in sec).
nwin: Input. The total number $K$ of data tapers used when forming the multitaper spectral estimate.
nwdt: Input. The (total sample time) $\times$ (frequency resolution bandwidth) product $N\Delta t\cdot W$.
ReImSpec12: Output. ReImSpec12[0...padded_length+1] is an array of double precision variables containing the values of the high resolution multitaper spectral estimate of the (complex-valued) cross-correlation spectrum $\tilde o_1^*(f)\tilde o_2(f)$. ReImSpec12[2*j] and ReImSpec12[2*j+1] contain, respectively, the values of the real and imaginary parts of
\begin{displaymath}
\Delta t^2\ {1\over K}\sum_{k=0}^{K-1}{1\over \lambda_k}
\le...
...t(\sum_{n=0}^{N-1} o_2(t_n) h_{k,n}\ e^{+i2\pi n j/N_p}\right)
\end{displaymath} (16.22.333)

where $h_{k,n}$ is the $n$th element of the $k$th order $NW \Delta t$ discrete prolate spheroidal sequence data taper, and $\lambda_k$ is its associated eigenvalue. (Here $j=0,1,\cdots,N_f-1=N_p/2$.) The data tapers are normalized so that $\sum_{n=0}^{N-1} h_{k,n}^2=N$.

If you want to obtain the same normalization as that used in the avg_spec() routine described by equation ([*]) for the case where $o_1(t)=o_2(t)$ then the output array ReImSpec12 should be multiplied by a factor of $2/ N \Delta t$.

Author:

Adapted from the original code (Lees and Park) by Bruce Allen (ballen@dirac.phys.uwm.edu), Adrian Ottewill (ottewill@relativity.ucd.ie), and Joseph Romano (romano@csd.uwm.edu).

Comments: None.


next up previous contents
Next: Structure: struct removed_lines Up: GRASP Routines: General purpose Previous: Function: multitaper_spectrum()   Contents
Bruce Allen 2000-11-19