- Email: [email protected]
Analysis of pulsations
Planet. Printed
SpaceSci., in Great
Vol. 30, No. Britain.
12, pp. 1249-1258,
0032-0633/82/12124~1ou)3.00/0 Pergamon Press Ltd.
1982
ANALYSIS
OF PULSATIONS F. GLANGEAUD
CEPHAG, BP 46, 38402 St Martin d’Heures, France LA 346 au CNRS (Received in fina/ form 7 June 1982) Abstract-There is no universal method of data analysis. The scientist must first know what he intends to measure; then, it is possible to select the most suitable methods and give practical applications. We limit our investigation here to three domains. In the first one, estimation of the parameters, the number of degrees of freedom is estimated, taking into account the a priori knowledge in magnetic pulsation models. The random and deterministic parts of the signal are evaluated and the meaning of the confidence interval of the measurements is discussed. The use and limitations of maximum entropy and autoregressive filters for analysing multicomponent signals are surveyed. In the second one, we show that the filtering of a multicomponent signal consists of defining the sub-space containing the expected part of the signal and projecting the recorded data on this sub space. With such a representation just a few coefficients are sufficient to ensure a good determination of the physical parameters of pulsation signals. In the last one, an example of detection in multichannel data is presented. The art of geophysics is to know what assumptions may be used to define the set of signals and to test this choice. 1. ESTIMATION
OF PARAMETERS
(i) Single channel A pulsation is of short
duration when the interesting part has less than about S-30 oscillations depending on the number of resonances, the selectivity, their relative amplitude and the signal noise ratio. This is very often the case for Pi2 and Pc3-4-5. For PC recorded on the ground, the stationarity time may be long and the pulsation may have more than 30 oscillations, while in space, because the satellite or the magnetosphere is moving, PC will be more likely to appear as successive distinct short duration signals. As indicated ih the following Table, classical analysis with a lag-window is best applied to long duration signals like Pcl: it is the simplest method and provides information on energy. For a discussion of the accuracy see Glangeaud (1981). Accuracy depends mainly on the value of B, T: B, T must be as large as possible.
T: Time duration of the signal (or stationarity) Bandwidth of the lag-window
B,:
There is a special case where classical analysis can be applied to a short duration signal: it is when the
An invited review paper presented at the 4th IAGA Scientific Assembly, Edinburgh, 14 August 1981.
signal can be considered to be deterministic so that no statistics or averaging have to be applied to get a result.
I
I
I
For short duration, non-deterministic signals the classical analysis cannot give valuable results; for low signal to noise ratio nothing can be done; for signal to noise ratio greater than 5 dB, parametric methods like autoregressive moving average filter (ARMA) and maximum entropy (MEM) can be used. The purpose of this paper is to encourage geophysicists to make more intensive use of these methods by providing some examples. For short duration, non-deterministic signals with high signal to noise ratio, the Pisarenko (1973) method is recommended. This interesting method (for the determination of the position in frequency of pure sinusoids when noise is added) will be only mentioned in this paper. Estimation parameters are presented in this section for single and multiple channels both in one aspect: classical analysis and parametric one. 1249
1250
F. GLANGEAUD
(a) Finite Fourier Transform. Figure 1 shows that the Fourier Transform of a Pcl, ITF(Pcl)J, has regular fluctuations. The BT value of this Pcl sample is 30 (B and T are bandwidth and duration of the signal but not of the data). This means that we need many parameters to completely describe the signal. In fact, to interpret spectral observations, fewer parameters are necessary. 4/S of the original data can be discarded. There are three equivalent ways to do it. in Fig. 1 first the autocorrelation is estimated using the complete data set (Td = 16 min) then a section of the central part of the signal is isolated by multiplication with a lag window O(T) and analysed. Multiplication in the time lag domain (7) is equivalent to convolution or smoothing in the frequency (v) domain so we obtain the spectral density y(v). Another way to do this is to divide the signal into 5 intervals each of length TJS and calculate an average of X(V). Y(V)* for each of the intervals where X(V) and Y(Y) are Fourier transforms of x(t) and y(t) and Y* is the conjugate of Y. For Pcl we see that using the autocorrelation function we lose the effect of echoes. For Pc3-4 smoothing is also necessary since the end of the correlation function is mainly statistical fluctuations. (b) AR methods. In classical analysis, spectral density is estimated by the Fourier Transform of w(T).~(T). The correlation function is well estimated only for small lag times 7; for large T values, there are only statistical fluctuations which are rejected by the window w(T). The frequency resolution is limited by the time width of w. In a way this procedure is equivalent to replacing r(r) by zero values for 7 5, TV. However, in this procedure, the physical signal lacks such a truncated correlation function. For instance, when the geophysicist says that he recorded a sine wave when only a few oscillations can be seen in the recording, he makes the assumption that the signal continues to oscillate for a much longer period of time than the interval of observation. He hypothesizes a long duration correlation function and, unfortunately, uses the classical method: that is a contradiction. The best way to obtain an estimation of the correlation function is to use the few well defined oscillations of r(T) and predict the other part. This prediction supposes that the part of the observed signal is the output of a filter with white noise at the input. For example, broad-band noise in interplanetary space is modified by resonances within the magnetosphere. In digital form this filter can be
written
with a0 = bO= 1 and nk the random series of the noise at the input. Using the z transform the output X(z) is related to the input by the transfer function T(z)== with
z
=
ez’m&
I t, = sampling duration.
The concept of parametric analysis differs completely from classical methods because of a direct estimation of the coefficients aj and bk. The autocorrelation of the signal is as close as possible to r(T) at the beginning and then differs. The interest of the t transform is to deliver a direct estimation of the power spectral density:
with P. the power spectral density of the noise. The transfer function may have three different forms: (i) only zeros (moving average MA) (ii) only poles (Autoregressive AR) (iii) poles and zeros (ARMA). Direct spectral analysis based on pole presentation provides a new analysis tool (Glangeaud et al., 1982). An example of time frequency analysis using the two ways least mean square autoregressive filter (Reid, 1979) appears in Glangeaud (1981). Two ways least mean square autoregressive filter (Reid, 1979). We want to find the variation of the
frequency of a Pcl vs time. The signal may be considered as stationary over M = 13 samples. We make the assumption that there are not more than three frequencies present at the same time so that the order of the AR filter will be 6. For a specific interval we keep a peak in the spectrum if the corresponding energy is significant defined through the use of a threshold. Then, with the technique of moving window, we analyse the event continuously. The possibility of adaptating the length of the memory M to the signal has been introduced so that M increases in the stationary part and provides the greatest accuracy of the frequency.
Analysis of pulsations
*1
-8
Time
mn
-1-a 0.7 Frequency
0.8
0.9 Hz
FIG.1. SPECTRALANALYSISOFPcl RECORDEDIN DELGOCHELIERAND SOCRA BYRASPOFW ANDTROTSKAYA. Tfw autospectrum of Sop f(r) is muItipIied by a window W(T)= (l/2(1 + cos 2&/T))"with II = 16. The corresponding spectrum u(v) is presented before and after multiplication. Bandwidth of the window B, is indicated. Coherency coefficient C(V) between the two stations is evaluated with three different windows: n = 4, 16 and 64.
F.
1252
GLANGEAIJD
Pet are generated by wave particle interactions: the time of this generation is short, but several successive interactions may occur. The results, presented in Fig. 2, give a measurement of the frequencies and the positions in time of the interactions, after a time deIay due to propagation. The results must also be used with care, for the following reasons: (i) Some doubling effects may occur so that a pure frequency may be represented by one frequency on either side of the expected value. (ii) When one component becomes too weak, this frequency disappears, but it changes the result of the next closest component. (iii) The result depends also on the initial condition (phase). (iv) When the signal has been sharply filtered, this method produces parasite frequencies at the cutoff of the band pass filter. Even with these imperfections, the method is still better than classical analysis, and it must be used for time frequency analysis when no assumption on the function y(t) are available. For classical frequency time analysis, see the review by Fraser (1979). (ii) Multichannel (a) Finite Fourier Transfom. Signals with small BT value like Pi2 may be considered as deterministic then a comparison between Pi2 recorded
in two stations may be performed. The comparison between two signals x(t) and y(t) can be performed, without smoothing, by the transfer function rX,l yXX : its modules gives the attenuation and its phase provides information about the phase shift of the time delay. When the value of BT is much different to unity, the coherency coefficient C(v) = yx,lfrxx - x,) is very useful to estimate the relation between two signals. If C(V) = 1, there is a linear and stationary relation between the two signals which can be expressed by the transfer function or complex gain ~_Jy_. Knowing x(t) and the value of this complex gain filter is enough to determine y(t) exactly. Before any interpretation of the value of C(V), we have to check that: -Any time delay between the two signals has been corrected. -The modulus of the three spectral densities yX.v before smoothing have random 3/V, %.yy.y? fluctuations such that the smoothing operation has a real statistical effect. When a pure frequency is present, results may not be significant within the bandwidth B, centered on this frequency. -The time duration product of the analyser (B,T) (or number of averages) must be smaller than the time duration product BT (or number of degrees of freedom of the signal). -For low values of B, T only results of C(V) close to one are significant.
b 16
0 Time,
rnn
FIG. 2. POSITIONOF THE PEAKS OF ENERGYAS A FUNCTIONOF FREQUENCYAND TIME FOR A Pcl EVENT RECORDEDAT HEARD BY GENDRININ 1971. The determination is made by autoregressive filter. The memory length is T. The order of the filter of
6 can provide the determination
of 3 simultaneous frequencies.
Analysis of pulsations
resolution of the maximum entrapy (or AR filter) is l/T = l/32 min. It is not limited by any window and no truncation effect appears because the cross-correlation function is close to zero at the limit of the interval. For this reason, we prefer the determination of the frequencies made by this method to the classical analysis. The signal analysed is supposed to be stationary. For a given length of the signal, the number of poles N to be determined is evaluated by the Akaike criterium but geophysicists should use their own assumptions to choose the best order (Glangeaud et al., 1982). Three dimensional coherency. In the preceeding section, we have used the coherency to determine the relation between two signals. When a multicomponent signal has more than two components, relations between two components can always be studied but generally they do not provide a description of the part of the signal which is the same (through various filters) for all components, that part which is called by Samson and Olson (1980) the total polarized wave. The multicomponent signal may be the various components of an electromagnetic wave recorded at the same point. Or it can be ULF recorded at different stations, in which case we assume that we have been able to compensate for ail the time delays. Samson and Olson have proposed a measurement of the degree of polarization P(Y) for an n-multi-
This last point is illustrated in the bottom part of Fig. 1 where the 90% confidence interval has been drawn for n value of 16 and B,T = 10. For values of C(Y) close to one a slight smoothing (h = 4, B,T = 5) is enough while for low values of C(V) (for v>O.83 Hz) a higher average is necessary (h = 64 and B,T = +20). This last lag-window does not give a good estimation of C(u) at v = 0.78 because C(V) changes rapidly at this frequency and too large a smoothing mixes the different values. C(V) expresses the ratio of the part of the signal y(t) which has a filter relation with ~J’.yrr/yJ~‘* and the total observed part ‘yrvwhich includes the uncorrelated added noise. When C(V) is tow, the uncorrelated part is large; C(V) can be also called the normalized cross spectrum. When two signals have some statistical relationship (good values of C(v)) it is possible to use this dependance in order to either reject the joint part (as in the next paragraph) or keep the common one (see Section 2). Autoregressive filter (Burg). Burg’s Algorithm has provided a two-dimensional program for pulsation studies for many years (Ionnidis, 1975). It is used to obtain spectra of Pe3 to PCS pulsations. Here we compare it with the rut? cross correlation function using Pc3 from Geos which is indicated by the thin line in Fig. 3. The thick line represents the same function estimated by Burg’s method. Both functions are identical for small values of the lag time (less than 3 min), then they differ. The
-16
-3
1253
variate
0
3
signal.
It
can
be
called
a generalized
16 mn
FIG. 3. CROSSCORRELA~ONFUNCTIONBETWEENPc3 RECORDED ABOARDGEOS 1 BY GENDRINAND ON THEGROUNDIN HUSAFELD(ICELAND) BY THE FRENCHMOBILESTATION(S. PERRAUT).
The estimation made using the m~mum
entropy method (Heavy curve) differs from the ruff estimation using a Fourier Transform for time delay bigger than 3 mn.
1254
F. GLANGEAUD
coherency coefficient. The important point is that P(v) is invariable regardless of the set of sensors used while C(Y) is not. One assumes that there is only one geophysical source with some added noise. Each component of the recorded signal has a part which can be related to this unique source by a filter. The study in the presence of several sources is more complicated and the spectral matrix S = (rij) must be diagonilized. In this case, there is more than one eigen value (Lacoume, 1979). The estimation of P(u) proposed by Samson and Olson is very direct uses only the traces of S and S2. For a three dimensional signal:
A good determination of P(V) requires the same 5 conditions as the C(Y) evaluation. The relation between C(v) and P(v) for a two dimensional signal is:
P(U) = C(Y) if which gives the foIlowing OGP(Y)C 1
yII
=
~22
properties:
JC(v)l= 13 IP(
P(v) real,
= 1
C(u)=O;r,P=ObutP(v)
resonance. The shapes of the classical analysis auto and cross spectra are mainly due to windowing effects while the maximum entropy method leads to a nice narrow peak. The coherency modulus at 0.02Hz obtained using the classical method is entirely wrong and the confidence interval has no meaning. The MEM gives a coherency value that is easier to interpret in terms of the relationship between satellite and ground; it is localized to the resonance region. Notice that, from other examples, when the peaks of the two autospectra are at the same frequency, the coherency may also be very low and in this case the corresponding filtered wave forms look very different. Unfortunately the MEM is not perfect. (i) We are not able to put any confidence interval on our results. (ii) Good values of C(v) obtained for small peaks in the autospectrum have to be rejected. (iii) It can be shown that the auto-spectrum, at one station, obtained by this method may differ, depending on the other station used for comparison. (iv) The spectrum is sensitive to any sharp variation of the power within the bandwidth. For example, if the recording instrument cut off has a large inthrence on the spectrum then the MEM will deliver a peak around the cut off to reproduce the spectrum. It is the responsibility of the geophysicist to interpret the results: for example, the coherency between magnetic pulsation and VLF amplitudes has been calculated by Sato and Kokubun (1980; Fig. 11). We may be confident of their result because, for several stations, they find good coherency only when using D components; values obtained with the H components are low; this contrast is a sort of proof of the validity of the conclusion. 2. FILTERHG (i) Single channel There are a lot of publications about filtering. It seems to me interesting to point to two aspects: first, for geophysicists, it is not only the effect of rejection in the frequency domain which is important (dB of the side lobe) but also the duration of the impulse response of the filter in order not to spoil all the signal when an accident is present. Notice that simple filter like those suggested by Oppenheim and Shafer (1975) with constant gain power in the bandwidth and just two transition points are very easy to put into operation.
Analysis of pulsations (ii)
(a) Rejection
of noise with a noise reference. Generally, there is a mixing of signals and noise on each sensor. The signal is what we want to study, the noise the undesirable part. Noise from one sensor to the next may be uncorrelated (in-
strument
noise) or correlated.
We present
here a
specific sations
case where electromagnetic ULF pulhave been recorded under the sea (Lacoume et al., 1981). In Fig. 5b, the main part of the recorded data M = (ULF + noise) is due to the
induction effect of the swell. There is a very good relation between the sea surface movement and this induction effect. So the record of the pressure (5a) in the sea is correlated to the noise N we want to reject.
1255
(b) Data adaptative polarization filter. Notice once again that the term “polarization” is not restricted to electromagnetism; it includes multicomponent signals having a stationary and linear part in all components. The geophysicist likes to see what the signal looks like. The wave form is the best representation. It is much clearer if it is possible to remove the noise. In this study the assumption is that at each time the expected geophysical signal is totally polarized and comes from only one source: then the use of P(U), r, and r2 (axes of the complex ellipse) gives the best filter which sorts out the components of the signal. A simpler (but not optimal) method has been proposed by Samson and Olson (1979): they use the same filter for all components and the gain of
77
26
.02 Frequency(Hz) FIG.~.
SPECTRALDENSITYOFSELECTIVESIGNALLIKE ESTIMATEDBYCLASSICALMETHOD(S)WHILEMAXIMUM
Pc4
.03
(DATAFROMGEOS 1 GENDRIN~ARENOTWELL ENTROPY ME EXHIBITCLEARPEAK.
Comparison of the two signals with coherency (S) give wrong results. In our opinion, the good value of the coherency ME obtained at 0.019 is more in agreement with the expected geophysical theory.
F. GLANGEAUD
1256
1
the bottom part of Fig. 6. This filter is an adaptative filter since P(V) is different for each stating time to. 3.DETECTIONINMULTICHANNELDATA Statevector pr~jecti~~~
FIG.5. ELECTROMAGNETICSIGNAL (b) RECORDED UNDER THE
SEA
AS
CORRELATEU
WELL
AS
PRESSURE
PART WHICH
VARIATIONS
IS EVALUATED
(a)
HAVE
A
(C) BY USE OF THE
FILTER F. The signal (d) has no coherency with (a). This part is related to magnetic variations.
this filter is simply P(Y) A single pass of the data through the filter leaves too much noise. So they have chosen a more selective filter P”(Y) which has the same effect as running g times through the filter P(V). The value of g determined by eye until the filtered components allow a significant geophysica1 interpretation, as shown in Fig. 6. In this example the stationarity time T was chosen to be 21 min about equal to the time duration of a Pi2, since the polarization can change from one Pi2 to the next. From looking at the unsmoothed spectrum (not shown here) one gets B,T = 6. The choice of T and B,T allows the calculation of P(V) for any starting interval (t,, toi- T). The filtered components through p”(v) with g = 9 are presented in
We shall restrict our presentation to a very good example. A theoretical calculation by Southwood and Hughes (1978) predicted the polarization pattern of a three component signal (Pcrl-5, H, D, Z) as a function of the latitudinal distance between the observation point and the source assumed to be on the same meridian. A poiarization is determined either by phase differences and attenuation between any two or by 2 complex numbers components, (attenuation +phase) for D and Z component referenced to H. The model obtained by Southwood and Hughes (1978) does not depend on frequency and can be roughly represented by two discrete states, either North or South of the source, the absolute latitudinal distance having a small effect. For instance, when the source is North, the complex numbers are 0.15-0.21 i for D and 0.25-0.43 i for Z. The equivalent complex vector is [fl); (0.15-0.21 i); (0.25-0.43 i)] with i = (- 1)“2. After normalization of this vector it is possible to calculate the spectral matrix S, of this model when the source is North. This matrix is of course deterministric. We want to know if the observed spectral matrix y(v) is closer to S, or S, in order to locate the source. The energy of a signal is represented by the trace of its spectral matrix, so we compare the trace of the projection of y on S, and S,. In order to work mainly on the signal and reject a part of the noise, Samson and Olson have proposed to use the detector D”: L),(Y) = P(v).Tr [~(Y).S,]. This method is very easy to operate. The comparison of D, and D, vs time and frequency is plotted in Fig. 4 of Samson and Olson (1980). The success of such a simple detector depends on choosing the right geophysical parameters: when the expected stationarity time is 30 min, B,T = 5; the time interval of signals analysed is selected by a moving window. A clear determination of the variation in frequency and time of the position of the source is obtained for two reasons: first, because the Southwood and Hughes model fits well and second, there is a drastic difference between the North and South filters as the angle between them is close to 90”. However this angle
1257
Analysis of pulsations
After
IO
nT
TIME (minutes) FIG. 6. Pi2 (H, D, Z) RECORDED SAMPLE
METHOD
OF SAMSON
AND
BY SAMSON ARE PRESENTED BEFORE AND AFTER FILTERING BY THE OLSON WHICH KEEP ONLY PART OF SIGNAL CONTAINED IN THE POLARIZATlON 3ANDWioTH.
(After Samson and Olson, 1979.)
is not exactly 90", thus the S, model does not have a zero projection on S,. Because of this, the method just described is not optimal and more sophisticated program are available from Samson and Olson which takes this fact into account. For a very short burst of pulsations all these methods can be still applied but parametric analysis must be used for the determination of the spectral matrix. Other methods using eigenvalues of the spectral matrix and additional assumptions on the sources
have been developed (Lacoume ef aI., 1982) but they take a long time on computer. CONCLUSION
Analysis of geophysical data needs the use of many different spectral methods. Depending on the value of BT either classical analysis using Fourier Transforms or parametric methods have to be used. Geophysicists have to transform their physical
1258
F. GLA.NGEAUD
assumptions about the signal into a clear determination of the parameters to be introduced in the analysis (duration of stationarity in time and in frequency.. .). This is the key to success, We hope, with the few examples presented in this paper, to convince geophysicists to use spectral analysis more and more. Acknawtedgements-I thank the referee of the original version of this paper for his helpful and constructive criticism. RRFERENC~ Akaike, (1974) A new look at the statistical model identification. IEEE Trans. of Automatic Control. Vol. AC-19, No. 6. Burg, J. P. (1967) Maximum entropy spectral analysis. Proceedings of the 37th meeting of the Society of Exploration Geophysicists. Fraser, B. J. (1979) The use of multistation networks in geomagnetic pulsation studies. Ann. Tkie’commun. 34, 220. Glangeaud, F. (1981) Signal Processing for Magnetic Pulsations _r. atmos. terr. Pftys. 43(9), 981. Glangeaud, F., Gharbi, M., Martin N. Lacoume, J. L. (1982) Use of multidimensional MEM spectral analysis in geophysics. Proceedings of ICASSPIZ, Vol. 3, pp. 1886-1889. Ioannidis, G. A. (1975) Application of multivariate autoregressive spectrum estimation to ULF waves. Radio Sci. 10, 1043.
Lacoume, J. L., Bouthemy, B., Glangeaud, F., Latombe, C. and Silvent, A. (1979) Caracterisation par analyse interspectrale du champ d’ondes recu sur un rCseau de capteurs. Applications. Colloque GRETSI. Nice. Lacoume, J. L., Bouthemy, B., Glangeaud, F. and Latombe. C. (1982) Use of snectral matrix for sources identification. Proceedings of the Inst. of Acoustics, Underwater Acoustics Group Conference, London, April 1982. Lacoume, J. L., Glangeaud, F., Lorenzino, P., Baudois, D. and Pretet G. (1981) Filtrage de signaux mufticomposantes utilisant les correlations intercomposantes. GRETSI, Huitieme Colloque sur le Traitment du Signal et ses Applications, Nice, 31 l-318. Oppenheim, A. V. and Schafer, R. W. (1975) Digital Signal Processing. Prentice Hail, New York. Pisarenko, V. F. (1973) The retrieval of Harmonies from a Covariance Function. Geophys. .J. R. astr. Sot. 33, 347-366. Reid, J. S. (1979) _r.geophys. Res. 84,5289. Samson, J. C. and Olson, J. V. (1979) Data-adaptive PolarizationFilters for Multichannel Geophysical Data. University of Alberta, Canada. Samson, J. C. and Olson, J. V. (1980) Generalized power spectra and the stokes vector representations of ultralow frequency micropulsation states. Can. J. Phys. 58, 123. Sate,. N. and Kokubun, S. (1980) Interaction between ELF-VLF emissions and magnetic pulsations. L geophys. Res. 85, 101. Southwood, D. J. and Hughes, W. J. (1978) Source induced vertical components in geomagnetic pulsation signals. Planet. Space Sci. 26, 715.
SpaceSci., in Great
Vol. 30, No. Britain.
12, pp. 1249-1258,
0032-0633/82/12124~1ou)3.00/0 Pergamon Press Ltd.
1982
ANALYSIS
OF PULSATIONS F. GLANGEAUD
CEPHAG, BP 46, 38402 St Martin d’Heures, France LA 346 au CNRS (Received in fina/ form 7 June 1982) Abstract-There is no universal method of data analysis. The scientist must first know what he intends to measure; then, it is possible to select the most suitable methods and give practical applications. We limit our investigation here to three domains. In the first one, estimation of the parameters, the number of degrees of freedom is estimated, taking into account the a priori knowledge in magnetic pulsation models. The random and deterministic parts of the signal are evaluated and the meaning of the confidence interval of the measurements is discussed. The use and limitations of maximum entropy and autoregressive filters for analysing multicomponent signals are surveyed. In the second one, we show that the filtering of a multicomponent signal consists of defining the sub-space containing the expected part of the signal and projecting the recorded data on this sub space. With such a representation just a few coefficients are sufficient to ensure a good determination of the physical parameters of pulsation signals. In the last one, an example of detection in multichannel data is presented. The art of geophysics is to know what assumptions may be used to define the set of signals and to test this choice. 1. ESTIMATION
OF PARAMETERS
(i) Single channel A pulsation is of short
duration when the interesting part has less than about S-30 oscillations depending on the number of resonances, the selectivity, their relative amplitude and the signal noise ratio. This is very often the case for Pi2 and Pc3-4-5. For PC recorded on the ground, the stationarity time may be long and the pulsation may have more than 30 oscillations, while in space, because the satellite or the magnetosphere is moving, PC will be more likely to appear as successive distinct short duration signals. As indicated ih the following Table, classical analysis with a lag-window is best applied to long duration signals like Pcl: it is the simplest method and provides information on energy. For a discussion of the accuracy see Glangeaud (1981). Accuracy depends mainly on the value of B, T: B, T must be as large as possible.
T: Time duration of the signal (or stationarity) Bandwidth of the lag-window
B,:
There is a special case where classical analysis can be applied to a short duration signal: it is when the
An invited review paper presented at the 4th IAGA Scientific Assembly, Edinburgh, 14 August 1981.
signal can be considered to be deterministic so that no statistics or averaging have to be applied to get a result.
I
I
I
For short duration, non-deterministic signals the classical analysis cannot give valuable results; for low signal to noise ratio nothing can be done; for signal to noise ratio greater than 5 dB, parametric methods like autoregressive moving average filter (ARMA) and maximum entropy (MEM) can be used. The purpose of this paper is to encourage geophysicists to make more intensive use of these methods by providing some examples. For short duration, non-deterministic signals with high signal to noise ratio, the Pisarenko (1973) method is recommended. This interesting method (for the determination of the position in frequency of pure sinusoids when noise is added) will be only mentioned in this paper. Estimation parameters are presented in this section for single and multiple channels both in one aspect: classical analysis and parametric one. 1249
1250
F. GLANGEAUD
(a) Finite Fourier Transform. Figure 1 shows that the Fourier Transform of a Pcl, ITF(Pcl)J, has regular fluctuations. The BT value of this Pcl sample is 30 (B and T are bandwidth and duration of the signal but not of the data). This means that we need many parameters to completely describe the signal. In fact, to interpret spectral observations, fewer parameters are necessary. 4/S of the original data can be discarded. There are three equivalent ways to do it. in Fig. 1 first the autocorrelation is estimated using the complete data set (Td = 16 min) then a section of the central part of the signal is isolated by multiplication with a lag window O(T) and analysed. Multiplication in the time lag domain (7) is equivalent to convolution or smoothing in the frequency (v) domain so we obtain the spectral density y(v). Another way to do this is to divide the signal into 5 intervals each of length TJS and calculate an average of X(V). Y(V)* for each of the intervals where X(V) and Y(Y) are Fourier transforms of x(t) and y(t) and Y* is the conjugate of Y. For Pcl we see that using the autocorrelation function we lose the effect of echoes. For Pc3-4 smoothing is also necessary since the end of the correlation function is mainly statistical fluctuations. (b) AR methods. In classical analysis, spectral density is estimated by the Fourier Transform of w(T).~(T). The correlation function is well estimated only for small lag times 7; for large T values, there are only statistical fluctuations which are rejected by the window w(T). The frequency resolution is limited by the time width of w. In a way this procedure is equivalent to replacing r(r) by zero values for 7 5, TV. However, in this procedure, the physical signal lacks such a truncated correlation function. For instance, when the geophysicist says that he recorded a sine wave when only a few oscillations can be seen in the recording, he makes the assumption that the signal continues to oscillate for a much longer period of time than the interval of observation. He hypothesizes a long duration correlation function and, unfortunately, uses the classical method: that is a contradiction. The best way to obtain an estimation of the correlation function is to use the few well defined oscillations of r(T) and predict the other part. This prediction supposes that the part of the observed signal is the output of a filter with white noise at the input. For example, broad-band noise in interplanetary space is modified by resonances within the magnetosphere. In digital form this filter can be
written
with a0 = bO= 1 and nk the random series of the noise at the input. Using the z transform the output X(z) is related to the input by the transfer function T(z)== with
z
=
ez’m&
I t, = sampling duration.
The concept of parametric analysis differs completely from classical methods because of a direct estimation of the coefficients aj and bk. The autocorrelation of the signal is as close as possible to r(T) at the beginning and then differs. The interest of the t transform is to deliver a direct estimation of the power spectral density:
with P. the power spectral density of the noise. The transfer function may have three different forms: (i) only zeros (moving average MA) (ii) only poles (Autoregressive AR) (iii) poles and zeros (ARMA). Direct spectral analysis based on pole presentation provides a new analysis tool (Glangeaud et al., 1982). An example of time frequency analysis using the two ways least mean square autoregressive filter (Reid, 1979) appears in Glangeaud (1981). Two ways least mean square autoregressive filter (Reid, 1979). We want to find the variation of the
frequency of a Pcl vs time. The signal may be considered as stationary over M = 13 samples. We make the assumption that there are not more than three frequencies present at the same time so that the order of the AR filter will be 6. For a specific interval we keep a peak in the spectrum if the corresponding energy is significant defined through the use of a threshold. Then, with the technique of moving window, we analyse the event continuously. The possibility of adaptating the length of the memory M to the signal has been introduced so that M increases in the stationary part and provides the greatest accuracy of the frequency.
Analysis of pulsations
*1
-8
Time
mn
-1-a 0.7 Frequency
0.8
0.9 Hz
FIG.1. SPECTRALANALYSISOFPcl RECORDEDIN DELGOCHELIERAND SOCRA BYRASPOFW ANDTROTSKAYA. Tfw autospectrum of Sop f(r) is muItipIied by a window W(T)= (l/2(1 + cos 2&/T))"with II = 16. The corresponding spectrum u(v) is presented before and after multiplication. Bandwidth of the window B, is indicated. Coherency coefficient C(V) between the two stations is evaluated with three different windows: n = 4, 16 and 64.
F.
1252
GLANGEAIJD
Pet are generated by wave particle interactions: the time of this generation is short, but several successive interactions may occur. The results, presented in Fig. 2, give a measurement of the frequencies and the positions in time of the interactions, after a time deIay due to propagation. The results must also be used with care, for the following reasons: (i) Some doubling effects may occur so that a pure frequency may be represented by one frequency on either side of the expected value. (ii) When one component becomes too weak, this frequency disappears, but it changes the result of the next closest component. (iii) The result depends also on the initial condition (phase). (iv) When the signal has been sharply filtered, this method produces parasite frequencies at the cutoff of the band pass filter. Even with these imperfections, the method is still better than classical analysis, and it must be used for time frequency analysis when no assumption on the function y(t) are available. For classical frequency time analysis, see the review by Fraser (1979). (ii) Multichannel (a) Finite Fourier Transfom. Signals with small BT value like Pi2 may be considered as deterministic then a comparison between Pi2 recorded
in two stations may be performed. The comparison between two signals x(t) and y(t) can be performed, without smoothing, by the transfer function rX,l yXX : its modules gives the attenuation and its phase provides information about the phase shift of the time delay. When the value of BT is much different to unity, the coherency coefficient C(v) = yx,lfrxx - x,) is very useful to estimate the relation between two signals. If C(V) = 1, there is a linear and stationary relation between the two signals which can be expressed by the transfer function or complex gain ~_Jy_. Knowing x(t) and the value of this complex gain filter is enough to determine y(t) exactly. Before any interpretation of the value of C(V), we have to check that: -Any time delay between the two signals has been corrected. -The modulus of the three spectral densities yX.v before smoothing have random 3/V, %.yy.y? fluctuations such that the smoothing operation has a real statistical effect. When a pure frequency is present, results may not be significant within the bandwidth B, centered on this frequency. -The time duration product of the analyser (B,T) (or number of averages) must be smaller than the time duration product BT (or number of degrees of freedom of the signal). -For low values of B, T only results of C(V) close to one are significant.
b 16
0 Time,
rnn
FIG. 2. POSITIONOF THE PEAKS OF ENERGYAS A FUNCTIONOF FREQUENCYAND TIME FOR A Pcl EVENT RECORDEDAT HEARD BY GENDRININ 1971. The determination is made by autoregressive filter. The memory length is T. The order of the filter of
6 can provide the determination
of 3 simultaneous frequencies.
Analysis of pulsations
resolution of the maximum entrapy (or AR filter) is l/T = l/32 min. It is not limited by any window and no truncation effect appears because the cross-correlation function is close to zero at the limit of the interval. For this reason, we prefer the determination of the frequencies made by this method to the classical analysis. The signal analysed is supposed to be stationary. For a given length of the signal, the number of poles N to be determined is evaluated by the Akaike criterium but geophysicists should use their own assumptions to choose the best order (Glangeaud et al., 1982). Three dimensional coherency. In the preceeding section, we have used the coherency to determine the relation between two signals. When a multicomponent signal has more than two components, relations between two components can always be studied but generally they do not provide a description of the part of the signal which is the same (through various filters) for all components, that part which is called by Samson and Olson (1980) the total polarized wave. The multicomponent signal may be the various components of an electromagnetic wave recorded at the same point. Or it can be ULF recorded at different stations, in which case we assume that we have been able to compensate for ail the time delays. Samson and Olson have proposed a measurement of the degree of polarization P(Y) for an n-multi-
This last point is illustrated in the bottom part of Fig. 1 where the 90% confidence interval has been drawn for n value of 16 and B,T = 10. For values of C(Y) close to one a slight smoothing (h = 4, B,T = 5) is enough while for low values of C(V) (for v>O.83 Hz) a higher average is necessary (h = 64 and B,T = +20). This last lag-window does not give a good estimation of C(u) at v = 0.78 because C(V) changes rapidly at this frequency and too large a smoothing mixes the different values. C(V) expresses the ratio of the part of the signal y(t) which has a filter relation with ~J’.yrr/yJ~‘* and the total observed part ‘yrvwhich includes the uncorrelated added noise. When C(V) is tow, the uncorrelated part is large; C(V) can be also called the normalized cross spectrum. When two signals have some statistical relationship (good values of C(v)) it is possible to use this dependance in order to either reject the joint part (as in the next paragraph) or keep the common one (see Section 2). Autoregressive filter (Burg). Burg’s Algorithm has provided a two-dimensional program for pulsation studies for many years (Ionnidis, 1975). It is used to obtain spectra of Pe3 to PCS pulsations. Here we compare it with the rut? cross correlation function using Pc3 from Geos which is indicated by the thin line in Fig. 3. The thick line represents the same function estimated by Burg’s method. Both functions are identical for small values of the lag time (less than 3 min), then they differ. The
-16
-3
1253
variate
0
3
signal.
It
can
be
called
a generalized
16 mn
FIG. 3. CROSSCORRELA~ONFUNCTIONBETWEENPc3 RECORDED ABOARDGEOS 1 BY GENDRINAND ON THEGROUNDIN HUSAFELD(ICELAND) BY THE FRENCHMOBILESTATION(S. PERRAUT).
The estimation made using the m~mum
entropy method (Heavy curve) differs from the ruff estimation using a Fourier Transform for time delay bigger than 3 mn.
1254
F. GLANGEAUD
coherency coefficient. The important point is that P(v) is invariable regardless of the set of sensors used while C(Y) is not. One assumes that there is only one geophysical source with some added noise. Each component of the recorded signal has a part which can be related to this unique source by a filter. The study in the presence of several sources is more complicated and the spectral matrix S = (rij) must be diagonilized. In this case, there is more than one eigen value (Lacoume, 1979). The estimation of P(u) proposed by Samson and Olson is very direct uses only the traces of S and S2. For a three dimensional signal:
A good determination of P(V) requires the same 5 conditions as the C(Y) evaluation. The relation between C(v) and P(v) for a two dimensional signal is:
P(U) = C(Y) if which gives the foIlowing OGP(Y)C 1
yII
=
~22
properties:
JC(v)l= 13 IP(
P(v) real,
= 1
C(u)=O;r,P=ObutP(v)
resonance. The shapes of the classical analysis auto and cross spectra are mainly due to windowing effects while the maximum entropy method leads to a nice narrow peak. The coherency modulus at 0.02Hz obtained using the classical method is entirely wrong and the confidence interval has no meaning. The MEM gives a coherency value that is easier to interpret in terms of the relationship between satellite and ground; it is localized to the resonance region. Notice that, from other examples, when the peaks of the two autospectra are at the same frequency, the coherency may also be very low and in this case the corresponding filtered wave forms look very different. Unfortunately the MEM is not perfect. (i) We are not able to put any confidence interval on our results. (ii) Good values of C(v) obtained for small peaks in the autospectrum have to be rejected. (iii) It can be shown that the auto-spectrum, at one station, obtained by this method may differ, depending on the other station used for comparison. (iv) The spectrum is sensitive to any sharp variation of the power within the bandwidth. For example, if the recording instrument cut off has a large inthrence on the spectrum then the MEM will deliver a peak around the cut off to reproduce the spectrum. It is the responsibility of the geophysicist to interpret the results: for example, the coherency between magnetic pulsation and VLF amplitudes has been calculated by Sato and Kokubun (1980; Fig. 11). We may be confident of their result because, for several stations, they find good coherency only when using D components; values obtained with the H components are low; this contrast is a sort of proof of the validity of the conclusion. 2. FILTERHG (i) Single channel There are a lot of publications about filtering. It seems to me interesting to point to two aspects: first, for geophysicists, it is not only the effect of rejection in the frequency domain which is important (dB of the side lobe) but also the duration of the impulse response of the filter in order not to spoil all the signal when an accident is present. Notice that simple filter like those suggested by Oppenheim and Shafer (1975) with constant gain power in the bandwidth and just two transition points are very easy to put into operation.
Analysis of pulsations (ii)
(a) Rejection
of noise with a noise reference. Generally, there is a mixing of signals and noise on each sensor. The signal is what we want to study, the noise the undesirable part. Noise from one sensor to the next may be uncorrelated (in-
strument
noise) or correlated.
We present
here a
specific sations
case where electromagnetic ULF pulhave been recorded under the sea (Lacoume et al., 1981). In Fig. 5b, the main part of the recorded data M = (ULF + noise) is due to the
induction effect of the swell. There is a very good relation between the sea surface movement and this induction effect. So the record of the pressure (5a) in the sea is correlated to the noise N we want to reject.
1255
(b) Data adaptative polarization filter. Notice once again that the term “polarization” is not restricted to electromagnetism; it includes multicomponent signals having a stationary and linear part in all components. The geophysicist likes to see what the signal looks like. The wave form is the best representation. It is much clearer if it is possible to remove the noise. In this study the assumption is that at each time the expected geophysical signal is totally polarized and comes from only one source: then the use of P(U), r, and r2 (axes of the complex ellipse) gives the best filter which sorts out the components of the signal. A simpler (but not optimal) method has been proposed by Samson and Olson (1979): they use the same filter for all components and the gain of
77
26
.02 Frequency(Hz) FIG.~.
SPECTRALDENSITYOFSELECTIVESIGNALLIKE ESTIMATEDBYCLASSICALMETHOD(S)WHILEMAXIMUM
Pc4
.03
(DATAFROMGEOS 1 GENDRIN~ARENOTWELL ENTROPY ME EXHIBITCLEARPEAK.
Comparison of the two signals with coherency (S) give wrong results. In our opinion, the good value of the coherency ME obtained at 0.019 is more in agreement with the expected geophysical theory.
F. GLANGEAUD
1256
1
the bottom part of Fig. 6. This filter is an adaptative filter since P(V) is different for each stating time to. 3.DETECTIONINMULTICHANNELDATA Statevector pr~jecti~~~
FIG.5. ELECTROMAGNETICSIGNAL (b) RECORDED UNDER THE
SEA
AS
CORRELATEU
WELL
AS
PRESSURE
PART WHICH
VARIATIONS
IS EVALUATED
(a)
HAVE
A
(C) BY USE OF THE
FILTER F. The signal (d) has no coherency with (a). This part is related to magnetic variations.
this filter is simply P(Y) A single pass of the data through the filter leaves too much noise. So they have chosen a more selective filter P”(Y) which has the same effect as running g times through the filter P(V). The value of g determined by eye until the filtered components allow a significant geophysica1 interpretation, as shown in Fig. 6. In this example the stationarity time T was chosen to be 21 min about equal to the time duration of a Pi2, since the polarization can change from one Pi2 to the next. From looking at the unsmoothed spectrum (not shown here) one gets B,T = 6. The choice of T and B,T allows the calculation of P(V) for any starting interval (t,, toi- T). The filtered components through p”(v) with g = 9 are presented in
We shall restrict our presentation to a very good example. A theoretical calculation by Southwood and Hughes (1978) predicted the polarization pattern of a three component signal (Pcrl-5, H, D, Z) as a function of the latitudinal distance between the observation point and the source assumed to be on the same meridian. A poiarization is determined either by phase differences and attenuation between any two or by 2 complex numbers components, (attenuation +phase) for D and Z component referenced to H. The model obtained by Southwood and Hughes (1978) does not depend on frequency and can be roughly represented by two discrete states, either North or South of the source, the absolute latitudinal distance having a small effect. For instance, when the source is North, the complex numbers are 0.15-0.21 i for D and 0.25-0.43 i for Z. The equivalent complex vector is [fl); (0.15-0.21 i); (0.25-0.43 i)] with i = (- 1)“2. After normalization of this vector it is possible to calculate the spectral matrix S, of this model when the source is North. This matrix is of course deterministric. We want to know if the observed spectral matrix y(v) is closer to S, or S, in order to locate the source. The energy of a signal is represented by the trace of its spectral matrix, so we compare the trace of the projection of y on S, and S,. In order to work mainly on the signal and reject a part of the noise, Samson and Olson have proposed to use the detector D”: L),(Y) = P(v).Tr [~(Y).S,]. This method is very easy to operate. The comparison of D, and D, vs time and frequency is plotted in Fig. 4 of Samson and Olson (1980). The success of such a simple detector depends on choosing the right geophysical parameters: when the expected stationarity time is 30 min, B,T = 5; the time interval of signals analysed is selected by a moving window. A clear determination of the variation in frequency and time of the position of the source is obtained for two reasons: first, because the Southwood and Hughes model fits well and second, there is a drastic difference between the North and South filters as the angle between them is close to 90”. However this angle
1257
Analysis of pulsations
After
IO
nT
TIME (minutes) FIG. 6. Pi2 (H, D, Z) RECORDED SAMPLE
METHOD
OF SAMSON
AND
BY SAMSON ARE PRESENTED BEFORE AND AFTER FILTERING BY THE OLSON WHICH KEEP ONLY PART OF SIGNAL CONTAINED IN THE POLARIZATlON 3ANDWioTH.
(After Samson and Olson, 1979.)
is not exactly 90", thus the S, model does not have a zero projection on S,. Because of this, the method just described is not optimal and more sophisticated program are available from Samson and Olson which takes this fact into account. For a very short burst of pulsations all these methods can be still applied but parametric analysis must be used for the determination of the spectral matrix. Other methods using eigenvalues of the spectral matrix and additional assumptions on the sources
have been developed (Lacoume ef aI., 1982) but they take a long time on computer. CONCLUSION
Analysis of geophysical data needs the use of many different spectral methods. Depending on the value of BT either classical analysis using Fourier Transforms or parametric methods have to be used. Geophysicists have to transform their physical
1258
F. GLA.NGEAUD
assumptions about the signal into a clear determination of the parameters to be introduced in the analysis (duration of stationarity in time and in frequency.. .). This is the key to success, We hope, with the few examples presented in this paper, to convince geophysicists to use spectral analysis more and more. Acknawtedgements-I thank the referee of the original version of this paper for his helpful and constructive criticism. RRFERENC~ Akaike, (1974) A new look at the statistical model identification. IEEE Trans. of Automatic Control. Vol. AC-19, No. 6. Burg, J. P. (1967) Maximum entropy spectral analysis. Proceedings of the 37th meeting of the Society of Exploration Geophysicists. Fraser, B. J. (1979) The use of multistation networks in geomagnetic pulsation studies. Ann. Tkie’commun. 34, 220. Glangeaud, F. (1981) Signal Processing for Magnetic Pulsations _r. atmos. terr. Pftys. 43(9), 981. Glangeaud, F., Gharbi, M., Martin N. Lacoume, J. L. (1982) Use of multidimensional MEM spectral analysis in geophysics. Proceedings of ICASSPIZ, Vol. 3, pp. 1886-1889. Ioannidis, G. A. (1975) Application of multivariate autoregressive spectrum estimation to ULF waves. Radio Sci. 10, 1043.
Lacoume, J. L., Bouthemy, B., Glangeaud, F., Latombe, C. and Silvent, A. (1979) Caracterisation par analyse interspectrale du champ d’ondes recu sur un rCseau de capteurs. Applications. Colloque GRETSI. Nice. Lacoume, J. L., Bouthemy, B., Glangeaud, F. and Latombe. C. (1982) Use of snectral matrix for sources identification. Proceedings of the Inst. of Acoustics, Underwater Acoustics Group Conference, London, April 1982. Lacoume, J. L., Glangeaud, F., Lorenzino, P., Baudois, D. and Pretet G. (1981) Filtrage de signaux mufticomposantes utilisant les correlations intercomposantes. GRETSI, Huitieme Colloque sur le Traitment du Signal et ses Applications, Nice, 31 l-318. Oppenheim, A. V. and Schafer, R. W. (1975) Digital Signal Processing. Prentice Hail, New York. Pisarenko, V. F. (1973) The retrieval of Harmonies from a Covariance Function. Geophys. .J. R. astr. Sot. 33, 347-366. Reid, J. S. (1979) _r.geophys. Res. 84,5289. Samson, J. C. and Olson, J. V. (1979) Data-adaptive PolarizationFilters for Multichannel Geophysical Data. University of Alberta, Canada. Samson, J. C. and Olson, J. V. (1980) Generalized power spectra and the stokes vector representations of ultralow frequency micropulsation states. Can. J. Phys. 58, 123. Sate,. N. and Kokubun, S. (1980) Interaction between ELF-VLF emissions and magnetic pulsations. L geophys. Res. 85, 101. Southwood, D. J. and Hughes, W. J. (1978) Source induced vertical components in geomagnetic pulsation signals. Planet. Space Sci. 26, 715.