Sheila Peacock
Dear Brian and Arthur and all,
<quote>
Sheila,
I guess there are a couple of options here. I think the easiest is
the use the freqlimits option in transfer to limit where the transfer
is applied. The other possibility is to implement an option for the
FAP transfer to do what you have done and extrapolate the FAP outside
it's prescribed range.
I think Arthur had some rationale for making the decision we did about
the FAP code and not extrapolating, but I cannot remember what it was.
Brian
__________________
Sheila and others,
When writing fap.c, we decided not to do any extrapolation beyond the given limits.
The idea is that those would be handled by other arguments on the TRANSFER line --
just as when doing a polezero deconvolution. (See the TRANSFER help file at
http://www.iris.edu/software/sac/commands/transfer.html.)
Applying an LP response directly to an impulse is asking for trouble. Abrupt
changes as one sees in that waveform mean very high frequencies that must be
gotten rid of before applying the response. When I did the following:
SAC> fg impulse
SAC> fft
DC level after DFT is 1
SAC> plotsp am
the amplitude was constant for all frequencies. One cannot apply an
instrument correction to a signal that does not go to zero at the Nyquist.
I have also worked on the interpolation routine in SAC, and the methods of
interpolation do not work well with such abrupt changes. SAC routines are
meant to be applied only to real data.
If I have misunderstood your points, please let me know.
Arthur
<unquote>
Thank you very much for your replies. I think we might have to agree
to differ here, but I'm sure we agree to deprecate the use of FAP files
where poles-and-zeros would allow sensible extrapolation. Here is
the rest of my train of thought, in the hope that it will be useful.
My problem arose from trying to duplicate the procedure used by
the CTBTO International Data Centre (IDC) to derive surface-wave
magnitudes. All seismograms there are converted to the response of a
KS36000 LP seismometer, the response for which is a bandlimited FAP
between 0.001 and 0.3 Hz (I sent it to you). I naturally
inspected its impulse response, by convolving it with "funcgen impulse..."
with the data-collection parameters I'd expect from an IDC broadband
installation, though severely at 100-Hz sampling rate (many IDC broadband
seismometers are sampled at 40 Hz, but lately we have been installing
100-Hz systems at our own sites). This was how I got
the alarming spike. I agree that a real seismogram shouldn't
look like an impulse, and I would probably have frequency-limited the
decon of the original instrument if I'd been doing it for real.
A purist approach to bandlimited FAP files like the KS36000LP
would be to avoid any extrapolation of the FAP at all. For the
KS36000LP, it would be necessary to decimate the
real seismogram to 0.6 Hz (with appropriate anti-aliasing) before
convolving it with the KS36000 LP FAP. If the response of the
original instrument is also a bandlimited FAP, like the other
example I sent you (top frequency 10 Hz) then a 40- or 100-Hz seismogram has
to be decimated before deconning that FAP, then decimated again
before convolving the KS36000 FAP, which prevents you from doing a
single "transfer".
You can avoid the decimation by using "freqlimits" in the "transfer" to
crack down on frequencies outside the band of the FAP, but for
bulk processing using the motley collection of FAPs provided
by the IDC, you would have to maintain a motley collection of "freqlimits".
So there is a case for tapering off the edges of FAP tables as a
substitute for unifying the frequency ranges of the world's variety
of FAP files.
The present SAC FAP routine is somewhere in between the last two paragraphs.
It doesn't force (or persuade with a warning message) the user to
decimate or filter the seismogram to fit the FAP bounds, and it doesn't try to
make the FAP file conform to the seismogram by extrapolating it
off to zero at/before the seismogram's Nyquist. So the user has to
be more intelligent than I was being at avoiding having any meaningful
signal outside the exact frequency bounds of the FAP table.
Best wishes,
Sheila Peacock.
-
Sheila,
Before I get to your question/problem, I have a "list" question: Did you
get a copy of the messages you sent to sac-dev? I ask because I did not
get a copy of the message I sent, but I see that my message is listed in
the archives. The sac-dev listserv is used much less than sac-help, so I
cannot remember if it has different rules. I will see if I get a copy of
this message.
I am not sure I follow your comments entirely. Specifically, I am not
sure on which issues "we might have to agree to differ" as I do not see
much disagreement between us on the basics -- except perhaps you might
vote for adding an option to fap that would allow it to extrapolate
outside the prescribed frequency range.
You are trying to duplicate results done previously with deconvolution
using FAP files. How did the earlier procedures handle extrapolation? If
you can find that out, I would suggest you write a program outside of SAC
that uses that procedure and can handle bulk processing. If you want to
use SAC, you may be able to write a script that can deal with your "motley
collection." By the way, the default for decimate is to have the
antialiasing turned on, so there would be one less thing to worry about.
Again, if I have missed something let me know.
Arthur
-
Dear Arthur,
Thanks for your email. To answer your list question first: no, I didn't
get a copy of my email from the list, but I had fortunately CCd it to
myself anyway and I got that. I did that because my first email to
the list on this matter, the one with the attached diagrams, was
held for moderation because it was too large. I assume that my latest
one was small enough not to fall into that trap.
As for the FAP decon, I have some code that I had bolted on to SAC-101.3
for doing FAP decon, which extrapolated the FAP table linearly from
the first two and the final two points. It came from a proprietary source
so I was grateful to find a FAP reader in SAC-101.4 that would allow me
to install a no-privacy-worries copy of SAC.
On being caught out by that spike, I doctored your new code in 101.4 to
do the same linear extrapolation as my proprietary code had done
(I also put in a few extra checks for log(0) and division by zero).
The doctored code shows a foolish user like me who convolves a 0.3-Hz
bandlimited FAP with an untreated 100-Hz spike more-or-less what they
think they want to see. With the flat continuation in the original
101.4 code rather than my linear extrapolation, the foolish
user gets the spike, or the Gibbs' phenomenon in the other example I sent.
The user should then see that some operation on the FAP file or the
seismogram prior to or during "transfer" is called for - if the user
is not too foolish.
(My code gives an impulse response indistinguishable from the original
SAC-101.4 code when given a sensible FAP file that goes down to zero
of its own accord within its bandwidth - "dinosaur-shaped spectrum", as
I recall the description in a signal-processing book, referring to
the textbook silhouette of a diplodocus.)
My colleagues aren't foolish, but they are very keen on backward
compatibility (and hence wary of upgrades). They had various
things running in SAC-2000, which need to go on working. They had got
used to the FAP code I'd attached to SAC-101.3, so my aim was to make
the 101.4 FAP "transfer" work as similarly to that one as possible.
I hope that this background is some help.
Regards,
Sheila Peacock.
-
Sheila,
We, along with IRIS, are attaching metadata to the polezero files. We might be able to do such a thing for the FAP files as well. Correct me if I am wrong, but it looks like you would like to associate a freqlimits with the FAP files. It that is the case we can add in a comment that tells SAC to extrapolate the response or truncate using freqlimits. See possible example below. It would impose more control over how the response is applied/removed from your data.
Thoughts ?
Brian
* First one (LP seismo)
* Freqlimits : 0.011 0.01 10.0 11.0
* Extrapolate Response : NO
* Comment Here
* Another Comment Here
1.00000E-02 0.394768 282.522
2.00000E-02 3.11046 146.286
3.00000E-02 6.05993 53.0681
4.00000E-02 7.18071 -21.5560
5.00000E-02 6.35654 -82.4200
6.00000E-02 4.83388 -130.514
......
On Jun 27, 2011, at 10:57 AM, Sheila Peacock wrote:
Dear Arthur,
Thanks for your email. To answer your list question first: no, I didn't
get a copy of my email from the list, but I had fortunately CCd it to
myself anyway and I got that. I did that because my first email to
the list on this matter, the one with the attached diagrams, was
held for moderation because it was too large. I assume that my latest
one was small enough not to fall into that trap.
As for the FAP decon, I have some code that I had bolted on to SAC-101.3
for doing FAP decon, which extrapolated the FAP table linearly from
the first two and the final two points. It came from a proprietary source
so I was grateful to find a FAP reader in SAC-101.4 that would allow me
to install a no-privacy-worries copy of SAC.
On being caught out by that spike, I doctored your new code in 101.4 to
do the same linear extrapolation as my proprietary code had done
(I also put in a few extra checks for log(0) and division by zero).
The doctored code shows a foolish user like me who convolves a 0.3-Hz
bandlimited FAP with an untreated 100-Hz spike more-or-less what they
think they want to see. With the flat continuation in the original
101.4 code rather than my linear extrapolation, the foolish
user gets the spike, or the Gibbs' phenomenon in the other example I sent.
The user should then see that some operation on the FAP file or the
seismogram prior to or during "transfer" is called for - if the user
is not too foolish.
(My code gives an impulse response indistinguishable from the original
SAC-101.4 code when given a sensible FAP file that goes down to zero
of its own accord within its bandwidth - "dinosaur-shaped spectrum", as
I recall the description in a signal-processing book, referring to
the textbook silhouette of a diplodocus.)
My colleagues aren't foolish, but they are very keen on backward
compatibility (and hence wary of upgrades). They had various
things running in SAC-2000, which need to go on working. They had got
used to the FAP code I'd attached to SAC-101.3, so my aim was to make
the 101.4 FAP "transfer" work as similarly to that one as possible.
I hope that this background is some help.
Regards,
Sheila Peacock.
_______________________________________________
sac-dev mailing list
sac-dev<at>iris.washington.edu
http://www.iris.washington.edu/mailman/listinfo/sac-dev
-
-
Arthur Snoke
Brian Savage
