Best-Basis Analysis of Broadband Tremor Signals

Active volcanoes usually generate highly non-stationary broadband tremor signals. Short-time shock events with a frequency content of several decades are superimposed on a stationary narrow band continuous tremor. Tremor signals of this type can be observed in the near field of many active volcanoes. In this paper we will demonstrate the analysis of such signals using a specific tremor signal of Mt. Stromboli (Sicily). We used the Best-Basis Algorithm (BBA) in order to compute a spectrogram which is adapted to signal properties on highly different scales. It turns out that the BBA can reveal better fitting properties of the tremor in the time-frequency plane compared to standard methods like Short-Time Fourier Transformation (STFT). Moreover, this very effective algorithm can be used for real time monitoring in the time-frequency plane, for data compression or for de-noising of the tremor signals.


Introduction
A conventional approach in order to reveal the properties of non-stationary signals consists of splitting up the signal into adjacent unit length time intervals.When quasi-stationary behaviour of the signal within an interval can be assumed a transformation for the interval can be applied to study the signal properties.The temporal variations are due to the subsequent time intervals.A standard method like STFT is based on this concept.
Problems arise e.g. when due to information on different time scales the signal can not be divided into quasi-stationary blocks.In this case a more exible concept which is provided by the adaptive waveform analysis as introduced e.g. in (Coifman and Wickerhauser, 1992), (Mallat and Zhang, 1993), (Wickerhauser, 1991) or (Wickerhauser, 1996) has to be used.
The Adapted waveform analysis describes how to decompose a signal with respect to a dictionary of basis functions.The properties of these basis functions in uence mainly the signal features which can be extracted.Assuming a dictionary with an inherent binary tree structure (see gure 1) e cient algorithms like the BBA can be applied.
Spectrograms related to adapted waveform analysis result in a more suitable resolution of time and frequency compared to STFT or Wavelet Transform (WT).Therefore, a frequency range of several decades can be resolved properly.Due to the wide frequency range of seismic broadband records these spectrograms proved to be very useful for obtaining further insight into the nature of tremor signals.

Wavelet Packets
In principle, basis dictionaries may consist of any type or family of orthogonal or non-orthogonal basis functions.Examples are impulses, harmonic signals, modulated windows at di erent time-scales, wavelets, or wavelet packets.A generalization of the discrete WT (Louis et al., 1994), (Chui, 1992), (Daubechies, 1992), (Strang and Nguyen, 1996) leads to the Wavelet Packet Transform (WPT) (Coifman and Wickerhauser, 1992), (Mallat and Zhang, 1993), (Wickerhauser, 1991), (Wickerhauser, 1996) which has been used in this approach since the resulting algorithm is very e cient.In the sequel, a short introductory description is given.For further details see (Wickerhauser, 1996).
To generate wavelets, the scheme of multi- Figure 1: A nested sequence of subspaces (approximation spaces) V0, V?1, V?2, ... and their related orthogonal complements (wavelet or detail spaces) W?1, W?2, ... span a MRA.The projection of a discrete signal f(n) onto these spaces is computed by ltering and down-sampling (operators G and H).By additionally decomposing all the projections onto wavelet spaces an over-complete collection of projections called a wavelet packet table is generated (gray arrows).
resolution analysis (MRA) can be used to provide a wide variety of di erent wavelets with di erent properties.A MRA consists of a sequence of nested subspaces V j , called approximation spaces, f0g V j ?1 V j V j+1 L 2 (R) (1) as depicted in gure 1.These spaces possess the following properties: which represents a ladder between spaces at di erent levels.(t) 2 V 0 can be decomposed by the basis functions (2t ?k) 2 V 1 as V 0 V 1 .The orthogonal complement of each approximation space V j with respect to V j+1 is called a detail space W j .It is spanned by translates of the function (t) termed wavelet.It can be stated: V j ?W j ; W i ?W j 8 i 6 = j V j+1 = V j W j : Therefore, it is obvious that (t) can be written as (4) Normalizing (t) and using the notation j;k (t) = 2 j=2 (2 j t ?k); (5a) j;k (t) = 2 j=2 (2 j t ?k); (5b) the sets of j;k (t); k 2 Z and j;k (t); k 2 Z constitute orthonormal bases of V j and W j , respectively.
In order to establish a MRA an appropriate set of coe cients a(k) has to be found by means of which a solution (t) of (3) can be obtained.(?1) k a(1 ?k), see (Daubechies, 1992).Further additional restrictions like e.g.orthogonality, biorthogonality, number of vanishing moments, compact support or smoothness of the basis functions (t) result in di erent types of wavelets.
The decomposition starts by projecting a signal (6) Then the coe cients of the projection onto the next coarser space V ?1 are given by c ?1 (k) = Z f (t) ?1;k (t) dt: Using (3,5a) and exploiting the orthogonality of the j;k (t) it can be shown that Obviously, the c ?1 (k) are obtained by convolution of a(?k) and c 0 (k) followed by a decimation of factor 2.
The operator which produces approximations on different levels acts as a low-pass lter (H) while the operator (G) can be recognized as a high-pass lter.The coe cients of a MRA can e ciently be calculated by successively ltering the coe cients of all approximations followed by a down-sampling operation.This procedure is called Mallat's Algorithm with asymptotically less computations compared with the FFT (Louis et al., 1994).
An appropriate choice of successive approximation and detail spaces results in the representation of signals by wavelet packets.The resulting expansion is called a wavelet packet table, see gure 1.An overcomplete collection of high-pass, low-pass and bandpass basis systems is provided.Because of the inherent binary tree structure of the basis systems e cient algorithms exist.E.g. the Best Basis (Coifman and Wickerhauser, 1992), Matching Pursuit (Mallat and Zhang, 1993) or Atomic Pursuit (Chen et al., 1996) algorithms select an `optimal' decomposition of the signal out of this dictionary.Because of less computational a ord we applied BBA.

Coifman Wavelets (Coi ets)
The probably most famous wavelets were introduced by I. Daubechies (Daubechies, 1988), (Daubechies, 1992).They constitute an orthonormal family and are obtained by the requirements

and
(t) should be as regular as possible.Similar to the Daubechies' wavelets the Coifman wavelets, often called `Coi ets', lead to a MRA.They are obtained by the additional requirement of vanishing moments for the scaling function Their support length (time localization) compared with other MRA wavelets increases.On the other hand an increasing number of vanishing moments provides increasing symmetry properties of the wavelet and the scaling function and results in a decreasing approximation error in (6) when the signal samples are used as the approximation coe cients of the initial space (Louis et al., 1994).In comparison to the Daubechies' lters the Coi et lters show approximately almost linear phase.For this reason they are very useful since they induce less phase distortions and therefore the center of the coe cients when interpreted in the time-frequency plane can be well estimated (Wickerhauser, 1996).Otherwise the stop band damping is slightly better in case of the Daubechies' lter which results in decreasing aliasing e ects due to down-sampling, see gure 3.
A Coi et lter with length 30 has been chosen in order to obtain a good trade-o between the support length of the wavelet and the localization properties in the frequency domain, see gure 2. The scaling function and the wavelet function are orthogonal for integer translates k; k 2 Z.The frequency localization is quite acceptable.By generating a wavelet packet dictionary out of Coi ets the maximum decomposition level has to be controlled because of the spreading of the packets in the frequency domain (Wickerhauser, 1994).Figure 4 shows a Coi et packet basis function at level 5 which is badly localized in the frequency domain.Apparently, an interpretation as a single tile in the time-frequency plane is not given.

