The shape of the 3-D prediction-error filter

Imagine a 3-D PEF and through the ``1.0'' put a plane that has an infinitesimal tilt so as to get all zeros on one side of the plane and all fitting coefficients are on the other. Figure 3 shows the volume of free coefficients of the three-dimensional prediction-error filter . The top plane is the 2-D filter seen earlier. The top plane can be visualized as the area around the end of a helix. Above the top plane are zero-valued anticausal filter coefficients.

 

3dpef Figure 3 Three-dimensional prediction-error filter.

Module lagtable constructs a linear subsript of filter lags along the helix from the 3-D cartesian description. To print out the cartesian form of a helix filter, you can use module format .

 

module lagtable {          # Find 3-D PEF filter coef locations on a helix.
contains
  subroutine table3( n1,n2, lag1,lag2, a1,a2,a3, aa, lag) {
    integer, intent (in)            :: n1,n2, lag1,lag2, a1,a2,a3
    real,    dimension (:), pointer :: aa	         # prepare for helicon.
    integer, dimension (:), pointer :: lag	         # prepare for helicon.
    integer                         :: i,j, i1,i2,i3, na
    na = (a3-1)*a1*a2 + (a2-lag2)*a1 + (a1-lag1)
    allocate ( aa( na), lag( na));  aa = 0.	    # a(lag1,lag2,1) = 1.
    do i3 = 1, a3 {
     do i2 = 1, a2 {    if( i3 == 1 .and. i2 < lag2                  ) cycle
      do i1 = 1, a1 {   if( i3 == 1 .and. i2 == lag2 .and. i1 <= lag1) cycle
	  i = i1-lag1 + (i2-lag2)*a1 + (i3-1)*a1*a2
          j = i1-lag1 + (i2-lag2)*n1 + (i3-1)*n1*n2	
          lag( i)  =  j                             # a(i1,i2,i3)
          }}}
    }
}

 

module format {
contains
  subroutine print3( n1, n2, lag1, lag2, a1, a2, a3, aa, lag) {
    integer,                intent( in) :: n1, n2, lag1, lag2,  a1, a2, a3
    real,    dimension( :), intent( in) :: aa
    integer, dimension( :), intent( in) :: lag
    integer                             :: ia, i1, i2, i3, i
    real, dimension( a1, a2, a3)        :: filt            # cartesian filter
    filt = 0.;    filt( lag1, lag2, 1) = 1.
    do ia = 1, size( lag) {
       i = lag( ia) + lag1 + ( lag2-1)*n1 - 1
       i1 = mod( i        , n1)
       i2 = mod( i/n1     , n2)
       i3 =      i/(n1*n2)
       filt( i1+1, i2+1, i3+1) = aa( ia)
    }
                    write( 0, *) "----------------------"
    do i3 = 1, a3 {                                        # loop over planes
       if( i3 > 1)  write( 0, *) "+++++++++++++++++++++++++++++++++++++++++++"
       do i1 = 1, a1
                    write( 0, "(10f7.3)") filt( i1, :, i3) # print filter row
    }
                    write( 0, *) "----------------------"
  }
}

A reasonable arrangement for a small filter is a1=5 a2=3 a3=2 lag1=1+a1/2=3 lag2=1+a2/2=2. I set all the filter coefficients to 2 except for a(3,2,1)=1. Then I invoked print3. The output is:

-------------------------------------------
  0.000  0.000  2.000
  0.000  0.000  2.000
  0.000  1.000  2.000
  0.000  2.000  2.000
  0.000  2.000  2.000
 +++++++++++++++++++++++++++++++++++++++++++
  2.000  2.000  2.000
  2.000  2.000  2.000
  2.000  2.000  2.000
  2.000  2.000  2.000
  2.000  2.000  2.000
 -------------------------------------------

where i1 runs vertical, i2 runs horizontal in a block, i3 jumps between blocks.