# Seismological Research Letters

- © 2012 by the Seismological Society of America

*Online Material:* Matlab codes and data example of instrument corrections.

## INTRODUCTION

Of all the filters applied to recordings of seismic waves, which include source, path, and site effects, the one we know most precisely is the instrument filter. Therefore, it behooves seismologists to accurately remove the effect of the instrument from raw seismograms. Applying instrument corrections allows analysis of the seismogram in terms of physical units (e.g., displacement or particle velocity of the Earth’s surface) instead of the output of the instrument (e.g., digital counts). The instrument correction can be considered the most fundamental processing step in seismology since it relates the raw data to an observable quantity of interest to seismologists. Complicating matters is the fact that, in practice, the term “instrument correction” refers to more than simply the seismometer. The instrument correction compensates for the complete recording system including the seismometer, telemetry, digitizer, and any anti‐alias filters. Knowledge of all these components is necessary to perform an accurate instrument correction.

The subject of instrument corrections has been covered extensively in the literature (Seidl, 1980; Scherbaum, 1996). However, the prospect of applying instrument corrections still evokes angst among many seismologists—the authors of this paper included. There may be several reasons for this. For instance, the seminal paper by Seidl (1980) exists in a journal that is not currently available in electronic format and cannot be accessed online. Also, a standard method for applying instrument corrections involves the programs TRANSFER and EVALRESP in the Seismic Analysis Code (SAC) package (Goldstein *et al.*, 2003). The exact mathematical methods implemented in these codes are not thoroughly described in the documentation accompanying SAC.

We describe a general method for causal instrument correction that is applicable to data from a wide range of seismometers and present a set of codes for implementing the method. In doing so, we provide an alternative to the SAC instrument correction codes and show the detailed mathematics that form the basis of our method. Special attention is paid throughout to preserving the causal property of the instrument response, in order to allow observation of first motions and maintain relative timing information among different frequency components. Furthermore, the use of acausal filters can compromise the accuracy of first arrival picking. The method is presented in detail since causal instrument corrections are based on the precise distinction between analog and digital filters.

We test the method with data from two stations in the network operated by the Alaska Volcano Observatory (AVO). The two stations consist of co‐located broadband and short‐period seismometers within the local network at Mt. Spurr Volcano, Alaska. In addition to these two seismometers having different frequency responses, the output of the short‐period seismometer is transmitted over a radio telemetry system. Therefore, the output of the co‐located seismometers differs in several aspects with regards to the entire recording system. Correcting two raw seismograms from these stations to match each other in physical units poses a considerable challenge. We compare the instrument‐corrected seismograms for the two stations over a frequency band where both instruments are able to resolve common ground motions above their respective noise floors (0.1–10 Hz) and find that the instrument‐corrected particle‐velocity seismograms are in excellent agreement. Furthermore, we compare instrument‐corrected seismograms obtained by both our method and the SAC codes. The method we develop compares favorably with the SAC instrument correction, although the SAC implementation must be performed in a specific and non‐standard way to preserve causality.

The method we propose has several distinguishing features. Most importantly, we place emphasis on the preservation of causality in the instrument correction. By causality, we mean that the instrument‐corrected seismograms should be identically zero prior to the arrival of seismic energy. The importance of causality for instrument correction has recently been demonstrated on data from explosive events at Kilauea Volcano by Patrick *et al.* (2011). In their study, Patrick *et al.* (2011) showed that the relative arrival times of different frequency components of the wave field depended on whether the processing was performed with a one‐pass (causal) or two‐pass (acausal) filter. The relative arrival times of the different frequencies are critical measurements since they provide strong constraints on the physics of the source process.

We implement the bilinear transform (Seidl, 1980; Scherbaum, 1996; Dost and Haak, 2006) to design the instrument correction filter and, in the process, ensure causality in the filtering. However, instead of warping the frequency axis (Ferber, 1989; Scherbaum, 1996), as is commonly done with the bilinear transform, we rely on interpolation and oversampling to achieve accuracy for the correction. This is a distinguishing feature of our method that we feel simplifies the approach compared to frequency warping. To facilitate implementation, Matlab scripts and functions are provided online (available in the electronic supplement to this paper). The data example involving the aforementioned co‐located broadband and short‐period seismometers is included along with the Matlab codes. The data example constitutes the ultimate test for any instrument correction method: match the seismograms in physical units between two vastly different co‐located instruments. We show that the method we describe passes this requirement.

The codes we provide for instrument correction have been designed for use with data from the AVO network although, in principle, they can be applied with minor modifications to data acquired from many types of systems. The majority of the instrument responses for stations in the AVO network fall into one of two categories: either (a) broadband instruments with on‐board digitizers (Guralp CMG‐6TD) or (b) short‐period instruments (Mark Products L4, L4‐3D, and L22 or Teledyne‐Geotech S13) whose signals are sent through a telemetry system prior to digitizing (Dixon *et al.*, 2008). The telemetry includes a voltage‐controlled oscillator, or McVCO (McChesney, 1999), that affects the response of the system. The co‐located broadband and short‐period stations addressed in our data example come from each of these two categories. As of 2007, these two categories include all but 4 of the 193 permanent seismograph stations in the AVO network. The 4 exceptions include a strong‐motion seismometer (station code AU22) at Augustine Volcano, a CMG‐40T seismometer telemetered prior to digitizing (AKT), a Trillium‐40 seismometer (AMKA) on Amchitka Island, and an L4 short‐period seismometer that is telemetered using gain‐ranging on Mt. Spurr (NCG). Therefore, although the Matlab codes we provide could, in principle, be applied to raw data from any seismometer, they can directly be used to perform instrument corrections on more than 97% of the stations in the AVO network.

