Multiple and changing cycles of active stars
I. Methods of analysis and application to the solar cycles
Z. Kolláth  K. Oláh
Konkoly Observatory of the Hungarian Academy of Sciences, PO Box 67, 1525, Hungary
Received 6 November 2008 / Accepted 1 April 2009
Abstract
Context. Longterm observational data have information on the magnetic cycles of active stars and that of the Sun. The changes in the activity of our central star have basic effects on Earth, such as variations in the global climate, so that understanding the nature of these variations is extremely important.
Aims. The observed variations related to magnetic activity cannot be treated as stationary periodic variations, therefore methods like Fourier transform or different versions of periodograms only give partial information on the nature of the light variability. We demonstrate that timefrequency distributions provide useful tools for analysing the observations of active stars.
Methods. We tested and used different methods, such as shortterm Fourier transform, wavelet, and generalised timefrequency distributions, for analysing temporal variations in timescales of observational data.
Results. With test data we demonstrate that the observational noise has practically no effect on the determination in the longterm changes of timeseries observations of active stars. The rotational signal may modify the determined cycles, therefore it is advisable to remove it from the data. Wavelets are less powerful in recovering complex longterm changes than other distributions that are discussed. By applying our technique to the sunspot data we find a complicated, multiscale evolution in the solar activity.
Key words: Sun: activity  methods: data analysis
1 Introduction
It is well known that solar cycles cannot be represented by single periodic variations. Besides the 11 year quasiperiodicity (Schwabe 1843) of solar activity, lots of different rhythms have been discovered or at least suggested by different studies (see e.g. Frick et al. 1997; ForgácsDajka & Borkovits 2006). One remarkable and wellstudied periodicity of the solar variation was found by Gleissberg (1966) about 70 years ago. He described this variability as a ``longperiodic fluctuation of the sunspot numbers, the period being about seven or eight sunspot cycles''. Frick et al. (1997) derived the length of the Gleissberg cycle as 100 years from wavelet analysis of a 384 years long dataset (Hoyt et al. 1993). The nonstationarity in the longterm variation of active stars has already been found in active stars as well (Oláh et al. 2007).
Timefrequency representation of a nonstationary signal yields information about characteristics of the dataset in the timefrequency plane. Then transient events or frequency/amplitude changes can be easily monitored with these tools. In the literature, a frequently used method is the waveletmap and some other techniques based on wavelets. However, wavelet is only one of the several timefrequency representations, and other methods may give better results than wavelets. The key problem is the balance between temporal and frequency resolution. The frequencytime uncertainty principle gives strict constraints for the time duration () and the bandwidth ( ) of a signal: . These measures of broadness are determined by the standard deviation of frequency and time, and the uncertainty principle with these quantities applies to any timefrequency distribution. The linear transforms (such as wavelet or windowed Fouriertransform) have a limit for joint resolution of time and frequency for local structures; for example, in general it is not possible to determine the frequency precisely at high temporal resolution. Bilinear timefrequency transformation, called the kernel method (Cohen 1995) can break the above principle locally, but its cost is paid by artifacts at other frequencies.
While wavelets and other methods based on Fourier transforms are exclusively used to analyse stellar activity, more sophisticated methods have been available for decades: the Wigner distribution (Wigner 1932) and its extensions by Cohen (1995). These methods have also been proven to be useful in analysing of variable star data of chaotic origin or nonstationary behaviour (see Buchler et al. 1995; Kolláth & Buchler 2001). In this paper we present a method that can be used to study changes of stellar activity cycles in time using data with yearly gaps and with shortterm signals (rotational modulation).
As example we reanalysed the available sunspot and radio records with modern timefrequency methods to give a complete picture of the beating of the Sun in the period range from a few years to 200 years. In the companion paper (Oláh et al. 2009, Paper II) we apply the method described in the present paper to 20 active stars, and discuss the time variability of their cycles.
2 Timefrequency distributions
A simple but still powerful method is the ``shortterm Fourier transform'' (STFT hereafter). In this case the dataset is multiplied by a Gaussian window centred on a given epoch t_{0}. If this windowed part of the signal is then Fouriertransformed, the frequency spectrum of the data around the time t_{0} can be generated. By shifting the centre of the Gaussian window in time, a sequence of Fourier spectra is obtained, which gives the timefrequency representation of the data. If the sampling of t_{0} is dense enough, then a twodimensional map of the energy distribution on the timefrequency plane can be plotted.
We demonstrate this method on the yearly sunspot numbers  a dataset that is well known and analysed by several authors. Figure 1 displays the STFT of the solar data. The Fourier spectra of the individually windowed subsets of the data are also plotted with thick lines on the timefrequency distribution map. By changing the width of the Gaussian window, the balance between time and frequency resolution can be modified. All the locations on the map, of course, satisfy the uncertainty principle , and the ratio of time and frequency resolution is constant on the map.
Figure 1: Shortterm Fourier transform ( bottom) of the sunspot number time series ( top). Fouriertransforms of the subsets of the dataset are also plotted. The colour scale of relative variations in the distributions is presented in Fig. 7. 

