SUMMARY AND COMPUTATION

Putting one polynomial into another or one infinite series into another is an onerous task, even if it does lead to a wavelet that is exactly causal. In practice we do operations that are conceptually the same, but for speed we do them with discrete Fourier transforms. The disadvantage is periodicity, i.e., negative times are represented computationally like negative frequencies. Negative times are the last half of the elements of a vector, so there can be some blurring of late times into negative ones. The subroutine we will use for Fourier transformation is fts()

subroutine fts( signi, nx, rr )
#   complex fourier transform.  if (signi>0) scale = 1; else scale=1/nx
#
#                      nx          signi*2*pi*i*(j-1)*(k-1)/nx
#   rr(k)  =  scale * sum rr(j) * e
#                     j=1             for k=1,2,...,nx=2**integer
#
integer nx, i, j, k, m, istep, pad2
real    signi, arg
complex rr(nx), cmplx, cw, cdel, ct
if( nx != pad2(nx) )    call erexit('fts: nx not a power of 2')
if( signi < 0.)
	do i= 1, nx
		rr(i) = rr(i) / nx
j = 1;  k = 1
do i= 1, nx {
        if (i<=j) { ct = rr(j); rr(j) = rr(i); rr(i) = ct }
        m = nx/2
        while (j>m && m>1) { j = j-m; m = m/2 }         # "&&" means .AND.
        j = j+m
        }
repeat {
        istep = 2*k;   cw = 1.;   arg = signi*3.14159265/k
        cdel = cmplx( cos(arg), sin(arg))
        do m= 1, k {
                do i= m, nx, istep
                        { ct=cw*rr(i+k);  rr(i+k)=rr(i)-ct;  rr(i)=rr(i)+ct }
                cw = cw * cdel
                }
        k = istep
        if(k>=nx) break
	}
return; end

integer function pad2( n )
integer n
pad2 = 1
while( pad2 < n )
        pad2 = pad2 * 2 
return; end