Station Analysis Tools

STATION_ANALYSIS_TOOLS
current tools:
compute_psd
compute_coherence
compute_incoherence
compute_pole_zero_response
general_pdf
get_mean_general_pdf
smooth_general_pdf
get_nth_percentile_from_pdf
interpolate_psd
spline_interpolate_psd
log10_linear_interpolate
compute_tri_incoherence

INSTALL:
This package is distributed via the normal linux method. Simply untar

tar -xvf station_analysis_tools.tar.gz

cd to the new directory

cd STATION_ANALYSIS_TOOLS

use the attached configure script to create a makefile for your system

./configure

compile with make and if you wish for a full install make install

make
sudo make install

To remove local executables, use make clean
make clean

EXECUTABLE DISCRIPTIONS:
compute_psd: A welch's average method in which we break up a seismogram into several overlapping windows, compute the spectral power (real * real + imag * imag) normalized by the number of points and sample rate (see paper by McNamara and Buland, 2004). This implementation returns the mean and standard deviation from those overlapping estimates. The function contains subroutines to remove a pole-zero response (using standard sac format) and convert the recorded signal to ground acceleration. Waveform data and header data can only be read in standard sac format.

compute_coherence: Similar in implementation to compute_psd. This code however returns the magnitude squared coherence which is a function of similarity between two waveforms in the frequency domain..

compute_incoherence: Similar implementation as compute_psd. Returns (1-coherence) * psd of a pair of waveforms. It becomes unstable when coherence = 1 therefore we put a large negative value (-999) to indicate a coherence of 1.

compute_pole_zero_response: This is essentially an educational tool. Given a sac formatted pole zero response file, it computes the spectral characteristics similar to evalresp. Unlike the evalresp library/function distributed by IRIS, this does not work on RESP.. files.

general_pdf: This reads a list of ascii files with x values in column 1 and y values in column 2. It stacks in a probabilistic sense giving the probability of a measurement landing in a grid cell which is defined by dx and dy given in the parameter file. The output returns the input x and y in columns 1 and 2 with the probability in column 3.

get_mean_general_pdf: This parses a file output from general_pdf and simply returns the y values with the highest probability. This is not robust to bi-modal distributions, and thus a more stationary median value can be obtained from get_nth_percentile_from_pdf with the percentile value of 50 given.

smooth_general_pdf: Applies a smoothing function onto a pdf. This helped stablize values returned by get_mean_general_pdf and makes the graphical representation of a pdf cleaner.

get_nth_percentile_from_pdf: Finds the y value at each x in a pdf where the values to the left have an N percent chance of occurring. To obtain the median, use an N of 50. If N is 10, then there is a 10% chance of a measurement being less than the line extracted.

interpolate_psd: Returns a regularly spaced linearly interpolated spectral measurement (psd, coherence, or incoherence). If the step given is small, large files result, but if its too big, then the result will be undersampled. Helpful when using the gmt 'surface' command to grid.

spline_interpolate_psd: Not very useful. Idea was to use a spline interpolate instead of a linear function. However, it creates a weird shape near the long period end and it does not really return a significant improvement. Currently still in a debugging mode. Included to give the splining functions which may be useful for future implementations.

log10_linear_interpolate: Similar to interpolate_psd in that it implements a linear interpolant function. However, the output array is hardwired to go from the minimum to the maximum (explicitly avoiding extrapolation) using a base 10 log scale.

compute_tri_incoherence: In the event of having three co-located and co-aligned sensors, this tool follows Sleeman et al., 2006 to compute a relative transfer function between two signals to find the noise of the third. It repeats this for all three channels.

RUN:
Programs are driven via a parameter file passed on the command line.
If run without argument or the -h option, then the parameter file format is displayed.
Parameter file format for compute_psd:
sac_file_name
pole_zero_file_name
n_sec n_skip unit_flag
output_file_name

sac_file_name is the name of a file formatted for SAC (binary 158 array header, then the data)
    pole_zero_file_name is the name of a pole zero response file relative to displacement (as output from rdseed)
    n_sec is the number of seconds in each window. This determines maximum period to compute
    n_skip is the offset between windows. The welch's average calls for several overlapping segments to make a reasonably robust estimate
    unit_flag is 0 for displacement, 1 for velocity, or 2 for acceleration. Other values will quit the program.
    output_file_name is the name of a file which will be produced. It has three columns as:
        period    mean    standard_deviation
        period    mean    standard_deviation
           .     .            .
           .     .            .
                   .     .            .

EXAMPLE:
in the example/ directory there is a data file from the Berkeley station YBH channel LHZ during a teleseism. (BK.YBH.LHZ.SAC)
Also included is the pole zero response (BK.YBH.LHZ.pz)
and a parameter file (param.psd) which can be used as an example.

KNOWN BUGS:
3.29.2011
No standard deviation computed in compute_tri_coherence yet and its unclear the standard deviations in compute_coherence and compute_incoherence are correct.

//FIXED on 8.29.2012 {
The compute_pole_zero_response.c program assumed all zeros at the origin when it removed a zero from the pole-zero structure. This has been updated to properly search for zeros located at the origin which typically determine the number of differentiations.
}

// FIXED on 3.29.2011 {
Welch's averaging must occur prior to the amplitude being calculated. Skipping this step leads to a coherence defined as 1. Further, it was found that to replicate other results I needed to change the cosine taper to a hann taper and correct with a standard value of sqrt(8/3) for the hann window.
3.25.2011
Coherence computation only considers the real part of the cross spectra. The averaging of the cross spectra has to be done over multiple parts before the amplitude is taken. This requires a change from the recursive averaging method currently implemented. Working on comparing output with GMT spectrum1d and Matlab mscohere functions.
}

//FIXED on 2.19.2011 {
Tests on Sun SPARC i86pc systems with gcc 3.4.3 cause segmentation faults during the ifft. A simple test of the fftw3 library on this system returns a non-sensical imaginary component which may be leading to this error. If you encounter this, upgrade your gcc and recompile fftw3 and this code.
Update on Sun SPARC tests (2.18.2011):
fixed the segmentation fault error, but imaginary component is still wonky on the SPARC system.
also, if an empty file is returned and the code segmentation faults, remove the empty file and retry. The files distributed were built on a 64bit system which might not be able to be overwritten by a 32 bit system.
}

BE AWARE:
Currently configured to have a maximum of 25 poles and 25 zeros in the response. This can be adjusted in station_analysis_tools.h on lines 18 and 19
Currently configured to read a maximum of 17,200,000 points from a sac file. Should you need more, adjust NPTS_SAC_FILE_MAX in the station_analysis_tools.h file and recompile
Spectral measurements (psd, coherence, incoherence) use a running mean smoothing algorithm for cleaner presentation. The calls are found in library functions "calculate_psd", "calculate_coherence", and "calculate_incoherence". A simple commenting out (// at the beginning of the line) will remove this.
Does not currently support miniseed, evalresp, or variable smoothing options. Libraries exist for miniseed and evalresp calls, but to reduce the dependence on external libraries, I chose to not code those into this package.
No plotting package is included. I suggest using gmt to create plots, but you can also use gnuplot, octave, matlab, or any other plotting package which can read ascii data.
the src/ and include/ directories contain copies of the source and headers. The actual compiled source is in the main directory
The sac file format may have differences due to 32 vs 64 bit and endian style differences. Those are not accounted for in the read_sac routine included in the lib_station_analysis_tools.c code.