Designing a separate filter for each patch

Recall the prediction-error filter subroutine find_pef() . Given a data plane, this subroutine finds a filter that tends to whiten the spectrum of that data plane. The output is white residual. Now suppose we have a data plane where the dip spectrum is changing from place to place. Here it is natural to apply subroutine find_pef() in local patches. This is done by subroutine find_lopef(). The output of this subroutine is an array of helix-type filters, which can be used, for example, in a local convolution operator loconvol .  

module lopef {            # Local PEF estimation in patches.
  use patch		  # Estimate a vector of filters, one for each patch.
  use misinput
  use pef
contains
  subroutine find_lopef( wall, aa, npatch, nwall, nwind, mask) {
    optional                                 :: mask
    integer,       dimension(:), pointer     :: npatch, nwall, nwind
    real,          dimension(:), intent( in) :: wall, mask
    type( filter), dimension(:)              :: aa
    real,          dimension(:), pointer     :: windata, winmask
    integer			             :: i, stat
                        allocate( windata( product( nwind)))   # a patch
    if( present( mask)) allocate( winmask( product( nwind)))   # missing inputs	
    call patch_init( npatch, nwall, nwind)
    do i = 1, product( npatch) {                	       # do all patches
       stat = patch_lop( .false., .false., wall, windata)      # get a patch
       if( present( mask)) {
	   stat = patch_lop( .false., .false., mask, winmask)
	   call find_mask( (winmask /= 0.), aa (i))             # missing data
           }
       if( count(.not.aa(i)%mis) > size(aa(i)%lag))             # enuf eqns?
       	  call find_pef( windata, aa(i), niter=size( aa(i)%lag)) # find PEF
       else if( i > 1) 
          aa(i)%flt = aa(i-1)%flt				# use last PEF
       call patch_close()
       }
    deallocate( windata)
    if( present( mask))  deallocate( winmask)
    }
}
#       if( size(aa(i)%mis) - count(aa(i)%mis) > size(aa(i)%lag)) # enuf eqns?

We notice that when a patch has fewer regression equations than the filter has coefficients, then the filter is taken to be that of the previous patch.

 

# successive invocations apply successive filters from a vector.
#		      (will fail a dot product test?  Oft used with patching.)
module loconvol {
use helicon
integer, private  :: i
type( filter), dimension(:), pointer :: aa
#% _init( aa)
i = 0
#% _lop( xx, yy)
 integer stat1;   i = i + 1
 call helicon_init( aa( i))
 stat1 = helicon_lop( adj, .false., xx, yy)
}