ABSTRACT
 
 

The author is interested in measuring and characterising the dynamical behaviour of materials and structures. New techniques for the interpretation of transient and intermittent data, using wavelets, allow system and excitation properties to be deduced from measured data with more precision and greater speed than before.

 The applications that are described in this paper come from the author's field of mechanical engineering. Four examples will be considered, three giving results that have not been published before. They involve the analysis of transient vibration data. The practical objective is to extract as much information as possible from measured results. The data may be short in duration because the phenomenon it represents happens quickly. Or the characteristics of the measured data may change with time because of changes in its underlying physical cause.

 One example is the analysis of vibration data recorded during the transmission of bending waves in a beam subjected to impact loading. This is an essentially intermittent phenomenon as wave reflections occur and energy is transmitted backwards and forwards along the beam. Another similar, but more complicated example uses data for pressure fluctuations recorded in an acoustic waveguide. Here there are many different waves interfering with each other. A third example uses data for ground vibration recorded near an underground train in London. Disturbance from the rumble of underground trains is becoming increasingly intrusive but it is very hard to predict. Finally, the fourth example computes time-varying cross-spectra for multi-channel measurements of soil vibration in a centrifuge test designed to model earthquake response. Simultaneously measured acceleration records at different points allow the changing soil properties that occur under dynamic loading to be explored. The first two examples are laboratory demonstrations which are used as student experiments in the author's department. The second two examples are taken from research in progress at Cambridge for which wavelet analysis now provides an investigative tool of considerable importance.
 
 

HARMONIC WAVELET THEORY

 The theory of harmonic wavelets has been described in previous papers by the author (1993, 1994a, b, 1998, 1999). In their simplest form, orthogonal harmonic wavelets provide a complete set of complex exponential functions whose spectrum is confined to adjacent (non-overlapping) bands of frequency. Their real part is an even function which is identical with the Shannon wavelet. Their imaginary part is a similar but odd function. Their equal spacing along the time axis is twice that of the corresponding set of Shannon wavelets. In their practical application, the boxcar spectrum of harmonic wavelets is smoothed (to improve localisation in the time domain) and the spectra of adjacent wavelet levels are overlapped to give over-sampling in order to improve time-frequency map definition.

 These wavelets have been found to be particularly suitable for vibration and acoustic analysis because their harmonic structure is similar to naturally-occurring signal structures and therefore they correlate well with experimental signals. They can also be computed by a numerically-efficient fft-based algorithm.

 For time-frequency mapping, there are similarities between the harmonic wavelet transform (HWT) and the short-time Fourier transform (STFT). The advantage of the HWT over the STFT is that the HWT is a computationally-efficient variable bandwidth transform. Therefore the time-frequency map it generates can have a variable bandwidth basis, with the analysing wavelet's bandwidth altered from one frequency to another to suit the problem being studied. In contrast, a time-frequency map constructed by the STFT always has a constant bandwidth basis, giving the same frequency resolution at high frequencies as it gives at low frequencies. This means that the STFT is less flexible and may lead to a requirement for (much) more computation than is required by the harmonic wavelet transform. A detailed discussion of the merits of the two related methods is given in Newland (1998).
 
 

NUMERICAL ALGORITHM

 A practical algorithm for time-frequency analysis is illustrated diagrammatically in fig. 1 (Newland 1998). The correlation calculation that is at the heart of the wavelet method is carried out in the frequency domain, where it becomes a multiplication operation rather than a convolution. The input signal f(t) is represented by the N-term series f₀, f₁, f₂, ... f_N-1 in the top box. It is transformed from time to frequency by the FFT to give the Fourier coefficients F₀, F₁, F₂, ... , F_N-1 in the second box. The record length is assumed to be unit time, so that the sampling interval is 1/N and the Nyquist frequency Np. By definition the

 Figure 1. FFT ALGORITHM TO COMPUTE HARMONIC WAVELET COEFFICIENTS FOR WAVELETS IN THE FREQUENCY BAND m2p£ w< n2p (FOR A RECORD OF UNIT LENGTH).

(1)

where m< n (Newland 1994b). Because the wavelets are complex, their Fourier transform is one-sided, so that W_k remains zero for all negative frequencies.

 On carrying out the multiplication operation, the F_k and W_k terms are multiplied to generate the new series
 
 
 

,

