SPECTRAL FACTORIZATION CODE AND THEORY

My book PVI gives a full account of spectral factorization. Here is subroutine lapfac() which was used in the previous section to factor the Laplacian operator, and I give a sketch of the underlying Kolmogoroff theory. To invoke this subroutine, you need to supply one side of an autocorrelation function. For example, let us specify the negative of the 2-D Laplacian in a vector n = $256\times 256$ points long.

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

The reason for the 4.00004 instead of 4.0 arises in the theory below where we are required to take the logarithm of the spectrum. The spectrum of the Laplacian operator is k_x²+k_y², which vanishes at the point (k_x,k_y)=(0,0) where the logarithm is minus infinity. The .00004 prevents that.

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

module lapfac_mod { # Factor 2-D Laplacian.  Prepare input tables for helicon()
use autofac_mod
contains
  subroutine lapfac(e, n1, aa, lag) {
    real,                    intent (in)  :: e
    integer,                 intent (in)  :: n1
    integer, dimension (:),  intent (out) :: lag
    real,    dimension (:),  intent (out) :: aa
    complex, dimension ( 1024*64)         :: rr
    real                                  :: scale
    integer                               :: i, na
    rr = 0.
    rr(1)      = 4.00004 + e
    rr(2)      = -1.
    rr(1+1024) = -1
    call autofac( rr);   scale= rr(1);   rr = rr/scale
    na = size (lag)
    do i=1,na/2 {
	lag(i) = i;         lag(na/2+i) =      n1    -(i-1)
	 aa(i) = rr(i+1);    aa(na/2+i) = rr( 1024+1 -(i-1))
    }
  }
}

Subroutine lapfac() has its main job done by subroutine autofac() which will be shown after we examine the basic theory.