The recursion step

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

$$\left[\matrix{ m_0 \cr w_0 \cr m_1 \cr w_1 \cr m_2 \cr ... ... B_2 \cr B_3 \cr B_4 \cr B_5 \cr B_6 \cr B_7 \cr }\right]$$ (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 m_t of the multiplication of each odd row with the vector B as a weighted local average of adjacent vector coefficients B_t. The multiplication of each even row with the vector B results in a number w_t which measures the variation of these adjacent vector coefficients. Each pair m_t, w_t samples the same location in vector B. I refer to m_t, w_t as average or lowpass and detail or highpass coefficients, respectively. The distance between locations sampled by adjacent output pairs is two B_t samples; for example the pair m₁, w₁ is the weighted sum of the B elements 2 through 5, and the pair m₂, w₂ is the weighted sum of the B elements 4 through 7.

Obviously, not every set of c_k 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:

 

$$\left[\matrix{ B_0 \cr B_1 \cr B_2 \cr B_3 \cr B_4 \cr ... ... m_1 \cr w_1 \cr m_2 \cr w_2 \cr m_3 \cr w_3 \cr }\right]$$ (2)

This orthogonality condition further constrains the choice of the wavelet filter coefficients c_k, 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