Madagascar development blog

Monday, November 21. 2011

RNA Interpolation

A new paper is added to the collection of reproducible documents:

Seismic data interpolation beyond aliasing using regularized nonstationary autoregression

Posted by Sergey Fomel in Documentation at 15:11 | Comments (0) | Trackbacks (0)

Wednesday, November 16. 2011

Fourier finite differences

A new paper is added to the collection of reproducible documents:

Fourier finite-difference wave propagation

Posted by Sergey Fomel in Documentation at 21:26 | Comments (0) | Trackbacks (0)

Saturday, November 5. 2011

Program of the month: sfenvelope

Complex trace attributes were introduced into geophysics by the paper

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

If

is the input seismic trace, then the analytical trace is defined as the complex-valued signal

where

is the Hilbert transform of

$h(t) =\displaystyle \frac{1}{\pi} \int \frac{s(\tau)}{t-\tau} d\tau\;.$

The signal envelope is the positive signal

$e(t)=\sqrt{s^2(t)+h^2(t)}\;.$

A phase-rotated seismic signal is

$p(t)=s(t)\,\cos{\phi} +h(t)\,\sin{\phi}$

where $\phi$ is the phase of rotation.

By default, sfenvelope computes the signal envelope. It can also produce a phase-rotated signal if given hilb=y and phase=. If phase=90 (the default value), the phase-rotated signal will be simply the Hilbert transform of the input.

The following figure from book/rsf/rsf/sfenvelope illustrate an application of sfenvelope:

Computing the discrete Hilbert transform is not a trivial task. In the Fourier domain, the continuous Hilbert transform is given by

$\displaystyle H(\omega) = i\,\operatorname{sgn}(\omega)\,S(\omega)$

where sgn is the sign function. The discontinuity of the sign function in the frequency domain at $\omega=0$ is related to the slow

decay of the filter impulse response in the time domain. The discontinuity at the Nyquist frequency creates additional oscillations. Different practical implementations shorten the filter impulse response by effectively smoothing the Fourier-domain discontinuities.

The Madagascar implementation of the discrete Hilbert transform follows the algorithm described in

Pei, S.-C., and P.-H. Wang, 2001, Closed-form design of maximally flat FIR Hilbert transformers, differentiators, and fractional delayers by power series expansion: IEEE Trans. on Circuits and Systems, v. 48, No. 4, 389-398.

The accuracy/cost trade-off is controlled by two parameters: order= and ref=. The following figures frombook/rsf/rsf/sfenvelope illustrate the effect of the order= parameter:

The Seismic Unix implementation (suhilb program) applies a Hamming window in the time domain. For some reason, it has the filter polarity reversed:

A multidimensional analog of the Hilbert transform is the Riesz transform. It is implemented in the sfriesz program.

Posted by Sergey Fomel in Programs at 09:00 | Comments (0) | Trackbacks (12)

Tutorial

New users often ask for a brief comprehensible introduction to Madagascar.

Jeff Gowdin contributes A brief introduction to Madagascar, a tutorial book designed to provide such an introduction.

The purpose of this document is to teach new users how to use the powerful tools in Madagascar to: process data, produce documents and build your own Madagascar programs.

You can access the tutorial in the following forms: