Finding a causal wavelet from a prescribed spectrum

Now for the task of going from a spectrum to a causal wavelet. Take as given the spectrum of the causal wavelet B(Z). This means that we are not given B(Z) itself, but we are given $S(Z) = \overline{B} (1/Z) B(Z)$.Assuming no zeros in the spectrum $S(\omega )$,it is easy to find the log of the spectrum $U=\ln S$.The spectrum may be specified as autocorrelation coefficients or values on the unit circle $Z=e^{i\omega}$.Thus,  

$$\overline{B}(1/Z) \, B(Z) = S(Z) = e^{\ln S(Z)} = e^{U(Z)} = e^{\overline{C}(1/Z) + C(Z)} = e^{\overline{C}(1/Z)} \, e^{C(Z)}$$ (16)

Given the spectrum S(Z) for each value on the unit circle, we could deduce the log spectrum $U(Z) = \ln S(Z)$ at each point on the unit circle:

 

$$U(Z) \eq \ln [ S(Z) ] \eq {\overline{C}(1/Z) + C(Z)}$$ (17)

This is the answer we have been looking for. Given U(Z) for all real values of $\omega$, we could inverse transform to the time domain, obtaining the two-sided function $u_t = \bar c_{-t} + c_t$.Setting to zero the coefficients at negative times eliminates $\bar c_{-t}$, leaving just c_t; hence C(Z). And we already know that the exponential of C(Z) gives B(Z) with a causal b_t. This method is known as ``Kolmogoroff spectral factorization," after the mathematician who discovered it.

The autocorrelation factorization subroutine is autofac()  

module autofac_mod { # Autocorrelation factorization.
use fft
contains
  subroutine autofac( rr) {    
    complex, dimension (:) :: rr # input:  rr = autocorrelation
    integer n                    # output: rr = min phase wavelet
    n = size (rr)
    rr(n:n/2:-1) = conjg( rr(2:n/2+2));	   call fts( +1, rr) # make spectrum
    rr = clog( rr);                        call fts( -1, rr)
    rr(1    ) = rr(1    ) / 2.
    rr(1+n/2) = rr(1+n/2) / 2.
    rr(1+n/2+1: n) = 0.;                   call fts( +1, rr)
    rr = cexp( rr);                        call fts( -1, rr)
  }
}

where the inputs and outputs were explained earlier with subroutine lapfac() . Subroutine autofac() is based on the classic 1-D discrete FT subroutine fts() .