An IRIS‐supported program, the Portable Data Collection Center (PDCC), has been used to build dataless seed volumes for the AVO network (network code AV). The codes for performing causal instrument corrections, described in a later section of this paper, have served to check the accuracy of the dataless seed files for the AVO network. The dataless seed files and waveform data for the AV network are available from IRIS. In providing these codes and a data example, we aim to demonstrate the high quality of the data within the AV network and the reliability of the associated instrument response information.

## BROADBAND SEISMOMETERS

We begin with a description of instrument correction for broadband seismometers in the AVO network, since it is simpler and more generic than the corrections for the telemetered short‐period instruments. An overview of the different recording system components for broadband and short‐period seismometers in the AVO network is shown in Figure 1. As stated before, the methods we use are general and can be applied to data from any seismometer described by poles and zeros. We assume the seismograph is a causal, linear time‐invariant (LTI) system whose response in the Laplace transform domain is of the form: (1)where *G* is the gain with dimensions of counts/m/s, *C* is a normalizing constant that renders *T*/*G* equal to unity at a reference frequency *f*_{R}, *z*_{j} stands for the *L* zeros, and *p*_{k} stands for the *N* poles. This notation is in agreement with standard conventions (Type A transfer function; Standard for the Exchange of Earthquake Data [SEED], 2010). The variable in the Laplace domain is given by *s* and it is related to angular frequency by *s*=*iω*. With this specification for the seismograph response, the input and output of the system are related in the *s*‐domain by multiplication (convolution in the time domain): (2)where the output is the raw seismogram *S* and the input is particle velocity *V*, the time‐derivative of displacement. The seismogram has dimensions of digital counts and the particle velocity has dimensions of length over time. Note that a common way of performing a “poor‐man’s” instrument correction would be to simply multiply the raw seismograms by the inverse gain, 1/*G*. This is sometimes called a calibration factor or the mid‐band sensitivity. However, this correction would only be exact for signals in the flat portion of the instrument response, which includes the reference frequency. The raw seismogram would not be instrument‐corrected over a wide frequency band. In contrast, the data example we show in a later section is accurately instrument‐corrected from 0.1 to 10 Hz.

The crux of instrument correction is that we must design a digital filter that emulates the analog filter *T* in equation (1). The digital filter is implemented in a computer on digitized data. Thus, instrument correction begins with digital filter design. For this purpose, we make extensive use of the bilinear transform (Seidl, 1980; Scherbaum, 1996) for creating a digital filter from equation (1). The bilinear transform is a conformal mapping relating the analog Laplace transform variable to the digital *Z*‐transform variable: (3)where *Z* is the *Z*‐transform variable and Δ*t* is the time sampling interval. The bilinear transform can be considered as a particular Pade approximation of the exponential function. The desirable aspect of the bilinear transform is that it uniquely maps the left‐hand side of the complex *s*‐plane to the inside of the unit circle in the complex *Z*‐plane (Scherbaum, 1996). As a result, it preserves causality in the transition from analog to digital. The non‐desirable aspect of the bilinear transform is that it is not exact: it is an approximation, similar to a finite‐difference approximation of a derivative. An approximation is needed, however, to make the relationship between the Laplace‐transform domain and the *Z*‐transform domain expressible as a polynomial. An advantage of such a polynomial relationship is that the resulting filter can be represented as an infinite impulse response (IIR) filter in the time domain. In contrast to the bilinear transform, the most intuitive method for designing a digital filter is to simply sample the analog filter *T*(*s*) at all possible discrete frequencies by employing the relationship *s*=*iω* in equation (1). The drawback to this process is that the portions of the analog filter *T*(*s*) beyond the Nyquist frequency become aliased (Scherbaum, 1996; fig. 7.3) and can appear as acausal artifacts in the time‐domain. This type of aliasing is avoided by using the bilinear transform, since the bilinear transform as stated above maps the entire left‐hand side of the complex *s*‐plane to the inside of the unit circle in the complex *Z*‐plane (Scherbaum, 1996; fig. 9.15).

By substituting equation (3) for *s* in equation (1), the digital filter describing the seismometer is (4)where the notation *T*_{D} indicates the digital approximation to the analog filter *T*. Through some manipulation, this expression can be rendered into a form describing a digital filter (Oppenheim and Schafer, 1975; Press *et al.*, 1992) in terms of its digital poles and zeros: (5)where the factor *F*(*Z*) is given by (6)

The digital poles and zeros are given in terms of the analog poles and zeros, and the time sampling interval by (7)and (8)Note that the digital poles and zeros are dimensionless.

Equation (5) represents the final digital approximation to the analog seismometer. Notice that all the analog poles and zeros have corresponding digital poles and zeros. Whereas the analog poles and zeros lie in the left‐hand side of the complex *s*‐plane, the digital poles and zeros lie inside the unit circle. Depending on the number of analog poles and zeros, the factor *F*(*Z*) can contribute additional poles and zeros on the unit circle. Thus, the bilinear transform ensures that causality is preserved in the digital approximation to the seismometer. Once the digital approximation to the instrument response has been built, we implement deconvolution through spectral division of the raw seismogram by the response in the frequency domain.

Although straightforward in principle, there are some important details concerning the implementation of equation (5) for instrument corrections. The first concerns the accuracy of the digital approximation to the true analog filter. The bilinear transform is similar in many ways to a finite difference approximation to a derivative. That is, it works well for data that are highly sampled, but not for data that approach the Nyquist frequency. As a rule of thumb, we have found that the bilinear transform provides an acceptable approximation up to a frequency equal to one‐fifth of the Nyquist frequency. This poses somewhat of a problem since we would ideally like to correct data up to the Nyquist frequency. Our solution to this problem is to interpolate the data to at least five times the sampling rate. Such an interpolation renders the original Nyquist frequency to be one‐fifth of the new Nyquist frequency. As a result, the approximate digital filter in equation (5) is highly accurate over the original frequency band on the interpolated data.

Given the need for interpolation, another issue is which interpolation method should be used. The most obvious method would be sinc, or Fourier domain, interpolation. However, this type of interpolation does not preserve the causality property; that is, a single impulse becomes acausal when sinc‐interpolated at five times the sample rate. To preserve causality, we have found that nearest‐neighbor interpolation is preferred. In a later section of this paper, we demonstrate the causality‐preserving property of nearest‐neighbor versus sinc interpolation for oversampling the original data. After interpolating our seismograms with the nearest‐neighbor method and applying the instrument correction, we simply decimate the result to return to our original sampling rate. There is the possibility that nearest‐neighbor interpolation introduces new frequencies that were not present in the original signal; however, we have found this problem to be negligible compared to the advantage of preserving causality.

Note that the traditional way of dealing with the approximate nature of the bilinear transform is to warp the digital frequency axis, a process called “pre‐warping.” In principle, this matches the digital approximation to the analog filter at certain specified control points, for example the poles of the filter. In our testing, we have found that this method results in unacceptable errors away from the control points and have instead opted for interpolation and oversampling. As will be seen in the data example, the method we describe results in highly accurate instrument‐corrected seismograms.

Although the bilinear transform preserves causality, there is the possibility that numerical errors may arise that compromise this property in practice. For this reason we exploit the properties of real causal signals when building the instrument response. Recall that, due to the Kramers–Kronig relation, the real and imaginary parts of a causal signal are related. For the instrument response this means that (9)where *H* is the Hilbert transform. Note that this condition applies to both the analog and digital versions of the instrument response. Additionally, the values of the instrument response at positive and negative frequencies are related by virtue of the instrument response being real (10)where the asterisk denotes complex conjugation. Therefore, since the instrument response is real and causal, we need only to build its real part for positive frequencies—its imaginary part and its values for negative frequencies can be computed subsequently. This forces the causal property to be obeyed to the level of intrinsic numerical precision.

Ideally, one would like to apply an instrument correction over the entire frequency band, from 0 to the Nyquist frequency. However, this would cause instability due to the amplification of noise outside of the valid frequency band, which depends on the signal‐to‐noise ratio and how steeply the instrument response falls off from its passband (Scherbaum, 1996; fig. 9.5). As a result, both low‐ and high‐frequency cut‐offs need to be specified in practice. Another approach to stabilizing the filtering is to use a “water level method” as discussed in Scherbaum (1996) and Aster *et al.* (2005). For our implementation of instrument correction, we use a cascade of both a causal Butterworth low‐pass and causal high‐pass filter to define the band over which the signals are rectified. Limitations must be set on the causal band‐pass filter, however, since a tradeoff exists between the narrowness and sharpness of the filter and its stability. Claerbout (1992) states that a causal band‐pass filter “is almost a contradiction in terms” and demonstrates the types of instabilities possible when performing causal band‐pass filtering. In the Matlab codes described in a later section, we require that the high‐frequency cut‐off be at least 10 discrete frequency intervals greater than the low‐frequency cut‐off. The discrete frequency interval Δ*f* is equal to the inverse of the total time duration of the seismogram *τ*, or Δ*f*=1/*τ*. To further avoid instabilities, Claerbout (1992) suggests using low filter orders for causal Butterworth filters. We require that the filter orders of the low‐pass and high‐pass filters fall in the ranges of 3–7 and 2–4, respectively. In the data example described in a later section, we have chosen low‐pass and high‐pass filters of order 5 and 3, respectively, to define the band‐pass filter. These filter orders represent a tradeoff between stability and sharpness of the band edge.

The final aspect of instrument correction for a broadband seismometer concerns the filtering action of the analog‐to‐digital (A/D) converter or digitizer. In the course of A/D conversion for a broadband seismometer, a series of decimation steps change the input signal to lower and lower sampling rates until the desired sampling rate is reached. An ideal A/D converter would act like an anti‐alias low‐pass filter equal to unity below the decimation frequency and zero otherwise. In practice, A/D converters filter the raw seismograms below the decimation frequency, notably resulting in precursory oscillations prior to earthquake phase arrivals (Scherbaum, 1996). In a later section of this paper, we show an example of these oscillations. The precursory oscillations are always composed of frequencies close to the Nyquist frequency. This is because the digital anti‐alias filters within the A/D converter do a good job emulating an ideal low‐pass filter away from the Nyquist frequency, but a poor job close to it. The precursory oscillations are troubling in the sense that they violate the desired causal property of the instrument response. Our practical solution to this complication is to define the high‐frequency cut‐off for the instrument‐correction well below the Nyquist frequency. This has the effect of filtering out the frequencies corrupted by the precursory oscillations. In addition, well below the Nyquist frequency there is almost no filtering action from the A/D converter, besides a constant multiplier converting the output of the seismometer (volts) to the digital output (counts or bits).

A more rigorous method for dealing with the precursory oscillations has been proposed by Scherbaum (1996). This method modifies the acausal response of the A/D converter (i.e., the source of the precursory oscillations) to be its causal equivalent. Although valid in principle when applied to data prior to decimation, this method encounters difficulties in practice when applied to data after decimation (Scherbaum, 1996; fig. 8.9). The origin of this problem lies in the generally unknown type of decimation scheme used within the A/D converter. Since the decimation schemes for instruments within the AVO network are not known, we have chosen the practical solution described previously to deal with the precursory oscillations; that is, we have chosen a high‐frequency cut‐off well below the Nyquist frequency. In the data example shown later, we use a high‐frequency cut‐off of 10 Hz for data with a 25 Hz Nyquist frequency. As a result, the problem of the precursory oscillations is solved, but at the cost of reduced high‐frequency content in the instrument‐corrected seismograms.

## SHORT‐PERIOD SEISMOMETERS

The treatment of short‐period seismometers includes all the steps previously described for broadband seismometers, except for the issue of decimation within A/D conversion, which does not apply. However, additional complications arise for the short‐period seismometers since they are transmitted through a telemetry system that further filters the data stream, as shown in Figure 1. The telemetry system adds additional poles and zeros to the instrument response, besides those related to the seismometer itself. In addition, the telemetry system causes a slight delay in the process of time‐stamping the data. We describe how to treat this pure delay present for the short‐period data. Finally, an additional anti‐alias filter is applied to some of the short‐period data in the AVO network (Fig. 1).

For the instrument corrections with short‐period seismometers, we use data directly measured by network engineers at AVO. These data include the delay due to the telemetry system and the frequency–amplitude pairs for an anti‐alias filter applied at some of the stations. The network engineers estimate the telemetry delay to be 55 ms for short‐period stations with an anti‐alias filter and 35 ms for stations without such a filter. For short‐period stations that have an anti‐alias filter, the network engineers have measured the filter’s amplitude at certain frequencies (i.e., frequency–amplitude pairs). Two anti‐alias filters are used within the AVO network. In Figure 2, we show the measured frequency–amplitude pairs for the two types of anti‐alias filters. Also plotted in Figure 2 are fits to those data points using a particular kind of digital filter. The fitting is necessary since we need the response of these filters at frequencies other than the frequencies where measurements were made in order to perform instrument corrections. Instead of simply interpolating the data, we choose to compute a best‐fitting finite‐impulse‐response (FIR) digital filter using least‐squares. We do so since we need to define a filter, as opposed to frequency–amplitude pairs, in order to deconvolve the effects from the instrument. However, since only a minimum phase filter has a stable inverse, we modify the FIR filter *A*(*ω*) to its minimum phase equivalent according to the relation (Claerbout, 1976) (11)which is a type of Kolmogoroff spectral factorization. This filter has the same amplitude spectrum as the FIR filter, but with the minimum delay property. The fact that it has a stable inverse means it can be used for instrument corrections. This *ad hoc* treatment of the anti‐alias filter could be avoided if we had the filter’s poles and zeroes instead of frequency–amplitude pairs. However, we feel this approach is an acceptable use of the frequency–amplitude data for the purpose of practical instrument correction.

The frequency–amplitude data only cover the frequency interval up to 50 Hz, which is the Nyquist frequency for the short‐period instruments in the AVO network. Because of this, we cannot apply the correction for the anti‐alias filter simultaneously with the instrument correction, since the instrument correction is applied to oversampled data with a higher Nyquist frequency. As a result, we apply the correction for the anti‐alias filter after applying the instrument correction, once the data has been decimated back to its original sampling rate. We apply the correction by convolving the instrument‐corrected seismogram with the inverse of the minimum phase filter given in equation (11).

Returning to the issue of telemetry delay, we correct for the delay after instrument correction, but prior to decimation of the oversampled data, by explicitly removing it. For the case when the delay is a multiple of the time sampling interval for the oversampled data, this simply involves shifting the samples of the instrument‐corrected seismogram. However, in the more general case where the delay is not a multiple of the time sampling interval, we face a problem. A Fourier‐domain time shift over a fractional number of samples induces acausal behavior, much in the same way as sinc interpolation breaks causality, as described earlier. To avoid this, we use a causal fractional delay function within Matlab. This ensures that we account for the fractional delay while at the same time preserving causality. More details about the implementation of the fractional delay filter are given in Appendix A.

The poles and zeros for a short‐period instrument include those related to the seismometer and those related to the telemetry system. For example, an underdamped seismometer with a natural frequency *f*_{0} and a damping coefficient *h* has two poles given by (Scherbaum, 1996) (12)and two zeros at 0 Hz. Note that for an underdamped seismometer *h*<1. In addition to these poles, there are other poles related to the telemetry system. In the AVO network, the telemetry is composed of a voltage‐controlled oscillator (McChesney, 1999) and a discriminator. These act as the send and receive components in the radio communication. Both of these components are modeled as low‐pass Butterworth filters with a corner frequency at 30 Hz. Therefore, the additional poles are those computed for a low‐pass Butterworth filter. Using these analog poles and zeros, we proceed by applying the bilinear transform to compute digital poles and zeros just as in the case of the broadband seismometers described earlier.

## DATA EXAMPLE AND MATLAB CODES

We now demonstrate our method of instrument correction with original Matlab codes and a data example from the AVO network. The instrument correction codes are described in Appendix B. In Figure 3, we show the result of the data example using the co‐located short‐period and broadband seismometers at Mt. Spurr. The broadband seismometer, known by its station code SPCR, consists of a Guralp CMG‐6TD. The short‐period instrument consists of a Mark Products L4 and sits at station CKT. Note that the seismograms in Figure 3 have been instrument‐corrected to particle velocity over the frequency band from 0.1 to 10 Hz. To compute displacement, the particle velocity seismograms would need to be numerically integrated. The Matlab codes produce a version of Figure 3 in which the two instrument‐corrected seismograms are plotted together. For improved visualization here, we have plotted the two seismograms in separate panels in Figure 3. The seismograms are 5 minutes (300 s) long, with the short‐period sampled at 100 Hz and the broadband at 50 Hz. The data are from 02:29 to 02:34 UTC on 5 November 2005. It is clear from Figure 3 that the instrument‐corrected seismograms are highly similar; in subsequent figures we illustrate different parts of the seismograms in detail and quantify the similarity. The three different parts of the seismograms that warrant further analysis are the microearthquake at ∼55 s, the large amplitude signal in the short‐period seismogram at ∼120 s, and the low‐frequency microseismic noise persisting throughout both recordings. We selected these data due to the variety of signals occurring within this short time interval.

We zoom in on the microearthquake in Figure 4. This is a high‐frequency volcano‐tectonic earthquake with an epicenter about 15 km southeast of Mt. Spurr. The *M*_{L} 1.5 earthquake occurred at 16.5 km depth below sea level. The origin time of the earthquake was estimated to be 52.2 s from the beginning of the seismograms (02:29:52.20 UTC on 5 November 2005). Clear *P*‐wave and *S*‐wave arrivals register at the co‐located short‐period and broadband seismometers. The *P*‐wave arrival has an upward first‐motion in the particle velocity seismogram. After numerical integration to displacement, this upward first‐motion would persist as a compressional arrival. Although some minor differences exist between the seismograms in Figure 4, the microearthquake signals are observed to be highly similar. The high degree of similarity can be appreciated in Figure 5, which plots the raw seismograms in counts. The two raw seismograms are substantially different and contain a significant amount of incoherent high‐frequency signal. The precursory oscillations discussed earlier are evident in the raw seismogram for the broadband seismometer.

We zoom in on the large amplitude signal occurring at ∼120 s time in the short‐period seismogram in Figure 6. Although we do not plot it here, the raw seismogram for the short‐period instrument contains a single large, negative amplitude data point at this same time. This is indicative of a noise spike due to the radio telemetry system that the short‐period instrument passes through. The noise spike is not related to the instrument response and so it results in high‐amplitude noise in the instrument‐corrected seismogram. Interestingly, the instrument‐corrected signal arising from the spike can be considered to be an approximation of the “impulse response” of the instrument‐correction filter. As seen in Figure 6, the approximate impulse response has a sharp onset and reverberates afterward, a consequence of the causal property of the instrument correction.

In Figures 7 through 9, we further explore this issue by examining the actual impulse response of the instrument‐correction filters for both the short‐period and broadband seismometers. A comparison of the impulse responses is shown in Figure 7 for an impulse at a time of 100 s. The impulse response for the broadband seismometer SPCR more closely approximates an ideal impulse than the impulse response for the short‐period seismometer CKT. Note the similarity between the impulse response for CKT and the instrument‐corrected noise spike in Figure 6. In Figure 8, we re‐plot Figure 7 but with a vertical axis exaggerated by a factor of 1000. We show the plot in Figure 8 in order to emphasize the causal property of the impulse response; that is, the instrument‐corrected seismogram is zero prior to the impulse at a time of 100 s. This brings up the issue we discussed earlier concerning which interpolation scheme should be used for oversampling: sinc or nearest‐neighbor. The plot in Figure 8 is the result for our preferred interpolation scheme, nearest‐neighbor. In Figure 9, we show the same plot as in Figure 8 expect that we use sinc interpolation for oversampling. As we argued earlier, sinc interpolation results in an acausal instrument‐correction, as seen by comparing Figures 8 and 9. From these figures, we conclude that the instrument corrections we have designed are indeed causal.

Returning to Figure 6, we note the high degree of similarity for the continuous microseismic noise between the seismograms. This similarity persists over the entire 300 s of the recording. This similarity is remarkable considering that the microseismic noise exists below the natural frequency of the short‐period seismometer (1 Hz). Thus, for these frequencies the amplitude spectrum of the short‐period instrument response rolls‐off sharply and the phase spectrum has considerable variation. The agreement between the two seismograms testifies to the precision of the instrument response information for the short‐period seismometer.

In Figure 10, we attempt to quantify the similarity observed in the two instrument‐corrected seismograms. To do so, we compute the cross‐spectrum between the seismograms over the first 100 s and analyze the cross‐spectral amplitude and phase. We limit the analysis to the first 100 s since the noise spike occurs after this time and obviously degrades the amount of similarity between the seismograms. In the top panel of Figure 10, we show the amplitude of the cross‐spectrum or coherence. The two seismograms are highly coherent in the band from 0.1 to 0.7 Hz and from 3 to 10 Hz. These two frequency bands represent the parts of the spectrum where actual signal exists, either microseismic noise (0.1–0.7 Hz) or earthquake arrivals (3–10 Hz). The middle panel of Figure 10 shows the relative timing error between the seismograms when the amplitude of the coherence exceeds 0.65. Note that this value for the coherence is exceeded in the same 0.1–0.7 Hz and 3–10 Hz frequency bands. To compute the relative timing error, we converted the phase of the cross‐spectrum to phase delay by dividing by angular frequency. We then divided the phase delay by the period to obtain relative timing error. The relative timing error does not exceed 0.2% in absolute value over the frequency bands in which a high degree of coherence exists. Finally, the logarithmic power ratio is given in the bottom panel when the amplitude of the coherence exceeds 0.65. Perfect matching of the power spectra between the two seismometers would result in a logarithmic power ratio equal to zero. As can be seen in the bottom panel, the logarithmic power ratio fluctuates around zero over the entire frequency band when the coherence exceeds 0.65. This confirms that the instrument correction has effectively equalized the amplitude and phase for the two seismograms over a wide frequency band, from 0.1 to 10 Hz.

We further investigate the coherency of instrument‐corrected signals between the co‐located seismometers in Figure 11 by analyzing a regional earthquake instead of a local microearthquake. We perform this test to check whether the coherency is high in the frequency band from 1 to 3 Hz that lacked signal in Figure 10. A regional earthquake has a broader spectrum than a local microearthquake and should have signal within this band. As seen in Figure 11, the coherency for the regional earthquake is high for all frequencies above 0.1 Hz. Since the frequency content of the signal is more continuous than that of the microearthquake considered previously, a subtle linear trend is evident in the relative timing between the two instruments, as shown in the middle panel. This small trend may represent the limits of the precision with which the timing, or the phase portion of the instrument response, is known for these instruments. The logarithmic power ratio in the bottom panel demonstrates that the power spectra for the two signals are highly similar.

Finally, in Figure 12 we compare instrument‐corrected seismograms for the broadband station SPCR using the method we have developed in this paper and the SAC command TRANSFER. An important detail is that to produce a causal instrument response, the SAC commands must be executed in a specific and non‐standard way. This issue is discussed in Appendix C. As shown in Figure 12, our method for instrument correction and SAC produce highly similar results for the microearthquake plotted in Figure 4. This validates our method against the widely used SAC instrument correction capabilities. Furthermore, we have found a suitable sequence of SAC commands for producing a causal instrument correction.

We wish to point out that an effort must be made on the part of the users of the provided codes to preserve the causal property of the instrument corrections in subsequent processing of the particle velocity seismograms. For instance, numerical integration to displacement is causal but often a low‐frequency noise or drift exists afterward. The usual approach to this problem is to apply a high‐pass filter after numerical integration. The high‐pass filter in this case would need to be causal itself. This also applies to any band‐pass filtering of the instrument‐corrected seismograms. We have limited the discussion of instrument correction in this paper to the central step of converting from counts to particle velocity. We leave subsequent integration or filtering to the discretion of the user.

## DISCUSSION

In this paper, the overriding goal has been to forge a method for instrument correction that preserves the property of causality. The data examples described above show that we have achieved this goal. Certainly, the preservation of causality makes sense for real seismic signals—the seismometer does not begin to shake prior to the arrival of seismic waves. However, one may wonder about the drawbacks of a causal instrument correction. Are there times when a causal instrument correction is not preferred, does not work well, or is not important? Regarding this issue, the main drawback to a causal instrument correction is that it can be prone to numerical instability, for the reasons discussed earlier regarding the low‐pass and high‐pass Butterworth filters. As discussed by Claerbout (1992), the inherent risk of instability with causal bandpass filters can be appreciated from equation (11). A bandpassed signal has an amplitude spectrum *A*(*ω*) close to zero outside of the passband; however, the logarithm of a number close to zero, as shown in equation (11), can become numerically unbounded. An example of causality not being important involves the instrument correction for surface waves. Surface waves, being a type of normal mode, do not, strictly speaking, satisfy causality (Aki and Richards, 1980; Box 6.7). This is consistent with the fact that there is no such thing as a “first motion” for a surface wave. Therefore, the issue of causality for instrument corrections of surface waves is not critical. However, when comparing surface wave recordings between different instruments, the amplitude and phase differences between the seismometers must be compensated for in a consistent way. Our point is that this compensation, or instrument correction, does not necessarily have to be causal for accurate surface wave analysis. Another example of causality not being important is for analysis based on the power spectrum. In this case, the phase spectrum is not utilized and whether causality is preserved or not is immaterial.

## CONCLUSIONS

We have presented a method to perform instrument corrections and provided Matlab codes to apply this method to raw seismograms. We have also conducted the strongest test of the instrument corrections by applying the method to a co‐located broadband and short‐period seismometer deployed in the field. We find that the algorithm accurately corrects the raw seismograms over a wide frequency band. We have further tested and proven the method by comparing it to the standard instrument correction provided within the SAC package. We encourage others to seek out instances of co‐located seismometers in seismic networks to test the accuracy of known instrument response information. Recently, Ringler and Hutt (2010) and Ringler *et al.* (2011) have demonstrated the utility of coherency tests between co‐located seismometers.

The fact that the oceanic microseism is accurately recorded on the short‐period seismometer has implications for the emerging technique called ambient noise tomography (ANT). The observation that the L4 seismometer can accurately record the oceanic microseism has been made previously by Riedesel *et al.* (1990). Although some of the earliest implementations of ANT used networks of exclusively broadband seismometers (Gerstoft *et al.*, 2006), short‐period instruments can also be used for frequencies above 0.1 Hz. The ability to perform accurate instrument corrections becomes important when applying ANT to networks consisting of both short‐period and broadband seismometers (Masterlark *et al.*, 2010). Without equalizing the phase between different types of seismometers, delay times obtained through cross‐correlation of ambient noise may not reflect only propagation delay between the stations. In addition, the codes we provide allow the analysis of low‐frequency signals recorded on a network made up of different types of instruments, since the causal instrument correction accurately removes the instrument response.

Although the codes we provide have been tested on data from the AVO network, we expect that they are applicable to a wide variety of seismometers as well as accelerometers. This is especially true of the code provided for broadband seismometers, since no AVO‐specific filters are involved. The instrument correction for short‐period seismometers involves many AVO‐specific filters, but we believe the description has value as a specific example of incorporating non‐ideal types of instrument metadata (e.g., frequency–amplitude pairs instead of poles and zeros). Since the codes are provided in Matlab, we envision these codes being widely used in current seismological investigations.

## ACKNOWLEDGMENTS

This work has been supported in part by the USGS Mendenhall Program. Further support has been provided by USGS/ARRA award G10AC00016. Matthew Haney was formerly with Boise State University. We thank Stephanie Prejean (USGS‐AVO) for assistance with the SAC instrument response comparison. This manuscript has benefitted from reviews by Dan McNamara (USGS) and David Wilson (USGS). The comments of an anonymous reviewer have also improved the manuscript. The efforts of the network engineers at the Alaska Volcano Observatory are gratefully acknowledged.

## APPENDIX A. FRACTIONAL DELAY

The causal fractional delay function we utilize within Matlab is not a standard Matlab function or a function within the popular Signal Processing Toolbox. All other functions included in the codes we have developed exist as either standard Matlab functions or as functions within the Signal Processing Toolbox. The fractional delay function exists in the Filter Design Toolbox for Matlab version 2010a and the DSP System Toolbox for Matlab version 2011a. Prior to implementing the fractional delay, we check on the availability of these toolboxes. In the case that the toolbox containing the fractional delay function is not available, the fractional part of the delay is not applied as part of the instrument correction.

In the worst case, when the necessary toolbox is not available, the omission of the fractional delay causes a small timing error in the instrument correction. The error is by definition small, though, since it is a fraction of the time sampling interval of the oversampled data. For the short‐period seismometer example in this paper, the data had a sample rate of 100 Hz and were subject to 55 ms of telemetry delay. For an oversampling rate of 5, the sampling rate of the data becomes 500 Hz, or 2 ms per sample. Thus, in terms of a multiple of the time sampling interval, the telemetry delay could be approximated as 54 ms. The additional 1 ms of delay in this case would be the resulting small timing error when the fractional delay function was unavailable. Note that this small error could be avoided for the short‐period data by choosing an oversampling rate of 10 instead of 5. An oversampling rate of 10 would result in a sampling rate of 1000 Hz, or 1 ms per sample. Thus, the 55 ms of telemetry delay would be a simple multiple of the time sampling interval and the fractional delay would not be necessary. Although these complications can be avoided if the fractional delay function is available, even when it is not available the error can in general be made to vanish by properly choosing the rate of the oversampled data.

When the necessary toolbox is not available, another option consists of shifting the sample times associated with the data by the necessary delay. In this case, the delay would not be a multiple of the sampling time. Thus, nothing is done to the data and changes are made to the associated time raster. For the comparison between the co‐located seismometers shown in this paper, such an approach would result in the instrument‐corrected seismograms existing on different time rasters. Although this process is not ideal, it does preserve causality.

## APPENDIX B. MATLAB EXAMPLE

The entire package is provided as a ZIP file at the following web address: http://volcanoes.usgs.gov/misc/haney/causal_instrument_corrections_final.zip.

The methodology and examples have also been built into the response tools in the GISMO Matlab Tools (Reyes and West, 2011). The ZIP file consists of the following files:

3 Matlab scripts:

example_instrum_correct.m

rm_instrum_resp_spcr.m

rm_instrum_resp_ckt.m

1 Matlab function:

rm_instrum_resp.m

2 ascii data files:

spcr_20051105022900_20051105023400.txt

ckt_20051105022900_20051105023400.txt

1 text README file:

README_instrum_corrections.txt

The data example is executed by putting all these files into a directory, going to that directory within Matlab, and typing:

>> example_instrum_correct

at the Matlab command prompt. This main script runs the other two scripts and produces a figure that overlaps the instrument‐corrected seismograms for the co‐located broadband (station code SPCR) and short‐period (CKT) seismometers. The three Matlab scripts and one Matlab function have a significant amount of inline comments that describe the different processes within the codes.

The Matlab function “rm_instrum_resp.m” applies the instrument correction to a raw seismogram. This function requires the 12 input values listed here:

Sample rate, in Hz

Low‐frequency cutoff, in Hz

Order of the high‐pass Butterworth filter

High‐frequency cutoff, in Hz

Order of the low‐pass Butterworth filter

Bad data value

Analog zeros, in rad/s

Analog poles, in rad/s

Inverse gain 1/

*G*, in m/s/countReference frequency

*f*_{R}, in HzOversampling rate (=5 by default)

Intrinsic system delay, in ms (=0 s for broadbands, 35 ms or 55 ms for AVO short‐periods)

The 6th input parameter, the bad data value, refers to the value given to samples in the seismogram when no data are received (e.g., a telemetry dropout). In the acquisition system used by the AVO network, this value is set to the most negative 32‐bit integer, −2^{∧}31. By forcing the user to specify this, we intend to keep the user aware that our Matlab instrument correction codes do not solve the “missing data” problem. It is up to the user to somehow fill in the missing data according to a certain methodology before using these codes. By specifying a bad data value and then attempting to run the codes on a raw seismogram with bad data, the codes will give an error message and stop. All other input parameters have been discussed already in this paper. Note that it is important to be aware that the analog poles and zeros are to be given in angular frequency (radians/s) instead of linear frequency (Hz).

## APPENDIX C. SAC INSTRUMENT CORRECTION

To apply a causal instrument correction in SAC, we execute the following series of commands on the raw seismogram, which for this example has the filename SPCR.BHZ.AV:

SAC> r SPCR.BHZ.AV

SAC> rmean

SAC> rtr

SAC> taper width 0.05

SAC> lowpass npoles 5 corner 10

SAC> rmean

SAC> rtr

SAC> taper width 0.05

SAC> highpass npoles 3 corner 0.1

SAC> rmean

SAC> rtr

SAC> taper width 0.05

SAC> trans from evalresp to vel

SAC> rmean

SAC> rtr

SAC> taper width 0.05

SAC> w SPCR.BHZ.AV.corr

These series of SAC commands can be combined in a single SAC macro file. Note that the trans command requires the presence of a RESP file in the directory for station SPCR. The RESP file, called RESP.AV.SPCR..BHZ, can be obtained from the AVO miniseed available on the IRIS website. The RESP file for SPCR is also available at: http://www.iris.edu/mda/AV/SPCR

In the above instance of the trans command, we deliberately avoid using the freqlimits option that is commonly invoked with this SAC command. The freqlimits option is intended to define a passband over which the instrument‐corrected seismogram is returned. However, as implemented in SAC, the freqlimits option uses a two‐pass, zero‐phase filter instead of a one‐pass, causal filter. This fact is pointed out in the SAC Users Manual. Thus, the instrument‐corrected seismogram returned after invoking the freqlimits option is acausal. To define a passband while also ensuring a causal response within SAC, we have elected to apply causal high‐pass and low‐pass filters prior to correcting with the trans command. The high‐ and low‐frequency corners of these filters, respectively 10 Hz and 0.1 Hz, and number of poles have been chosen to match the example in this paper. Finally, we generously utilize tapering, mean removal, and trend removal prior to and following all filtering steps. Employing these additional steps ensures that instabilities and errors due to the finite length of the seismogram are avoided.