An automatic method to determine the frequency scale of the ionospheric Alfv´en resonator using data from Hylaty station, Poland

. ULF/ELF magnetic ﬁeld data recorded at the “Hylaty” station in Poland (49 ◦ 19’ N, 22 ◦ 56’ E; L ’ 2) are analysed to ﬁnd the characteristics of spectral resonance structures (SRS) in the frequency range 1–5 Hz, related to the ionospheric Alfv´en resonator (IAR). An automatic procedure is employed to SRS events observed at “Hylaty” during the nighttime in 2001–2003, to calculate the parameter which determines the separation between the harmonics of the resonator, termed the frequency scale. Diurnal and sea-sonal variations of the frequency scale within the range of 0.4–0.8 Hz have been found. The usefulness and disadvantages of this particular method of SRS analysis, and of other methods, are discussed.


Introduction
Since 1985, when spectral resonance structures (SRS) were discovered in the frequency range from 1-10 Hz, in the spectra of the magnetic field observed at the mid-latitude observatory "Novaya Zhizn" (Belyaev et al., 1987), several long-term observations have been made in different parts of the globe and statistical analyses of the phenomenon performed. The structures were observed first at high latitudes in Kilpisjärvi (Belyaev et al., 1999) and then in Sodankylä (Yahnin et al., 2003;Hebden et al., 2005), Finland. Recent reports show that the SRS is observed at a very high-latitude site, Barentsburg in Svalbard, in the polar cap region (Semenova et al., 2005). At low latitudes the effect was observed in Crete (Bösinger et al., 2002) and fine spectral resonance structures (FSRS) have been discovered Correspondence to: A. Odzimek (ao64@ion.le.ac.uk) there (Bösinger et al., 2004). The observations were continued at the Nizhny Novgorod observational station "Novaya Zhizn" (Belyaev et al., 1997) and more events were observed at "Borok" and "Mondy" (Pokhotelov et al., 2003), and Karimshino (Molchanov et al., 2004), Russia. Single observations of the SRS were also made at low and middle latitudes, and noted by Hickey et al. (1996), Kułak et al. (1999), and Odzimek (2004a). The structures are believed to be the effect of the ionospheric Alfvén resonator (IAR); both observations and modelling confirm this (see, for example, Belyaev et al., 1990;Demekhov et al., 2000;Yahnin et al., 2003).
The spectral resonance structures observed at different latitudes have common characteristics, originating from the same phenomenon of the ionospheric Alfvén resonance, but there are also differences arising from variations on a regional scale. The main parameters revealing this diversity are the occurrence rates and the frequency scales f of the structures. The frequency scale is defined as the separation between consecutive resonance frequencies (Polyakov and Rapoport, 1981). This depends on the parameters of the height profile of the Alfvén speed in the ionospheric Fregion, where the resonator is located: where v AF is the minimal Alfvén speed near the maximum ion density in the F-region and l AF is the geometric size of the resonator, determined by the geometric parameters of the Alfvén speed profile in the F-region. The SRS observed on the Earth's surface at a particular observing station consists of IAR lines of the resonator above the point of observation, and, consequently, the SRS is characterised by the same frequency scale (Belyaev et al., 1990). In this paper some results of an analysis of the spectral resonance structures derived from ground-based measurements of ULF/ELF magnetic field made in Poland are presented.
An automatic numerical algorithm has been developed to calculate the SRS frequency scale, and diurnal variations of this parameter have been obtained using this method; both are discussed. In a separate section, the relative advantages and disadvantages of the algorithm are considered, as well as those of other methods which has been recently published.

