Next: 1-D seismic signal decomposition
Up: INTRODUCTION TO WAVELETS
Previous: INTRODUCTION TO WAVELETS
The following FORTRAN routine performs wavelet decomposition and reconstruction.
It has been written based on hints from Strang's article.
The half-cycle square-wave wavelet requires no trigonometric functions.
Wavelet basis functions are recursively computed from previous iterations.
The progression in this code is from Nyquist wavelengths to DC.
Each pass is half the the length of the previous pass, totaling order n
computations.
C 1-D WAVELET TRANSFORM
SUBROUTINE wavelet (x,n,direction)
INTEGER n, direction ! n is length of x, power of two
REAL x(n) ! x() is input and output
INTEGER m, i
REAL y(2048) ! y() is permutation buffer
C FORWARD
IF (direction .GT. 0) THEN
C POWER OF 2 PASSES, EACH HALF AS LONG
m = n
C BUTTERFLY
10 continue
DO i=1,m,2
y(i) = (x(i) + x(i+1)) / 2
y(i+1) = (x(i) - x(i+1)) / 2
ENDDO
C PERMUTATION
m = m / 2
DO i=1,m
x(i) = y(i+i-1)
x(m+i) = y(i+i)
ENDDO
if (m .GT. 1) GOTO 10
C INVERSE
ELSE
C POWER OF 2 PASSES, EACH TWICE AS LONG
m = 1
20 continue
C PERMUTATION
DO i=1,m
y(i+i-1) = x(i)
y(i+i) = x(m+i)
ENDDO
m = m + m
C BUTTERFLY
DO i=1,m,2
x(i) = (y(i) + y(i+1))
x(i+1) = (y(i) - y(i+1))
ENDDO
if (m .LT. n) GOTO 20
ENDIF
END
Figure 1:
Wavelet decomposition and reconstruction of a 1-D seismic trace.
Upper-left: Unweighted basis functions--half-cycle square-waves.
Upper-right: Decomposition: signal, transform, weighted basis functions.
Lower-left: Reconstruction--adding rougher functions bottom to top.
The top and bottom traces match.
Lower-right: Reconstruction--adding smoother functions bottom to top.
The bottom and second from bottom traces match.
|
The transform output overlays the original data.
They are in order from DC to Nyquist like a Fourier transform.
There is one DC coefficient, one full wavelength coefficient, two half wavelength coefficients,
to n/2 Nyquist wavelength coefficients.
Within each scale group the wavelets cover the whole signal.
One can make an overall estimate of the relative frequency distribution
in a signal by viewing a plot of the coefficients, but the intervals are
not as smooth as a Fourier transform.
Figure 2:
The right image are the decomposition coefficients of the two dimensional
signal on the left.
The DC coefficient is in the upper left corner.
Wavelengths decrease to the lower right.
|
Next: 1-D seismic signal decomposition
Up: INTRODUCTION TO WAVELETS
Previous: INTRODUCTION TO WAVELETS
Stanford Exploration Project
1/13/1998