Seismic Analysis Code 4

The final two labs on SAC will cover how to do two common types of signal processing in SAC. This lab focuses on stacking, a common method of reducing noise in signals. SAC provides many handy tools for stacking seismic waveforms.

Stacking

Stacking involves adding signals together so that the signals constructively interfere, while the noise adds incoherently, improving the signal to noise ratio. Seismologists often use stacking when recorded waveforms are relatively weak so that the coherent signal adds, while the random noise cancels out. Stacking can be done in a rather rudimentary fashion, or using a more sophisticated set of commands in SAC known as the “signal stacking subprocess.”

What is a subprocess? A subprocess is like a module in Python – it adds some extended functionality to SAC, but we have to tell SAC that we intend to use it. When you run a subprocess, you have additional commands available that are not normally part of SAC, though you have a restricted set of basic SAC commands, usually just commands that can be used to manipulate data. To learn more about the commands available in a particular subprocess, consult the SAC manual.

Basic Stacking

To perform a basic stacking operation in SAC, we need the ADDF command. ADDF, short for “add file,” is used to add together a signal in memory with a signal loaded from a file. The two signals must have the same number of data points and the same sampling rate, or SAC will complain and not let you add them together. If you aren’t sure whether or not your signals can be added, check the header files.

addf works in a different manner depending on the number of files in memory and the number of files listed in the addf command. If the number of files in memory matches the list of files provided with addf, then the first file listed in the addf command will be added to the first file in memory, the second file in addf will be added to the second file in memory, etc. So for instance,:

read file1.sac file2.sac
addf file3.sac file4.sac

will add file1.sac and file3.sac and place it as the first signal in memory, and add file2.sac and file4.sac and place it as the second signal in memory.

If the number of files in addf is less than the number of files in memory, then the last file provided with addf will be added to all remaining files in memory that do not match a file in addf. Thus, I can add a single file to many different files using:

read file1.sac file2.sac file3.sac
addf file4.sac

Exercise: Use the addf function to combine an impulse function with a boxcar function. Try changing the number of points and time spacing on the functions to verify that you need two signals with the same time sampling to use addf.

Once we have addf working, it is simple to stack many files: we simply need to repeatedly call addf with all of the files that we wish to stack. This can be done most easily using a macro with a do loop. So for instance, if I have a series of files that I want to stack, I could stack them as follows:

read station1.sac
do file list station2.sac station3.sac station4.sac station5.sac
    addf $file
enddo
plot

This will first read in the data for the first station, then loop over the remaining stations and stack them on the signal already in memory. If I have a large number of files with complicated names and did not want to type all of them in, I could use a wildcard loop to read them in:

# first generate a signal of all zeros using funcgen
funcgen line 0 0
do file wild station[1-5].sac
    addf $file
enddo
plot

The first line is to initialize the signal in memory as all zeros before stacking the first file on top of that; this is simply to avoid the awkward case of reading in the first file but then excluding it from the loop. Note that this assumes that each station has time data that matches the default for funcgen. If the default signal length and spacing differs from the signals to be stacked, you will need to specify npts and delta to match your data when generating the initial string of zeros. Also note that addf does not accept wildcards, so if you want to use a wildcard to stack a number of files, you will need to use a do loop to accomplish this.

Exercise: I have provided a series of 100 noisy signals (“stack.tar.gz”) to practice stacking, which are available on the course website and in my public folder. Read in a few of them to verify that they appear quite noisy – while you can see the strongest arrivals at around 8 seconds, the rest of the signal is not particularly useful (see one example below). Write a SAC macro to stack them, plotting the result (see the second plot below for the stacked version). Try only using a subset of them to see how the signal improves as you add more data.

A noisy time signal. There are 100 files like this, where noise has been added that drowns out much of the signal.

The stacked version signal, which combines 100 noisy waveforms. This signal a much better signal to noise ratio, making the coda decay much clearer.

Exercise: Write a SAC macro using addf to stack a given number of sine functions that have had high amplitude noise added to them (the number of records to stack should be input to the macro as an argument). Your macro should first use a do loop to generate the records using funcgen sine and funcgen random, add them together using addf, and write to disk. Then, use a do loop to read the files and stack them, divide by the number of records (use div <number> to divide the signal in memory by a constant number), and plot the data.

Note: Random number generators do not produce truly random numbers. Instead, they produce a deterministic sequence of numbers that follow the same distribution as a set of random numbers. This means that calling fg random repeatedly will always produce the same sequence of numbers, which is not what we really want here (note however that this gives us a way to reproduce any calculation using random numbers, so this is in many ways a good thing). To produce a different sequence of numbers each time, we need to provide a different “seed” to the random number generator, which will ensure that the deterministic sequence starts at a different place. This can be done in SAC by calling funcgen random 1 <seed>, where <seed> is an integer that should be different each time through the loop.

Signal Stacking Subprocess

You can do more complicated things with stacking using the Stacking Subprocess. This gives you more sophisticated tools for aligning signals with different origin times so that they will produce coherent signals when stacked. We will not cover these time alignment tools in this class, but once you know the basics of how to stack signals using the subprocess it is not difficult to figure them out. To start a signal stacking subprocess, enter sss into SAC. Once you enter a subprocess, you will not have access to the full extent of SAC commands, so see the manual to know what commands are allowed.

Once you have started the stacking subprocess, you can stack waveforms in the same way as before, the only difference is the commands that you will use. Useful commands include:

• First, to do a stack you are required to define the time window of interest by entering timewindow <start> <end>, which defines the window to be used for the stack. Only time values within this window will be read into memory, when you add additional files to be stacked.
• To zero out the stack, enter zerostack to start with a clean signal.
• Once the time window is defined and the stack is zeroed, you can use addstack <file> to add a file to be stacked. This is the equivalent of addf within the stacking subprocess, but is different in that the original file in memory is not overwritten. This allows you to visually inspect the signals before performing the stacking operation.
• Once you have added files to be stacked, you can plot them by entering plotstack. This will show all of the files that have been added to the stack with their time axes aligned. If you have chosen a particular velocity model to align wave arrivals, the time axis will account for this correction, but since we are doing a simple stack here, it will plot all of the signals over the same time range.
• When you have added all of the desired files, stack them by entering sumstack. To see all of the signals, plus the stacked signal, you can enter plotstack to see the results.
• To write the stacked signal to disk, use writestack <filename>.
• Finally, you can exit the stacking subprocess by entering quitsub into SAC, which will return to normal SAC functionality.

This does essentially the same thing as the simple procedure outlined above for stacking files. However, the additional tools for defining time windows, the simple way for zeroing out the data in memory, and the ability to see the signals separately along with the stacked sum make the subprocess a better way to carry out stacking operations.

Exercise: Modify the macro you wrote above to stack the noisy signals that I have provided to use the signal stacking subprocess. Write your result to disk so that you can retrieve the stacked result in the main part of SAC.