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, 337–344.
Shampine, L. F., and H. A.
Watts, FZERO, a root-solving code, Report SC-TM-70-631, Sandia Laboratories,
1970.