This function calculates the multi-taper estimate of the
amplitude spectrum :

(16.21.332) |

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 .

The arguments are:

`data:`Input. Pointer to the time-domain data array,`data[0..npoints-1]`.`npoints:`Input. The total number 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 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) (frequency resolution bandwidth) product.`inorm:`Input. Determines choice of normalization. Possible values are- 1: Divide spectrum by .
- 2: Divide spectrum by .
- 3: Divide spectrum by .
- 4: Divide spectrum by .

`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 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 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 .`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.