Open with DEXTER 
While in this demonstration we used a timecentred Gaussian window to plot the curves of Fourier transforms, in the implementation of STFT, we use windowing in Fourier space:
(1) 
where denotes the Fourier transform of the signal. The inverse Fourier transform is only calculated for the positive frequencies. This method has several advantages. At low frequencies it reduces the distortion due to the components at negative frequencies; moreover, when is large, the transform reduces to the analytic function (see Buchler & Kolláth 2001). Depending on the functional form of , Eq. (1) reduces to STFT or a wavelet. With , where f_{0} is a constant frequency, the equivalent temporal window function of STFT is given by
(2) 
Usually it is a good choice to set f_{0} to the maximum or the central frequency of the studied period domain. Then provides a naturally well balanced temporal and frequency resolution. It can be seen from Eq. (2) that the half width of the timecentred windowing function is . Then, if depends on frequency, the width of the window also varies with frequency. Especially if then T(t,f) is equal to wavelet, with the Morlet kernel, since the width of the temporal windowing function decreases with 1/f.
Figure 2: Wavelet ( top and middle) and shortterm Fourier transform (STFT) distribution ( bottom) of the yearly sunspot numbers. 

Open with DEXTER 
In the standard Morlet wavelet is fixed to 1. The introduction of a free parameter () that controls the balance between frequency and temporal resolution has several advantages as already described by Kolláth & Szeidl (1993) and Soon et al. (1999).
As can be seen from the above definitions, a wavelet with Morletkernel is very similar to STFT; however, for wavelets the balance between temporal and frequency resolution depends on the frequency. This feature is demonstrated in the top panels of Fig. 2. The parameter (f_{0}) is defined in such a way, that for a given , the wavelet and STFT have the same frequency/time resolution at the frequency f_{0}. The top and middle panels of Fig. 2 show wavelets where the lower frequency part (top) and the highfrequency part (middle) was matched to the STFT, respectively. Resolution in time increases with frequency, since the same number of cycles is covered on a shorter timebase, at higher frequencies. The disadvantage of this representation is that harmonic structures, which display the same characteristics in STFT, are distorted. Figure 2 thus demonstrates that, while the timefrequency distribution in the lower panel clearly show the synchronised variation of the 11 year cycle with its double frequency harmonic, the wavelet destroys this parallel structure.
We used a different scaling of the amplitudes at different frequencies. In timefrequency distributions the lower amplitude parts are easily washed out by the contributions from power of other components. Far from our topics but similar in mathematical tools, it is well known in speech analysis that the obscured parts of the frequency distribution can be recovered by applying a preemphasis filter. A multilevel preemphasis Fourier filter is applied to visualise the whole studied period range equally. For the solar data relative to the highest amplitude 11 year periodicity (0.070.14 c/y), the signal in different frequency ranges were amplified by the following factors (): for 0.0 < f < 0.05, and above f=0.16 (f is given in c/y). A smooth change in is used in the transition regions. With this preemphasis filtering there is a factor of 3.5 difference in the amplitudes to get similar grey levels around 11 and 5.5 year periods.
The optimal balance of temporal and frequency resolution depends on the signal itself. A question arises as to whether it is possible to construct a distribution, where the resolution in both coordinates is optimal, but some other disadvantage is allowed. The answer is more general group of methods for signal processing introduced by Cohen (see e.g. Cohen 1995). The generalised timefrequency distribution is given by
C(t,f)  =  
(3) 
where s(t) is the analysed time series and the kernel of the distribution that determines the specific properties of the distribution. The simplest timefrequency distribution is the WignerVille transformation with . This distribution only works well for special datasets, otherwise the resulting map is heavily contaminated by cross terms of the different components. A biGaussian kernel provides a distribution that, in the limit of the kernel width, is equivalent to STFT. This socalled pseudoWigner distribution (PWD hereafter) is widely used in sound processing. For comparison purposes, we also use the ChoiWilliams distribution (CWD) (see Choi & Williams 1989), which applies an exponential kernel.
The top panel of Fig. 6 shows the PWD of the yearly sunspot numbers. It can be seen that the resolution in both coordinates is improved compared to the result with STFT displayed in Fig. 2 (bottom panel); however, there are artificial features visible on the PWD distribution that have no real meaning in the observed signal. The reason for these cross terms is the nonlinearity of the method; e.g., if two sinusoidal waves with different frequencies coexist in the signal, then false power arises between the two ridges representing the cycles.
3 Preprocessing observational data
Stellar observations in reality have lots of deficiencies, like observational noise and imperfect sampling. Then timefrequency methods cannot be applied directly to the observational data, and the analyses should be preceded by a sequence of preprocessing steps.
In the case of active stars, the situation is even more complicated, because different timescales of variability coexist, which interfere with the sampling in different ways. The rotational periods fall into the range of days to weeks, very close to the order of the sampling time. In contrast, longterm changes have a characteristic time above a year, so averaging and resampling on this level is needed to obtain information on these components. Thus, the preprocessing of the data requires extreme care.
To avoid false averages due to the undersampling of the rotational component, it is advisable to whiten all observations with the rotational period. Because the rotational component is modulated by e.g. the differential rotation, prewhitening cannot be performed simultaneously for the whole dataset.
The main problem is that the lightcurves are unevenly sampled, and they contain seasonal gaps. Although methods exist to calculate e.g. wavelets of this kind of data directly (like the adaptive wavelet), according to our previous experiments we prefer to interpolate the data when it is possible. When the sampling is denser than the dominant periodicities in the data, then a simple moving averaging can provide a continuous signal. Next, the averaged data can be smoothed with a cubic smoothing spline (Reinsch 1967). This operation provides a continuous function S(t) based on the discrete time series (
i=1,n) by minimising the integral of S''(t)^{2}, with the following constraint:
where , the smoothing factor, is a free parameter, which controls the smoothness of the spline. In theory its value should be chosen in the confidence interval corresponding to the left side of Eq. (4), thus should be of the same order of amplitude as the standard deviation of the observational noise. In practice, we usually manually check the resulting splinesmoothed curves to fine tune of the averaging and spline smoothing. As a final step, S(t) is resampled to generate an equally spaced time series.
4 Testing the methods
4.1 Artificial test data sets
We created artificial datasets to test the methods for a general signal representing the cycles of active stars. Signal ``I'' consists of a periodic signal with a constant frequency of 0.002 c/d, a component with linearly increasing frequency (chirp) with a modulation between f=0.0007 and f=0.0011 c/d, and a wave packet at f=0.0003 c/d with a Gaussian amplitude modulation: ). The amplitudes of the constant f=0.002 c/d signal, the chirp, and the wave packet are 0.02, 0.04, and 0.04, repectively.
In addition to the components of signal ``I'', signal ``II'' has a higher frequency variation, demonstrating the rotational period of active stars, and a uniformly distributed noise with amplitude comparable to the rotational component. To test the preprocessing and the timefrequency methods together, the test signal ``II'' was finalised by resampling the data with the time series of the observations of EI Eri (see Paper II) The data were averaged with a time base of 90 days, then the averages were spline smoothed with and resampled with days.
4.2 Observational noise and timefrequency distributions
It is well known that the Fourier transform has a very good noise rejection property in the case of periodic or quasiperiodic signals. The same is true for timefrequency distributions. If the observational noise and the signal are not correlated, then the method provides a good result even in the case of a large amount of noise. A nice example of this nature of timefrequency distributions is the analysis of the semiregular star V CVn (Buchler et al. 2004), in which two datasets are available for the same period of time. One of the records consists of visual magnitude estimates by amateur astronomers, while the other data are professional measurements by photoelectric photometry. It is clear that even the statistical properties of the noise sources are different in the two cases. Figure 2 in Buchler et al. (2004) clearly demonstrates that the timefrequency plots are almost fully independent of the choice of the dataset, so the effect of observational noise is negligible for signals with good SNR and sampling.
Figure 3: Timefrequency distribution of test data. Panel A): PWD of signal ``I''; panel B): PWD of signal ``I'' with white noise; panel C): PWD of signal ``II'' with averaging and spline interpolation; panel D): PWD of signal ``II'' after Fourier filtering of the rotational period in the independent observational seasons and then averaging and spline smoothing; panel E): STFT of the same signal as in panel D); panel F): CWD of the same signal as in panel D). 

