COMPUTER PROCESSING OF NEW Z E A L A ND STRONG-MOTION ACCELEROGRAMS

SYNOPSIS: Recently implemented modifications to the computer processing of earthquake accelerograms are described. The main extensions to the system include dynamic instrument correction, robust band-pass filtering and accurate integration, all implemented in the frequency domain. The processed records now provide accurate representations of strong ground motion within the frequency band from 0.25 Hz to 24.5 Hz: full accuracy is maintained in the computed velocity and displacement signals. The use of the new processing system is demonstrated using two previously recorded accelerograms.


INTRODUCTION;
The Physics and Engineering Laboratory of DSIR has maintained and operated the New Zealand strong motion earthquake recording network for the last 16 years: there are currently 160 time-base accelerographs in operation providing accelerograms on 35 mm photographic film.
An important aspect of the operation of the complete network is the processing of the analogue film records into a digital representation which is readily useable by the engineering and scientific community.
Firstly, the accelerogram traces of records of sufficient intensity are digitised using a semi-automatic computer-based digitiser system: this process is described in detail elsewhere (1).
In this paper, the procedures for routinely processing these raw digital film trace-measurements are described.The processing performed is essentially of three types: (a) to correct the data for instrumental response and known imperfections; (b) to eliminate or reduce in effect the errors introduced in the digitisation process; and (c) to calculate other quantities of interest from the corrected acceleration data, for example, displacement records or engineering response spectra.
Work on this project began in 1978 (2), and had reached a state in 1981 that enabled the publication of corrected accelerograms and response spectra for twenty of the largest recorded New Zealand accelerograms (3).However, the experience gained from the preparation of this report, and the appearance of the data itself, suggested that further changes to the system of processing the records were required.
In particular, a new digitisation system was purchased to enable more rapid and accurate conversion of the analogue data to computer readable form and an investigation was begun into the use of more sophisticated digital signal processing techniques which would maintain signal integrity at all stages of processing.
The major enhancements recently applied to the computer processing system for New Zealand accelerograms can be summarised as follows: (a) correction for dynamic frequency *Physics and Engineering Laboratory, DSIR, Private Bag, Lower Hutt.

BULLETIN OF THE NEW ZEALAND NATIONAL SOCIETY FOR
response of the recorder sensors (in addition to the previouslyapplied static corrections); (b) higher initial sampling rates, followed by digital filtering to reduce potential aliasing effects; (c) removal of low-frequency components of the signals, which are potentially contaminated by digitisation and other noise; (d) use of frequency domain filtering and integration to ensure that signal integrity is maintained through all stages of processing; (e) appropriate handling of end effects to obtain zero initial and final velocities and displacements.
The signal processing system now used ensures that the computed ground velocities and displacements contain all the information present in the original accelerogram: previously this data had been calculated in a manner which suppressed the higher frequencies.This means that the signals calculated may be used in many diverse applications: for example, the generated displacement signals may be used directly as input to a displacement-controlled shaking table and need not be modified prior to use.
The digital fi incorporated in the system aim to provi extremely accurate band between 0.25 H components outside are attenuated in s from the processed ltering techniques routine processing de signals which are in the frequency z and 24.5 Hz; data this frequency range trength or eliminated signals.

