Multiple and changing cycles of active stars
Key Words.:
Sun: activity – methods: data analysisAbstract
Context:
Aims: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, like variations in the global climate. Therefore understanding the nature of these variations is extremely important.
Methods:The observed variations related to magnetic activity cannot be treated as stationary periodic variations, therefore methods like Fourier transform or different versions of periodograms give only partial information on the nature of the light variability. We demonstrate that timefrequency distributions provide useful tools for analysing the observations of active stars.
Results:With test data we demonstrate that the observational noise has practically no effect on the determination in the 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 which are discussed. Applying our technique to the sunspot data we find a complicated, multiscale evolution in the solar activity.
Conclusions:
1 Introduction
It is well known that solar cycles cannot be represented by single periodic variations. Besides the 11 year quasiperiodicity (Schwabe schwabe ()) of solar activity lots of different rhythms were discovered or at least suggested by different studies (see e.g. Frick et al. frick (), ForgácsDajka & Borkovits emese ()). One remarkable and wellstudied periodicity of the solar variation was found by Gleissberg (gleissberg ()) 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. (frick ()) derived the length of the Gleissberg cycle as 100 years from wavelet analysis of a 384 years long dataset (Hoyt et al. hoyt ()). The nonstationarity in the longterm variation of active stars has already been found in active stars as well (Oláh et al. cycles3 ()).
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 monitored easily 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 for 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 kernel method (Cohen cohen ()) 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, wigner ()) and its extensions by Cohen (cohen ()). These methods have also been proven to be useful in the analyses of variable star data of chaotic origin or nonstationary behaviour (see Buchler et al. buchler1 (), Kolláth & Buchler kollath1 ()). In this paper we present a method which 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 reanalyzed 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. paper2 ()) 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 around a given epoch . If this windowed part of the signal is then Fourier transformed, the frequency spectrum of the data around the time 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 is dense enough, then a twodimensional map of the energy distribution on the time frequency plane can be plotted.
We demonstrate this method on the yearly sunspot numbers – a dataset which is well known and analysed by several authors. Fig. 1 displays the STFT of the solar data. On the timefrequency distribution map individual windowed Fourier spectra are also plotted with thick lines. 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.
While in this demonstration we used a time centred Gaussian window to plot the curves of Fourier transforms, in the implementation of STFT we use windowing in the Fourier space:
(1) 
where denotes the Fourier transform of the signal. Note, that the inverse Fourier transform is calculated only 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 and Kolláth buchler2 ()). Depending on the functional form of , Eq. 1 reduces to STFT or wavelet. With , where is a constant frequency, the equivalent temporal window function of STFT is given by
(2) 
Usually it is a good choice to set 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 time centered windowing function is . Then, if depends on frequency, the width of the window also varies with frequency. Especially if then is equal to wavelet, with the Morlet kernel, since the width of the temporal windowing function decreases with .
Note that in the standard Morlet wavelet is fixed to . The introduction of a free parameter () which controls the balance between frequency and temporal resolution has several advantages as already described by Kolláth & Szeidl (kollath2 ()) and Soon, Frick & Baliunas (soon ()).
As can be seen from the above definitions, 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 () is defined in such a way that for a given the wavelet and STFT have the same frequency/time resolution at the frequency . The top and middle panels of Fig. 2 show wavelets where the lower frequency part (top) and the high frequency part (middle) was matched to the STFT, respectively. Resolution in time increases with frequency, since the same number of cycles is covered in a shorter timebase, at higher frequencies. The disadvantage of this representation is that harmonic structures, which display the same characteristics in STFT, are distorted. Fig. 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.
Note that 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, in speech analysis it is well known that by the application of a preemphasis filter the obscured parts of the frequency distribution can be recovered. 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.07 – 0.14 ) the signal in different frequency ranges were amplified by the following factors (): for , and above ( is given in ). In the transition regions a smooth change in is used. 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 years periods.
The optimal balance of temporal and frequency resolution depends on the signal itself. A question arises whether it is possible to construct a distribution, where the resolution in both coordinates is optimal, but some other disadvantage is allowed as its cost. The answer is a more general group of methods for signal processing introduced by Cohen (see e.g. Cohen cohen ()). The generalised timefrequency distribution is given by
(3) 
where is the analysed time series and is the kernel of the distribution that determines the specific properties of the distribution. The simplest timefrequency distribution is the WignerVille transformation with . This distribution works well only for special datasets, otherwise the resulting map is heavily contaminated by cross terms of the different components. A biGaussian kernel provides a distribution which, 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 and Williams cwd ()), 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 which 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 deficiency, like observational noise and imperfect sampling. Then timefrequency methods cannot be applied directly to the observational data, the analyses should be preceded by a sequence of preprocessing steps.
In the case of active stars the situation is even more complicated: 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. On the contrary, long term 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 due to 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 in addition, they contain seasonal gaps. Although there exist methods 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 reinsch ()). This operation provides a continuous function based on the discrete time series ( ) by minimising the integral of , with the following constraint:
(4) 
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 as the standard deviation of the observational noise. In practice, we usually check manually the resulting spline smoothed curves to make a fine tuning of the averaging and spline smoothing. As a final step is resampled to generate an equally spaced time series.
4 Testing the methods
4.1 Artificial test data sets
To test the methods for a general signal representing the cycles of active stars, we created artificial datasets. Signal ’I’ consists of a periodic signal with a constant frequency of , a component with linearly increasing frequency (chirp) with a modulation between and and a wave packet at with a Gaussian amplitude modulation:). The amplitudes of the constant 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 in case of periodic or quasiperiodic signals the Fourier transform has a very good noise rejection property. The same is true for timefrequency distributions. If the observational noise and the signal are not correlated, then the method provides good result even in the case of large noise. A nice example of this nature of timefrequency distributions is the analysis of the semiregular star V CVn (Buchler, Kolláth & Cadmus buchler3 ()), 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, Kolláth & Cadmus (buchler3 ()) 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.
We demonstrate this noise rejection property of timefrequency distributions on our first test timeseries. 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 on the top panel of the figure. While the noise has practically no effect on the lower frequencies, the 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 tested 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 distorts 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, 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 to, while on panels D, E and F it is filtered out from the data. Averaging and smoothing spline interpolation introduce a lowpass filtering of the signal, which can be compensated by the preemphasis filtering. In these cases we selected for . With these settings and processing, all the components above about 2 years of characteristic time are recovered well, especially in the case when the rotation period is filtered out. Traces of the c/d (yrs) component can be found, but it is strongly distorted. The reason for this distortion is, that due to the sampling, these periodicities cannot be exactly recovered. Except 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, if at all, recovered.
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 2000 and 4000 days at low frequencies. These are well defined signatures of structures with two different frequencies, even short living 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, in some cases PWD is superior 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 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. Fig. 4 shows the False Alarm Probability , i.e., the probability that a structure with amplitude larger than 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. On the figure from left to right the curves belong to = 0.0025, 0.005, 0.01, 0.02, 0.04 and 0.08 mags, representing a wide range of observational noise. The different distributions behave similarly in this test, slight differences have been found only. The results with STFT (black) and CWD (grey) are displayed on the figure, the curves obtained with PWD are located between the two plotted ones. It is well seen, 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 high amount 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 have performed our tests with a broad band (white) noise. For correlated, e.g. color 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), similar shift is expected in the curves. However, as our experience suggests, in most applications the effect of a realistic noise turns out to be negligible.
During the MonteCarlo simulation the most probable timefrequency distribution due to observational noise can be calculated. 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 low frequency filter, so depends on . 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. Fig. 5 displays this transfer function for the added noise of for STFT and PWD. With PWD the most probably range of spurious components are narrower in time than that with STFT. This test confirms the above findings, that in the high frequency part of the distribution (around 1 year) the amplitude of the signal is decreased. For times where the sampling is very poor, in addition, an amplification is possible at low frequencies. This amplification is located on the timefrequency plane at periods which are comparable to the length of the longest gaps. However, according to our experience, with a thorough interactive control of the parameters, this effect can be drastically reduced. In the MonteCarlo test of course it is not possible to do the same fine tuning, and as results the structures are visible in Fig. 5.
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 (usoskin ()), with an impressive list of references on the subject. The longterm behaviour of the sunspot group numbers were analysed using wavelet technique by Frick et al. (frick ()) who plotted the changes of the Schwabe cycle (length and strength) and studied the grand minima. However, there are still remaining features which 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
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 timespan. 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 independent measures from that 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 about the errors in the timefrequency distributions. Fig. 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 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 followings 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 of solar activity fluctuates between 9 and 13 years which is well seen in our results (Fig. 2 and Fig. 6). However, the previous investigations neglected to check the harmonics of the frequency . Since the solar cycle is not a sinusoidal variation one should find some power at and possibly at . Since the harmonic frequencies are synchronised to , by the investigation of these harmonics one gets an independent estimate of the variation of 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 after the minimal activity it increased to 11 years.
Till the middle of the XXth 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 have only odd harmonics) – this is another indication of the Waldmeier effect and related phenomena (see e.g. Cameron and Schüssler cameron ()).
The ’bubble’ like structures in timefrequency plots indicate that for a finite time interval two frequencies coexist close to each other. The power is clearly split into two frequency parts after 1950, and it is best 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 it is the most prominent in PWD as already mentioned in Section 4.3. As an additional test we have calculated the timefrequency distribution from the 10.7cm 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 DRAO ()), three times a day, and are given in solar flux units or Janskys. The result is plotted in Fig. 7, and it is well seen that the same structure is found in this diagram, than that of Fig. 6 from sunspot numbers; however, we can not 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 Fig. 6 and Fig. 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 well known ones ( and , i.e., the Gleissberg and the Schwabe frequencies, respectively).
5.2 Long term variations in the solar cycle
The temporal evolution of the Gleissberg cycle can also be seen on the timefrequency distribution of the solar data. The Gleissberg cycle has been found to be variable as the Schwabe cycle. It has two higher amplitude occurrences: first around 1800 (during the Dalton minimum), and then around 1950. A very interesting behaviour is the continuous decrease of 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 of the Gleissberg cycle. This systematic change in the frequency does not support the conclusion of Lawrence, Cadavid and Ruzmaikin (lawrence ()) that the power at long periods is a fingerprint of a period doubling cascade. The regular pattern of frequency splitting after 1950 gives the possibility to estimate the value of the Gleissberg period for the last decades. can be precisely determined from the splitting at and . Similarly is well defined. From these two ranges (= 0.0250.03 and = 0.0310.032, so ) we suggest that the period of Gleissberg cycle has already reached a level above 140 years. This result agrees with the the rate of frequency decrease from 1800 to 1950.
We extended the Schove (schove ()) series which 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 spline smoothed 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 low frequency part of the PWD of the recent data (Fig. 6). Fig. 8 exhibits the smoothed Schove series together with the variation of 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 of the Gleissberg period (after 1700) coincides with its termination. An increase of 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 indicating 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 shorter period, then we can expect faster changes in the long term activity level, perhaps a sudden termination of the high activity period. If the 30 year cycle (i.e.,, yr) is just a temporary one, then a very slow decrease of 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 which 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 and Tobias (bushby ()) discussed 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: 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 long scale 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 long term changes in our climate and solar activity (see Velasco & Mendoza (velasco_mendoza ()) 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 are due to the referee, Dr. F. Baudin, whose suggestions helped to improve this paper as well as the companion paper. K.O. acknowledges supports from the Hungarian Research Grants OTKA T048961 and T068626.Footnotes
 offprints: Z. Kolláth
References
 Buchler, J.R., Kolláth, Z., Serre, T. & Mattei, J., 1995, ApJ, 462, 489
 Buchler, J.R. & Kolláth, Z., 2001, in Nonlinear Studies of Stellar Pulsation, M. Takeuti, D. Sasselov, Eds., Atrophysics and Space Science Library, 47, 185
 Buchler, J.R. Kolláth, Z. & Cadmus, R.R., 2004, ApJ, 613, 532
 Bushby, P.J., Tobias, S.M. 2007, ApJ 661, 1289
 Cameron, R., & Schüssler, M. 2008, ApJ 685, 1291
 Choi, H., Williams, W.J. 1989. IEEE. Trans. Acoustics, Speech, Signal Processing, 37, 862.
 Cohen, L. 1995, TimeFrequency Analysis, (PrenticeHall, Englewood Cliffs).
 ForgácsDajka, E., Borkovits, T. 2006, MNRAS 374, 282
 Frick, P., Galyagin, D., Hoyt, D.V. et al. 1997, A&A 328, 670
 Gleissberg, W. 1939, The Observatory 62, 158
 Hoyt, D.V., Schatten, K.H., NesmeRibes, E. 1994, Geophys. Res. let. 21, 2067
 Kolláth, Z., & Buchler, J.R. 2001, in Nonlinear Studies of Stellar Pulsation, M. Takeuti, D. Sasselov, Eds., Astrophysics and Space Science Library, 47, 29
 Kolláth, Z. & Szeidl, B. 1993, A&A, 277, 62
 Lawrence, J. K., Cadavid, A. C. & Ruzmaikin, A. A., 1995 ApJ, 455, 366
 Oláh, K., Strassmeier, K.G., Granzer, T., Soon, W., & Baliunas, S.L. 2007, AN 328, 1072
 Oláh, K. Kolláth, Z., Granzer, T. et al. 2008, A&A, submitted
 Reinsch, C. H., 1967, Numer. Math., 10, 177
 Schove, D. J., 1955, J. Geophys Research, 60, 127
 Schwabe, H., 1843, reprinted in Early Solar Physics, A.J. Meadows, Ed., (Pergamon, London, 1970), pp.9597.
 Soon, W., Frick, P. & Baliunas, S., 1999, Apj 510, L35
 Tapping, K.F., Charrois, D.P. 1994, Sol. Phys. 150, 305
 Usoskin, I.G. & Mursula, K. 2003, Sol. Phys. 218, 319
 Velasco, V.M., Mendoza, B. 2007, Adv. Space Research 42, 866
 Wigner, J. 1932, Phys. Rev. 40, 749