Next: The butterfly Up: WAVELET TRANSFORM OF A Previous: WAVELET TRANSFORM OF A

## The recursion step

In each recursion step we filter an input vector by multiplication with a quadratic matrix C:
 (1)

The matrix is very similar to a convolution matrix. A wraparound boundary condition is implemented in the lower right corner.

Let us imagine the result mt of the multiplication of each odd row with the vector B as a weighted local average of adjacent vector coefficients Bt. The multiplication of each even row with the vector B results in a number wt which measures the variation of these adjacent vector coefficients. Each pair mt, wt samples the same location in vector B. I refer to mt, wt as average or lowpass and detail or highpass coefficients, respectively. The distance between locations sampled by adjacent output pairs is two Bt samples; for example the pair m1, w1 is the weighted sum of the B elements 2 through 5, and the pair m2, w2 is the weighted sum of the B elements 4 through 7.

Obviously, not every set of ck filter coefficients has such a highpass/lowpass characteristic. The coefficients are correspondingly chosen (see the section The lowpass filter coefficients''). In order for the DWT to be easily invertible, the matrix C must be orthogonal. The inverse matrix is then the transpose of C:

 (2)

This orthogonality condition further constrains the choice of the wavelet filter coefficients ck, as I discuss in the section The low pass filter coefficients.'' The subroutine dwt.rt executes the recursion step.

#
#%
subroutine dwt(conj, b,n, c,nop)
integer        conj,   n,   nop, i,j, no2
real 	             b(n),c(nop)
temporary real tmp(n)

no2 = int(n/2)
if(conj != 0)
Do j = 1, no2{
tmp(2*j-1) =  b(j    )
tmp(2*j  ) =  b(j+no2)
}
call conjnull(conj, 0, b,n, tmp,n)

Do j = 1, n - 1, 2
Do i = 1, nop
if(conj == 0){
tmp(j  ) = tmp(j  ) +           c(      i    ) * b(mod(j+i-2,n) + 1)
tmp(j+1) = tmp(j+1) - (-1)**i * c(nop - i + 1) * b(mod(j+i-2,n) + 1)
}
else{
b(mod(j+i-2,n) + 1) = b(mod(j+i-2,n) + 1) +          c(    i  ) * tmp(j  )
b(mod(j+i-2,n) + 1) = b(mod(j+i-2,n) + 1) -(-1)**i * c(nop-i+1) * tmp(j+1)
}
if(conj == 0)
Do j = 1, no2{
b(j    ) = tmp(2*j-1)
b(j+no2) = tmp(2*j  )
}
return; end


Next: The butterfly Up: WAVELET TRANSFORM OF A Previous: WAVELET TRANSFORM OF A
Stanford Exploration Project
11/18/1997