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 .  

#$=head1 NAME
#$
#$lopef - estimate pef in patches
#$
#$=head1 SYNOPSIS
#$
#$C<call find_lopef(wall,aa,npatch,nwall,nwind,mask)>
#$
#$=head1 PARAMETERS
#$
#$=over 4
#$
#$=item  wall - C<real(:)>
#$
#$=item  aa   - type(filter)
#$
#$       helix filter
#$
#$=item  npatch -C<integer(:)>
#$
#$       size of patch
#$
#$=item  nwall  -C<integer(:)>
#$
#$=item nwind  -C<integer(:)>
#$
#$=item mask   -C<real(:) >
#$
#$      Mask of known and unknown data
#$
#$=back
#$
#$=head1 DESCRIPTION
#$
#$Estiamte a PEF in patches
#$
#$=head1 SEE ALSO
#$
#$L<pef>,L<misinput>,L<patch>
#$
#$=head1 LIBRARY
#$
#$B<geef90>
#$
#$=cut

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)
}