PROCESSING TECHNIQUES: (a) Instrument correction
The strong motion accelerograph used exclusively throughout the New Zealand strong-motion network serviced by PEL is the MO-series instrument (4).This instrument records on 35 mm photographic film the response of three nominally orthogonal accelerometer pendulums to earthquake excitation.
The raw trace readings obtained from the digitiser are initially converted to units of acceleration (mm/sec 2 ) using the appropriate static sensitivity, with corrections being made for sensing-axis direction errors and cross-axis effects: this processing is performed using data obtained from a static calibration of the individual accelerometer pendulums (5).
The crossaxis effects are of greatest importance during the strongest phases of ground motion, but the correction is applied without regard to signal amplitude.
The coordinate system used for later presentation of results may at this stage be rotated about the vertical axis to any desired direction.
Since the nominal sensing directions of the accelerometer axes are at 45° to the accelerometer case, a default counter-clockwise rotation of 45° is applied to align the trace directions with the longitudinal and transverse axes of the case.
Usually the accelerometer case itself is aligned along dominant structural axes, however, recently some accelerometers (for example in Bowen Street Bridge, Wellington) have been installed so that the nominal sensing directions are aligned with major structural features, to permit simpler interpretations to be made from a visual inspection of the original film record.
Here an axis rotation would generally be inappropriate.
In addition to the static corrections outlined above, a correction is made for the dynamic frequency response of the accelerometer pendulums.
The transducers in the MO-series accelerograph are stiff and heavily damped, giving a raw trace reading which is essentially equal to the instantaneous acceleration in the opposite direction.
The equation of motion relating the true acceleration a(nAt) to the trace reading y(nAt) at time t = nAt is: a (nAt) = " % (y (nAt) + ^iy(nAt) + n -K y( nAt )K (i) n where o> n is the natural undamped circular frequency of oscillation of the transducer, and g is the critical damping ratio.Here the term outside the brackets in (1), 2 -uj n , represents the static sensitivity already applied; the dynamic correction enclosed in the brackets is imposed through the use of Fourier transforms.
Denoting the discrete Fourier transform of the appropriately rotated and statically corrected trace displacement record y(nAt) as Y, we define (6) the forward transform and inverse transform where At is the sampling interval (sec), N is the total number of samples in the signal and Af = 1/(NAt)is the frequency resolution (Hz) of the transform.The numerical considerations involved with the use of the Fourier transform are discussed later.
Substituting (3) into (1) and performing the differentiations we obtain the transformed acceleration A(kAf) as 2 A(kAf) = -a> H A (kAf) Y(kAf) , k = 0 ... N-l (4) where the dynamic instrument correction is performed by a complex multiplication of the transformed trace reading with the transfer function of instrument response, H A(kAf). Here, where f = GO /2TT is the natural frequency n n " of the transducer, measured in hertz.- The nominal dynamic properties of the MO-series accelerometer block are currently used for the dynamic instrument correction; these are a natural frequency f n of 30 Hz and a critical damping ratio £ n of 0.60.
As shown in Figure 1, an oscillator with these characteristics has a transfer function whose amplitude is very close to unity throughout the frequency band of interest (0-25 Hz), while the phase of the higher frequency components is significantly affected.This phase shift can produce significant changes to the peak acceleration applicable to a particular trace.

