5. Program Description and Algorithm

Program Description

\BIDO-win.bat         A batch file to activate MSYS on Windows. Not used on Linux.

\bin                       Executable files for use on Windows are stored in \bin\winbin, while executable files for use on Linux are stored in \bin\linbin. Executable files for use on MSYS are stored just beneath \bin.

\demo                   Contains demo data (used both Windows and Linux).

\etc                       Contains scripts to activate MSYS. Not used on Linux.

\run.sh                  A B shell script to activate the program (used both Windows and Linux).

\script                    Contains B shell scripts to link Fortran codes (used both Windows and Linux).

\src                       Contains Fortran codes (used both on Windows and Linux).

 

Algorithm

Description of the general flow and individual procedures. In parentheses are the file names of relevant B shell scripts and Fortran codes.

 

General Flow (\script\circle.sh)

1) Select portions of the data that are good to use.

2) Azimuthally average data around the circle (output from step 1) is not used here)

3) Estimate spectral densities by using the output from steps 1) and 2).

4) Estimate spectral ratios, phase velocities, NS ratios etc.

5) Repeat steps 3) and 4) as many times as there are segment clusters.

6) Calculate means and standard deviations using the output from step 4).

7) Plot the output from step 6).

 

[1] Automatic Selection of Segments (\script\mksegment.sh, \src\evalrms.F, \src\segment.F)

Segments are selected as follows:

1) The following procedure is performed on all components of all sensors.

- Subtract a linear trend from the original waveform data and calculate the RMS (let this be called RMSall).

- For every portion of the data with a prescribed segment duration into which the original waveform is divided (the portions are extracted so that they mutually overlap by half), subtract a linear trend and calculate the RMS. Normalize the RMS values by RMSall (let these be called RMSseg).

2) Make a histogram of all RMSseg, for all components of all sensors, at intervals of 0.1. Identify the interval of the largest frequency in the histogram.

3) Pick up data segments in which all RMSseg of all components and all sensors fall into the interval of the largest frequency simultaneously. Mark them as the segments good to use, and catalog them in the segment file.

 

[2] Azimuthal Averaging of Data around the Circle (\script\mkcrcle*.sh, \src\mkcrcl_*.F)

All methods adopted in this program start with taking weighted azimuthal averages of records around the circumference. Weighted azimuthal averaging corresponds to calculating Fourier coefficients in the Fourier series expansion around the circle. Our program calculates Fourier coefficients of the zeroth and first orders by default. Our theory is adaptable to unevenly spaced sensors around the circle, and adaptability to practical cases has been investigated for the CCA method (Reference [1]). We have not, however, closely investigated the adaptability for all methods, so we recommend the use of equidistant arrays to the extent that that is possible (we are particularly uncertain about the adaptability to methods that use cross-spectral densities).

 

[3] Estimating Spectral Densities (\script\estspec*.sh, \src\estspec.F)

This program eliminates the trend from each segment, applies a cosine-type taper with a 50-percent overlap (Carter et al., 1973), and estimates spectral densities by using both techniques of segment averaging and spectral smoothing for the raw FFT spectral densities (Bendat & Piersol, 1971). The raw spectral densities are smoothed with a Parzen window before they are segment-averaged. The segment duration and the number of segments in the averaging over multiple segments are given in the dialogue on activating the program: "Duration of data segments for the evaluation of spectra" and "Number of data segments over which averages are taken. Enter 0 or a very large number if you wish to use all segments simultaneously." The bandwidth of smoothing with a Parzen window is given in the dialogue: "Band width of the Parzen spectral window."

 

[4] Calculating Spectral Ratios (\script\specratio*.sh)

Ratios are taken, with no frills, between spectral densities estimated in the above-described procedure, except when the denominator is zero. Different types of spectral ratios are linked to the phase velocities via formulae described in References [2, 5]. Note that the autocorrelation coefficient of the SPAC method is defined here by a spectral ratio according to the formulation of Reference [2] (in usual practice, the SPAC coefficient is calculated as an azimuthal average of complex coherences between the central point and a peripheral point).

Calculating the H/V spectrum

Starting with Version 1.2.2, I added a feature that estimates horizontal-to-vertical (H/V) spectral ratios, provided that the data have three components, at the one station that is indicated at the top of the seism file (H/V refers to the ratio of the power-spectral density of the horizontal components to that of the vertical component. The power of horizontal motion is defined as the sum of the NS and EW component powers). Accordingly, even when the seism file describes a single measurement station alone (even if this does not constitute an array), H/V spectral ratios are calculated as long as there are three-component records. Once the calculation is over, the logarithmic mean and standard deviation are plotted.

Calculating Phase Velocities (\script\spec2pv*.sh, \src\sctr2pv.F)

Spectral ratios are equated to Bessel functions according to the formulae, and a root-solving method that combines bisection and the secant method (Shampine & Watts, 1970) is used to search for rk (radius times wavenumber), the argument of the Bessel functions, in the range [0, rkmax]. rkmax corresponds to the first maximum or minimum of the function value. The rk obtained is used to calculate the phase velocity, c=2p f/k (f stands for frequency).

 

[5] Repeat Steps [3] and [4] as Many Times as There Are Segment Clusters

 

[6] Calculating Means and Standard Deviations (\script\mkave.sh, \src\calave.F)

If the number of segments (number of segments over which averages are taken when estimating spectral densities with the segment averaging method. In other words, the integer value that you enter in the dialogue "Number of data segments over which averages are taken. Enter 0 or a very large number if you wish to use all segments simultaneously" on activating the program) satisfies

         (number of all segments catalogued in the segment file) > (number of segments),

then more than one spectral density estimates are obtained from the given waveform data. If we define

         (number of segment clusters) = (number of all segments catalogued in the segment file) / (number of segments),

there will be as many estimates of spectral densities, spectral ratios derived from them and phase velocities as the number of segment clusters (the remainder of division is discarded). This program calculates means and standard deviations on the basis of those estimates.

When the program is executed, folders with names \RESULT\(a number) are generated beneath the data folder. The number here represents that of a segment cluster, and these numerically named folders contain the corresponding analysis results. The \RESULT\ave folder contains output of statistical processing of the analysis results stored in those numerically named folders.

To eliminate "outliers," the maximum and minimum values are excluded from the statistical processing if and only if there are more segment clusters than the number NROBUST4AVERAGE_INC indicated in \src\PARAM.h. NROBUST4AVERAGE_INC is set at 8 by default. Some of the analysis results are averaged arithmetically and others logarithmically. Logarithmic averaging is used when averaging ratios like AmpRV_R.d, nsr.d, nsrlim_cca.d, nsrlim_cca.lwapx.d and powratio_R2L.d.

 

Bendat, J. S., and A. G. Piersol, Random Data: Analysis and Measurement Procedures, John Wiley & Sons, 1971.

Carter, G. C., C. H. Knapp, and A. H. Nuttall, 1973, Estimation of the magnitude-squared coherence function via overlapped Fast Fourier Transform processing: IEEE Transactions on Audio Electroacoustics, AU-21, 337344.

Shampine, L. F., and H. A. Watts, FZERO, a root-solving code, Report SC-TM-70-631, Sandia Laboratories, 1970.

 

[RETURN] [BIDO TOP]