SUB FFT(g(,),p)
    ! fast Fourier transform of complex input data set g
    ! transform returned in g
    DIM Wp(2),factor(2),temp(2)
    LET N = 2^p                   ! number of data points
    LET N2 = N/2
    ! rearrange input data according to bit reversal
    LET j = 1
    FOR i = 1 to N-1
        ! g(i,1) is real part of f((i-1)*del_t)
        ! g(i,2) is imaginary part of f((i-1)*del_t)
        ! set g(i,2) = 0 if data real
        IF i < j then             ! swap values
           LET temp(1) = g(j,1)
           LET temp(2) = g(j,2)
           LET g(j,1) = g(i,1)
           LET g(j,2) = g(i,2)
           LET g(i,1) = temp(1)
           LET g(i,2) = temp(2)
        END IF
        LET k = N2
        DO while k < j
           LET j = j-k
           LET k = k/2
        LOOP
        LET j = j + k
    NEXT i
    ! begin Danielson-Lanczos algorithm
    LET jmax = 1
    FOR L = 1 to p
        LET del_i = 2*jmax
        LET Wp(1) = 1             !  Wp initialized at W^0
        LET Wp(2) = 0
        LET angle = pi/jmax
        LET factor(1) = cos(angle)     ! ratio of new to old W^p    
        LET factor(2) = -sin(angle)
        FOR j = 1 to jmax
            FOR i = j to N step del_i
                ! calculate transforms of length 2^L
                LET ip = i + jmax
                LET temp(1) = g(ip,1)*Wp(1) - g(ip,2)*Wp(2)
                LET temp(2) = g(ip,1)*Wp(2) + g(ip,2)*Wp(1)
                LET g(ip,1) = g(i,1) - temp(1)
                LET g(ip,2) = g(i,2) - temp(2)
                LET g(i,1) = g(i,1) + temp(1)
                LET g(i,2) = g(i,2) + temp(2)
            NEXT i
            ! find new W^p
            LET temp(1) = Wp(1)*factor(1) - Wp(2)*factor(2)
            LET temp(2) = Wp(1)*factor(2) + Wp(2)*factor(1)
            MAT Wp = temp
        NEXT j
        LET jmax = del_i
    NEXT L
END SUB