Function: remove_spectral_lines()

0 void remove_spectral_lines(float *data, int npoints, int padded_length, float nwdt, int nwin, int max_lines, int maxpass, int *num_removed, struct removed_lines *line_list,
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 $\hat C_i$ 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 $N_p$ of zero-padded data that will be analyzed. Note that $N_p$ 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) $\times$ (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, ${\tt fimin} < {\tt fimax}$.

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, ${\tt fimin} < {\tt fimax}$.

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:

1. The mean value is subtracted from the data-set, and it is zero padded to the specified length.
2. The set of Fourier coefficients for the tapered data sets are determined.
3. >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 $1-1/{\tt npoints}$ (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.
4. The working list is now sorted into order of decreasing F-values.
5. 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.
6. 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.