FT by Z-transform

The program slowft() is unnecessarily slow, requiring us to compute a complex exponential at each step. By reorganizing easily using the Z-transform, the computational load can be reduced by about a factor of five (from a complex exponential to a complex multiply) at every step.

For simplicity we consider a signal that is only four points long:  

$$B(\omega) \eq b_0 + b_1 Z + b_2 Z^2 + b_3 Z^3$$ (10)

Reorganizing the polynomial (10) by nesting gives

$$B(\omega) \eq b_0 + Z(b_1 + Z(b_2 + Z(b_3 )))$$ (11)

A subroutine for evaluating $B(\omega)$ in this way is polyft(). 

# Fourier transform by polynomial evaluation.
subroutine polyft( nt,tt, nw,cww )
integer nt      # number of points in the time domain
integer nw      # number of points in the fourier transform
real    tt(nt)  # sampled function of time
complex cww(nw) # sampled fourier transform
        integer it, iw
        real omega
        complex cz, cw
        do iw= 1, nw {
                omega =  3.14159265 * (iw-1.) / ( nw-1.)
                cz = cexp( cmplx( 0., omega ) )
                cw = tt(nt)
                do it= nt-1, 1, -1      # loop runs backwards
                        cw = cw * cz + tt(it)
                cww(iw) = cw
                }
        return; end