Acquisition and analysis of data
Observations in the Polish part of the Bieszczady Mountains, at a site named "Hylaty" (49.19 • N, 22.56 • E; L 2), were performed to study Schumann resonances (∼ 8,14,20... Hz) in the ELF magnetic field. Observations of the magnetic N-S horizontal component have been made there sporadically since May 1996, andup to November 2003, the number of 24-h measurements amounted to about eighty. The usual schedule of observations was twelve 5-min measurements each hour, except for the earlier period to September 2001 (when the signal was measured four times each hour, for 5 min). The antenna at "Hylaty" was a solenoidal magnetic antenna of length 1 m, consisting of 240 000 turns of copper wire distributed over five coils (to reduce the capacitance of the antenna and provide a wide transmission band). Passing the signal from the antenna through a lownoise pre-amplifier limited the low-pass frequency of the signal to 0.03 Hz. Moreover, the signal was integrated before being converted to digital, so that the output was proportional to the magnetic field itself. The antialiasing filter cut off at 55 Hz and the signal was recorded with a sampling frequency of 178 Hz. Since the frequency range where the SRS is usually observed (1-5 Hz) lies within the observed range, it was possible to explore the occurrences of the SRS.
There are not many fast and effective methods to detect SRS events and to estimate parameters of the structures, such as the frequency scale. Different methods of analysis were tested with the data recorded at "Hylaty". In the end, a special procedure was designed for these data and applied to find and analyse the spectral resonance structures. An important aim was to make the analysis automatic as much as possible. The algorithm consisted of: I. Calculating averaged spectra. A single 5-min time series was analysed using the standard spectral analysis technique known as the Fast Fourier Transform method. The series was divided into smaller series, each with a length of 2048 points (equivalent to 11.5 s), and Fourier-transformed. A power spectrum was calculated as an incoherent average of the transforms (Lyons, 2000). A 11-point Hamming window with 50% overlap of the data was used in the frequency domain to obtain a smoother spectrum for a "pre-detection" algorithm, described in the next step of the procedure.
II. Searching for SRS cases. For every calculated power spectrum the frequencies of its local maxima were found by the condition that the amplitude of the maximum was higher than the amplitudes at the two adjacent frequencies (at 0.0056-Hz intervals). This condition was checked for each power spectrum up to 7 Hz. The frequencies of the maxima were plotted on a two-dimensional frequency versus time plane. In the plot of the positions of the maxima in the frequency-time domain, when there is no SRS, the maxima should occur at random, spread out over the entire area. When the SRS is present, the points corresponding to the maxima tend to be grouped along several stripes, the separations of which change with time. The stripes are related to IAR frequencies and the separation between the stripes reflects the frequency scale of the structure -see left windows in Figs. 1a-f, for six different nights during 2002 and 2003, when the SRS has been revealed. Such an approach, which at first sight seems to be primitive, usually gives a clearer picture of an SRS than a standard spectrogram, because the dependence on amplitude is suppressed. Although this method defines the SRS events, it still requires visual examination and is not fully automatic. This procedure is referred to as "pre-detection".
III. Calculation of the SRS frequency scale. The SRS cases were analysed to estimate the frequency scale f by employing an automatic procedure (the algorithm will be referred later as the "fitting method"). The frequencies of the local maxima of a spectrum serve as the input to the procedure, chosen automatically in the procedure of pre-detection (step II), in the range of the spectrum in which the SRS structure appears. The range can be specified earlier on the basis of the pre-detection plot. Then the calculation of the SRS frequency scale was made according to the algorithm of Odzimek (2004a), with some modifications that seemed to improve the method. The steps of the procedure can be summarised as follows: -The experimental frequencies for the n-th resonance are assumed to obey the theoretical result from Trakhtengerts and Feldstein (1981), which is valid for the condition in the nighttime ionosphere: -The parameter f in the argument of the Bessel function J 1 on the left side of the equation may be estimated by a least-squares fitting to the predicted formula, i.e. the minimum of the function should be found with respect to the parameter f in the range of possible values of f : 0.2-1.2 Hz was allowed in our calculations. Heref m is one of the experimental frequencies assumed to be a resonance one (they are found in step II), and M is the number of these frequencies in the considered frequency range.  The results of the minimisation are usually not unique in the chosen range, and an additional criterion is added to choose the optimal value of f by reviewing the correspondence between the theoretical resonance frequencies and the experimental ones. First, for each f found by the minimisation, the set of theoretical resonance frequencies {f n } is calculated; they are related to the zeroes of J 1 by f/π (Eq. 2). The identification of an experimental frequencyf m with a theoretical frequency f n is assumed to occur when |f m −f n |<δ, where δ depends on the current value of f and is set to 30% of f .
-While the sets of theoretical and experimental frequencies are compared, a so-called "fitting parameter" χ is derived, which tells how these two sets vary from each other. It is calculated by adding |f m −f n | if there is a theoretical frequency f n in {f n } n=1...N identified with one of the frequenciesf m from {f m } m=1...M , or it is increased by δ every time there is no corresponding experimental frequencyf m to a theoretical frequency f n , or if there is a frequencyf m in the set without a corresponding frequency f n . The final result is that value of f , obtained by the minimisation of the function Eq. (3), for which the fitting parameter χ is minimal.
The calculation of f with the fitting method, for spectra of consecutive time series, is now fully automatic. The variations of the calculated f are expected to follow the variation of f during an SRS event; in the case of no SRS, rather random values appear all the time. In the latter case, the method was tested with a white noise signal that simulates a time derivative of the magnetic signal. The data were created by a random number generator and then the analysis described in steps I and II was done. Consecutive values of " f " calculated for these data in step III did not seem to show a regular pattern of variations; however, most of them were concentrated in the range of 0.3-0.5 Hz, which is possibly connected with local maxima in a power spectrum calculated with the parameters of step I. It means that the algorithm applied to the spectra calculated with those parameters may not be efficient with structures of smaller f .
For 39 days of measurements at "Hylaty", from October 2001 to November 2003, employing the pre-detection procedure, we distinguished ten events when the SRS occurred and lasted for several hours. This seems to be rather a small number of events; however, some observations were made on days of enhanced solar or geomagnetic activity, for Schumann resonance studies, and these conditions might be not favourable for SRS occurrence (Yahnin et al., 2003). Figure 1 presents some examples of the diurnal variations of f obtained with this fitting method. Occasionally the results form a regular pattern of f variations, as in Fig. 1f, and sometimes only in some parts of the SRS events, as in Figs. 1a, b, e. Random values accompany these variationsthis effect depends on the efficiency of the fitting method, to a large extent, which will be discussed further in Sect. 3.
In Fig. 2 variations of the respective hourly means of the values obtained with the fitting method are shown. Higher values of the standard deviation of the mean correspond to times when a structure was changing quickly or when it was not very clear, so that the fitting method was not effective -both effects usually arise at the beginning and at the end of the structure. There is a larger dispersion of subsequent values of f and, moreover, the estimated values of f are smaller -and so is the mean (compare, for example, Figs. 1f and 2f at about 00:00 UT). Both the single value evaluations and the mean values indicate that the frequency scale increases during the night. A larger spread of f values depends on the dynamics of the SRS, and also on the efficiency of the fitting method.
Beside the effects depending on the efficiency of the method, the results of the f calculations are in agreement with the visual estimates of the pre-detection figures. They give estimates for variations of the local f parameter of the SRS, observed at "Hylaty" in 2001-2003, of approximately 0.4-0.8 Hz, depending on the hour, and the season (higher values in September, October than in August). This is in agreement with the results of other observations at middle latitudes (Belyaev et al., 1990(Belyaev et al., , 1997, especially the observations in Karimshino, at the site of similar L-shell: L 2 (Molchanov et al., 2004).

Comparison of the fitting method with other methods
The purpose of this section is to summarise the advantages and disadvantages of the method developed compared with the manual method most often employed to study spectral resonance structures and with the method developed by Molchanov et al. (2004). Most SRS observations have been analysed manually thus far on the basis of a visual study of the spectrograms. The frequency scale f of an SRS had been estimated by measuring the intervals between the maxima or minima of the structure in a spectrogram (A. Yahnin and N. Semenova, personal communication). Molchanov et al. (2004) presented a numerical algorithm for calculating the frequency scale and the quality parameter for an SRS, but this method can be applied only for longer time series, to obtain estimates with sufficient accuracy. Hebden et al. (2005) used a manual cursor-clicking technique, and also an automatic method very similar to this introduced by Molchanov et al. (2004). The continuous, automated search for the structures and determination of their parameters (at least the frequency scale) is essential with the increasing amount of experimental data available. The development of fast and more fully automatic algorithms, like the fitting method discussed here, is needed. The method includes a particular model of the IAR to perform the calculation of the SRS frequency scale. In the model used f is a parameter of the Bessel function J 1 that determines the positions of the eigenfrequencies (Trakhtengerts and Feldstein, 1981;Odzimek, 2004a). The model assumes that the reflection coefficient of Alfvén waves from the bottom ionosphere is close to 1, which is thought to be correct for the nighttime mid-latitude ionosphere. In this model and method the positions of the frequencies are important. The positions predicted by the model are slightly different from the model, which considers both the first eigenfrequency and the spacing between neighbouring eigenfrequencies to be equal (Polyakov and Rapoport, 1981). However, if f is estimated as a mean interval between subsequent resonance frequencies and the second or the third frequency is taken into account (which, in reality, is more frequent because sometimes the first could not be observed or is hard to recognise), or if f is estimated by the method of Molchanov et al. (2004), a similar value of f should be obtained in these approaches.
The main positive features of the fitting method are: -It is automatic, not too time-consuming and not subjective in interpreting the structures, as are the manual methods. It can be used to analyse large amounts of data.
-The reaction of the method to data with developed and clear SRS structures, when it calculates the values of f for subsequent time series, is peculiar in comparison to pure noise. This gives the prospect that it could be used both to estimate f and to search for SRS at the same time.
-It does not require enumeration of the resonance frequencies, and the sequence of the resonance frequencies is obtained simultaneously with the result of the fit.
-With other models of the IAR, other formulas determining the resonance frequencies might be used instead of J 1 and the sense of calculation could be maintained. But if it has to include more parameters the problem that can happen is that the minimisation no longer remains a one-dimensional problem.
The drawbacks of the method are: -The method is not so good when it is applied to data with fine or unclear structures, and where there is more noise in the signal.
-The important reason for false estimates of f could be a wrong or imperfect choice of input frequencies, assumed to be the resonance frequencies and constituting an SRS. The method of manual estimate is probably able to do this step better, because it relies on the human eye, which is able to resolve the structure quickly and accurately.
-The accuracy of a single evaluation of f is not determined. Another useful value that may be related to an estimate of the frequency scale of an SRS would be an index or a parameter that could indicate a level of quality of the structure, an example of one is already provided by Molchanov et al. (2004). There is no quality parameter introduced for the fitting method at present. Also, a quality parameter of the whole SRS event might be suggested.
-The model used is an ideal one and the observed resonance frequencies should have smaller values than an ideal model predicts, as usually happens with damped resonators. When we refer f to ionospheric parameters of the ion density profile from the theoretical model, the f we observe is probably a deformed value of f given by Eq.
(1). This issue should be remembered when referring the observed values of f to ionospheric parameters. Now a comparison of the fitting method with the method developed recently by Molchanov et al. (2004) will be discussed; the latter method will be called the "double transform method". To explain the details of the analysis with the double transform method, we mention the main steps in this method. It is based on a model with equally spaced resonance frequencies. A second Fourier transform is performed on the part of a spectrum where an SRS appears, after subtracting the trend of the background in the spectrum. The main period of the frequency repeat is obtained by this transform and the frequency scale is the period (in Hz) corresponding to the maximum of the amplitudes of the secondary spectrum. Though it is improper for the data collected at "Hylaty" because of the limited accuracy of the results, some examples of how the two methods work with these data are shown in Fig. 3.
The left-hand column of Fig. 3 corresponds to the results of the fitting method and the right-hand second column shows the analysis and results of the double transform method, for the same spectra and structures. The comparison of the two methods, including the results from the cases presented in Fig. 3, can be summarised as follows: -Both methods are automatic and can be applied to large amounts of data in a short amount of time.
-The important advantage of the double transform method is that it contains known standard methods, such as FFT, and the "quality parameter" can be calculated for a single f estimate (but only discrete results are available, and long time series are needed to obtain the result with good accuracy).
-Both methods give very similar results when the structure is clear -see Figs. 3a and b or Figs. 3e and f.
-Both methods fail with unclear or fine structures, at least for the data from "Hylaty" -compare the different results for Figs. 3c and d or 3g and h.
-In both methods the result can depend on the input parameters of the method, for example, the range of a spectrum chosen for further analysis (not showed in Fig. 3). This more often happens with unclear structures, when some resonance peaks are not sufficiently distinguished in the background.
The best way to validate such methods, especially in questionable cases, would be a comparison with real ionospheric parameters obtained from an ionospheric diagnosis above the station at the same time. This is not easy since the IAR parameters depend on upper ionosphere characteristics, such as the topside ion density profile, obtainable from sounding by satellites or with incoherent scatter radars on the ground. As for high-latitude studies in Europe, the well developed EIS-CAT system is available, for mid-latitudes the only possibility to use such a facility might be the Kharkov incoherent scatter radar in the Ukraine. This would require an essential undertaking before we could use any of these methods to monitor the SRS parameters and consequently, have a possibility to use them as an additional method of diagnosis of the ionosphere, as was suggested by Belyaev et al. (1990). The additional methods would broaden our knowledge of this part of space, and improve topside ionospheric models. Early attempts to compare the observed values of some ionospheric parameters inferred from the observed f (or directly from f ) and respective modelled values indicate quantitative differences in these two approaches (Yahnin et al., 2003;Odzimek, 2004b). Moreover, the important advantage of diagnosis by observations of the SRS is that this method is cheap and does not require powerful equipment.

Conclusion and future prospects
At the observational site "Hylaty" in Poland, ten strong SRS events have been discovered among forty 24-h observation runs from 2001 to 2003. A new automatic procedure, the fitting method, has been described and used to obtain the diurnal variations of the SRS frequency scale of the events. Spec-tral resonance structures need to be analysed with faster and automatic methods, such as the fitting method, to study the phenomenon on a large scale. Connected with the possibility to use it to diagnose the upper ionosphere, first the efficiency of these methods should be improved and the results should be checked with the results of standard ionospheric diagno-sis. The main objectives to improve the fitting method can be summarised as follows: -Improving the efficiency of the method when analysing noisy signals or those SRS events with a smaller frequency scale.
-Deriving a parameter specifying the quality of the SRS when using the method.
-Deriving the accuracy of a single evaluation of the frequency scale, possibly connected with a "quality" parameter.
-Using the diagnosis of the ionosphere above the point of observation simultaneously with the automatic analysis of the observed SRS events.
-Checking how the method works with other data, recorded at other observational sites at different latitudes.
-Comparing the results of the fitting method with other methods available with larger amounts of observational data.
A final challenge would be using an effective analysis of observations of the IAR spectral resonance structures as an additional method of investigating the topside ionosphere.