float *mtap_spec_init, float *mtap_spec_final, int dospecs, int fimin, int fimax)

This routine automatically identifies and removes ``spectral lines"
from a time series. The procedure followed is described in Percival
and Walden Chapter 10. A worked example may be found in Section 10.13
of that book, and the next subsection of the GRASP manual includes two
example programs which use `remove_spectral_lines()`. Upon return,
`remove_spectral_lines()` provides both an ``initial" multi-taper
amplitude spectrum, of the original data, and a ``final" multi-taper
amplitude spectrum, after line removal.
Upon return, the data set has the spectral lines
subtracted. This routine also returns a list of the lines removed.
For each line it provides the frequency bin (for the padded data set)
in which the line falls, the value of the F-test for that line, and the
complex coefficient defined by Percival and Walden eqn (499a)
which defines the line.

The arguments are:

`data:`Input. A pointer to the time-series array`data[0..npoints-1]`.`npoints:`Input. The number of points in the previous array.`padded_length:`Input. The number of points of zero-padded data that will be analyzed. Note that must be an integer power of two greater than or equal to`npoints`. We recommend that you use at least a factor of four greater, to obtain sufficient frequency resolution to accurately identify/remove spectral lines.`nwdt:`Input. The (total sample time) (frequency resolution bandwidth) product.`nwin:`Input. Number of Slepian tapers. See previous sections.`max_lines:`Input. The maximum number of spectral lines that you want removed. The array`line_list[0..max_lines-1]`must have at least this length.`maxpass:`Input. The maximum number of iterations or passes through the line-removal loop described below. Set to a large number to make as many passes as needed to remove all the spectral lines.`num_removed:`Output. The actual number of spectral lines subtracted from the data.`line_list:`Ouput. A list of structures`line_list[0..num_removed-1]`containing the frequency bin, real and imaginary parts of the removed line, and the F-test significance value associated with the*first*removal of the line. Upon return from this function, the elements of`line_list[]`are sorted into increasing frequency-bin order.`mtap_spec_init:`Output. The multi-taper estimate of the amplitude spectrum of the*initial*`data[]`array, including both DC and Nyquist frequency bins. The array range is`mtap_spec_init[0..padded_length/2]`. Warning -this is an*odd*number of entries. The user must provide a pointer to sufficient storage space.`mtap_spec_final:`Output. The multi-taper estimate of the amplitude spectrum of the*final*`data[]`array, with the spectral lines subtracted, including both DC and Nyquist frequency bins. The array range is`mtap_spec_final[0..padded_length/2]`. Warning -this is an*odd*number of entries. The user must provide a pointer to sufficient storage space.`dospecs:`Input. If set non-zero, then the initial/final amplitude spectra (pointed to by`mtap_spec_init`and`mtap_spec_final`) are calculated. If set to zero, then to save time in situations where all that is needed a list of spectral lines and their amplitudes and phases, then neither of the amplitude spectra are calculated.`fimin:`Input. In situations where all that is needed is a list of spectral lines and their amplitudes, and it is desired to limit the search to a restricted range of frequencies, then`fimin`defines the lower bound of the range of (padded) frequency bins which are searched for spectral lines. The range of`fimin`is 0..`padded_length/2`. Also, .`fimax:`Input. In situations where all that is needed is a list of spectral lines and their amplitudes, and it is desired to limit the search to a restricted range of frequencies, then`fimax`defines the upper bound of the range of (padded) frequency bins which are searched for spectral lines. The range of`fimax`is 0..`padded_length/2`. Also, .

The algorithm used by `remove_spectral_lines()` is an automated
version of the procedure illustrated in Percival and Walden Section 10.13.
The steps followed are:

- The mean value is subtracted from the data-set, and it is zero padded to the specified length.
- The set of Fourier coefficients for the tapered data sets are determined.
- >From these coefficients the F-statistic is determined for each frequency bin (Percival and Walden eqn (499c)). If the confidence level (that the frequency bin contains a spectral line) exceeds (Percival and Walden pg 513), an estimator of the spectral line coefficients is constructed, and the line is placed onto a working list. If no frequency bins exceed this level of confidence, we are finished.
- The working list is now sorted into order of decreasing F-values.
- To ensure that we do not remove the same line twice, the spectral line associated with each spectral line on the working list is subtracted from the data-set, provided that it does not lie within a frequency width of of a stronger (larger F-value) line.
- We return to step 1 above, iterating this procedure, provided that the
number of times that we have passed by step 1 is less than or equal to
`maxpass`.

- Author: Bruce Allen (ballen@dirac.phys.uwm.edu) and Adrian Ottewill (ottewill@relativity.ucd.ie).
- Comments:
If
`max_lines`is not large enough, then the`line_list[]`array may not contain all of the possible spectral lines, which exceed the confidence level above. This may even be the case if`num_removed`is less than`max_lines`. We suggest that you make`max_lines`somewhat larger than`num_removed`. One ought to be able to improve on this routine, by using the array of F-values generated internally and interpolating to find the frequency of the lines more precisely. One might also be able to fit to a model of two closely separated lines to better remove certain ``split" features, or to fit to an exponentially-decaying model to remove other broadened features.