Symmetric range (–T/2 to T/2)  (-F/2 to F/2)

            This document, the transform document and code described here are in SYMFFT.ZIP. – A more recent version of fftsr along with tests and a welcome describing it is in fftsrHth.zip.  The version in this zip allows for NIN and NOUT symmetrical time and frequency points and the power of 2 size N2P > max(NIN,NOUT).  This N2P –NIN points are symmetrically zero padded.  Only NOUT points are in the output file. 

6/6/2011 –

The Fortran file fftsr.for  SUBROUTINE FFTSR(IIN,IIOUT,NIN,NOUT,N2P,ISIGN) starts with NIN complex*16 data points in the file IIN, then uses N2P (a power of 2 > max(NIN,NOUT)) to transform these into NOUT compex*16 points in the file IIOUT. 

The Fortran file fftsrnr.for  SUBROUTINE FFTSNR(DIIN,DIIOUT,NIN,NOUT,N2P,ISIGN)  starts with NIN pairs of real*8 data points in DIIN arranged as (real,imag) with the odd numbers real, and converts these into NOUT, DIIOUT pairs using N2P (a power of 2 > max(NIN,NOUT)).

The points from –NP2/2 to –NIN/2-1 and from NIN –NIN/2 to N2P-1-N2p/2 are all zero.  This is discussed in ZeroPadding.doc htm ZeroPaddingF.doc .htm.  It is appropriate for functions that are zero, with zero derivatives, at –NIN/2 and NIN-1-NIN/2 but introduces problems for functions that are not.  

 

            As defined in Discrete Fourier Transform.doc#SymmetricRange .htm the Fourier sum is

        (1)

        (2)

      (3)

            The simplest implementation of (3) uses t = i Δf and f = m Δf with i  and m having both negative and positive values.  The functions h and H are stored in unformatted files and the evaluations are made directly.  dftsy.zip.  This method is not as fast as the FFT described below and may have some truncation problems for large values of f ´ t.

 

Equation (3) approximates

           (4)

The values of ti and fm are fundamental to the theory.  In the case that N is a power of 2, the Cooley-Tukey Fast Fourier Transform allows a very fast evaluation of

(5)

The sums in (5) are not quite the ones desired owing to the different ranges for i and m.  By subtracting N/2-1 from each value of i in the exponent
      (6)

This has values for m=-N/2, -N/2+1, …, N/2-1.  By using m’ = m+N/2+1, this exists for m’ values from 1 to N. 

(7)

Define

         (8)

(9)

(10)

As a practical matter to make a Fourier sum from –T/2 to T/2

  1. Find h(tk) at N = 2M evenly spaced points from –T/2 to T/2 – The point at T/2 is excluded
  2. Change the sign of every other value of h.
  3. call fft(h,N,1) to convert h hat into H hat.
  4. Change the sign of every other value of H hat to convert H hat into H.

 

The transform described here is FFTSR.FOR which handles the overhead above in approximating the integral.

      SUBROUTINE FFTSR(AIN complex*16 input data array of N values

,ADOUB real*8 temporary file array of 2N values, with an equivalence can be the input value

,AOUT, complex*16 array of N values output array

TIME,  real*8 – used to find deltat

N,  integer – number of data points

ISIGN  1 for hà H -1 for Hà h )

This is tested by in fftsr.zip using TransformOfGaussian.doc .htm as a test function.

           (11)

Setting w =0 .7, the range should be ~ -3.5 to 3.5 to have zero at both ends.  Set t0 = 0.3 to avoid an accidental symmetry.

Figure 1  The offset Gaussian to be transformed.  (4096 data points between t = -3.5 sec and 3.5 seconds)

        (12)

          (13)

Figure 2  Transform of Gaussian using 4096 pts. (frequencies from -2048/7 Hz and 2048/7 Hz)

Figure 3  Detail of peak

            An interesting feature of this is that only about 32 points are required to produce this detail.  The extra points are at frequencies too high or too low to be other than zero.  – Note 64  points produced a visually identical picture.

Figure 4  The absolute value of  (theoretical transform -  numerical transform)

            This one is subject to errors from using a finite number of integration points and from using a finite time.  The curve with 64 points peaks at about 9.5e-13.

Figure 5  h(t) – (back transform of (H(f) = forward transform of h(t))

            This is a measure of the extra error introduced by the precision of the transform.

Figure 6  The forward transform - backward transform using 64 points.

            This indicates that the sum of truncation errors contributes the bulk of the random error in the transform.