(b)
Band-pass filter The digital readings obtained from the analogue 35 mm film traces are subject to at least two potential forms of error: experience has shown that these errors are most severe at the extreme ends of the frequency spectrum (6).
Firstly, there is an unknown amount of aliasing caused by the act of discretely sampling the continuous waveform.
The highest frequency present in the sampled record (the Nyquist frequency) is equal to one-half of the sampling frequency.Any higher frequency content present in the original record is "folded" about • the Nyquist frequency, appearing in the sampled version as a lower frequency component.
This effect is illustrated in Figure 2 for a 40 Hz sinusoid sampled at 50 Hz, which appears as a 10 Hz component in the sampled record.
The effects of aliasing are minimised in the early stages of digitisation by sampling the traces at extremely high rates (>130 Hz), and linearly interpolating to give equally spaced samples at a sampling rate of 100 Hz.
An overlay graph is prepared, and from a visual examination of this it is clear that aliasing can only be a problem at frequencies close to the Nyquist frequency at 50 Hz, where signal strengths are expected to be extremely small.The potential for aliasing within the frequency band of interest (0-25Hz) is eliminated by low-pass filtering prior to decimation of the original trace readings to give a standard sampling frequency of 50 Hz (20 m-sec sampling interval) with Nyquist The second source of error is noise generated during the digitisation process; this digitisation noise contains components at all frequencies in the band of interest.However, at the lower frequency end the true signal is generally small in amplitude, leading to a relatively low signal-to-noise ratio.
The consequences of the lowfrequency components introduced into the signal are not evident until the velocity and displacement signals are computed: here the low-frequency noise components are magnified (components of frequency w in the accelerations are enhanced by factors of 1 and 1 in the velocities and displacements respectively) and can completely dominate the displacement record, as in (3).The low-frequency noise components present in the digitised traces, along with any valid signal components in this range, are removed through high-pass filtering.
The magnitudes of the velocities and displacements which are generated from acceleration components having an amplitude equal to the approximate resolution of the digitisation system (0.001 g) are plotted in Figure 3 against the frequency of the acceleration signal.
Here the large displacement errors resulting from the limited resolution of the system are apparent, for instance, a sinusoidal acceleration signal at 0.001 g amplitude and 0.1 Hz produces a displacement amplitude of 25 mm.
As the basic resolution of the digitisation system is a lower bound on the magnitudes of the errors, it is possible that the errors in the displacement record associated with low-frequency errors in the digitisation may be somewhat larger than this example indicates.
Ideally, the band-pass filtering resulting from the above considerations would allow all frequency components between a given low-frequency limit f and a given high-frequency limit f to be passed without alteration, and for frequency components outside this band to be completely eliminated.
When dealing with finite duration data, this approach is impractical since the impulse response of the filter (the inverse Fourier transform of its transfer function) is infinite in length and decays very slowly.Some compromise in filter specification must therefore be made: usually a transition band is provided near the frequency cutoff points to allow the filter transfer function to change smoothly from unit in the pass band to zero in the stop band.
The approach taken in the processing of United States records at the California Institute of Technology ( 7) was to provide a transition band in which the magnitude of the transfer function varied linearly (as shown in Figure 4): this is called an Ormsby filter.
However, the filter transition bandwidth selected at the lowfrequency end is very narrow, resulting again in a filter with an extremely long impulse response.
This effect is particularly noticeable in the calculation of displacements: the impulse response for the standard Caltech Ormsby high-pass filter contained with double integration is shown in Figure 5.
The "ringing" of the response is indicative of the narrowness and the form of the transition band.The original implementation of this system used multiple implementations of a filter with truncated impulse response operating on a decimated (and strongly aliased) version of the accelerogram, and had the potential for introducing significant distortions into the processed records.
In addition, the filter transients at each end of the record, legitimately arising through the use of the non-causal Ormsby filter, were ignored.This results in velocity and displacement signals with non-zero initial and final values, and may be responsible for long-period drifts present in the records.
A recently proposed modification to the processing methods used for United States accelerograms (8) includes a recommendation for the use of a zero-phase infinite impulse response (IIR) elliptical filter to achieve the band-pass filtering required.
However, this approach requires a lengthy filter design process to be used for each different set of filter specifications , and is generally a complicated method to implement in a routine processing system.
The procedure that is adopted here is a frequency sampling approach to filter design where the filter transfer function is specified directly in the frequency domain; a sampled form of this is used to window the discrete Fourier transform of the input sequence.
For the band-pass filter, a transition band is selected such that the transfer function varies between zero and unity in such a way that the function and its first derivative are continuous (c.f. the Ormsby filter in which only the function itself was continuous).
This gives a substantially more rapid decay in the impulse response of the filter, as shown in Figure 6 for a filter using the same transition frequencies as the US standard Ormsby • filter, but having a transfer function of the form T where f L and f L are the roll off initiation and termination frequencies of the high-pass filter, fJ and fT are the corresponding lowpass filter frequencies and f >T is the N original Nyquist frequency.The actual transition roll off initiation and termination frequencies selected are of extreme importance in determining the quality of the filter used.

STANDARD ORMSBY FILTER IMPULSE RESPONSE
For the routine accelerogram processing conducted at PEL, the default frequencies used are as follows: (i) 'Low-pass filter.
Consistent with the desire to produce accelerograms with a 20 msec sampling interval, we select a Nyquist frequency of 25 Hz.
The transition band is of width 1 Hz centred on this frequency, giving f* = 24.5 Hz and f? = 25.5.Hz.
Here we do not have the luxury of specifying a transition bandwidth of 1 Hz as a significant amount of valuable information would be lost along with the removal of the noise.Following the general approach adopted by Basili and Brady (9) we choose the cutoff frequency, fT as the greater of 0.10 Hz and 3/T T , where T is the duration of L L the digitised accelerogram in seconds.The default rolloff initiation frequency is selected to provide a transition bandwidth of 0.15 Hz in the sampled transfer function of the filter.
Thus for an accelerogram of digitised duration 30 sec, we obtain fT = 0.10 Hz and fl = 0.25 Hz.
The band-pass filter is applied to a signal with a sampling rate of 100 Hz.The filtered signal is then decimated in the frequency domain to obtain the final sampling rate of 50 Hz: the transform of the data is halved in length and amplitude to account for the decimation.The combination of band-pass filtering selected above, along with this decimation, provides for accelerograms of normal duration, a signal which is unaliased with a frequency content which is unattenuated between limits of 0.25 Hz and 24.5 Hz.