Best-Basis Algorithm
By decomposing a signal in a wavelet packet table we obtain the representations v d;b (k) of the signal in the over-complete collection of subspaces (sub-bands).d denotes the decomposition level and b counts the basis blocks from left to right ( gure 1).The table v d;b (k) contains a lot of redundancy.Therefore, an algorithm is needed to select an orthonormal basis which describes the whole L 2 and which nds a basis which is best adapted to the signal with respect to a cost functional.Because of the rst two requirements the sum of adjacent sub-bands must not have  gaps or overlaps.See gure 5 for two valid selections.Many well-known criteria can be used as information cost functional (Wickerhauser, 1996).In the present approach Shannon's entropy cost functional has been used.For each basis block the entropy cost is calculated.The algorithm starts at the bottom (see gure 5) and compares the sum of entropy of the two children basis systems with the entropy cost of the parent basis.A minimal cost factor indicates a more peaked and therefore better tting basis.

Time-Frequency Plane
To compute a time-frequency representation the bandwidth and time duration of a rectangular tile related to a coe cient has to be known.At a constant decomposition level each basis block (sub-band) is spanned by basis functions with constant bandwidth.Therefore, all coe cients in a basis block at level d and block number b 2 f0; ; 2 ?d ?1g are related  to tiles which approximately cover the frequency interval b 2 1?d f s , (b+1) 2 1?d f s ] (f s sampling-frequency).The initial space is spanned by orthonormal integer shifted Dirac impulses.To get a non-overlapping tiling along the time axis an interval of 1=f s is chosen.
Using (5) we get the time intervals k ?2 ?d?1 =f s ; k + 2 ?d?1 =f s ] for a tile which corresponds to a coe cient at position k in a basis block at level d.The position k has to be corrected according to the phase distortion which occurred to the v d;b (k) after applying a cascade of operators H and G (Wickerhauser, 1996).
Each tile is mapped with a color related to its amplitude (black is highest) in the mentioned intervals.The resulting representation is called idealized timefrequency plane as we did not calculated the dimensions of the rectangular tiles belonging to the used Coi et (Wickerhauser, 1996).
However, this interpretation of the coe cients does not take into account the misbehavior of wavelet packets at higher levels in the frequency domain.Therefore we have to limit the maximum decomposition level to 5-6 for the used lters.