Open with DEXTER 
We demonstrate this noise rejection property of timefrequency distributions on our first test time series. The top panels of Fig. 3 display the PWD of signal ``I'' (A) and the same signal modified with noise, with large relative temporal variance, added to it (B). The amplitude of the noise is seen in the top panel of the figure. While the noise has practically no effect on the lower frequencies, the f=0.002 c/d component is slightly distorted. The results with signal ``I'' are presented only for PWD, but we note that the other timefrequency methods (STFT, wavelet, CWD) have the same properties.
We have the sunspot observations for the sensitivity to noise, too. White or Gaussian noise with the same amplitude as the maximum sunspot number does not significantly distort the timefrequency plots.
4.3 Effect of uneven sampling and noise
Because of the averaging and spline interpolation, the effect of observational noise is mixed with the effects of sampling, so we performed our tests with illsampled data. The resulting timefrequency distributions obtained from the gapped signal ``II'' are displayed in the lower panels of Fig. 3, where, on panel C the rotational signal is included, while it is filtered out from the data in panels D, E, and F. Averaging and smoothing spline interpolation introduces a lowpass filtering of the signal, which can be compensated for by the preemphasis filtering. In these cases we selected for f> 0.0015. All the components above about 2 years of characteristic time are recovered well with these settings and processing, especially in the case where the rotation period is filtered out. Traces of the F=0.002 c/d (P=1.37 yrs) component can be found, but it is strongly distorted. The reason for this distortion is that these periodicities cannot be exactly recovered because of the sampling. Except for this deficiency, the main time scales of the illsampled data can be recovered by the methods described above. This result gives a warning that, because of the gapped nature of stellar data, activity cycles around 12 years cannot be easily recovered, if at all.
The above conclusions are valid for all the timefrequency distributions we tested. We display the results for the filtered signal ``II'' for three different transformations. The main components are recovered with all distributions. The cross terms are clearly visible in PWD and CWD (especially in CWD). Closely spaced frequencies, together with their cross terms, tend to form loops (or ``bubbles''), like the ones visible in all PWD figures between and 4000 days at low frequencies. These are welldefined signatures of structures with two different frequencies, even shortlived ones, and their appearance remains mostly invariant against changes in the preprocessing (e.g. smoothing). STFT with our test data clearly displays the same double frequency structure, however, according to our experience, PWD is superior in some cases to STFT in this property.
To give a quantitative test on the combined effect of sampling and noise in a statistical sense, we performed MonteCarlo simulation of the effect of observational noise on signal ``II''. We calculated 10^{4} realisations of the added Gaussian noise sequence with different amplitude (standard deviation). Then for each realisation we calculated the same procedure as for the original test data. The timefrequency distribution of the difference of the noisy and the noiseless smoothed data was obtained and the statistics of the amplitude of the largest feature (peak or ridge) was estimated. Figure 4 shows the falsealarm probability FAP(A), i.e., the probability that a structure with amplitude larger than A appears in the transform due to observational noise. This test was performed with different standard deviations () of the added Gaussian noise and for all the timefrequency distributions we use. In the figure from left to right the curves belong to , 0.005, 0.01, 0.02, 0.04, and 0.08 mag, representing a wide range of observational noise. The different distributions behave similarly in this test, with only slight differences. The results with STFT (black) and CWD (grey) are displayed in the figure, the curves obtained with PWD are located between the two plotted ones. It is clear that the location of the curves are scaled linearly with the amplitude of the added noise. Only very high noise, above the displayed ones, distorts this trend, as the spline interpolation becomes unstable. The amplitude of the observed features of the distributions of the original signal is displayed by a black rectangle in Fig. 4. It is clearly demonstrated that even a lot of noise is negligible compared to the real signal in the timefrequency distribution. We have to note that this result depends on the sampling and the signal content of the observations, so it should be repeated for datasets with different characteristics. The nature of noise also influences the above results. We performed our tests with a broad band (white) noise. For correlated, e.g. colour noise, with the same power, false structures appear with higher probability within the frequency range of the noise. Since the power density increases linearly with the inverse of the bandwidth of the noise (with the same total power), a similar shift is expected in the FAP(a) curves. However, as our experience suggests, the effect of a realistic noise turns out to be negligible in most applications.
The most probable timefrequency distribution due to observational noise can be calculated during the MonteCarlo simulation. For pure Gaussian noise, the expected value of the timefrequency distribution is the same for all timefrequency pairs. The averaging and spline smoothing interpolation, however, acts as a lowfrequency filter, so depends on f. Because of the uneven sampling, this filtering also depends on the time. Therefore, of the noise components can be used as a combined timefrequency transfer function of the preprocessing of the data. Figure 5 displays this transfer function for the added noise of for STFT and PWD. With PWD the most probable range of spurious components is narrower in time than that with STFT. This test confirms the above findings that the amplitude of the signal is decreased in the high frequency part of the distribution (around 1 year). For times when the sampling is very poor, in addition, an amplification is possible at low frequencies. This amplification is located on the timefrequency plane at periods that are comparable to the length of the longest gaps. However, according to our experience, this effect can be drastically reduced with a thorough interactive control of the parameters. In the MonteCarlo test, of course, it is not possible to do the same fine tuning, and the structures are visible as result in Fig. 5.
Figure 4: Falsealarm probability that a structure appears in STFT (black) and CWD (grey) due to observational noise. The amplitude of the added Gaussian noise ( from left to right) are 0.0025, 0.005, 0.01, 0.02, 0.04, and 0.08. The black rectangle represents the amplitudes of the visible real structures in Fig. 3. 