(2)

 in the next box in fig. 1. These A_k are the Fourier transforms of the wavelet coefficients. Computing their IFFT reverts to the time domain to give the series of wavelet coefficients a_r, r = 0 to N-1, which are shown in the bottom box.

 The time scale runs from andgives the result of calculating the wavelet coefficient for the wavelet centered at the chosen position t=r/N on the time axis. The usual circularity property of the discrete Fourier transform method applies and when a wavelet runs off the end of the unit time scale, it wraps round and reappears at the opposite end.

 The computation in fig.1 therefore gives N wavelet coefficients for reference wavelets in all the possible N positions along the time axis. Only n-m of these are needed to form an orthogonal set (Newland, 1994b) and usually less than the (large) number N is needed to produce adequate resolution along the time axis. This is achieved by selecting N1 equally-spaced values from the total available. If N1 is not a factor of N, appropriate methods of interpolation can be used. An efficient method of doing this is very important. The method used here is described in Newland (1999).

 Instead of computing the IFFT of the N-term series A_k, k = 0 to N-1, in the lowest but one box in fig. 1, this interpolation method computes the IFFT of a shorter N1-term series B_k, k = 0 to N1-1, whose first n-m terms are the non-zero A_k, k = m to n-1, and whose remaining terms are all zeros. It is shown in the reference that this generates a set of coefficients b_s, s = 0 to N1-1, which correspond to selected terms in the longer series a_r, r = 0 to N-1, provided that sN/N1 is an integer. If it is not, then the b_s interpolate between the nearest two values of a_r. The magnitudes of corresponding terms are the same. Therefore a time-frequency amplitude map drawn by computing the shortened N1-term series b_s defined above will faithfully represent an amplitude map computed from the full-length N-term series a_r. The phase angles of corresponding terms will generally be different according to (Newland 1999)
 
 
 

,

(3)

 but allowance can be made for these differences.

 The centre frequency of the wavelet Fourier transform in fig. 1 is (m+n)p and its bandwidth is . By changing the centre frequency, or bandwidth, or both, and repeating the calculation, a new series of wavelet coefficients, a_j, j = 0 to N-1, is generated. If this process is carried out for N4 different centre frequencies and each output series a_j downsampled to give N1 terms, the resulting N4´N1 array A(N4,N1) is generated. This array is used to draw time-frequency maps to show how the amplitude and phase of the wavelet coefficients change over time and frequency. In the author's programs using these principles, the parameters N1 and N4 have the above meaning. Wavelet bandwidth is allowed to change linearly from n-m=N2 to n-m=N3 over the full frequency range of the calculation.

 The algorithm in fig. 1 applies for all harmonic wavelets, namely wavelets defined in the frequency domain with a compact spectrum such that W(w ) = 0 outside a defined (generally narrow) band of frequencies. This no longer defines an orthogonal family of wavelets but since reconstruction of the signal being analysed is not required, that does not matter. For the results given below, the boxcar spectrum of orthogonal harmonic wavelets has been windowed by a Hanning function, so that the function in the third box in fig. 1 is given by
 
 
 

.

(4)

 This has been found to gives good localization in the time domain.
 
 

PHASE INTERPRETATION

 This calculation procedure generates complex wavelet coefficients, a_r. Their phase depends on the relative position of the signal and its analysing wavelet. This defines the ratio of the imaginary part of a_r (correlation with the odd part of the harmonic wavelet) to that of the real part of a_r (correlation with the even harmonic wavelet). When, for a constant harmonic signal, the wavelet is moved to a new position, its phase will be different. Therefore absolute phase is not a useful indicator because it depends on wavelet location. But phase gradient, defined as the rate of change of phase with time for wavelets in the same frequency band, is an interesting parameter because it is constant when f(t) is a harmonic of fixed frequency and phase. It is shown in Newland (1999) that, for a single harmonic of frequency w ₀, the phase gradient
 
 
 

(5)

 where W is the centre frequency and 2p B the bandwidth of the analysing wavelets.

 The essential property is that the rate of change of phase with time is constant for a harmonic of fixed frequency and phase so that a two-dimensional map of phase gradient, with on a frequency-time base, is sensitive to phase changes in the signal being analysed. This will be illustrated in one of the examples of wave propagation given below for which sudden changes in phase gradient occur in between successive reflections of energy in local frequency bands.

 A different example in which absolute phase can be used helpfully is the last example below, which is the analysis of simultaneous multi-channel recordings of ground movement after shock loading. Corresponding time records are analysed by the same wavelet arrays, with the wavelets in the same time location for the different channels. This enables differences in the phase of the channels to be detected and mapped as a function of time and position. It will be shown how changes in system properties (caused for example by soil slippage) can be detected by harmonic wavelet analysis from corresponding changes in the measured phase response.

