Wiki » History » Version 3

Robert Porritt, 08/07/2013 01:11 PM

1 1 Robert Porritt
2 1 Robert Porritt
  current tools:
3 1 Robert Porritt
4 1 Robert Porritt
5 1 Robert Porritt
6 1 Robert Porritt
7 1 Robert Porritt
8 1 Robert Porritt
9 1 Robert Porritt
10 1 Robert Porritt
11 1 Robert Porritt
12 1 Robert Porritt
13 1 Robert Porritt
14 1 Robert Porritt
15 1 Robert Porritt
16 1 Robert Porritt
17 2 Robert Porritt
This package is distributed via the normal linux method. Simply untar
18 1 Robert Porritt
19 2 Robert Porritt
@tar -xvf station_analysis_tools.tar.gz@
20 1 Robert Porritt
21 2 Robert Porritt
cd to the new directory
22 1 Robert Porritt
23 2 Robert Porritt
24 1 Robert Porritt
25 2 Robert Porritt
use the attached configure script to create a makefile for your system
26 1 Robert Porritt
27 2 Robert Porritt
28 1 Robert Porritt
29 2 Robert Porritt
compile with make and if you wish for a full install make install
30 1 Robert Porritt
31 2 Robert Porritt
32 2 Robert Porritt
sudo make install@
33 1 Robert Porritt
34 2 Robert Porritt
To remove local executables, use make clean
35 2 Robert Porritt
@make clean@
36 1 Robert Porritt
37 1 Robert Porritt
38 2 Robert Porritt
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.
39 1 Robert Porritt
40 2 Robert Porritt
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..
41 1 Robert Porritt
42 2 Robert Porritt
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.
43 1 Robert Porritt
44 2 Robert Porritt
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.
45 1 Robert Porritt
46 2 Robert Porritt
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.
47 1 Robert Porritt
48 2 Robert Porritt
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.
49 1 Robert Porritt
50 2 Robert Porritt
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.
51 1 Robert Porritt
52 2 Robert Porritt
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.
53 1 Robert Porritt
54 2 Robert Porritt
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. 
55 1 Robert Porritt
56 2 Robert Porritt
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.
57 1 Robert Porritt
58 2 Robert Porritt
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.
59 1 Robert Porritt
60 2 Robert Porritt
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.
61 1 Robert Porritt
62 1 Robert Porritt
63 2 Robert Porritt
Programs are driven via a parameter file passed on the command line.
64 2 Robert Porritt
If run without argument or the -h option, then the parameter file format is displayed.
65 2 Robert Porritt
Parameter file format for compute_psd:
66 3 Robert Porritt
67 1 Robert Porritt
68 1 Robert Porritt
		n_sec n_skip unit_flag
69 1 Robert Porritt
70 1 Robert Porritt
71 1 Robert Porritt
	sac_file_name is the name of a file formatted for SAC (binary 158 array header, then the data)
72 1 Robert Porritt
	pole_zero_file_name is the name of a pole zero response file relative to displacement (as output from rdseed)
73 1 Robert Porritt
	n_sec is the number of seconds in each window. This determines maximum period to compute
74 1 Robert Porritt
	n_skip is the offset between windows. The welch's average calls for several overlapping segments to make a reasonably robust estimate
75 1 Robert Porritt
	unit_flag is 0 for displacement, 1 for velocity, or 2 for acceleration. Other values will quit the program.
76 1 Robert Porritt
	output_file_name is the name of a file which will be produced. It has three columns as:
77 1 Robert Porritt
		period	mean	standard_deviation
78 1 Robert Porritt
		period	mean	standard_deviation
79 1 Robert Porritt
		   .     .            .
80 1 Robert Porritt
		   .     .            .
81 1 Robert Porritt
                   .     .            .
82 3 Robert Porritt
83 1 Robert Porritt
84 1 Robert Porritt
85 2 Robert Porritt
in the example/ directory there is a data file from the Berkeley station YBH channel LHZ during a teleseism. (BK.YBH.LHZ.SAC)
86 2 Robert Porritt
Also included is the pole zero response (BK.YBH.LHZ.pz)
87 2 Robert Porritt
and a parameter file (param.psd) which can be used as an example.
88 1 Robert Porritt
89 1 Robert Porritt
90 1 Robert Porritt
91 1 Robert Porritt
No standard deviation computed in compute_tri_coherence yet and its unclear the standard deviations in compute_coherence and compute_incoherence are correct.
92 1 Robert Porritt
93 1 Robert Porritt
//FIXED on 8.29.2012 {
94 1 Robert Porritt
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. 
95 1 Robert Porritt
96 1 Robert Porritt
97 1 Robert Porritt
98 1 Robert Porritt
99 1 Robert Porritt
100 1 Robert Porritt
// FIXED on 3.29.2011 {
101 1 Robert Porritt
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.
102 1 Robert Porritt
103 1 Robert Porritt
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.
104 1 Robert Porritt
105 1 Robert Porritt
106 1 Robert Porritt
//FIXED on 2.19.2011 {
107 1 Robert Porritt
	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.
108 1 Robert Porritt
	Update on Sun SPARC tests (2.18.2011):
109 1 Robert Porritt
		fixed the segmentation fault error, but imaginary component is still wonky on the SPARC system.
110 1 Robert Porritt
		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.
111 1 Robert Porritt
112 1 Robert Porritt
113 1 Robert Porritt
114 1 Robert Porritt
	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
115 1 Robert Porritt
	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
116 1 Robert Porritt
	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.
117 1 Robert Porritt
	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.
118 1 Robert Porritt
	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.
119 1 Robert Porritt
	the src/ and include/ directories contain copies of the source and headers. The actual compiled source is in the main directory
120 1 Robert Porritt
	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.