next previous
Up: Techniques used in

5. Power spectrum analysis

5.1. Power spectrum scaling

One must first choose a suitable scaling scheme to use for a spectrum that consists neither of incoherent noise nor of completely coherent sine waves, but of modes of finite lifetime, superimposed on a non-negligible noise background. The situation is further complicated because, in practice, a substantial fraction of the data in the time series is missing: our best 2-month spectrum so far has 94 per cent fill, but in some early years the best 2-month fills were around 25 per cent (early to mid 1980 s). For frequency measurements, of course, the scaling can be arbitrary, but in order to compare mode amplitudes obtained from different spectra the scaling needs to be carefully considered. There are two main methods of scaling power spectra: "power per unit frequency'' and "equivalent sine wave''. In the former method, the total power in the spectrum is equal to the sum of the squares of the velocity residuals - or, in our case, to the sum of squares of the residuals divided by the fill to allow for the missing data. This method is appropriate when comparing incoherent noise levels.

When there are zero-valued residuals in the time series - owing to the presence of gaps - power will be re-distributed into sidebands characteristic of the Fourier transform of the resulting window function. The power in the central peak will therefore be reduced. The equivalent sine-wave scaling scheme compensates for this, by dividing again by the fill, to give an estimate of the power in the oscillatory signal that is independent of the window function. This method is appropriate for estimating the power in a given signal.

5.2. Window function

The window function can usefully be characterized by generating an artificial spectrum in which the real data in a time series are replaced by a single sinusoid of a frequency such that it appears in a single bin of the power spectrum. There is generally a dominant 24-hour period due to the regular interruption of the data by night, which gives rise to a series of sidelobes associated with each peak, spaced at 1/day or approximately 11.57 tex2html_wrap_inline1306Hz. In the analysis of the low-tex2html_wrap_inline1216 p-mode spectrum, where the modes occur in pairs separated by about 5 to 15 tex2html_wrap_inline1306Hz, this can be a serious problem, as one peak may overlap with the sidelobe of its neighbour. With a multi-station network, although some remnant of the twenty-four hour period persists, the breaks are more random and the 1/day sidelobes are greatly reduced in size. Figure 3 (click here) shows typical window function transforms for a one-station and a six-station time series. In our best six-station spectra, the sidelobes have all but disappeared into the window-function noise.

  figure371
Figure 3: Window function power spectrum for May through June 1996. Here, a sine wave with period 320 s and amplitude tex2html_wrap_inline1318 has been forced through each window function. Top: Tenerife only, fill 32 per cent. Bottom: Six stations, fill 78 per cent. The ordinate is equivalent sine-wave scaled

In addition to the discrete sidelobes, the window function contains background noise and some power in the wings of the central peak. The broadening of the central peak, which is particularly severe if there are long gaps in the data, can corrupt the measurement of mode strengths and must be corrected for when comparing mode amplitudes between spectra. This correction has been discussed in detail by Elsworth et al. (1993). The correction applied for 1/day sidelobes is discussed below under the heading of "sidelobe fitting''.

5.3. Power spectrum statistics and smoothing

The ragged appearance of the power spectrum of the modes arises from the stochastically excited nature of the oscillations and the fact that each Fourier spectrum is only an estimate, or realization, of the true spectrum, where the power in each bin consists of the underlying spectrum multiplied by a random number distributed as chi-squared with two degrees of freedom (Duvall & Harvey 1986; Anderson, Duvall & Jefferies 1990). If a large number of such realizations were averaged together, a better estimate of the true spectrum would be obtained, with errors tending to a Gaussian distribution. However, unlike observers studying the high-degree modes by resolving the image of the Sun, we have only a few dozen modes and two to three dozen spectra of reasonable quality to work with. In very early work the repeating pattern of the spectrum made it useful to average together corresponding portions of the spectrum at different frequencies. However, more detailed examination of well-resolved spectra reveal that the properties (including the repeat interval) of the spectrum vary considerably with frequency: the averaging of different frequencies is therefore of limited use for accurate work. Moreover, the properties of the spectrum vary with solar activity (e.g., Libbrecht & Woodard 1990, 1991; Rhodes et al. 1993; Elsworth et al. 1993, 1994b; Regulo et al. 1994): averages of many spectra from different epochs obscure this information and must be treated with caution. We therefore work mainly with individual spectra, or occasionally with the mean of two spectra close together in time.