BENDING WAVE TRANSMISSION IN A BEAM

 The first example comes from a laboratory experiment in the Cambridge University Engineering Department. The experiment illustrates bending wave propagation in a thin steel beam. The beam is suspended on light cords with its long axis vertical and is hit gently at one end by a soft-ended impulse hammer. This generates lateral bending waves which travel to the other (far) end of the beam, where they are reflected and return to the point of impact, before undergoing successive reflections until eventually they are dissipated by damping after several seconds. A small accelerometer is mounted on the beam so as to detect lateral vibration. For the results shown below the accelerometer is positioned close to the first end of the beam, near its point of impact. The beam is 7.2 m long and has a rectangular cross-section 32.1×6.3 mm. The impulse hammer had a soft tip designed so that only low frequency vibrations were generated (up to about 1 kHz). The sampling frequency was 4096 Hz.

 Because the group velocity of bending waves depends on frequency (velocity proportional to frequency^1/2), groups of high frequency waves travel faster than low frequency waves. Therefore a time-frequency map should show more frequent reflections for high frequencies than for low frequencies. This behaviour is, of course, not at all evident from the recorded time-domain response, which is shown for one second duration in the top inset view in fig. 2, or from the spectrum which is drawn at the side of fig. 2 with the frequency scale running from 0 to 750 Hz approximately. A short length of this data has been included in Newland (1998) but the full record has not been considered previously.
 
 

 As time passes, the regular pattern of curved ridges in fig. 2 is interrupted by some transverse "valleys" that run from left to right in the figure. These appear to be caused by non-bending modes into which vibrational energy "leaks" as the wave propagation process continues. Within the frequency range of fig. 2, there are about 50 bending modes which are excited. There are also about 9 twisting modes and two longitudinal modes whose frequencies lie in range. Some of these may be unintentionally excited by the impulsive input being slightly away from the geometical centre of the beam or slightly off-line in direction, or they may be coupled to the bending modes, for example by the action of the supporting elastic cords or the mass loading of the accelerometer.

 When a harmonic wavelet with a narrower bandwidth is used for the analysis, reduced definition along the time axis is achieved; for a harmonic wavelet with a wider bandwidth, reduced definition along the frequency axis is obtained (specific
 
 

 The unavoidable smearing of spectral features that occurs in fig. 2 can be reduced by plotting only the ridges of the three-dimensional surface whose contours generate the figure. The exact identification of ridges is difficult (Eberly 1996) and identifying their precise position is complicated. The approach used by the author (Newland 1999) is to seek the height maxima of sections cut in the direction of the (smoothed) surface's greatest curvature. When this strategy is applied to the surface plotted in fig. 2, the result is shown in fig. 3. Each ridge marks the arrival at the measurement point of successive groups of bending waves. At high frequencies the group velocity is higher, so successive reflections arrive more quickly than at low frequencies, when the ridges are further apart. Knowing the length of the beam, by measuring the time between successive reflections, the group velocity can be estimated as a function of frequency.

RESPONSE OF AN ACOUSTIC WAVEGUIDE

A similar, but more complicated, example is provided by the reflection of pressure waves within an acoustic waveguide. This is also a laboratory experiment at Cambridge. Internal air pressure perturbations are generated in a closed circular duct of approximately 12 m length and 0.75 m diameter. These perturbations are caused by a pulse-like electrical input to a small loudspeaker mounted near the edge of one of the rigid ends of the duct. This excites several different families of acoustic waves which travel backwards and forwards within the duct. A microphone mounted at the centre of the end with the loudspeaker

 Fig. 4 shows the ridges of an amplitude time-frequency map, computed as described before. In addition to the main ridges, there are numerous small, generally horizontal ridges which arise from local fluctuations in surface height. They can be eliminated by introducing more smoothing before ridge detection and it is a matter of judgement to generate the ridge map which is the "best" for a required purpose. As for fig. 3, fig. 4 has the input time history shown for comparison along the top, and the modulus of this signal's Fourier transform plotted along the left-hand side (using arbitrary units). For convenient scaling, the square root of the Fourier transform is plotted.

 The underlying physical processes represented in this map are quite complicated. Axial, plane waves travel at constant velocity (independent of frequency) and bursts of energy from these plane waves arrive periodically at the microphone. They show as equally-spaced vertical lines in fig. 4. Knowing the dimensions of the duct and the acoustic properties of air, the position of these lines can be calculated and their theoretical position is superimposed on their experimental position in fig. 4. Within the frequency range of these maps there are two other families of non-axial waves which are dispersive (their group velocity depends on frequency). Their passage time between reflections is given by where L =12.16 m is the length of the duct, c=334 m/s is the speed of sound, w is the wave frequency and W₀ is the cut-off frequency. For plane waves, the cut-off frequency is zero, and for the first two families of non-axial waves which are detected by a microphone at the center of the duct it is W₀=3.83c/a and 7.02c/a where a=.386 m is the duct's radius (see, for example, Skudrzyk 1971, p. 431). These results have been used to plot the theoretical lines on fig. 4, measuring time forward in steps of T(w ) from the instant of impulsive excitation.

 The calculated cut-off frequencies are plotted as the horizontal lines in fig. 4. Successive reflections of the dispersive waves appear as the two families of curved lines that are asymptotic to the cut-off lines. The horizontal ridge at about 200 Hz is due to ringing of the loudspeaker's diaphragm and is not associated with the travelling wave acoustic phenomena.

 Fig. 5 shows a differential phase map of the same data. It is a contour plot of the three-dimensional surface obtained by plotting the modulus of phase gradient on a base of frequency versus time. It can be seen that phase perturbations occur generally between ridge positions indicating phase changes at every reflection of the travelling wave energy. A characteristic of this presentation is that the vertical distribution (distribution over frequency) of the phase perturbations is correlated with the position of the peaks in the Fourier transform of the input signal (plotted along the left-hand edge). By placing a straight edge across a similar diagram drawn to larger scale, it can be shown that the phase perturbations align quite closely with the positions of the troughs in the spectral data plotted on the left-hand side (Newland 1999). This is not evident so clearly in corresponding graphs of wavelet amplitude.
 
 

frequency / Hz (vertically), time / seconds (horizontally)

 Recently ground acceleration data has been measured near to a curved section of the Piccadilly underground line in London. For the results given here, an accelerometer was secured to a stone step in an adjacent building and 20 seconds of the passage of a train recorded. The time history of this process is shown in the upper view in fig. 6, the units being g's. The recording begins with the train already passing, and continues until it has passed out of hearing in 20 seconds.

 Fig. 6 shows a harmonic wavelet amplitude map for this vibration, covering the frequency range from zero to 250 Hz (half the Nyquist frequency of 500 Hz). It is evident that there is a broad band response as the train passes, with ground vibrational energy in a wide range of frequencies as a result of wheel and rail surface irregularities, wheel flange-to-rail contact, mechanical train noise and electrical collector noise.
 
 

GEOTECHNICAL CENTRIFUGE TESTING

 Fundamental knowledge of the (large amplitude) dynamic behaviour of soil under earthquake excitation is meagre. Studies at Cambridge in our Geotechnical Centrifuge Centre and elsewhere (Lee & Schofield 1988, Taylor (ed.) 1995, Butler 1999 and others) have obtained good data on the transient vibration of soil models under earthquake conditions. The levels of excitation cause large deflection inter-granular movements which lead to so-called soil liquefaction effects when the soil's response is closer to that of a fluid than a solid. Because excitation lasts only for a second or two with excitation frequencies ranging up to about 200 Hz, data analysis can only be done if there are good methods of transient vibration analysis. Wavelet methods make this possible and good preliminary results have already been achieved using harmonic wavelets (Newland & Butler, 1998). New research is concentrating on developing these methods to estimate time-varying cross-spectra between adjacent measuring points. This is seen as a very important area of further theoretical and experimental research.

EXPERIMENTAL DATA

 The test data used below is that published in Newland & Butler (1998) and was obtained from geotechnical centrifuge tests. The experimental system represented a saturated soil model poured at two relative densities and mounted within a flexible container. The container is shown in fig. 8. It is rectangular in shape with its side walls made of a series of flat rings, each mounted to the next by a rubber gasket. The intention is that the loaded container functions as an equivalent shear beam whose shear modulus matches approximately that of the enclosed soil medium. Sand was poured at a density 1576 kg/m³ in the lower 160mm of the container and at 1670 kg/m³ for the remaining 365mm to the top of the container.

 The container and its contents were centrifuged to apply an acceleration vertically downwards (in fig. 8) of 50g in order to simulate the response of a large ground volume in a small model. Horizontal shear excitation to the base of the container was supplied by a device called a stored angular momentum actuator. This consists of a flywheel which is connected to a reciprocating rack by a clutch assembly. When the clutch is engaged, there is a sudden burst of oscillatory energy which shakes the container and its contents while this is being centrifuged.

 Sand movement was detected by miniature piezo-electric accelerometers. As seen in fig.8, these were stacked vertically within the test specimen with four in the bottom layer of sand and ten in the top layer. Previous tests have shown that they have an accuracy of 5% within the frequency range of 20Hz to 2kHz. Their natural frequency when embedded in sand is estimated to be about 4 kHz compared with frequencies of interest up to about 400 Hz. Each transducer was carefully oriented in the sand to record the resulting horizontal motion (the x-direction in Figure 1) within the saturated model.

 Data was stored in a digital data acquisition system developed in Cambridge as part of the centrifuge's instrumentation, from which it is retrieved for detailed computer analysis. The input motion to the base of the model container had a fundamental frequency of 27Hz with a displacement amplitude of 1.5mm. The duration of the shaking excitation was set to 1.2 seconds.

 For purely harmonic movement, these displacements correspond to a lateral acceleration amplitude of about 4.4g. However loose-play and nonlinearities in the mechanism introduce a harmonic content to the excitation, as will be apparent from the measured results below. The measurement points are shown in fig. 8. The signals f₁ to f₆ were recorded at the following six locations: f₁ at 7726, f₂ at 7828, f₃ at 7319, f₄ at 7709, f₅ at 12611 and f₆ at 12612, but only two of these measurements f₁ and f₄ are used for the results given below.

 All signal processing computations have been done by the harmonic wavelet method using the algorithm described above. Results are shown as before as two-dimensional maps of three-dimensional surfaces plotted for the relevant parameter. For ease of identification for a multi-channel system, it is convenient to refer to (i) the power spectral density of a measurement, (ii) the amplitude of the cross-spectral density between two measurements, and (iii) the phase of the cross-spectral density. These terms are not strictly correct because they are defined for stationary random processes, whereas we are concerned with transient and non-stationary processes. However the amplitude squared of the wavelet coefficient is called a power spectral density since, for an orthogonal set of harmonic wavelets, the mean-square signal is equal to the sum of the (weighted) wavelet amplitudes squared (see e.g. Newland 1993). When the signal is oversampled to generate extra wavelet coefficients, the same analogy may be used. Similarly, the product of two wavelet amplitudes, when computed for the same wavelet at the same instant of time for two signals, represents the amplitude of the cross-spectral density between these signals (for that time instant and frequency band). Also the phase difference between the same two wavelet coefficients gives the phase of the cross-spectral density between these signals.

 Time-varying auto-spectral densities calculated this way are plotted for records f₁ and f₄ (positions 7726 and 7709 in fig. 8) in figs. 9 and 10.
 
 

 Figure 9. EXPERIMENTAL DATA: POWER SPECTRAL DENSITY FOR SIGNAL f₁ MEASURED AT POSITION 7726 (see fig. 8).

The amplitude map in fig. 9 shows that the vibration close to the bottom of the box remains approximately constant as shaking continues because the horizontal stripes have approximately constant width and continue for the full duration of the shaking process. In contrast, all the other power spectral densities, for example f₄ in fig. 10, show obvious changes with time. This must be due to the changing physical properties of the soil model as a result of its changing dynamic properties with time.

CROSS-SPECTRAL DENSITIES

 Power spectral density data indicates the total energy in a signal and its distribution over frequency and time. Relative changes in two signals are described by the cross-spectral density. This provides a measure of the local correlation between two signals. The cross-spectral density's amplitude is measured by the product of the amplitudes of both signals; its phase depends on the relative phase of the two signals. As for all time-frequency analysis, all spectral calculations are estimates for the chosen frequency band and time window considered. For the harmonic wavelet method, this is defined by the bandwidth and time duration of the chosen wavelet, which is under the control of
 
 

 For the experimental data, fig. 11 shows the phase of the time-varying cross-spectral density for f₄ with f₁ These results have not been published before, and some explanation of their interpretation is needed. The density legend is shown on the right side of fig. 11. It runs from p = 180° at one extreme to -p = -180° at the other extreme. Before the forced motion has been applied, and after it has finished, there are residual noise signals and these give rise to the haphazard phase represented by "marbling" on the right and left-hand sides of the map. During the forced motion, the map has obvious horizontal stripes, each corresponding to one harmonic of the motion (identified on the left side spectrum, which is for f₁). Around the fundamental frequency of 27 Hz, the density of the stripe indicates that f₄ and f₁ are approximately in phase with each other. The same is true for the second harmonic initially, but as motion continues there is a phase change from approximately zero through minus 90° and then through -180° to approach zero again. For higher harmonics, the density of each stripe changes as time passes, with the transition sometimes being quite sudden.

 In between each horizontal stripe there is a thinner marbled stripe where there is uncorrelated response between two harmonics and phase coherence is not maintained.
 
 

SIMULATED RESULTS