SUGABOR - Outputs a time-frequency representation of seismic data via the Gabor transform-like multifilter analysis technique presented by Dziewonski, Bloch and Landisman, 1969.

sugabor <stdin >stdout [optional parameters]

Required parameters:
if dt is not set in header, then dt is mandatory

Optional parameters:
dt=(from header) time sampling interval (sec)
fmin=0 minimum frequency of filter array (hz)
fmax=NYQUIST maximum frequency of filter array (hz)
beta=3.0 ln[filter peak amp/filter endpoint amp]
band=.05*NYQUIST filter bandwidth (hz)
alpha=beta/band^2 filter width parameter
verbose=0 =1 supply additional info

Notes: This program produces a muiltifilter (as opposed to moving window)
representation of the instantaneous amplitude of seismic data in the
time-frequency domain. (With Gaussian filters, moving window and multi-
filter analysis can be shown to be equivalent.)

An input trace is passed through a collection of Gaussian filters
to produce a collection of traces, each representing a discrete frequency
range in the input data. For each of these narrow bandwidth traces, a
quadrature trace is computed via the Hilbert transform. Treating the narrow
bandwidth trace and its quadrature trace as the real and imaginary parts
of a "complex" trace permits the "instantaneous" amplitude of each
narrow bandwidth trace to be compute. The output is thus a representation
of instantaneous amplitude as a function of time and frequency.

Some experimentation with the "band" parameter may necessary to produce
the desired time-frequency resolution. A good rule of thumb is to run
sugabor with the default value for band and view the image. If band is
too big, then the t-f plot will consist of stripes parallel to the frequency
axis. Conversely, if band is too small, then the stripes will be parallel
to the time axis.

Examples:
suvibro | sugabor | suximage
suvibro | sugabor | suxmovie n1= n2= n3=
(because suxmovie scales it's amplitudes off of the first panel,
may have to experiment with the wclip and bclip parameters
suvibro | sugabor | supsimage | ... ( your local PostScript utility)


Credits:

CWP: John Stockwell, Oct 1994

Algorithm:

This programs takes an input seismic trace and passes it
through a collection of truncated Gaussian filters in the frequency
domain.

The bandwidth of each filter is given by the parameter "band". The
decay of these filters is given by "alpha", and the number of filters
is given by nfilt = (fmax - fmin)/band. The result, upon inverse
Fourier transforming, is that nfilt traces are created, with each
trace representing a different frequency band in the original data.

For each of the resulting bandlimited traces, a quadrature (i.e. pi/2
phase shifted) trace is computed via the Hilbert transform. The
bandlimited trace constitutes a "complex trace", with the bandlimited
trace being the "real part" and the quadrature trace being the
"imaginary part". The instantaneous amplitude of each bandlimited
trace is then computed by computing the modulus of each complex trace.
(See Taner, Koehler, and Sheriff, 1979, for a discussion of complex
trace analysis.

The final output for a given input trace is a map of instantaneous
amplitude as a function of time and frequency.

This is not a wavelet transform, but rather a redundant frame
representation.

References: Dziewonski, Bloch, and Landisman, 1969, A technique
for the analysis of transient seismic signals,
Bull. Seism. Soc. Am., 1969, vol. 59, no.1, pp.427-444.

Taner, M., T., Koehler, F., and Sheriff, R., E., 1979,
Complex seismic trace analysis, Geophysics, vol. 44,
pp.1041-1063.

Chui, C., K.,1992, Introduction to Wavelets, Academic
Press, New York.

Trace header fields accessed: ns, dt, trid, ntr
Trace header fields modified: tracl, tracr, d1, f2, d2, trid, ntr