## Wiki » History » Version 3

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

1 | 1 | Robert Porritt | STATION_ANALYSIS_TOOLS |
---|---|---|---|

2 | 1 | Robert Porritt | current tools: |

3 | 1 | Robert Porritt | compute_psd |

4 | 1 | Robert Porritt | compute_coherence |

5 | 1 | Robert Porritt | compute_incoherence |

6 | 1 | Robert Porritt | compute_pole_zero_response |

7 | 1 | Robert Porritt | general_pdf |

8 | 1 | Robert Porritt | get_mean_general_pdf |

9 | 1 | Robert Porritt | smooth_general_pdf |

10 | 1 | Robert Porritt | get_nth_percentile_from_pdf |

11 | 1 | Robert Porritt | interpolate_psd |

12 | 1 | Robert Porritt | spline_interpolate_psd |

13 | 1 | Robert Porritt | log10_linear_interpolate |

14 | 1 | Robert Porritt | compute_tri_incoherence |

15 | 1 | Robert Porritt | |

16 | 1 | Robert Porritt | INSTALL: |

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 | @cd STATION_ANALYSIS_TOOLS@ |

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 | @./configure@ |

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 | @make |

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 | EXECUTABLE DISCRIPTIONS: |

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 | RUN: |

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 | sac_file_name |

67 | 1 | Robert Porritt | pole_zero_file_name |

68 | 1 | Robert Porritt | n_sec n_skip unit_flag |

69 | 1 | Robert Porritt | output_file_name |

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 | EXAMPLE: |

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 | KNOWN BUGS: |

90 | 1 | Robert Porritt | 3.29.2011 |

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 | 3.25.2011 |

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 | BE AWARE: |

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. |