(c) Integration
The integration algorithm used in the standard processing of United States accelerograms (7), and that used previously for New Zealand accelerograms (3), is the trapezoidal rule.
This integration method may be shown to suppress the higher frequency components (10) , with errors larger than five percent at frequencies greater than one-quarter of the Nyquist frequency, and complete elimination of the Nyquist frequency component.
Although most earthquake displacement signals are dominated by the components at low frequency, and the trapezoidal rule reproduces the overall pattern reasonably well, a more accurate digital integrator enables the displacement signal to retain the same level of information as the original accelerogram.
The recently proposed modification to the United States routine processing system (8) recommended the use of the Schuessler-Ibler digital integration filter (11) for calculation of earthquake velocities and displacements.This system, based on Newton*s 3/8 rule, is only valid for signals which are bandlimited at one-half of the Nyquist frequency, and may introduce errors leading to incorrect high frequency components of the integrated signal.
In common with the approaches adopted in the previous sections, the integration of the accelerogram traces to produce velocity and displacement signals is conducted here directly in the frequency domain using frequency sampling techniques.Thus the transfer function for single integration is: and for double integration:

IMPLEMENTATION
The product of the instrument correction, band-pass filter and integration transfer functions (from ( 5)-( 8)) operates directly on the discrete" Fourier transform of each component of the accelerogram.
Since the Fourier transform assumes that the original function is periodic (with period equal to the original sample length), precautions are necessary to ensure that linear (rather than circular) filtering is performed, with no interaction between the beginning and end of the record occurring (11).
Thus, prior to the taking of the Fourier transform, the statically-corrected zero-mean accelerogram is extended with zeroes so that the sample length is an integer power of two, and this sequence is doubled in length with additional zeroes.
The discrete Fourier transform algorithm used is a radix 4-2 FFT implementation (12).
This method is ideal for use in a routine scheme, being extremely fast and relatively economical in storage.
As an example, an 8192 point DFT as required in processing a record of 20-40 seconds duration, takes less than 2 seconds processor time on the DSIR VAX 11/780 using subroutine FAST in (12).
It is imperative that the impulse response of the combined filter (which is convolved with the extended accelerogram record) be zero for over one-half of its length.
Under these conditions, the filter transients generated at each end of the original accelerogram will not interact and will be found in the passed portion of the record.
The transient at the commencement of the record will, as a result of the periodicity implied by the use of the DFT, be found at the end of the extended record.
The combined filter response functions, which are calculated as the inverse DFT of the combined transfer function, are shown in Figure 7 for the cases of zero, one and two integrations.
Here it may be seen that the non-zero impulse response lengths are longest for the displacement calculations with two integrations, however, with the parameters used this length is substantially less than one-half of the total record length.
Although this filter would be quite expensive to implement in the time domain, the implementation in the frequency domain is relatively simple and the transient merely provides a transition, consistent with the filtering employed, between the assumed quiescence before and after the earthquake and the portion of the earthquake recorded and/or digitized.
In this sense, the filter transient represents the lack of knowledge about conditions existing before triggering and after the recorder stopped recording.
The complete processing package is very similar in concept to that adopted for the processing of ground motion signals in Japan (13).
Although the instrument correction uses different parameters, and the filter transfer function shape is different, the use of frequency sampling filtering and integration is the same as is adopted here.
The Japanese implementation in (13) does not include the filter transients (resulting from the non-causal band-pass filter) in the published records; the approach adopted here ensures that these transients are included.

1.
Read three components of raw accelerogram trace readings given at 10 msec intervals.

2.
Remove trace mean from each component.