Open with DEXTER 
Figure 5: Timefrequency transfer function of the processing. Top: STFT; bottom: PWD. 

Open with DEXTER 
5 Application to solar cycles
Longterm solar records have been analysed by several authors to date. A review of solar cycle evolution is presented by Usoskin & Mursula (2003), with an impressive list of references on the subject. The longterm behaviour of the sunspot group numbers have been analysed using wavelet technique by Frick et al. (1997) who plotted the changes of the Schwabe cycle (length and strength) and studied the grand minima. However, there are still features remaining that can be revealed by the use of timefrequency methods. Here we present some of these characteristics from sunspot numbers and solar radio data.
5.1 Variations related to the Schwabe cycle
Figure 6: PWD of the yearly sunspot number time series ( top) and the same of the datasets generated from the January ( middle) and July ( bottom) monthly averages. 

Open with DEXTER 
It is well known that the sunspot number is not the best physical indicator of solar activity; however, it is the only direct measure we have for a relatively long time span. Moreover, the sunspot number has an intrinsic fluctuation because of the arbitrariness of the number of visible spots and groups. If we assume that the averages of e.g. the January data are of independent measures those of the July data (i.e., the ``noise'' for the averages of a given month are mostly independent from that of the month half a year before and after), then it is possible to construct two independent datasets: one for the January and one for the July data. The comparison of the results for the two sets gives an idea of the errors in the timefrequency distributions. Figure 6 (middle and bottom panels) exhibits the timefrequency distributions for periods longer than 2.5yr for both data subsets. The PWD of the two subsets of the data is hardly distinguishable for periods longer than 4 years.
In Fig. 6 (top) we present the PWD of the whole yearly averaged sunspot data. The STFT (Fig. 2, bottom) shows basically the same structure as the PWD, but with less visibility and resolution. Also, because the STFT does not contain any false cross terms, it validates the higher resolution results of the PWD.
Both the Schwabe and the Gleissberg cycles can be easily identified in Fig. 6. In the following we denote the instantaneous frequencies of these cycles by and , respectively. The modulation of the Schwabe cycle is well known. The period of the main variation in solar activity fluctuates between 9 and 13 years, which is clear in our results (Figs. 2 and 6). However, the previous investigations neglected to check the harmonics of frequency . Since the solar cycle is not a sinusoidal variation, one should find any power at and possibly at . Since the harmonic frequencies are synchronised to , investigating these harmonics leads to an independent estimate of the variation in the cycle length. The inspection of the history of the component confirms that around 1780, the cycle had a period of about 9 years, then it increased to 11 years after the minimal activity.
Till the middle of the 20th century, the period decreased to 10 years, then it made a swing to a longer period. It is noticeable that the amplitude of the component is very small when the Schwabe period changes. It looks as if there were some locked phases, with almost constant period, when the power at is high. Another important feature related to the nature of the component is that its amplitude is increased for the solar cycles with high amplitude more than expected from the ratio of the maximum levels. Since the even harmonics indicate the asymmetry of the form of a solar cycle (symmetric waveforms only have odd harmonics), another indication of the Waldmeier effect and related phenomena (see e.g. Cameron & Schüssler 2008).
The ``bubble'' like structures in timefrequency plots indicate that two frequencies coexist close to each other for a finite time interval. The power is clearly split into two frequency parts after 1950, and it is most visible around 1970. The frequency splitting is the same at all three places ( ), which indicates a nonlinear connection between different periodicities. The splitting is also visible in the other timefrequency distributions (STFT, CWD), but most prominent in PWD as already mentioned in Sect. 4.3. As an additional test we calculated the timefrequency distribution from the 10.7 cm radio flux of the Sun, which is available from 1947, and is a good proxy for its activity. The solar radio dataset was recorded at the Dominion Radio Astrophysical Observatory (DRAO) at Penticton, Canada (Tapping & Charrois 1994), three times a day, and are given in solar flux units or Janskys. The result is plotted in Fig. 7, and it is clear that the same structure is found in this diagram like in Fig. 6 from sunspot numbers; however, we cannot resolve the Gleissberg cycle from the radio data because of the short time base. The possible frequencies appearing in the splitting of the power are indicated by horizontal dotted lines in Figs. 6 and 7. This test strongly confirms the physical origin of the structure that we have uncovered. To unfold all the regularly spaced frequencies (, , and , k=1, 2), at least one more frequency is needed in addition to the wellknown ones ( and , i.e., the Gleissberg and the Schwabe frequencies, respectively).
Figure 7: PWD of the monthly solar radio flux data. 

