WILSON-BURG SPECTRAL FACTORIZATION

(If you are new to this material, you should pass over this section.) Spectral factorization is the job of taking a power spectrum and from it finding a causal (zero before zero time) filter with that spectrum. Methods for this task (there are many) not only produce a causal wavelet, but they typically produce one whose convolutional inverse is also causal. (This is called the ``minimum phase'' property.) In other words, with such a filter we can do stable deconvolution. Here I introduce a new method of spectral factorization that looks particularly suitable for the task at hand. I learned this new method from John Parker Burg who attributes it to an old paper by Wilson (I find Burg's explanation, below, much clearer than Wilson's.)

Below find subroutine lapfac() which was used in the previous section to factor the Laplacian operator. To invoke the factorization subroutine, you need to supply one side of an autocorrelation function. For example, let us specify the negative of the 2-D Laplacian (an autocorrelation) in a vector n = $256\times 256$ points long.

        rr(1)     =  4.
        rr(2)     = -1.
        rr(1+256) = -1.

Subroutine lapfac() finds the helical derivative (factored negative Laplacian) and then prepares the required filter coefficient tables for the helix convolution and deconvolution subroutines.  

#$=head1 NAME
#$
#$lapfac -  factor a 2-D Laplacian
#$
#$=head1 SYNOPSIS
#$
#$C<aa= lapfac2(eps,n1,na)>
#$  
#$=head1 PARAMETERS
#$  
#$=over 4    
#$  
#$=item eps - real 
#$  
#$      small value for approximation
#$  
#$=item n1  - integer 
#$  
#$      length of n1 axis
#$
#$=item na - integer 
#$   
#$      1/2 number of filter coefs for representation
#$  
#$=back
#$  
#$=head1 RETURN VALUE
#$  
#$=over 4
#$  
#$=item aa    - type(filter)
#$
#$      Helix filter
#$  
#$=back
#$     
#$=head1  DESCRIPTION
#$              
#$Create a one sided, minimum phase filter that approximates a laplaciaan
#$     
#$=head1 SEE ALSO
#$     
#$L<wilson>
#$  
#$=head1 LIBRARY
#$
#$B<geef90>
#$  
#$=cut
module lapfac {                             # Factor 2-D Laplacian.  
  use wilson
contains
  function lapfac2( eps, n1, na)   result (aa) {
    type( filter)                       :: aa, lap
    real,                   intent( in) :: eps
    integer,                intent( in) :: n1, na
    integer                             :: i
    real                                :: a0, lap0
    call allocatehelix( lap, 2)         # laplacian filter
    lap0 = 4. + eps                     # zero lag coeff.
    lap%lag = (/ 1, n1 /)		# lag(1)= 1; lag(2)=n1  # one side only
    lap%flt = -1.			# flt(1)=-1; flt(2)=-1
    call allocatehelix( aa, 2*na)       # laplacian derivative
    aa%flt = 0.;			# probably done already in allocation. 
    do i = 1, na {
       aa%lag( i   ) =      i		# early lags (first row)
       aa%lag( na+i) = n1 + i - na	# late lags (second row)
       }
    call wilson_init( 10 * n1 )
    call wilson_factor( 20, lap0, lap, a0, aa)
    call wilson_close()
    call deallocatehelix( lap)
    }
}

Subroutine lapfacn() has its main job done by subroutine wilson_factor() shown after the Wilson-Burg theory.