3.
Convert to units of acceleration using static sensitivity values.

4.
Perform sensing direction and crossaxis corrections; rotate coordinates 'about vertical axis if required; compute peak statically corrected accelerations.
Perform steps 5-10 for each component:

5.
Remove trace mean; extend with zeroes to next higher integral power of two; double in length with additional zeroes.

6.
Take DFT of extended sequence.save windowed DFT on disk.

8.
Take inverse transform to give dynamically corrected and filtered accelerations; save with filter transients present in extended portion of record.

9.
Restore DFT; multiply by transfer function of single integration (7), take inverse transform to give velocity signal; save with transients.

10.
Restore DFT; multiply by transfer function of double integration (8); take inverse transform to give displacement signal; save with transients.

VERIFICATION OF THE PROCESSING SYSTEM
A test of the high-pass filter and double integration algorithms adopted for the processing of New Zealand accelerograms was made for the El Centro 1940 May 18, Component S00E accelerogram (14).For this test a comparison is made with the displacement signal (presented in Figure 5  does have an effect on trace displacements, while for the processing system adopted here the complete transient is calculated in all cases, and the displacement signal is independent of the duration of the transient retained.

EXAMPLE
The new processing programs have been applied to a number of recently recorded accelerograms and have been found to be convenient to use and economical, and retain sufficient flexibility that parameters may be readily altered to suit particular requirements.
An This low frequency component dominates the displacement signal (making interpretation of the record very difficult) and its retention cannot be justified with regard to the resolution of the digitisation system employed.
Using the new digitisation and computer processing systems, these low frequency errors are eliminated, giving a displacement trace which is more amenable to interpretation and which provides a more accurate representation of ground motion over the complete frequency band of interest.
As shown in Figure 10(b), the peak acceleration has changed insignificantly, the peak velocity has been increased by 6% while the peak displacement has been reduced by 76%.
In this case, the use of accurate digital signal processing techniques has produced a significant reduction in the peak displacements calculated from the accelerogram.

CONCLUDING REMARKS
The newly implemented computer processing system for New Zealand strongmotion accelerograms has the objective of producing digital records of strong ground shaking which are accurate over the widefrequency band of interest to engineers and scientists.
New Zealand strong-motion records are now processed using a scheme with the following features: (a) Full correction for dynamic frequency response of the recorder, in addition to the static sensitivity, sensing axis direction and crossaxis corrections applied previously.

(b)
A robust band-pass filtering system is employed to minimise the errors, at the extreme ends of the frequency band, which are introduced during the digitisation process.The implementation provides an extremely flexible and convenient method of tailoring the filtering frequencies to particular earthquake records or site conditions.

(c)
The calculation of velocities and displacements from the recorded 1 accelerograms uses an integration scheme which retains full accuracy in the calculated quantities across the complete frequency band.

(d)
The filter transients introduced at each end of the records due to the use of the non-causal band-pass filter are retained in the published data.These transients provide a consistent transition between the record and the assumed quiescence before and after the digitised motion, and ensure that the initial and final values of velocity and displacement are close to zero.
The normally applied filtering frequencies provide an accurate representation of the earthquake induced motion throughout the frequency band from 0.25 Hz to 24.5 Hz.All filtering and integration is conducted in the frequency domain, however care is taken to ensure that the combined impulse response is of limited duration.

Figure
Figure 5: Impulse response of Ormsby filter combined with double integration fj" = 0.07 Hz, = 0.05 Hz ij L , ,SINUSOIDAL TRANSITION Figure 7 Double integration Impulse response functions for standard sinusoidal transition filter fl = 0.25 Hz, fl = 0.10 Hz Figure 10(a) as presented in (3) , shows a large drift in the ground displacement record with a period of approximately 12 seconds and amplitude of approximately 40 mm (corresponding to a peak acceleration amplitude of 0.001 g).This low frequency component dominates the displacement signal (making interpretation of the record very difficult) and its retention cannot be justified with regard to the resolution of the digitisation system employed.

Figure 9 :
Figure 9: Comparison between MIT and New Zealand routine processing ELL Centro 1940 May 18 Component SOOE f L -0.32 Hz, ff = 0.224 Hz