32.1 Spec program

The SPEC program is used for making spectra of many seismic signals in a semiautomatic manner. It can be used for several investigations:

A: Making a large series of signal spectra, which can be corrected for instrument and path.

Average spectra are calculated. There are two options for further processing the calculated spectra:
Option (1): Calculate acceleration density spectra which are plotted compared to the Peterson noise model.
Option(2): Using the slope of the flat part of the displacement spectra to calculate the near surface attenuation kappa (see 8.12.

B: Making relative spectra of seismic events or background noise in order to determine the soil response. When using relative spectra of horizontal versus vertical components, this is referred to as the Nakamura method (Nakamura, 1989).

C: Making relative spectra of signals from two stations in order to determine Q. The program makes output files for generating GMT plots in addition to standard SEISAN plots.

Note: Parameter file has changed between SEISAN 7.2 and 8.0 (number of windows and overlap has been added).

The program can technically operate in two ways: (1) Making relative spectra of a series of pairs of stations terminated by the average spectra, (2) Making a series of spectra for a number of stations and events. The spectra can be corrected for distance, Q, and instrument response. In addition, the spectral levels can be expressed in moment or moment magnitude calculated in the same way and with the same units as in MULPLT. All relevant parameters are taken from the CAT files, the CAL files and the input parameter file for SPEC. Window selection for the spectra can be specified to be related to the P, S arrival times or the earthquake origin time and it is thus possible to automatically make e.g. S-wave spectra of a large set of stations and events. Optionally, noise spectra, can be calculated together with the signal spectra. The noise window is selected at the start of the waveform file.
Before the program is started up, the input files must be prepared. The program need two input files. The parameter file (default spec.par) gives the parameters to use and the list of stations to process. The event file (default spec.inp) is a CAT file with events to use or a filenr.lis type file with waveform file names (can only be used if no readings are needed, like for Nakamura studies). An example of a spec.par and spec.inp file is found in DAT. These files can be used immediately with the test data set.

The program produces several output files. The main output is in spec.out with the parameters used, the station event combinations used and error messages. The other files are giving output of most graphs shown. These ASCII output files can be used in other plotting programs, however they have been specifically formatted for the SEISAN GMTXY plotting script. Note that the numerical values given for the spectral output given in those files is the same that appear on the plots and the values are linear. So if a spctrum has been instrument corrected and smoothed, that is what is given. Similarly if a relative spectrum is used. The number of files depends on number of stations used. Examples of files could be

spec_all_ASK__S__Z.out All spectra from ASK, S Z
spec_all_BER__S__Z.out All spectra from BER, S Z
spec_all_gmt.out All spectra from ASK and BER
spec_ave_ASK__S__Z.out Average spectrum from ASK, S Z
spec_ave_BER__S__Z.out Average spectrum from BER, Z Z

In order to plot these files with GMTXY (only Unix), give e.g. command

gmtxy spec_all_ASK__S__Z.out
There is one more output file, spec_amp.out, which gives the log log spectra of of all traces calculated (possibly instrument corrected) before smoothing takes place.
Limitations of amount of data: The program is set up to handle 300 spectra of up to 30000 points each for one run. The dimensions can be increased in spec.for, however the program must then be recompiled. The spectral windows are 10% tapered. The analyzed signals will be checked for clipping and rejected if clipped. A message is then given in spec.out

The program is tarted with command spec. it can also be started with arguments:

-par file: give parameter file

-cat file: give cat file

-batch: no plot

if no par and cat files are given, the default names are assumed so only the argument -batch can be used.

The spec.par file: The file contains alternate lines of parameter names and parameter values, and must contain the number of lines shown in the example below.

selection criteria 1: P, 2: S, 3: S from P, 4: abs
window length, # of windows,overlap
number of times to smooth
gain factor of channel 1
noise spectrum 0: n 1: y
make relative spectras 1: y, 0; n
plot pics
frequency band to use
1.0 7.0 
response removal: 0: none, 1: displ., 2: vel., 3: accel. 4. noise pow. 5. kappa
rotate: 0: no, 1: yes
q0, qalpha and kappa
distance correction
minimum correlation and minimum sn for kappa
0.5  2.0
velocity and density
0.0 0.0
magnitude spectrum
stations and components, format a5,1x,a4,1x,a5,1x,a4
FOO   S  Z SUE   S  Z

The parameters:
Selection criteria: Determines how the start of the time window is selected. 1: Start with the P-arrival time, 2: Start with the S-arrival time, 3: Start with the S-arrival time calculated from the P-arrival time assuming a P to S velocity ratio of 1.78, 4: Start with 'start' (see next parameter) seconds after the origin time as given in the CAT file header. This option can be used if no readings are available in the CAT file. When using a P or S-time for start of window, the program uses the first P or S phase found in the CAT file for a given station. Component is of no importance here, so there is only a need for e.g. one P-time for the station being processed if 3 component data is used. This is also the case when rotating the signal, see below. However, on the trace plots, only readings on those components shown will be seen on the plots.

Start: If the selection criterion is 1,2 or 3, this is the number of P or S travel times (from the origin) used to find start time of window. Use 1.0 if the window shall start exactly at the phase time picked. If selection criteria is 4, start is the number of seconds after the origin time. Optionally, give a second input, number of seconds to start before the window.

Window length, #of windows, overlap:

- Window length: Window length in secs for both signal and noise (if selected).
- # of windows: If more than 1, spectra will be made in several windows following the first window and average spectra will be made. This option can only be used if selection criteria is 4. Used for noise studies or Nakamura studies.
- Overlap: Windows can overlap (factor $<$ 1.0) exactly follow each other (factor=1.0) or have gaps (factor $>$ 1.0). E.g. 0.9 is equal to 10 % overlap.

Number of times to smooth: Number of times to smooth, 0 means no smoothing.

Gain factor of channel 1: Factor that the spectral level for channel 1 is multiplied with. This can be used if the response shape is the same for the two channels and only the levels are different. If the shape is also different, set factor to 1 and use response removal below.

Noise spectrum: If 0, no noise spectrum, if 1, make noise spectrum. The noise window is taken from the beginning of the trace and the window length is the same as given above.

Make relative spectra: If zero, no relative spectra, if 1, make relative spectra. The relative spectra will appear one on each page, and the average relative spectrum on the last plot (see Figure 32.1 ). If no relative spectra are chosen, only one trace and one spectrum is shown per page and the average spectrum is shown on the final plot. MUST BE SET to 1 to calculate Q, see below.

Plot pics: If 1, the phase pics in the CAT file spec.inp will be plotted.

Frequency band used: Lower and upper frequency bands for the spectral plots.

Response removal:

If 0, no response is removed.

1: displacement

2: velocity

3: acceleration (units is nm, nm/s and nm/s*s)

4: Power spectral density in dB relative to ((1m/s**2)**2)/Hz. This option is used for seismic background noise studies.

5: Determine kappa. The flat part of the displacement spectrum (frequency below corner frequency) is approximated by a straight line and kappa calculated for each event and the average at the end (in spec.out file and on final plot). The spectrum will normally be corrected for Q, BUT NOT kappa. For more details, see (Havskov and Ottem├Âller, 2010).

6: Determine kappa. The flat part of the acceleration spectrum (frequency above corner frequency) is approximated by a straight line and kappa calculated for each event and the average at the end (in spec.out file and on final plot). The spectrum will normally be corrected for Q, BUT NOT kappa.

Make sure to set appropriate frequency limits and correct distance corrections. Can be used for both P and S-spectra.

A cal file for each channel must be available in the CAL directory (see section 4.6). For relative spectra, the response removal has no importance if the response is the same for the channels compared. A simple correction can be made with "Gain factor of channel 1" parameter above. NOTE: If moment or magnitude spectrum is made, response removal MUST be 1.

Rotate components: If 1, the horizontal components are rotated. This means that if the user has specified N or E, radial or transverse respectively will be used instead. The original data remain unchanged. If start time of spectra are chosen by using P or S, there must be a reading from those components if the pics are to be plotted. If the parameter is zero, no rotation is done. See also MULPLT for more details of rotation.

Q0, qalpha and kappa:

Parameters in Q-relation Q = Q0**qalpha used for spectral correction (see also section on MULPLT for standard attenuation relations). Only used if response is removed. If first 2 parameters are 0,0, no Q-correction. New from SEISAN7.2 is that a kappa correction also can be used (see MULPLT spectral section). Note that from version 10.1, Q=Q0 for f < 1 Hz.

Calculation of Q:
If Q0 and qalpha is set to -1,0, the relative spectra will be used to calculate q as a function of f (see standard relations in MULPLT section) and the plots will show q as a function of f. This can be used for both P-waves and S-waves. The distance correction MUST be set, S-velocity must be given (see below) and it is recommended to assume body wave spreading (amplitude proportional to 1/distance, factor is 1.0 below). If the response of the 2 stations is not identical, correction for response must also be made. There must be an origin time and phase readings for components used must be available in order to calculate Q. Q is calculated as pi * f * (t2-t1) / (ln(A2(f)/A1(f)) + alpha * ln(t2/t2)), where A1 and A2 are spectral levels at frequency f for the two stations, t1 and t2 are travel times and alpha is geometrical spreading exponent (1.0 is body wave spreading). Q values lower than 1 and higher than 5000 are not used, the Q(f) plot might then display a long straight line. The Q=Q0*f**qalpha is calculated from the 'good' values'.

Distance correction alpha: The spectral amplitudes are multiplied by R**(distance correction) if different from zero. This option MUST be set if moment or moment magnitude options (see below) are selected as well as calculation of Q. However, it can be used without instrumental correction. For body waves, use 1. Note that the geometrical spreading use here is simpler than used in MULPLT.

Minimum correlation coefficient and minimum signal to noise ratio for kappa and optional maximum hypocentral distance: The minimum correlation coefficient and signal to noise ratio for an event to be included in average kappa. The coefficients are from the linear fit to the flat part of the spectrum. The optional maximum hypocentral distance can be used to limit the data to e.g. 50 km. if used here, then the data in the input file does not have to be distance limitied.

Velocity and density: Velocity (km/sec) and density (g/cm*cm) used for calculating moment spectra. If set to 0,0, no moment spectra are calculated. See section on MULPLT for details of calculation.

Magnitude spectrum: If 1, the spectral level is converted to moment magnitude, see MULPLT for details of calculation.

Stations and components: Station-component pairs used, one pair per line, format (a5,1x,a4,1x,a5,1x,a4). If no relative spectrum is used, the first station-component on the line is used.

NOTE: Note for component, only the orientation need to be used like Z instead of S Z. This can be useful if the data contains different component codes for same data set like HH Z and S Z.

Calculation is done separately for each station so the the program stops between each and the average plots is shown unless option -batch is used.

Averaging in spec:
Q: For each frequency, the average linear 1/Q and corresponding sd is calculated. The upper and lower bounds are calculated by subtracting and adding the sd. These values are then converted back to Q and finally the log is taken. Only the `good' individual values are used. There is a possibility that the lower bound becomes negative. In that case, the log Q is set to zero. Because the average is made in 1/Q, the upper and lower bounds curves will not be symmetric around the average Q-curve.

Power spectrum: For each frequency the dB values are averaged and upper and lower curves should be symmetric.

Kappa: Same as for Power spectrum.

Other spectra: The linear spectra or relative spectra are averaged. The sd used in the log spectra are calculated by subtracting log average spectrum from log(average spectrum + sd).

If many data sets are to be processed at the same time one can have one CAT file for each station or station pairs and optionally also a location for the corresponding waveform files. On the station line, the CAT file name is then given starting from column 25 and the directory of the waveform files starting at column 105. In the example below, the names of the CAT files all have the same name but they are in different directories but they could also have had different names and placed in the working directory. The specification of the waveform directory is needed if the waveform files are not in the working directory or the in WAV.

If cat files are given on the channel line, there is no question later of name for the cat file.



STAT  COMP STAT COMP    Abs location CAT file, start col 25                                             Directory for WAV files, start col 105
KTK1     Z              /home/seismo/WOR/jens/kappa/KTK1/spec.inp.q                                     /home/seismo/WOR/jens/kappa/KTK1
KTK1     N              /home/seismo/WOR/jens/kappa/KTK1/spec.inp.q                                     /home/seismo/WOR/jens/kappa/KTK1
KTK1     E              /home/seismo/WOR/jens/kappa/KTK1/spec.inp.q                                     /home/seismo/WOR/jens/kappa/KTK1
ARA0     Z              /home/seismo/WOR/jens/kappa/ARA0/spec.inp.q                                     /home/seismo/WOR/jens/kappa/ARA0
ARA0     N              /home/seismo/WOR/jens/kappa/ARA0/spec.inp.q                                     /home/seismo/WOR/jens/kappa/ARA0
ARA0     E              /home/seismo/WOR/jens/kappa/ARA0/spec.inp.q                                     /home/seismo/WOR/jens/kappa/ARA0
SKAR     Z              /home/seismo/WOR/jens/kappa/SKAR/spec.inp.q                                     /home/seismo/WOR/jens/kappa/SKAR
SKAR     N              /home/seismo/WOR/jens/kappa/SKAR/spec.inp.q                                     /home/seismo/WOR/jens/kappa/SKAR
SKAR     E              /home/seismo/WOR/jens/kappa/SKAR/spec.inp.q                                     /home/seismo/WOR/jens/kappa/SKAR
LOF      Z              /home/seismo/WOR/jens/kappa/LOF/spec.inp.q                                      /home/seismo/WOR/jens/kappa/LOF
LOF      N              /home/seismo/WOR/jens/kappa/LOF/spec.inp.q                                      /home/seismo/WOR/jens/kappa/LOF
LOF      E              /home/seismo/WOR/jens/kappa/LOF/spec.inp.q                                      /home/seismo/WOR/jens/kappa/LOF

In this example, it was desired to calculate kappa for many different stations and since only stations in a given distance were used, a separate CAT files was needed for each station. Since a complete specification is given to where the waveform files are located, there is no requirement of having the waveform data in the the currem directory, WAV or the defualt WAV data base.

If all data is in the default data base, a simpler way to calcualte kappa for many stations is setting a distance limit in spec.par

minimum correlation, minimum signal to noise ratio, max distance


and having a spec.inp with the station in many different distances. Then the channel part of the par file could be:

KTK1     Z
KTK1     N
KTK1     E
ARA0     Z
ARA0     N
ARA0     E
LOF      Z
LOF      N
LOF      E
STEI     Z
STEI     N

Normally, at the end of calculation, there is a plot of all the results. If calculation for many stations, the plot can be skipped if spec is started with option -batch. There is a summary output file for the kappa value for each station and component called spec_kappa_average.out. The content is station, component, kappa, sd of kappa and number of values.

KTK1    Z   0.031  0.013   31
KTK1    N   0.037  0.014   38
KTK1    E   0.025  0.013   43
ARA0    Z   0.001  0.015   16
ARA0    N   0.009  0.022   15
ARA0    E   0.014  0.016   15
SKAR    Z   0.040  0.014    2
SKAR    N   0.025  0.003    3
SKAR    E   0.027  0.003    3
LOF     Z   0.038  0.018   22
LOF     N   0.035  0.022   26
LOF     E   0.033  0.018   26

For kappa, there are also ouput files in Nordic format with average kappa at a particular station. The epicenter is the station coordinates and the magnitude is proportinal with kappa. The file can be plotted with map or mapg and also used for furhter ananlysis, see example next section.

Running the program:
The program gets the first pair of stations (or one station) from spec.par, calculates the spectra using the list of events in spec.inp and at the end of the station list, calculates the average spectral ratios for all pairs (max 100). All spectra are then shown on one plot together with averages and standard deviation. Then the next pair of stations is processed in the same way and the program continues until the end of file spec.par. Each pair of stations with signals and spectra is plotted on one page. If no relative spectra are made, the plots look similar except that only one station is shown. Hard copy plots are made for each page and sent to the printer if specified (see below). The hard copy postscript file is called spec.eps and when the program finishes, a file with the last plot is available on the disk. For each spectrum (relative or single), the average spectrum (or Q) is calculated both as an average of the log spectrum and as an average of the linear spectrum. There is no frequency weighting and since all values shown on the plot are used, the average value will be more representative of the high frequency part of the spectrum since there are more values. This can be regulated by choosing another frequency range. The average spectra shown on the last plot are log-averages. If option to calculate Q is used, the plots show 1/Q as a function of frequency instead of relative spectra (proportional to relative spectra). For each event, Q0 and qalpha are calculated.

When calculating kappa, the average spectrum do not have much physical meaning since the averages are made from absolute spectra of events that might have very different moments. So the kappa calculated from the average spectrum is not to be used.

Interactive output of level and frequency: With a spectral ratio (or Q) plot on the screen, position the cursor at the point of interest on the spectrum and click. The level and frequency will now be displayed on the right side of the plot.

The output file spec.out gives details of the run like averages and missing data. The output file spec_ave.out gives the x and y-values of the average spectrum IF IT HAS BEEN PLOTTED ON THE SCREEN. File spec_rel.out gives the values of the relative spectra.

There are 4 interactive input options:

0: All spectra are calculated but not sent to the plotter or screen except the last plot with the average spectra (sent to both screen and printer). Used for checking the files or making a final run. If no relative spectrum is chosen, no final plot is made. For each station and event combination, check lines are written out on the screen.

1: All plots are shown on the screen, but not sent to the laser printer.

2: All plots are shown on the screen and at the same time sent to the laser printer.

3: No plots are shown on the screen, all are sent to the laser printer. For each station event combination, check lines are written to the screen.

4: Only final plot on laser.

5: No plotting so no graphics window come up anytime. Can be used e.g. for batch mode.

How to run the program with only waveform files available: Two options:

(A) Using S-files
Step 1: Generate S-files in your local directory with AUTOREG,
Step 2: Make the spec.inp file with COLLECT.

(B) Using filenr.lis
Step 1 : Make a dirf of waveform files to use
Then use filenr.lis as input file name

With only waveform files and no readings in the spec.inp file, it is only possible to use option 4 (absolute time) for start criteria. Since the events have not been located, the "origin time" read from the S-files will be identical to the waveform file start time, so the parameter "start" can then be set to number of seconds after waveform file start time. Figure 32.1 shows an example.

Figure 32.1: An example of using the SPEC program. On top the original traces are shown with windows chosen, in the middle the spectra of each channel and at the bottom, the relative spectrum. Lower right shows the input parameters used. In some cases (kappa and Q) the values calculated for this case are also shown.

Figure 32.2: An example of a GMT plot. The figure shows an example of making noise spectra of several traces.