In today's activity we will revisit some ideas about frequency content and filtering using some real seismic data. To begin with, we need to examine how continuous ground vibrations are digitized with a particular sampling rate when they are recorded into computer memory. To begin this activity, let's login to the linux server with our VNC viewer, and then move into the directory called
[username@linux username] cd source
[username@linux source] sac
SAC> r *SUR*SAC
SAC> p
You should notice that the data shows the earthquake prominently, but there is also some "noise" that is present before the earthquake begins, making it more difficult to see the initiation of the earthquake on the seismogram. In this activity we will work to remove this noise, by learning how the seismgoram is sampled to figure out the range of recorded frequencies, what the frequency content of the earthquake is, and then what the frequency content of the noise is. With all of that information, we will able to devise a filter to remove the noise and preserve the earthquake signal.
We will start with learning how the seismogram is sampled. Since sac data files have a variety of information stored in a "header" at the beginning of the file, we can use the LISTHDR command to see this information, which can be abbreviated to lh. Let's run this command on the data file you just read into sac. Note that since there is a pretty long list of information about this data file, you will need to hit the enter key when sac gives you the "Waiting" prompt.
SAC> lh
The first six header variables listed in this output will help you to see how this data file was recorded:
NPTS = number of point
B = begin time of data
E = end time of data
IFTYPE = type of data
LEVEN = whether the data is evenly spaced in time
DELTA = the time difference between each sample (delta-t)
Using the DELTA value, calculate what the sampling rate is (number of samples per second). Then calculate the maximum frequency one can observed in the data. Store the answers to these questions in a file called sample.txt.
Next we will examine the frequency content of this earthquake recording. Recall from our earlier SAC tutorial on frequency content that we use FFT command to perform the fast-fourier transform. Then the time series will be converted to the frequency domain, so the regular PLOT command will not work. Instead, we will use the PLOTSP to plot the frequency spectrum (abbreviation PSP). For this activity, we can use the AM and LOGLOG options to just plot the amplitude part of the frequency component and do a logarithmic X-axis and Y-axis.
SAC> fft
SAC> psp am loglog
What is the general trend in the amplitude vs. frequency plot. At what frequencies do we see the largest amplitudes? How might this result differ from a smaller earthquake? Store the answers to these questions in a file called eq-freq.txt.
Our next task is to figure our the frequency content of the noise that is recorded before the earthquake on the seismogram. To isolate this part of the signal, we will use the CUT command to restrict the data to 30 seconds before the earthquake (T1 marks the first arrival of the earthquake waves) and then re-read the seismogram.
SAC> cut t1 -30 0
SAC> r *SUR*SAC
SAC> p
You may notice that since we've shortened the seismogram, there is now a downward trend to the data. This can unfairly influence our frequency determination, so we can use the RTREND command to remove this trend.
SAC> rtrend
SAC> p
Looking at the waves in the time domain of this seismogram, make an estimate of the dominant frequency of the noise. You may find it easiest to measure the period from peak to peak, and then convert the period to a freqency. After making a time domain estimate, convert the time series to the frequency domain using FFT, and then plot the results with PSP. At what frequencies do you see the largest amplitudes? How do these two estimates of the dominant noise frequency compare? Store the answers to these questions in a file called noise-freq.txt.
Now that you have an idea what the frequency content of the earthquake and the noise are, we should be able to construct a filter that can remove the noise and preserve the earthquake signal. We will use the BANDPASS command to perform the filtering, which can be abbreviated bp. We then need to specify the frequency corners for this filter, which means the beginning and the end of the frequency band that we want to pass through the filter. Frequencies less the first corner and greater than the second corner will be removed. Here are the commands you should use to turn off the CUT command (we want the whole seismogram now), re-read the seismogram,apply the filter (but you need to specify numbers for FREQ1 and FREQ2), and then read the original seismogram back in for comparison. You can try several filters if you repeat this sequence of commands with different FREQ1 and FREQ2 values.
SAC> cut off
SAC> r *SUR*SAC
SAC> bp corner FREQ1 FREQ2
SAC> r more *SUR*SAC
SAC> p1
What are the FREQ1 and FREQ2 that best remove the noise and preserve the earthquake signal? Why did you choose thos values? How does it change the original seismogram? Store the answers to these questions in a file called noise-filt.txt.
lh |
sac command that lists the header variables for all files read into sac. |
rtrend |
sac command that removes the linear trend for all files read into sac. |
brudzimr@muohio.edu, 30th October 2007