As discussed by Anderson, Duvall & Jefferies, one approach to the treatment of individual spectra is to fit the model parameters to the peaks using a "maximum likelihood'' method with the appropriate statistics, rather than normal least-squares fitting, which assumes a Gaussian distribution of the data points about the model. We have recently started to use this method in studies of high-resolution spectra. In our earlier work (Elsworth et al. 1990a,b, 1991, 1993, 1994), however, we have adopted a different approach. The power spectra are first smoothed using a weighted moving mean (Pollard 1977). Pairs of peaks are then fitted to a model function using non-linear least-squares methods. A frequency range of about 50 tex2html_wrap_inline1306Hz (about 250 frequency bins) is used to fit an tex2html_wrap_inline1148, tex2html_wrap_inline1328 pair, and a range of 70 tex2html_wrap_inline1306Hz for an tex2html_wrap_inline1332, tex2html_wrap_inline1334 pair. As was discussed in our publication on mode strengths (Elsworth et al. 1993), the smoothing does introduce systematic errors in the height and width of the fitted peaks, but these can be corrected for. We have attempted to deconvolve the smoothing function by using the same smoothing on the model function as was used on the data, but this was found not to be effective. We have found empirically that the smoothing function chosen has no effect on a Lorentzian peak profile until the frequency range over which the smoothing is carried out exceeds the full width at half maximum of the peak. Thereafter the width of the peak is increased by the addition in quadrature of the effective width of the smoothing function, which is slightly less than one third of the full range of the function, while its height decreases.

We have carried out numerous tests in which artificial spectra, generated from known parameters, were smoothed and fitted in the same way as the real data. Although individual realizations behave in very different ways, the effect of the smoothing on the average parameters can be predicted from the behaviour of a simple Lorentzian peak under the same amount of smoothing, provided that the smoothing is not performed over too wide a range compared with the peak width. In practice, for two-month spectra, we use between 7 and 49 bins (about 1.5 to 10 tex2html_wrap_inline1306Hz), using the heavier smoothing for the wider modes at the high-frequency end of the spectrum. The fitted peak parameters are corrected for the effect of smoothing as follows. The width obtained from the fit is first corrected for the effect of smoothing by subtracting in quadrature the effective width of the smoothing function. A Lorentzian peak of this corrected width and known height is then smoothed and fitted in the same way as the real data, in order to estimate the correction to be applied to the peak height. By comparing the results from the same spectrum fitted with different amounts of smoothing, we have verified that this procedure satisfactorily removes the systematic effects of smoothing.

5.4. Sidelobe fitting

The model function used to fit a pair of neighbouring peaks consists of two Lorentzian profiles, for which the adjustable parameters are the height h, the full width at half-maximum tex2html_wrap_inline1340, and the central frequency tex2html_wrap_inline1342, on a flat background of which the size b is also an adjustable parameter. The power, tex2html_wrap_inline1346, in each single Lorentzian component, is given by


equation389

Associated with each peak is a pair of symmetrical sidelobes at 11.57 tex2html_wrap_inline1306Hz on either side of the central peak. The sidelobes are assumed to have the same width and Lorentzian shape as the central peak, but their height relative to the main peak is left as an adjustable parameter for each peak. The relative height of the sidelobes can vary significantly from one peak to another. This is partly because of interference with neighbouring modes and noise, but we believe that a real effect is also involved: the amplitude of individual modes can vary dramatically from one week to the next, so that if a particular mode happens to have been strong for a period when only one or two stations were producing good data, it will appear to have stronger sidelobes than one that was strong when the network was running well. A more rigorous method of allowing for the sidelobes would be to convolve the model function with the known window function spectrum at each step of the minimization. This form of deconvolution would be much more expensive in computing time, and makes no allowance for the sidelobe variability discussed above.

5.5. Maximum-likelihood estimation and rotational splitting

A commonly used method for extracting mode parameters from power spectra is to use a maximum-likelihood estimator that takes into account the exponential distribution of the data about the underlying spectrum. As discussed above, we have found our own smoothing and least-squares fitting methods adequate for obtaining frequency measurements from 2-month power spectra. However, in most of this analysis we have neglected the rotational splitting of the modes of tex2html_wrap_inline1350 and treated all modes as single Lorentzian peaks. As the rotational splitting is not clearly resolved - particularly in 2-month spectra - for frequencies above about 2.2 mHz, and there is neither strong empirical evidence for, nor any theoretical reason to expect (at least from a simple model) asymmetries in the amplitudes of splitting components, this seemed to us a reasonable approximation. However, when we wish to analyse higher-resolution spectra to derive the rotational splitting of lower-frequency modes, the smoothing method becomes less appropriate. As discussed in Elsworth et al. (1995a), we have therefore chosen to use maximum-likelihood estimation on power spectra covering periods up to, for example, 16 months, with model functions that explicitly include the splitting components at the theoretically expected amplitude ratios calculated by Christensen-Dalsgaard (1989). To satisfy the statistical requirements of some aspects of the fitting procedure, we vary the logarithms of some of the parameters, rather than the parameters themselves.

The calculation of formal errors for this kind of fit, where neither the parameters nor (because of incomplete fill) the data points are completely independent, is complicated. We have in the past chosen to use the scatter of many determinations from independent spectra to estimate the errors of the peak parameters. This becomes more difficult with long spectra, of which we have only a few examples.

A thorough discussion of the errors in BiSON data can be found in Chaplin et al. (1996e), in which we outline our reservations regarding the use of formal uncertainties.


next previous
Up: Techniques used in

Copyright by the European Southern Observatory (ESO)
web@ed-phys.fr