The fast wavelet transform

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.