Open with DEXTER 
5.2 Longterm variations in the solar cycle
Figure 8: PWD of the long term envelope of the sunpot data. 

Open with DEXTER 
The temporal evolution of the Gleissberg cycle can also be seen in the timefrequency distribution of the solar data. The Gleissberg cycle has been found to be as variable as the Schwabe cycle. It has two higher amplitude occurrences: first around 1800 (during the Dalton minimum), and then around 1950. Very interesting behaviour is the continuous decrease in the frequency (increase of period). While near 1750 the cycle length was about 50yr, it lengthened to approximately 130yr by 1950. This period variation explains why different works give different periods for the Gleissberg cycle. This systematic change in the frequency does not support the conclusion of Lawrence et al. (1995) that the power at long periods is a fingerprint of a perioddoubling cascade. The regular pattern of frequency splitting after 1950 gives the possibility of estimating the value of the Gleissberg period for the last decades, and can be precisely determined from the splitting at and . Similarly, is welldefined. From these two ranges ( and , so ), we suggest that the period of Gleissberg cycle has already reached a level above 140 years. This result agrees with the rate of frequency decrease from 1800 to 1950.
We extended the Schove (1955) series that lists activity maxima, minima, and amplitudes, based mostly on aurora records before 1610 and using optical records afterwards, by the envelope of the recent sunspot data. The dates and values of the maxima were splinesmoothed and interpolated to get an equally sampled envelope curve. The last part of the timefrequency distribution of the envelope displays the same structure as the lowfrequency part of the PWD of the recent data (Fig. 6). Figure 8 exhibits the smoothed Schove series together with the variation in the Gleissberg cycle for 500 years.
The Gleissberg cycle had a long period during the Maunder minimum. After 1700 the period jumped to a lower value (50 years) and started a slow increase. If we accept that the Gleissberg cycle gives the primary indication of the long term behaviour of sunspot numbers, we can interpret the 500yearlong history of solar activity as follows. The long period of the Gleissberg cycle is a good indicator of an extended Maunder minimum, while the subsequent, sudden increase in the Gleissberg period (after 1700) coincides with its termination. An increase in the amplitude of the Gleissberg cycle around 1800, when its period was relatively short, resulted in a shorter (Dalton) minimum. During the last century, the Gleissberg period became long, but now indicates a long lasting active stage of the Sun. Can we derive any prediction about the possible future of solar cycles from the extension of this history? The answer depends on the nature of the newly found splitting of the frequencies. If it is a sign of the restart of the Gleissberg cycle at a shorter period, then we can expect faster changes in the longterm activity level, perhaps a sudden termination of the high activity period. If the 30 year cycle (i.e., , 30yr) is just a temporary one, then a very slow decrease in activity level is expected. To rule out or confirm any of the two possibilities we have to wait for further observations or construct a physical model that is able to produce all the observed behaviour and unfold the nature and connection of all the different periodicities.
Unfortunately, at present we are very far from the construction of a proper physical model of the Sun, if it is at all possible. Bushby & Tobias (2007) discuss the possibilities of predicting the solar cycle through meanfield models, in which the cycle modulations originated from stochastic or deterministic processes. They found that good prediction even of the nearest cycles is impossible in both cases. However, timeseries analysis techniques, which find a mapping function for the data, may give better results at least for a shortterm prediction.
6 Conclusions
It is demonstrated that timefrequency distributions are useful tools for analysing nonstationary time series, because the shortterm Fourier transform, the Wigner distribution, and its extensions are more useful than wavelets. Longterm changes (cycles)of active stars can be correctly identified with the illsampled observational data after careful preprocessing. We found an extremely complicated multiscale evolution in the solar activity. All the observed features in the timefrequency history of the Sun should provide strong constraints on modelling the solar magnetism.
Whether it is possible to predict the longscale future of solar activity from the similarities of the structures that have emerged in the last decades and those in the early history of the Sun is uncertain. Its verification would be an important future project both in theory and data analysis, because of the probable connection between the longterm changes in our climate and solar activity (see Velasco & Mendoza 2007, and references therein, for more).
Acknowledgements
We are grateful to A. F. Lanza and W. Soon for useful discussions during the preparation of the manuscript. Our thanks go to the referee, Dr. F. Baudin, whose suggestions helped to improve this paper, as well as the companion paper. K.O. acknowledges support from the Hungarian Research Grants OTKA T048961 and T068626.
References
 Buchler, J. R., & Kolláth, Z. 2001, in Nonlinear Studies of Stellar Pulsation, ed. M. Takeuti, & D. Sasselov, Atrophysics and Space Science Library, 47, 185 (In the text)
 Buchler, J. R., Kolláth, Z., Serre, T., & Mattei, J. 1995, ApJ, 462, 489 [NASA ADS] [CrossRef] (In the text)
 Buchler, J. R., Kolláth, Z., & Cadmus, R. R. 2004, ApJ, 613, 532 [NASA ADS] [CrossRef] (In the text)
 Bushby, P. J., & Tobias, S. M. 2007, ApJ, 661, 1289 [NASA ADS] [CrossRef] (In the text)
 Cameron, R., & Schüssler, M. 2008, ApJ, 685, 1291 [NASA ADS] [CrossRef] (In the text)
 Choi, H., & Williams, W. J. 1989, IEEE. Trans. Acoustics, Speech, Signal Processing, 37, 862 (In the text)
 Cohen, L. 1995, TimeFrequency Analysis (Englewood Cliffs: PrenticeHall) (In the text)
 ForgácsDajka, E., & Borkovits, T. 2006, MNRAS, 374, 282 [NASA ADS] [CrossRef] (In the text)
 Frick, P., Galyagin, D., Hoyt, D. V., et al. 1997, A&A, 328, 670 [NASA ADS] (In the text)
 Gleissberg, W. 1939, The Observatory, 62, 158 [NASA ADS] (In the text)
 Hoyt, D. V., Schatten, K. H., & NesmeRibes, E. 1994, Geophys. Res. Lett., 21, 2067 [NASA ADS] [CrossRef] (In the text)
 Kolláth, Z., & Buchler, J. R. 2001, in Nonlinear Studies of Stellar Pulsation, ed. M. Takeuti, & D. Sasselov, Astrophysics and Space Science Library, 47, 29 (In the text)
 Kolláth, Z., & Szeidl, B. 1993, A&A, 277, 62 [NASA ADS] (In the text)
 Lawrence, J. K., Cadavid, A. C., & Ruzmaikin, A. A. 1995, ApJ, 455, 366 [NASA ADS] [CrossRef] (In the text)
 Oláh, K., Strassmeier, K. G., Granzer, T., Soon, W., & Baliunas, S. L. 2007, AN, 328, 1072 (In the text)
 Oláh, K. Kolláth, Z., Granzer, T., et al. 2009, A&A, 501, 703 [NASA ADS] [CrossRef] [EDP Sciences] (Paper II)
 Reinsch, C. H. 1967, Numer. Math., 10, 177 [CrossRef] (In the text)
 Schove, D. J. 1955, J. Geophys Research, 60, 127 [NASA ADS] [CrossRef] (In the text)
 Schwabe, H. 1843, reprinted in Early Solar Physics, ed. A. J. Meadows (London: Pergamon), 95 (In the text)
 Soon, W., Frick, P., & Baliunas, S. 1999, ApJ, 510, L35 [CrossRef] (In the text)
 Tapping, K. F., & Charrois, D. P. 1994, Sol. Phys., 150, 305 [NASA ADS] [CrossRef] (In the text)
 Usoskin, I. G., & Mursula, K. 2003, Sol. Phys., 218, 319 [NASA ADS] [CrossRef] (In the text)
 Velasco, V. M., & Mendoza, B. 2007, Adv. Space Res., 42, 866 [NASA ADS] [CrossRef] (In the text)
 Wigner, J. 1932, Phys. Rev., 40, 749 [NASA ADS] [CrossRef] (In the text)
