APPLICATIONS OF THE DFT/CLEAN TECHNIQUE TO SOLAR TIME SERIES PATRICK C. CRANE Interferometrics Inc., 14120 Parke Long Court, Suite 103, Chantilly, VA 20151, U.S.A. (e-mail:
[email protected])
(Received 11 December 2000; accepted 10 August 2001)
Abstract. A new technique of Fourier analysis, DFT/CLEAN, has been adapted for the study of solar time series. The technique was developed by Roberts and his collaborators (1987, 1994) to address the limitations of other techniques of Fourier analysis and the shortcomings of many astronomical time series. The utility of the technique is illustrated with several applications to solar time series.
1. Introduction Fourier analysis of solar time series has been an active area of research at least since Ward and Shapiro (1962), and new applications are still being found – for example, the study of the center-to-limb behavior of the solar contrast (Crane et al., 2001). The state of our knowledge of solar variability has been reviewed by Lean (1997). The time series for a measure of solar activity, such as the international sunspot number, R (Figure 1), exhibits obvious systematic variations on decadal and daily time scales. The former are dominated by a pseudocycle with a period of ≈ 11 years with individual solar cycles varying greatly in length and amplitude. The latter are dominated by ≈ 27-day quasi-periodicities imposed by solar rotation but the amplitudes and phases are not stable because of the complications introduced by the birth, evolution, and death of solar active regions. These quasi-periodicities are well established but variations also occur on intermediate time scales – e.g., 100–1000 days – which are difficult to characterize and to associate with a physical origin (Lean and Brueckner, 1989; Pap, Tobiska, and Bouwer, 1990; Oliver, Ballester, and Baudin, 1998). The difficulties in analysis arise both from the characteristics of the solar time series themselves and from the limitations of the techniques of Fourier analysis used, as discussed in Sections 2 and 3, respectively. Section 4 describes the DFT/CLEAN technique of Fourier analysis. Section 5 illustrates the utility of the technique with examples of its application to solar time series. Mathematical details of the technique are presented in the Appendix.
Solar Physics 203: 381–408, 2001. © 2001 Kluwer Academic Publishers. Printed in the Netherlands.
382
P. C. CRANE
Figure 1. Daily international sunspot numbers, R, during solar cycles 21 and 22.
TABLE I Selected measurements of solar activity Measurement
Abbreviation
Start
URL
International Sunspot Number Fe XIV Coronal Index Ca II-K Plage Index Solar Radio Flux Density Mount Wilson Sunspot Index Mount Wilson Magnetic Plage Strength Index He I 1083-nm Equivalent-Width Index SBUV Mg II Core-to-Wing Index UARS/SUSIM Ultraviolet Irradiances Global Oscillation Network Group
R CI PI F10.7 cm MWSI MPSI He I EW Mg II SUSIM GONG
1818 1939 1942 1947 1970 1970 1974 1978 1991 1995
1 2 3 4 5 5 6 7 8 9
URLs: 1. http://web.ngdc.noaa.gov/stp/SOLAR/SSN/ssn.html. 2. ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/CORONA/INDEX/. 3. ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/CALCIUM/INDEX/. 4. http://web.ngdc.noaa.gov/stp/SOLAR/FLUX/flux.html. 5. http://www.astro.ucla.edu/obs/150_data.html. 6. http://www.nso.noao.edu/nsokp/dataarch.htm. 7. http://ssbuv.gsfc.nasa.gov/solar.html. 8. http://wwwsolar.nrl.navy.mil/susim_uars.html. 9. http://www.gong.noao.edu/.
DFT/CLEAN TECHNIQUE FOR SOLAR TIME SERIES
383
2. Solar Data Sets Table I lists several long-running measures of solar activity, including the GONG helioseismological experiment for comparison. The corresponding data sets share several characteristics that complicate the application of Fourier analysis to their study: Nonstationarity. Solar activity is not a stationary or temporally homogeneous process, perhaps not even stochastic; it is best described as quasi-stationary and quasi-periodic. The behavior depends upon the time interval under study, with no repetition seen for the longest intervals yet studied; for example, variations in the amplitude and period of the solar cycle; the disappearance of ≈ 27-day variations near solar minimum when no sunspots are visible; and the great reduction in activity during the Maunder Minimum (circa 1645–1715; Eddy, 1976). Analysis of a section of a time series only applies locally, but much study is devoted to the differences between sections of solar time series. Nonnegativity. Most solar indices are measures of brightness, irradiance, area, or equivalent width – all constrained to be greater than or equal to zero. Most indices, fortunately, do not approach zero. For those that do, notably R, analysis such as interpolation or extrapolation may produce negative values. Sampling. The measures listed are nominally available on a daily basis, so the sampling interval is one day. Therefore, the highest-frequency component that can be measured properly has a frequency of 0.5 day−1 which is the Nyquist critical frequency (Press et al., 1992). The GONG data are uniformly sampled at 1-min intervals, and the Nyquist frequency is 8.3 mHz. Size. Since the various measures nominally are obtained daily, the largest data sets include tens of thousands of entries. Except for GONG, the longest time series is for R which comprises approximately 67 000 days (1818–2001). Among the listed measures, even the shortest interval, 1991–2001, includes more than 3000 days. The GONG time series includes 1 814 400 samples. Missing data. Many measurements are or have been based upon observations made with a single instrument, and weather or mechanical problems may prevent a useful measurement on a given day or even on several consecutive days. For example, during 1818, the first year of R, the index is available for only 213 days. Also, after the first measurement of F10.7cm on 14 February 1947 during solar cycle 18, the index is available for only 80% of the subsequent days during that cycle. Even the GONG time series is about 18% incomplete. Irregularly sampled data. The case of missing data is a subset of the more general case of irregularly sampled data. Irregular sampling may arise, as for measurements
384
P. C. CRANE
from spacecraft, because of the limits placed on solar visibility by orbital dynamics or avoidance of the South-Atlantic Anomaly. Another example is F10.7cm , the time of whose measurement changed from 17:00 UT to 19:00 UT when the site of the measurements moved from Ottawa, Ontario to Penticton, British Columbia on 1 June 1991, and whose values are replaced by adjacent measurements in the case of bad weather or solar flares. The analysis of such data long has been of concern, starting with Deeming (1975). Variations on all time scales. The time series are dominated by variations on decadal (≈ 11-year solar pseudocycle) and daily time scales (≈ 27-day solar rotation). But longer, intermediate, and shorter (≈ 13.5 days) variations all are present. Especially when comparing short-term variations during different parts of the solar cycle, the longer-term variations may introduce artifacts. Instantaneous (non-bandwidth-limited) measurements. The measurements upon which most indices are based are effectively instantaneous; the minute required for the image used to count sunspots or to measure a solar radio flux density is very much shorter than the daily sampling interval. Some indices are based upon multiple instantaneous measurements throughout the day, as from a spacecraft (UARS/SOLSTICE) or from multiple observatories around the globe; daily averages minimize the contribution from impulsive behavior such as solar flares. Or measurements that exhibit impulsive behavior have been replaced with nearby or interpolated values, as is done for F10.7cm . The GONG samples are integrated for, and sampled at, one-minute intervals and are bandwidth-limited. Aliasing. Instantaneous measurements may see variations with frequencies greater than the nominal Nyquist frequency. If present, such features will be aliased (falsely translated) to frequencies of interest below the Nyquist frequency. Irregularly sampled data will include some measurements more closely spaced than the nominal sampling interval and are less susceptible to aliasing. The bandwidth-limited GONG data are not subject to aliasing. Errors. The measurement errors in most of the solar time series listed in Table I – R, for example – simply have not been characterized. As discussed in the Appendix, sometimes they may be characterized as constant or proportional or both: for F10.7cm , the estimated error is the maximum of one solar flux unit and 1% of the measurement (Tapping and Charrois, 1994; Tapping, 2000). Or they may be very complicated, especially when examined in detail, as for the SUSIM and SOLSTICE instruments on the Upper Atmosphere Research Satellite (Woods et al., 1996).
DFT/CLEAN TECHNIQUE FOR SOLAR TIME SERIES
385
3. Traditional Analyses Studies of the long-running time series of measures of solar activity (e.g., R and F10.7cm ) generally utilize three techniques summarized in Section 3.1: correlation analysis and power-spectrum analysis (cf., Heath, Repoff, and Donnelly, 1984; Lean and Repoff, 1987; Bouwer, 1992) and modified-periodogram analysis (Horne and Baliunas, 1986). General references include Blackman and Tukey (1958), Davenport and Root (1958), Bracewell (1965), and Press et al. (1992). As discussed in Section 3.2, these techniques impose certain requirements or structures on the problem that may limit the analysis. 3.1. T ECHNIQUES The autocorrelation function and the power spectrum are a Fourier-transform pair – i.e., they are related through the Fourier transform and its inverse. As such, while containing exactly the same information, they display it from complementary perspectives that together may reveal complex behavior that otherwise might be overlooked. Correlation analysis. A correlation function measures the degree of linear correlation between two time series as a function of time delay between the two: Rxy (τ ) = x(t)y(t + τ )/x(t)y(t).
(1)
In other words, it measures the degree of resemblance between the two as time passes. If the time series are identical, the result is the autocorrelation function which is symmetric about a maximum at zero delay. Otherwise, the result is the cross-correlation function which is not necessarily symmetric. Correlation functions can be calculated only when the time series are regularly sampled but the definition can be generalized to handle missing data (cf., Heath, Repoff, and Donnelly, 1984) or calculated from power spectra (Scargle, 1989). Applications to solar time series can be found in Donnelly (1987). Power spectrum. The power spectrum gives the distribution in frequency of the power in a time series and routinely is calculated as the Fourier transform of the corresponding autocorrelation function: ∞ Pxx (f ) =
Rxx (τ ) exp(−i2πf τ ) dt.
(2)
−∞
As discussed below, the calculated power spectrum is the convolution of the true power spectrum and a spectral window.
386
P. C. CRANE
Modified periodogram. An irregularly sampled time series can be studied using an alternative approach to estimating the power spectrum. The modified periodogram (Lomb, 1976; Scargle, 1982; Horne and Baliunas, 1986) has been used extensively to study periodicities in irregularly sampled time series (e.g., Lean and Brueckner, 1989). Again, as discussed below and like the power spectrum, the periodogram is the convolution of the true periodogram and a spectral window (Rinehart et al., 2000). 3.2. L IMITATIONS Nonstationarity. A correlation function, power spectrum, or modified periodogram for nonstationary solar activity is a valid representation of the time interval analyzed, but it does not describe other time intervals. But, as noted above, the study of solar time series is devoted as much to the differences, as to the similarities, between time intervals. Deconvolution. Because of the finite length and discrete nature of any practical time series, the calculated power spectrum (or periodogram) is the convolution of the true power spectrum and a spectral window. The spectral window is given by the Fourier transform of the sampling function in which the individual elements of the autocorrelation function are replaced with unity. For uniform regular sampling, the spectral window is a sinc function whose width is inversely proportional to the length of the autocorrelation function. An additional windowing function (Blackman and Tukey, 1958) can be applied to the autocorrelation function to reduce the sidelobes of the spectral window with a loss of spectral resolution, relative to the sinc function. However, the sidelobes cannot be reduced to zero and still will extend to infinity; consequently, adjacent features will interfere with each other and weak features may be obscured by the sidelobes of strong distant features. Furthermore, because the power spectrum, although symmetric, is intrinsically two-sided, features at negative frequencies can interfere with those at positive frequencies and vice versa. In some cases a window function may be useful as an alternative to deconvolution; see, for example, the discussion of multitaper techniques by Percival and Walden (1993). The tradeoff between sidelobe level and spectral resolution cannot be avoided. Deconvolution is the process of undoing the effects of the convolution with the spectral window to estimate the true power spectrum; little work on the deconvolution of solar power spectra appears in the literature. The modified periodogram also suffers from spectral leakage; this has been treated by (repeatedly if necessary) fitting and subtracting a sinusoid in the time domain and recalculating the periodogram (Horne and Baliunas, 1986), but more recently, versions of the CLEAN algorithm have been used (Rinehart et al., 2000; Zhao, Bower, and Goss, 2001). The number of operations involved in the deconvolution of N data points typically scales as N 2 per iteration.
DFT/CLEAN TECHNIQUE FOR SOLAR TIME SERIES
387
Fast Fourier transform (FFT) or direct Fourier transform (DFT). Calculation of a Fourier transform (as in calculating the power spectrum from the autocorrelation function) is computationally intensive. The simple, brute-force implementation of the algorithm is the direct Fourier transform (DFT) and the number of operations scales as N 2 . However, the algorithm can be rewritten as a fast Fourier transform or FFT (Press et al., 1992) which for values of N that are powers of 2 requires only of order N log2 N operations; slower FFT implementations exist when N has small prime factors greater than 2. The FFT is greatly superior in speed to the DFT and, consequently, almost always has been the implementation of choice. But the FFT imposes a certain structure on the calculation: the input time series must be regularly sampled; the input time series and the output Fourier transform have the same number of values; implicit in the algorithm is the assumption that both input and output are infinitely periodic. The input time series can be padded with extra zeros to reach an optimum value of N and missing values can be replaced with zeros. Irregular samples can be assigned to a regular grid by interpolation, which has been done in a few cases (Pap, Tobiska, and Bouwer, 1990; Bouwer, 1992) but is not recommended by Press et al. (1992), or by convolution, with a loss of coherence and signal at higher frequencies, as is done in synthesis imaging in radio astronomy (Briggs, Schwab, and Sramek, 1999). Binning, smoothing, and detrending. Binning of large solar time series once was necessary to reduce the number of values to be transformed, even when using an FFT, to a manageable number like 512 or 1024 (cf., Hughes and Kesteven, 1981). Smoothed versions of solar time series are used for many purposes, such as determining the time of solar minimum, but one common use is for detrending time series. Since such smoothing (or filtering) could not be done in the frequency domain (because calculating the necessary large Fourier transforms was impractical), astronomers turned to alternative, ad hoc approaches in the time domain: (1) using running means of perhaps 365 days or 27 days, depending upon the time scales of interest (Pap, Tobiska, and Bouwer, 1990); or (2) subtracting an nth-order polynomial fit to the time series where values of n such as 4 or 12 have been used, again depending upon the time scales of interest (Pap, Tobiska, and Bouwer, 1990; Bouwer, 1992). Both approaches have shortcomings: An n-day (boxcar) running mean only partially discriminates against variations on time scales shorter than n days since its Fourier transform is a sinc function whose sidelobes extend to infinite frequency, introduces biases in the presence of nonzero second derivatives (Press et al., 1992), and does not handle irregularly sampled data. Subtracting a polynomial baseline results in undesirable artifacts such as overshooting and ringing (Gershenfeld, 1999) and erroneous spectral features, especially at low frequencies (Bouwer, 1992). Dynamic studies. Studies of the short-term, dynamic behavior of solar activity are necessary and of interest because of the quasi-stationary and quasi-periodic nature
388
P. C. CRANE
of such activity. Detrending the time series of interest generally is a prerequisite for such studies, with the limitations described above. Two approaches typically have been used to study the time dependence of short-term variations: complex demodulation (Lean and Repoff, 1987) which becomes cumbersome when used to study behavior at more than a few selected periods; and dynamic power spectra (Bouwer, 1992) or periodograms (Cebula and Deland, 1998) which examine behavior for a full range of periods. Both approaches are susceptible to spectral leakage.
4. The DFT/CLEAN Technique There is another approach to the calculation of a power spectrum. Instead of calculating the autocorrelation function and its Fourier transform or calculating the periodogram, the Fourier transform of the time series itself is calculated (Deeming, 1975). While the power spectrum can be calculated by squaring the amplitude, the Fourier transform has greater utility than the power spectrum because it is complex. The Fourier transform is intrinsically two-sided – i.e., extends to both negative and positive frequencies (but for a real time series is Hermitian). Because of the retention of the phases (unlike the power spectrum), the transformation from the time domain to the frequency domain is reversible, and manipulations (e.g., convolution and multiplication) that are difficult in the one domain are simple in the other. The most obvious application is to filter the Fourier transform in the frequency domain by applying a low-pass or bandpass filter, for example, to produce the desired smoothed or detrended time series. Such filters can be ideal rectangles in the frequency domain. The Fourier transform, whether calculated as a DFT or FFT, is the convolution of the true Fourier transform with a spectral window. As discussed above, deconvolution is necessary to remove the effects of spectral leakage. In 1987 Roberts, Lehár, and Dreher proposed that the CLEAN algorithm could be adapted for the deconvolution of the one-dimensional Fourier transform of a time series. CLEAN (Högbom, 1974; Schwarz, 1978, 1979) was developed by radio astronomers to deconvolve the two-dimensional images of the radio sky produced by synthesisimaging arrays like the Very Large Array radio telescope; the theoretical understanding of CLEAN is still relatively poor (Cornwell, Braun, and Briggs, 1999). The initial implementation combined an FFT and CLEAN (hence FFT/CLEAN), which could be used to analyze time series of regularly sampled data with missing samples. As faster computer hardware became available, Roberts and Lehár (1994) recognized that the technique could be made more general, capable of handling time series of irregularly sampled data, by replacing the FFT with a brute-force DFT (DFT/CLEAN). Depending upon the size of a time series and whether it is irregularly sampled, either FFT/CLEAN or DFT/CLEAN can be used. The implementation of the CLEAN algorithm performs an approximate deconvolution of the true Fourier transform from the calculated dirty Fourier transform.
DFT/CLEAN TECHNIQUE FOR SOLAR TIME SERIES
389
(Note that quantities calculated prior to CLEANing are referred to as ‘dirty’ and those subsequent, as ‘clean’). The final product is a set of discrete clean Fourier components which when convolved with the dirty spectral window reproduces the dirty transform and whose inverse Fourier transform reproduces the original time series. For a real time series, the Fourier transform and, consequently, the clean components are complex and Hermitian. The procedure used by CLEAN is to locate the peak amplitude in the dirty Fourier transform, scale the dirty spectral window by the (complex) value of the peak multiplied by a clean gain, shift the scaled dirty spectral window to the position of the peak, and subtract it from the dirty Fourier transform, producing a residual Fourier transform. The Hermitian property of the Fourier transform of a real time series can be used to limit the search window for peaks to, for example, positive frequencies; the positive-frequency clean component and its negative-frequency Hermitian conjugate are then subtracted simultaneously from the dirty Fourier transform. The procedure is repeated upon successive residual Fourier transforms until only noise remains. The clean component for each iteration is specified by the frequency of the peak and the complex value used to scale the dirty spectral window. For improved stability usually only a fraction of the peak value is used to scale the dirty spectral window; a clean gain of 0.1 often is used. And the number of points across the width of the dirty spectral window is chosen to exceed the Nyquist criterion by a factor of 3–5 (see Perley, 1999). Criteria for deciding when only noise remains – i.e., the stopping condition – are not definitive; the one used here is to stop when the product of the number of clean components and of the clean gain equals the number of points in the dirty Fourier transform (consequently, the total number of operations will scale as N 3 and the time required will be correspondingly greater than for the FFT or DFT). The goal, though, is to clean sufficiently deeply into the noise that all significant features are transferred to the clean components. But if pushed deeply enough, CLEAN eventually may diverge because of computational problems such as finite precision arithmetic (Cornwell, Braun, and Briggs, 1999). One particular instability of CLEAN is well known, but rare, in the analysis of synthesis images (Cornwell, Braun, and Briggs, 1999), but whether there is a counterpart in the analysis of time series is unknown. For cosmetic purposes the clean components may be convolved with a gaussian clean spectral window. Usually, the width of the gaussian is chosen to match that of the dirty spectral window. Finally, to preserve the noise level and retain any features not well represented by the clean components, the last residual Fourier transform is added to produce a clean Fourier transform, but the two terms will be characterized by different sidelobes. Unlike the inverse Fourier transforms of the dirty Fourier transform and of the clean components, the inverse transform of the clean Fourier transform will not reproduce the original time series. The CLEAN algorithm also may be applied to the deconvolution of power spectra and modified periodograms (Rinehart et al., 2000; Zhao, Bower, and Goss, 2001), although no phase information will be obtained.
390
P. C. CRANE
5. Applications to Solar Time Series The utility of the DFT/CLEAN technique for the study of solar time series is illustrated in this section with several applications. The primary application of the DFT/CLEAN technique is the Fourier analysis of an irregularly sampled time series, which is illustrated in Section 5.1. Other applications include resampling (Section 5.2), filtering (Section 5.3), dynamical studies (Section 5.4), super-resolution (Section 5.5), error estimation (Section 5.6), and, finally, in Section 5.6 the use of DFT/CLEAN for functional optimization. A mathematical discussion of the DFT/CLEAN technique is given in the Appendix. 5.1. DFT/CLEAN The steps in such an analysis are illustrated in Figure 2 for the Mg II core-to-wing index (Floyd et al., 1998) produced by the Solar Ultraviolet Spectral Irradiance Monitor (SUSIM, Brueckner et al., 1993) on the Upper Atmosphere Research Satellite (UARS). The various versions of the Mg II index have been found to be useful indices of solar activity (Heath and Schlesinger, 1986; Viereck and Puga, 1999). The SUSIM version is a good example of an irregularly sampled time series because the times of the daily medium-resolution (1.1 nm) observations from which it is derived vary throughout the day, for example, because of the exigencies of orbital dynamics or avoiding the South Atlantic Anomaly. The original time series, x(t), is shown in Figure 2(a) covering the time range 12 October 1991 to 16 September 1998 (UARS days 31–2562). The effects of the irregular sampling can be seen clearly in the spectral window, W , in Figure 2(c), shown for low frequencies. If the time series were regularly sampled, W would decrease as sin(x)/x (sinc function) but, as can be seen, the amplitudes and spacings of the sidelobe peaks of W vary irregularly. The most important test, however, is to do the deconvolution: The dirty Fourier transform (shown at low frequencies as X in Figure 2(e)) has been CLEANed (deconvolved) using the spectral window, W , in Figure 2(c) to produce the set of clean components, c, in Figure 2(d) for periods longer than 12.5 days. The differences, shown in Figure 2(b), between the original time series and that reconstructed from the clean components are smaller than one part in a million. The dirty and clean Fourier transforms are compared in Figure 2(e), which is restricted to periods longer than 12.5 days and which shows significant features at the longest periods and at periods of ≈ 27 and ≈ 13.5 days, which arise from real episodes of ≈ 13.5-day variations (Floyd et al., 1997; Crane, 1998). The most significant differences are at low frequencies where CLEAN has corrected for the near-in sidelobes of the strong long-period variations (which are a factor of three greater than displayed or ≈ 30 × 10−4 ) and between the significant features where CLEAN has reduced the level of the ‘fuzz’ by correcting for the distant sidelobes of the stronger features. The comparison of the dirty and clean Fourier transforms at high frequencies is shown in Figure 2(f), for the frequency
DFT/CLEAN TECHNIQUE FOR SOLAR TIME SERIES
391
Figure 2. Illustration of DFT/CLEAN technique applied to UARS/SUSIM Mg II core-to-wing index (version V19R3): (a) original time series, x, for 12 October 1991 to 16 September 1998 (UARS days 31-2562); (b) differences, x, between original time series and that reconstructed from clean components; (c) amplitudes of dirty spectral window, W (dotted line), and clean spectral window, C (solid line); (d) amplitudes of clean components, c, for periods > 12.5 days; (e) amplitudes of dirty Fourier transform, X (dotted line), and clean Fourier transform, Xc (solid line), for periods > 12.5 days; and (f) amplitudes of dirty Fourier transform, X (dotted line), and clean Fourier transform, Xc (solid line), at high frequencies. The technique was applied using f = 0.00005 days−1 , g = 0.1, and niter = 100 000.
392
P. C. CRANE
Figure 3. Examples of resampling: (a) UARS/SUSIM Mg II core-to-wing index (version V19R3) for 8 May 1992 to 16 August 1992 (UARS days 240 – 340) as measured (solid circles), calculated from NOAA-9 Mg II core-to-wing index (solid line), and missing values calculated from clean components (open diamonds); and (b) international sunspot number, R, for the same interval as measured (solid line), with same gaps as Mg II index (solid circles), and missing values calculated from clean components (open diamonds).
range [0.42, 0.5] days−1 and shows how the spectral leakage from the distant strong long-period variations dominated before CLEAN was applied – afterwards the average amplitude is a factor of three smaller, from 1.04 × 10−5 to 0.35 × 10−5 . 5.2. R ESAMPLING The DFT/CLEAN technique describes the time series in terms of an equivalent set of discrete complex clean (Fourier) components which, when added together and sampled at the same times, reproduce the original time series, as in Figure 2(b). When sampled at other times, this sum of clean components can be used to calculate a regularly sampled version of the time series. More significantly, similar calculations can be done to estimate missing samples, For example, not only is the time series for the SUSIM Mg II index irregularly sampled on a daily basis, it also includes intervals when the SUSIM was turned off and no measurements were obtained. The longest of these gaps covers the interval UARS days 266–312 (3 June–19 July 1992), except for days 304 and
DFT/CLEAN TECHNIQUE FOR SOLAR TIME SERIES
393
305. A longer interval covering UARS days 240–340 is shown in Figure 3(a) with the actual and estimated (from the clean components derived in Section 5.1) SUSIM measurements. The reliability of the reconstruction of the missing data can be judged by comparison with the similar measurements of the Mg II index by the NOAA-9 SBUV/2 instrument (Deland and Cebula, 1998) which overlap the SUSIM time series. The NOAA-9 values shown are based upon a linear regression between 1070 common measurements in the two time series (linear correlation coefficient of 0.99). While the agreement between the estimated SUSIM and the NOAA-9 measurements is not exact, the linear correlation coefficient is still high (0.89 for 41 measurements). If there were no independent measurements available, the reconstructed data would well represent the behavior of the Mg II index in this interval. For comparison, a time series of R was selected to correspond to the SUSIM Mg II time series, and has been processed identically and is presented in Figure 3(b). The actual values of R are presented for comparison with the estimated values. Surprisingly, the degree of agreement between the actual and estimated sunspot numbers is considerably less (linear correlation coefficient of 0.62 for 45 measurements) than it was between the Mg II indices. But the reconstructed data represent the actual data reasonably well and would be better than none at all. The DFT/CLEAN technique works well for the interpolation of individual or a few data; for regularly sampled data it resembles the numerical application of the sampling theorem described by Press et al. (1992). In the examples above, the technique successfully reconstructed long sequences of missing data. Because of the non-stationary nature of solar activity, such success is not assured – e.g., during transitions between ≈ 27-day and ≈ 13-day variations such as occurred near 27 September 1994 and 4 June 1995 (Crane, 1998). 5.3. F ILTERING As discussed in Section 3.2, the standard approaches to smoothing and detrending time series have serious shortcomings. The alternatives include digital filtering in the time or frequency domains. Usually, except for real-time applications, filtering in the frequency domain is favored (Press et al., 1992). Furthermore, digital filtering in the time domain is not simply applicable to irregularly sampled time series. The natural choice in the context of Fourier analysis is digital filtering in the frequency domain. Having decomposed a time series into a set of discrete complex Fourier components, as represented by a set of clean components, one can apply ideal filters in the frequency domain to obtain time series that have been smoothed or detrended, as desired, without the artifacts associated with the standard approaches. An ideal low-pass, high-pass, or bandpass filter is applied simply by selecting the clean components at the desired frequencies and calculating the inverse Fourier transform. The application of filtering to the SUSIM Mg II index is illustrated in Figure 4
394
P. C. CRANE
Figure 4. UARS/SUSIM Mg II core-to-wing index (version V19R3): (a) unfiltered (dots) and low-pass-filtered (periods > 90 days, solid line) time series and (b) bandpass-filtered (periods = 2 – 90 days) time series.
which shows three versions of the time series: unfiltered, smoothed with a lowpass filter (periods > 90 days), and detrended with a bandpass filter (periods = 2–90 days). 5.4. DYNAMIC STUDIES Because of the non-stationary nature of solar activity, one wants to examine how the short-term, ‘dynamical’ behavior of a solar time series changes from time interval to time interval, and a useful approach is to follow Bouwer (1992) but utilizing Fourier transforms rather than power spectra. Having obtained a detrended time series without artifacts, Fourier transforms are calculated for shorter intervals (smaller subsets) of a time series selected using a small sliding window. The series of Fourier transforms then exhibits the dynamical behavior of the original time series. An example based on Crane (1998) is shown in Figure 5. In this case
DFT/CLEAN TECHNIQUE FOR SOLAR TIME SERIES
395
Figure 5. Dynamic Fourier transforms for the international sunspot number, R, during solar cycle 20, calculated from the bandpass-filtered (periods = 2 – 90 days) time series. The transforms are presented as the amplitudes of (a) the dirty Fourier transforms, (b) the clean Fourier transforms, and (c) the clean components and show the behavior of the persistent periodicity at ≈ 27 days and an unusual episode of ≈ 13.5-day variations near mJD 40500. The values of the grey-scale contours increase from white to black, and the maxima are (a) 27.0, (b) 27.2, and (c) 19.9.
396
P. C. CRANE
Figure 6. Amplitudes of dynamic clean Fourier transforms for the international sunspot number, R, for intermediate periods, calculated from the bandpass-filtered (periods = 90 – 900 days) time series, with the low-pass-filtered time series (periods > 900 days) shown for comparison. The values of the grey-scale contours increase from white to black.
dynamic Fourier transforms have been calculated for a sliding 100-day window applied to a detrended time series (filtered to include periods = 2–90 days) of R during solar cycle 20. The three panels of Figure 5 present the amplitudes of (a) the dirty Fourier transforms, (b) the clean Fourier transforms, and (c) the clean components. The most prominent feature in each panel is a wavy structure at periods near 27 days running vertically from bottom to top. White and lighter shades of gray are more common in Figure 5(b) than in Figure 5(a) which shows the deleterious effects of the low-level clutter from the sidelobes of the dominant features. The clean components in Figure 5(c) present the clearest picture of what is happening, although no additional significant features are visible. The wavy structure, especially near solar maximum, shows beating between active regions appearing and disappearing at different longitudes. Also of interest is the extended interval of ≈ 13-day variations near mJD 40500 which corresponds to real peaks and valleys in the time series (Crane, 1998). A second example, again for R, during solar cycles 21 and 22, is shown in Figure 6. The dynamic Fourier transforms have been calculated using a sliding 1000-day window applied to a bandpass-filtered time series (including periods = 100–900 days). For comparison the corresponding low-pass-filtered time series (including periods > 900 days) is plotted alongside. The evanescent non-stationary
DFT/CLEAN TECHNIQUE FOR SOLAR TIME SERIES
397
nature of solar activity is evident: activity is low near solar minimum. Transient features with varying periods near 110, 160, and 300 days appear during both solar cycles although they do not repeat in detail. Similar results averaged over solar cycle 21 were reported by Lean and Brueckner (1989) but their analysis lacked the fine time resolution to see the dynamical behavior. 5.5. S UPER - RESOLUTION As described in Section 4, one of the final steps in the CLEAN algorithm is to convolve the clean components with a gaussian clean spectral window. As noted, the width of the gaussian is usually chosen to match that of the dirty spectral window. But other choices are possible. If a larger width is used, the clean Fourier transform will be smoothed, emphasizing larger scale features. If a smaller width is used, it may be possible to achieve super-resolution, to see detail finer than the inherent resolution limit. The results can be misleading and are accurate only for strong, isolated features (Fomalont, 1999); in Figure 5(c), for example, no smoothing was done and while no additional features are visible, the presentation is much less cluttered. For additional discussion see Cornwell, Brown, and Briggs (1999). 5.6. E RROR ESTIMATION As noted above, the measurement errors of solar time series are often complicated and/or poorly characterized, and the determination of the errors of the Fourier transform is problematic. Because solar time series are dominated by ≈ 27-day quasi-periodicities imposed by solar rotation and longer periodicities, their Fourier transforms at higher frequencies (> 0.1 days−1 , for example) are dominated by noise and an alternative ad hoc approach is feasible. As radio astronomers do for two-dimensional synthesis images of the radio sky, the errors can be estimated by examining signal-free parts of the Fourier transform. The root-mean-squared noise, σ , can be estimated at higher frequencies where day-to-day measurement errors and noise dominate; the amplitudes follow a Rayleigh distribution (see the Appendix) for which σ = (2/π )1/2A.
(3)
For example, for the Mg II index in Figure 2, the average amplitude of the clean Fourier transform at high frequencies, [0.42, 0.5] days−1 , was 0.35 × 10−5 ; the estimate of the root-mean-squared noise, σ , is 0.28 × 10−5 . This example also shows the importance of applying CLEAN before estimating the noise. 5.7. F UNCTIONAL OPTIMIZATION Fourier transforms can also be used to provide a figure-of-merit in the optimization of the definition of a function which is used to calculate the elements of a time
398
P. C. CRANE
Figure 7. Comparison of mid-resolution (1.1-nm) and high-resolution (0.24-nm) UARS/SUSIM spectra (version V19) in the vicinity of the Mg II feature at a wavelength of ≈ 280 nm based on 709-day averages in the interval UARS days 259 – 1036. The high-resolution spectrum shows the emission cores at the centers of the broad, overlapping k and h absorption features. The minimum of the mid-resolution Mg II absorption feature has a wavelength, λmin , as indicated; an integration window of width 2W is indicated by the solid horizontal line.
series. The time series must exhibit quasi-periodic or periodic behavior, and the figure-of-merit is the ‘signal-to-noise ratio’ (SNR) of the corresponding feature in the Fourier transform of the time series. Depending upon the sampling of the time series, the FFT or the DFT may be used; i.e., the latter will be necessary if the sampling is irregular. Deconvolution with CLEAN will be necessary if the Fourier transform is complicated and spectral leakage is important. But the important step is that Section 5.6 has shown how the noise can be estimated. An example is the optimization of the definition of a Mg II core-to-wing index (a similar but more extensive optimization has been performed for version V19R3 of the SUSIM Mg II core-to-wing index by Floyd et al. (1998).) This index is calculated from the SUSIM daily mid-resolution (1.1 nm) measurements of the Mg II absorption feature near 280 nm. At high resolution (0.01 nm, Hall and Anderson, 1988) the feature consists of broad overlapping absorption in the Mg II h
DFT/CLEAN TECHNIQUE FOR SOLAR TIME SERIES
399
and k lines, narrow emission cores, and even narrower absorption features on top of the cores. The background continuum originates in the solar photosphere and the hierarchy of absorption, emission, and absorption features originates ever higher in the chromosphere. At the best resolutions available with the SUSIM (e.g., 1.1 nm and 0.24 nm as shown in Figure 7), these features are blended together. One can calculate the ‘core’ value used in the core-to-wing ratio by integrating or averaging over a window of variable width centered on the minimum of the SUSIM spectrum; the size of the window is specified by its halfwidth, W . The greater the halfwidth, the more the noise will average out but the more the photospheric continuum will dilute the chromospheric signal of interest. In the case of the Mg II index, the signal of interest is the quasi-periodic variations associated with solar rotation (≈ 27 days). To find the optimum halfwidth, the time series is calculated for the core value corresponding to each value of the halfwidth. The DFT (since the sampling is irregular) is calculated for each value of the halfwidth (since this is only an example no deconvolution is done) and the maximum amplitude, A27 , in the frequency interval [0.035, 0.040] days−1 (periods ≈ 27 days) is found. The root-mean-squared noise, σ , is estimated using Equation (3) in the frequency interval [0.1, 0.5] days−1 which avoids periods of ≈ 13.5 days. The value of the halfwidth for which SNR = A27 /σ
(4)
is maximum is then the optimum choice. The SNR is shown in Figure 8 as a function of W for the interval UARS days 259–1036; a distinct maximum (≈ 43) occurs for W ≈ 0.3 nm. The calculation of the hundreds of DFTs involved in this example would have been impractical only a few years ago.
6. Summary Despite the nonstationary nature of solar time series, Fourier analysis is a powerful tool for their study, mostly of how their behaviors change with time. The most powerful tool is the Fourier transform; yet regardless of its many advantages, with a few exceptions (e.g., Pap, Tobiska, and Bouwer, 1990; Bouwer, 1992), most recent studies of solar time series have used, not the FFT or DFT, but the modified periodogram. The primary reason has been concern about the handling of missing and irregularly sampled data (Lean and Brueckner, 1989). At the same time only primitive methods have been used to deal with spectral leakage; it has not been treated as a deconvolution problem. Recent advances in the speed of computer hardware allow practical implementations of the DFT and CLEAN algorithms. The DFT handles irregularly sampled data and CLEAN addresses the deconvolution problem. And since irregular sampling leads to worse spectral leakage, the combined DFT/CLEAN technique is especially useful when studying such time
400
P. C. CRANE
Figure 8. The signal-to-noise ratio as a function of halfwidth, W , for the mid-resolution (1.1-nm) UARS/SUSIM spectra included in Figure 7. The peak signal-to-noise ratio of ≈ 43 occurs for a halfwidth, W , of ≈ 0.3 nm.
series, while the FFT/CLEAN technique is advantageous for the study of regularly sampled time series. The DFT/CLEAN technique has been shown to be quite useful for the study of solar time series. Because it retains phases the technique has advantages over other techniques of Fourier analysis used in earlier studies and overcomes many of the shortcomings of the time series themselves. It allows unified and consistent analyses of the behaviors of the time series on all time scales without the shortcomings of the standard approaches to smoothing and detrending. And for the many time series whose measurement errors are complicated and/or poorly characterized, it can estimate the errors. The ability to use the DFT for irregularly sampled time series opens new avenues of study for such time series, as in functional optimization. The most important lesson is that spectral leakage is unavoidable in the Fourier analysis of time series. For some problems that may be satisfactory but, in general, deconvolution is recommended. The deconvolution algorithm CLEAN discussed
DFT/CLEAN TECHNIQUE FOR SOLAR TIME SERIES
401
here has been adapted from radio astronomy, and one can look there for other algorithms as well, including hybrids (Cornwell, Brown, and Briggs, 1999).
Acknowledgement This work was supported by NASA-Defense Purchase Requests S14789D and S10108X.
Appendix. Implementation of the DFT/CLEAN Technique The DFT/CLEAN technique combines the DFT and CLEAN algorithms described in Sections A.1 and A.2. The combination of the two algorithms is described in Section A.3. Probability distributions and the propagation of errors are discussed in Sections A.4 and A.5. A.1. T HE DIRECT F OURIER TRANSFORM (DFT) The DFT/CLEAN technique is applied to a time series of N measurements: xi = x(ti ), i = 0, . . . , N − 1.
(A1)
The simple, brute-force direct Fourier transform (DFT) of the time series is given by X(f ) =
1 N−1 xi exp(−2π if ti ), N i=0
(A2)
which can be written as the sum of a real and an imaginary term, X(f ) = R(f ) − iI (f ),
(A3)
or in terms of amplitude and phase, X(f ) = A(f ) exp[−2π iφ(f )],
(A4)
where A(f ) = [R 2 (f ) + I 2 (f )]1/2 ,
(A5a)
φ(f ) = tan−1 [I (f )/R(f )].
(A5b)
and
For a real time series, the real and imaginary terms are given by the cosine and sine transforms of xi :
402
P. C. CRANE
1 N−1 R(f ) = xi cos(2πf ti ), N i=0
(A6a)
and I (f ) =
1 N−1 xi sin(2πf ti ). N i=0
(A6b)
In practice, the two-sided DFT can be calculated on a uniform grid of 2n + 1 frequencies fj = j f, j = −n, . . . , −1, 0, 1, . . . , n,
(A7)
where the maximum frequency satisfies fmax = nf ≥ fc .
(A8)
For a regularly sampled time series (interval t), the Nyquist (critical) frequency, fc , is given by fc =
1 . 2t
(A9)
For an irregularly sampled time series, the appropriate time interval that determines the Nyquist frequency is the minimum separation between samples: fc =
1 . 2tmin
(A10)
For a bandwidth-limited signal the Nyquist frequency is the maximum frequency for which the time series satisfies the sampling theorem. The length of the time series determines the longest period (lowest non-zero frequency) that can be measured so that the frequency interval should be f ≤
1 2(tN−1 − t0 )
.
(A11)
The two-sided Fourier transform is the most convenient form to use for the deconvolution because the sidelobes of features at negative frequencies extend to positive frequencies (and vice versa) and the two sides must be deconvolved jointly. For a real time series the Hermitian property of the transform can be used to calculate simply and quickly the negative-frequency side of the transform from the positive-frequency side (i.e., X(−f ) = X ∗ (f )). The last step before proceeding with the deconvolution is to calculate the spectral window: 1 N−1 exp(−2π if ti ). W (f ) = N i=0
(A12)
DFT/CLEAN TECHNIQUE FOR SOLAR TIME SERIES
403
The spectral window, however, must be calculated to ± 2fmax in order that the effects of a feature in the transform at −fmax can be removed from a feature at +fmax (and vice versa); i.e., W is calculated for 4n + 1 frequencies fk = kf, k = −2n, . . . , −1, 0, 1, . . . , 2n.
(A13)
If the time series is sampled symmetrically about the mean time, ti , subtracting the mean time will result in a real W . For a regularly sampled time series for which the fast Fourier transform can be used, the inequalities in Equations (A8) and (A11) are equalities. A.2. A CLEAN ITERATION Starting with the (i − 1)th residual Fourier transform, Ri−1 , and the spectral window, W , the CLEAN algorithm calculates the ith residual Fourier transform, Ri , according to the following procedure: (1) Search Ri−1 (f ) for the maximum value of |Ri−1 |; determine the corresponding values (Ri−1 )p and fp . (2) Calculate the ith residual Fourier transform, Ri , by shifting, scaling, and subtracting a copy of the spectral window, W , from Ri−1 . In the present case, X is Hermitian, and one need search only non-negative frequencies and subtract simultaneously at positive and negative frequencies: Ri (f ) = Ri−1 (f ) − g(Ri−1 )p W (f − fp ) − g(Ri−1 )∗p W (f + fp ), Ri (f ) = Ri−1 (f ) − g(Ri−1 )p W (f ),
fp = 0.
fp > 0, (A14a) (A14b)
The gain, g, must be in the range (0, 2) for convergence (Schwarz, 1979). Smaller values (≈ 0.1) allow a smoother approach to convergence in the presence of noise and for deeper deconvolutions. (4) Save the ith clean component, at both negative and positive frequencies (initially c ≡ 0): c(fp ) = c(fp ) + g(Ri−1 )p ,
fp > 0,
c(−fp ) = c(−fp ) + g(Ri−1 )∗p , c(0) = c(0) + g(Ri−1 )p ,
fp > 0,
fp = 0.
(A15a) (A15b) (A15c)
A.3. T HE DFT/CLEAN TECHNIQUE As noted earlier, X and W are referred to as the ‘dirty’ Fourier transform and spectral window, respectively, and the results of the deconvolution are identified as the ‘clean’ components (c), spectral window (C), and Fourier transform (Xc ).
404
P. C. CRANE
(1) Compute X and W for the input time series using the DFT. (Subtraction of the mean value, xi , from the data prior to calculating X will speed up convergence of the deconvolution if xi is large.) (2) The initial residual Fourier transform, R0 , is given by X: R0 = X.
(A16)
(3) Perform a CLEAN iteration. (4) Test against stopping criteria. If satisfied, proceed to step 5; otherwise return to step 3. Typical stopping criteria that are used, perhaps singly or in combination, are to stop after a specified number of interations or when the peak amplitude of the residual Fourier transform reaches a specified value. (5) Calculate the clean spectral window C(f ) = exp[−4 ln 2 (f/w)2 + 2π i(tsub − ti )],
(A17)
where the full width at half maximum, w, has been determined from the dirty spectral window. (6) Calculate the clean Fourier transform by adding the convolution of the clean components with the clean spectral window to the final residual Fourier transform: Xc (f ) = c(f ) ∗ C(f ) + Ri (f ).
(A18)
A.4. P ROBABILITY DISTRIBUTIONS If the errors in R and I have independent gaussian probability distributions with zero means and identical σ , the probability distribution for the amplitude is given by the Rice distribution (Rice, 1944, 1945; Vinokur, 1965; Papoulis, 1965); the last two references also derive the probability distribution for the phase. The properties of these distributions are discussed extensively in the radio-astronomical literature (cf., Moran, 1976; Thompson, Moran, and Swenson, 1986). In the presence of a signal the distributions are complicated, but in the absence of one they simplify to the Rayleigh and uniform distributions, respectively: p(A) = Aσ −2 exp[−A2 /(2σ 2 )],
A ≥ 0,
(A19a)
and p(φ) = (2π )−1 ,
0 ≤ φ < 2π.
(A19b)
A useful property of the Rayleigh distribution is that A = σ (π/2)1/2,
(A20)
which can be used to estimate σ from multiple measurements of A. The Rayleigh distribution is also the basis for significance testing.
DFT/CLEAN TECHNIQUE FOR SOLAR TIME SERIES
405
A.5. P ROPAGATION OF ERRORS The propagation of errors through the DFT/CLEAN technique must be considered separately for the DFT and CLEAN algorithms. Since the DFT is a linear operation, propagation of errors for the DFT is relatively straightforward and will be discussed below in more detail. The CLEAN algorithm, on the other hand, is highly nonlinear and the propagation of errors is unknown (Cornwell, Brown, and Briggs, 1999). In practice, for the deconvolution of two-dimensional synthesis images of the radio sky, the ad hoc assumption that CLEAN has little affect on the errors has proven satisfactory and, furthermore, the errors often are estimated by examining signalfree areas of the images. A.5.1. Propagation of Errors and the DFT Discussion of the propagation of errors most simply starts with the real and imaginary terms, R(f ) and I (f ). For a real time series, xi , with independent, randomly distributed errors with dispersions σi , the errors in R(f ) and I (f ) are given by σR2 (f )
1 N−1 = 2 σ 2 cos2 (2πf ti ), N i=0 i
(A21a)
1 N−1 σ 2 sin2 (2πf ti ). N 2 i=0 i
(A21b)
and σI2 (f ) =
If the σi are known, σR (f )and σI (f ) can be calculated and, in general, will be functions of frequency. If not, they can be estimated from the measurements themselves. A general discussion of measurement errors is beyond the scope of this paper (cf., Rabinovitch, 1995). Discussions of constant and proportional errors will be used to provide guidance. A.5.2. Weights and Windows If the dispersions σi , are known and not identical, the Fourier transform can be calculated using weights wi = (σi )−2 to minimize the squared errors: X(f ) =
1 N−1 wi xi exp(−2π if ti ), W
(A22a)
W (f ) =
1 N−1 wi exp(−2π if ti ), W i=0
(A22b)
σR2 (f ) =
1 N−1 wi cos2 (2πf ti ), W 2 i=0
(A22c)
406
P. C. CRANE
and σI2 (f )
1 N−1 = 2 wi sin2 (2πf ti ), W I =0
(A22d)
where W =
N−1
(A23)
wi .
i=0
A window function also can be applied by replacing wi with wi ti , where ti represents the window function. As described earlier, such functions are applied to shape the spectral window – i.e., to reduce the amplitudes of the sidelobes which worsens the spectral resolution – but obviously will degrade the errors. Degraded errors and resolution will reduce the ability of CLEAN to identify and discriminate between discrete Fourier components. In other circumstances, applying window functions is a powerful tool of spectral analysis (Percival and Walden, 1993), even for irregularly sampled time series. A.5.3. Constant Errors Constant errors are described by σi ≡ σ . For large N and non-pathological distributions of ti , N−1 i=0
cos2 (2πf ti ),
N−1
sin2 (2πf ti ) → N/2,
(A24)
i=0
and σR (f ), σI (f ) → (2N)−1/2 σ.
(A25)
A.5.4. Proportional Errors Proportional errors are described by σi = axi , where a is a constant. For large N and non-pathological distributions of ti , and at high frequencies, cos(2πf ti ), sin(2πf ti ) change more rapidly with ti than the xi , σR2 (f ), σI2 (f ) → a 2 (2N 2 )−1
N−1
xi 2 ,
(A26)
xi2 .
(A27)
i=0
which can be estimated by σR2 (f ), σI2 (f ) ≈ a 2 (2N 2 )−1
N−1 i=0
References Blackman, R. B. and Tukey, J. W.: 1958, The Measurement of Power Spectra, Dover, New York.
DFT/CLEAN TECHNIQUE FOR SOLAR TIME SERIES
407
Bouwer, S. D.: 1992, Solar Phys. 142, 365. Bracewell, R.: 1965, The Fourier Transform and Its Applications, McGraw-Hill, New York. Briggs, D. S., Schwab, F. R., and Sramek, R. A.: 1999, in G. B. Taylor, C. L. Carilli, and R. A. Perley (eds.), Synthesis Imaging in Radio Astronomy II, Astronomical Society of the Pacific, San Francisco, p. 127. Brueckner, G. E., Edlow, K. L., Floyd, L. E., Lean, J. L., and VanHoosier, M. E.: 1993, J. Geophys. Res. 98, 10695. Cebula, R. P. and DeLand, M. T.: 1998, Solar Phys. 177, 117. Cornwell, T., Braun, R., and Briggs, D. S.: 1999, in G. B. Taylor, C. L. Carilli, and R. A. Perley (eds.), Synthesis Imaging in Radio Astronomy II, Astronomical Society of the Pacific, San Francisco, p. 151. Crane, P. C.: 1998, Solar Phys. 177, 243. Crane, P. C., Floyd, L. E., Herring, L. C, Cook, J. W., and Prinz, D. K.: 2001, Astrophys. J. (submitted). Davenport, W. B. and Root, W. L.: 1958, An Introduction to the Theory of Random Signals and Noise, McGraw-Hill, New York. Deeming, T. J.: 1975, Astrophys. Space Sci. 36, 137. Deland, M. T. and Cebula, R. P.: 1998, Solar Phys. 177, 105. Donnelly, R. F.: 1987, Solar Phys. 109, 37. Eddy, J. A.: 1976, Science 192, 1189. Floyd, L., Brueckner, G., Crane, P., Prinz, P., and Herring, L.: 1997, in A. Wilson and B. Fleck (eds.), Proc. 31st ESLAB Symp., Correlated Phenomena at the Sun, in the Heliosphere and in Geospace, ESA SP-415, ESTEC, Noordwijk, The Netherlands, pp. 235 – 242. Floyd, L. E., Crane, P. C., Herring, L. C., and Prinz, D. K.: 1998, EOS Transactions 79, F10. Fomalont, E. B.: 1999, in G. B. Taylor, C. L. Carilli, and R. A. Perley (eds.), Synthesis Imaging in Radio Astronomy II, Astronomical Society of the Pacific, San Francisco, p. 301. Gershenfeld, N.: 1999, The Nature of Mathematical Modeling, Cambridge University Press, New York. Hall, L. A. and Anderson, G. P.: 1991, J. Geophys. Res. 96, 12927. Heath, D. F. and Schlesinger, B. M.:1986, J. Geophys. Res. 91, 8672. Heath, D. F., Repoff, T. P., and Donnelly, R. F.: 1984, NOAA Tech. Memo. ERL ARL-129, NOAA ERL, Boulder, Colorado. Högbom, J. A.: 1974, Astron. Astrophys. Suppl. 15, 417. Horne, J. H. and Baliunas, S. L.: 1986, Astrophys. J. 302, 757. Hughes, V. A. and Kesteven, M. J. L.: 1981, Solar Phys. 71, 259. Lean, J.: 1997, Ann. Rev. Astron. Astrophys. 35, 33. Lean, J. L. and Brueckner, G. E.: 1989, Astrophys. J. 337, 568. Lean, J. L. and Repoff, T. P.: 1987, J. Geophys. Res. 92, 5555. Lomb, N. R.: 1976, Astrophys. Space Sci. 39, 447. Moran, J. M.: 1976, in M. L. Meeks (ed.), Methods of Experimental Physics, Vol. 12, Astrophysics, Part C: Radio Observations, Academic Press, New York, p. 228. Oliver, R., Ballester, J. L., and Baudin, F.: 1998, Nature 394, 552. Pap, J., Tobiska, W. K., and Bouwer, S. D.: 1990, Solar Phys. 129, 165. Papoulis, A.: 1965, Probability, Random Variables and Stochastic Processes, McGraw-Hill, New York. Percival, D. B. and Walden, A. T.: 1993, Spectral Analysis for Physical Applications: Multitaper and Conventional Univariate Techniques, Cambridge University Press, New York. Perley, R. A.: 1999, in G. B. Taylor, C. L. Carilli, and R. A. Perley (eds.), Synthesis Imaging in Radio Astronomy II, Astronomical Society of the Pacific, San Francisco, p. 275. Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P.: 1992, Numerical Recipes in Fortran, 2nd edn., Cambridge University Press, New York.
408
P. C. CRANE
Rabinovitch, S.: 1995, Measurement Errors: Theory and Practice, AIP Press, New York. Rice, S. O.: 1944, Bell Syst. Tech. J. 23, 282. Rice, S. O.: 1945, Bell Syst. Tech. J. 24, 46. Rinehart, S. A., Hajian, A. R., Houck, J. R., and Terzian, Y.: 2000, Pub. Astron. Soc. Pacific 112, 977. Roberts, D. H. and Lehár, J.: 1994, private communication. Roberts, D. J., Lehár, J., and Dreher, J. W.: 1987, Astron. J. 93, 968. Scargle, J. D.: 1982, Astrophys. J. 263, 835. Scargle, J. D.: 1989, Astrophys. J. 343, 874. Schwarz, U. J.: 1978, Astron. Astrophys. 65, 345. Schwarz, U. J.: 1979, in C. van Schooneveld (ed.), Image Formation from Covariance Functions in Astronomy, D. Reidel, Dordrecht, The Netherlands, p. 261. Tapping, K. F.: 2000, private communication. Tapping, K. F. and Charrois, D. P.: 1994, Solar Phys. 150, 305. Thompson, A. R., Moran, J. M., and Swenson, G. W.: 1986, Interferometry and Synthesis in Radio Astronomy, John Wiley & Sons, New York. Viereck, R. A. and Puga, L. C.: 1999, J. Geophys. Res. 104, 9995. Vinokur, M.: 1965, Ann. Astrophys. 28, 412. Ward, F. and Shapiro, R.: 1962, J. Geophys. Res. 67, 541. Woods, T. N., Prinz, D. K., Rottman, G. J., London, J., Crane, P. C., Cebula, R. P., Hilsenrath, E., Brueckner, G. E., Andrews, M. D., White, O. R., VanHoosier, M. E., Floyd, L. E., Herring, L. C., Knapp, B. G., Pankratz, C. K., and Reiser, P. A.: 1996, J. Geophys. Res. 101, 9541. Zhao, J.-H., Bower, G. C., and Goss, W. M.: 2001, Astrophys. J. 547, L29.