Best-Basis Analysis of Tremor Signals
Mt. Stromboli is one of the best investigated volcanos.Its seismic signals possess a typical continuous tremor with spectral energy in several frequency bands within 1 < f < 10Hz and superimposed shock events with spectral content from f << 0:5Hz up to the Nyquist frequency.
For our investigations we used velocity data recorded by J. Neuberg during an array measurement with Guralp broadband seismometers on the Stromboli volcano in 1995 Figure 6 shows a spectrogram computed with the BBA and gure 7 shows two spectrograms computed with the STFT of a seismic signal of a duration of 33sec (2048 samples) which contains a shock event of Mt.Stromboli.The partitioning of the frequency axis in the case of the STFT is necessary since a short window is needed in the upper frequency band in order to get a good time resolution and a long window is needed in order to catch the deepest 10sec periodicity.The BBA is inherently computing this kind of adaptation to the signal.Therefore, a better tting partitioning of the frequency axis in adjacent frequency intervals is obtained.When the signal changes its spectral properties in time signi cantly di erent algorithms like Matching Pursuit (Mallat and Zhang, 1993) or Atomic Pursuit (Chen et al., 1996) should be used instead of BBA.They decompose the signal by recursively searching for the best tting basis function instead of choosing a whole sub-band basis system.Consequently they do not x the frequency resolution along the whole time axis.On the other hand one can restrict the window length so that only one shock event is under consideration.
The BBA obviously results in a better decorrelated representation.The coarse time-frequency distribution is equal in both representations but the BBA gives a clearer pattern.Above 5Hz the STFT reveals several spectral lines which appear quite less apparent in case of the BBA.
For studying the dominant sub-bands, individual frequency bands out of the BB spectrogram have been reconstructed in gure 8.The low-pass (LP) trace on the bottom shows a bump shaped event with its maximum just on the very beginning of the shock.The band-pass (BP) signal with center frequency 0.35Hz consists of a precursor to the shock.The correspon-  ding tile can clearly be recognized at the shock onset in gure 6.In the higher frequency bands some short duration wave groups occur temporally during the shock (BP 1Hz) while others exist for the whole duration (BP 5Hz).
The analysis of several shock events revealed that most shocks consist of a deep-frequency bump shaped event during the onset of the shock and several dominant sub-bands above 1Hz which consist of randomly distributed wave groups.The center frequencies of these sub-bands are stable during the occurrence of the shock but vary between di erent shock events.
As the BB basis is better adapted to the signal properties compared with the STFT basis it results in good compression ratios when only dominant coefcients are used for signal storage.Clipping all coecients below a threshold such that 96% of the original signal energy is preserved we achieved the following results: original: 2048 coe cients (100%) reduced: 310 coe cients (b =15.14%)In gure 8 the top three traces contain from top to bottom the original signal, the de-noised signal after thresholding the coe cients and inverse transforming, and the di erence signal.

Conclusion
The time-frequency analysis of broadband tremor by the BBA resulted in a characteristic and more clear representation compared with the standard STFT method.The computation time is similar to the STFT and therefore real-time computation for volcano monitoring and simultaneous compression and de-noising for data storage is possible.An e cient reconstruction algorithm allows the use of this algorithm as an adaptive bandwidth paraunitary lter bank.
One limitation is given by the maximum achievable frequency resolution because of frequency spreading of the wavelet packets.Another limitation is thenite time window for the investigation of the signal where upon strong changes in the spectral content of the signal cannot be well measured by the BBA.Last but not least the BB spectrograms are obtained by time-variant operations because of the downsampling when computing the wavelet packet table.Therefore the sequence of minimaand maximawithin a sub-band might not correspond to the real energy distribution of the signal.

Figure 2 :
Figure 2: Coi et scaling function and wavelet of lter order 30 in time and Fourier domain.(!: frequency)

Figure 3 :
Figure 3: Frequency responses of the magnitude and phase of a Coi et (solid) M = 3 and a Daubechies (dotted) M = 5 lters.

Figure 4 :
Figure 4: Time and frequency domain of a Coi et (30) packet basis function out of a badly frequency localized basis block at level 5. ( : normalized frequency)

Figure 5 :
Figure 5: The cost functional of each basis block is calculated (number).The BBA starts at the bottom of the basis tree a).From bottom to top the cost of the parent basis is compared with the sum of costs of the associated two children bases.b) The basis with the lower cost is chosen.

Figure 6 :
Figure 6: BBA spectrogram of a shock event from Stromboli volcano.Most of the shock events consist of a deep frequent bump occurring during the onset of the signal.

Figure 7 :
Figure 7: STFT spectrograms computed with two di erent window lengths in order to get a good time resolution and to measure long period events.

Figure 8 :
Figure 8: From top to bottom the rst three traces (identically scaled) show the original shock event, the reconstructed signal based on 310 coe cients, and the di erence signal.The following four traces (scaled to maximum) show the reconstructed signal (without compression) based on coe cients which correspond to the center frequency of one band.