All Figures
Figure 1: Shortterm Fourier transform ( bottom) of the sunspot number time series ( top). Fouriertransforms of the subsets of the dataset are also plotted. The colour scale of relative variations in the distributions is presented in Fig. 7. 

Open with DEXTER  
In the text 
Figure 2: Wavelet ( top and middle) and shortterm Fourier transform (STFT) distribution ( bottom) of the yearly sunspot numbers. 

Open with DEXTER  
In the text 
Figure 3: Timefrequency distribution of test data. Panel A): PWD of signal ``I''; panel B): PWD of signal ``I'' with white noise; panel C): PWD of signal ``II'' with averaging and spline interpolation; panel D): PWD of signal ``II'' after Fourier filtering of the rotational period in the independent observational seasons and then averaging and spline smoothing; panel E): STFT of the same signal as in panel D); panel F): CWD of the same signal as in panel D). 

Open with DEXTER  
In the text 
Figure 4: Falsealarm probability that a structure appears in STFT (black) and CWD (grey) due to observational noise. The amplitude of the added Gaussian noise ( from left to right) are 0.0025, 0.005, 0.01, 0.02, 0.04, and 0.08. The black rectangle represents the amplitudes of the visible real structures in Fig. 3. 

Open with DEXTER  
In the text 
Figure 5: Timefrequency transfer function of the processing. Top: STFT; bottom: PWD. 

Open with DEXTER  
In the text 
Figure 6: PWD of the yearly sunspot number time series ( top) and the same of the datasets generated from the January ( middle) and July ( bottom) monthly averages. 

Open with DEXTER  
In the text 
Figure 7: PWD of the monthly solar radio flux data. 

Open with DEXTER  
In the text 
Figure 8: PWD of the long term envelope of the sunpot data. 

Open with DEXTER  
In the text 
Copyright ESO 2009