Function: multitaper_spectrum()

0 multitaper_spectrum(float *data, int npoints,int kind, int nwin, float nwdt, int inorm, float dt, float *ospec, float *dof, float *fvalues, int padded_length, float *cest, int dospec)

This function calculates the multi-taper estimate of the amplitude spectrum $\sqrt{S(f)}$:

$$\sqrt{\hat S^{(mt)}(f_j)}:= \sqrt{{1\over K}\sum_{k=0}^{K-1}... ..._{m=0}^{N-1}x_m h_{k,m}\ e^{-i2\pi mj/N_p}\right\vert^2}\ .$$ (16.21.332)

Here $x_m$ is the $m$th element of the input data array, $h_{k,m}$ is the $m$th element of the $k$th order $NW \Delta t$ discrete prolate spheroidal sequence data taper, $\lambda_k$ is its associated eigenvalue, and $f_j$ is the discrete frequency $f_j:=j/N_p\Delta t$, where $j=0,1,\cdots,N_f-1=N_p/2$. The above multi-taper estimate differs slightly from Equation (333) in Percival and Walden: (i) there is the square root; (ii) the data tapers are normalized so that $\sum_{n=0}^{N-1} h_{k,n}^2=N$; (iii) there is no factor of $\Delta t$; (iv) the individual estimates are weighted by $1/\lambda_k$. (But for $K< 2NW \Delta t$, $1/\lambda_k\approx 1$).

For the sake of efficiency, multitaper_spectrum() computes, and then stores internally, the Slepian taper functions, so that if it is called a second time (and needs the same tapers) they do not need to be re-computed. If called with different parameters it recomputes the Slepian tapers for the new parameters.

If you want to obtain the same normalization as that used in the avg_spec() routine described by equation (), the output array ospec[] should first be squared, and then multiplied by a factor of $2\Delta t/N$.

The arguments are:

data: Input. Pointer to the time-domain data array, data[0..npoints-1].

npoints: Input. The total number $N$ of data points in the data array.

kind: Input. If set to 1, compute the normal (i.e., ``high resolution'') multi-taper estimate of the amplitude spectrum. If set to 2, compute the ``adaptive" multi-taper estimate $\sqrt{\hat S^{(amt)}(f_j)}$ of the amplitude spectrum, defined by Percival and Walden Equation (370a).

nwin: Input. The number of tapers to use.

nwdt: Input. The (total sample time) $\times$ (frequency resolution bandwidth) product.

inorm: Input. Determines choice of normalization. Possible values are

1: Divide spectrum by $N^2$.

2: Divide spectrum by $\Delta t^2$.

3: Divide spectrum by $N$.

4: Divide spectrum by $1$.

dt: Input. Sample interval (only used for normalization).

ospec: Output. The output spectrum, including both DC and Nyquist frequency bins. The array range is ospec[0..padded_length/2]. Warning -this is an odd number of entries. The user must provide a pointer to sufficient storage space.

dof: Output. The effective number of degrees of freedom of the spectral estimator at a given frequency, defined by Percival and Walden eqn (370b). The number of degrees of freedom is the constant nwin-1 for kind=1 above, and only useful in the adaptive case where kind=2. The array range is dof[0..padded_length/2]. Warning - this is an odd number of entries. The user must provide a pointer to sufficient storage space.

fvalues: Output. The value of the F-statistic in each frequency bin spectrum, including both DC and Nyquist. This is defined by Percival and Walden equation (499c), and roughly speaking is the ratio of the energy explained by the hypothesis that one has a fixed-amplitude spectral line at that frequency to the energy not explained by this hypothesis. The array range is fvalues[0..padded_length/2]. Warning -this is an odd number of entries. The user must provide a pointer to sufficient storage space.

padded_length: Input. The padded length $N_p$ is an integer power of 2, greater than (or equal to) npoints. The (tapered) data is zero-padded out to this length. You generally want $N_p$ to be around four to eight times larger than the length of your data set, to get decent frequency resolution. The number of frequency bins (including DC and Nyquist) in the output spectrum is $N_f=1+N_p/2$.

cest: Output. The estimated Fourier coefficients of any spectral lines in the data. The real and imaginary parts at DC are contained in cest[0],cest[1]. The next higher frequency bin has its real/imaginary parts contained in cest[2],cest[3], and so on. This pattern continues up to and including the Nyquist frequency. The length of the array is cest[0..padded_length+1]. The normalization/sign conventions are identical to Percival and Walden eqns (499a) and the example on line 20 of page 513, except that the sign of the imaginary part is reversed, because the Percival and Walden FFT conventions eqn (65ab) are opposite to Numerical Recipes. The user must provide a pointer to sufficient storage space.

dospec: Input. If set non-zero, then the power spectrum (pointed to by ospec) is calculated. If set to zero, then to save time in situations where all that is needed is cest, the power spectrum ospec is not calculated.

Author: Adapted from the original code (Lees and Park) by Bruce Allen (ballen@dirac.phys.uwm.edu) and Adrian Ottewill (ottewill@relativity.ucd.ie).

Comments: None.