share/numeric/fft.usg

   1 FAST FOURIER TRANSFORMS
   2
   3 To load the routines do LOADFILE(FFT,FASL,DSK,SHARE);
   4
   5 The basic functions are:
   6         FFT --> fast fourier transform
   7         IFT --> inverse fast fourier transform
   8
   9 These functions perform a (complex) fast fourier transform on either
  10 1 or 2 dimensional FLOATING-POINT arrays, obtained by:
  11 ARRAY(<ary>,FLOAT,<dim1>); or ARRAY(<ary>,FLOAT,<dim1>,<dim2>);
  12 For 1D arrays <dim1> must equal 2^n-1, and for 2D arrays
  13 <dim1>=<dim2>=2^n-1 (i.e. the array is square).  (Recall that
  14 MACSYMA arrays are indexed from a 0 origin so that there will be
  15 2^n and (2^n)^2 arrays elements in the above two cases.)
  16
  17 Calling sequence is:
  18 FFT(<real-array>,<imag-array>); or IFT(<real-array>,<imag-array>);
  19
  20 The real and imaginary arrays must of course be the same size.
  21 The transforms are done in place so that <real-array> and <imag-
  22 array> will contain the real and imaginary parts of the transform.
  23 (If you want to keep the transformed and un-transfromed arrays
  24 separate copy the arrays before calling FFT or IFT using the
  25 FILLARRAY function - see SHARE;ARRAY USAGE.)
  26
  27 The definitions of the Fast Fourier Transform and its inverse
  28 are given here.  Here A is the array to be transformed and AT is
  29 its transform.  Both A and AT are complex arrays, although as noted
  30 above FFT and IFT can only deal with separate real arrays for
  31 the real and imaginary parts of A and AT.  N (or N^2) is the number
  32 of elements in A in the 1D (or 2D) case.  (In fact these definitions
  33 are not of the FFTs but of the discrete Fourier transforms.  The
  34 FFT and IFT functions merely provided efficient algorithms for the
  35 implementation of these definitions.)
  36
  37 1D case:
  38
  39       N - 1
  40       ====                     - 1
  41       \          2 %I %PI I K N
  42 AT  =  >    A  %E
  43   K   /      I
  44       ====
  45       I = 0
  46
  47           N - 1
  48           ====                        - 1
  49       - 1 \           - 2 %I %PI I K N
  50 A  = N     >    AT  %E
  51  I        /       K
  52           ====
  53           K = 0
  54
  55 2D case:
  56
  57          N - 1 N - 1
  58          ====  ====                                - 1
  59          \     \             2 %I %PI (I K + J L) N
  60 AT     =  >     >    A     %E
  61   K, L   /     /      I, J
  62          ====  ====
  63          I = 0 J = 0
  64
  65              N - 1 N - 1
  66              ====  ====                                   - 1
  67          - 2 \     \              - 2 %I %PI (I K + J L) N
  68 A     = N     >     >    AT     %E
  69  I, J        /     /       K, L
  70              ====  ====
  71              K = 0 L = 0
  72
  73 Other functions included in this file are:
  74
  75 POLARTORECT(<magnitude-array>,<phase-array>);
  76 converts from magnitude/phase form into real/imaginary form
  77 putting the real part in the magnitude array and the imaginary part
  78 into the phase array.  (<real>=<magnitude>*COS(<phase>)
  79 and <imaginary>=<magnitude>*SIN(<phase>).)
  80
  81 RECTTOPOLAR(<real-array>,<imag-array>);
  82 undoes POLARTORECT.  The phase is given in the range
  83 from -%PI to %PI.
  84
  85 Like FFT and IFT these function accept 1 or 2 dimensional
  86 arrays.  However, the array dimensions need not be a power of 2,
  87 nor need the 2D arrays be square.
  88
  89 (The above 4 functions return a list of their arguments)
  90
  91 EXAMPLE - do LOADFILE(FFT,FASL,DSK,SHARE); DEMO(FFT,DEMO,DSK,SHARE);
  92
  93 DISCLAIMER - I didn't write these routines, but copied them from
  94 AI:LIBDOC;FFT BWOOD1.  However you can report bugs to me.
  95
  96                                                 CFFK@MC
  97 GJC@MIT-MC 05/15/81 20:40:29 Re: FFT
  98 To: JPLJMR at MIT-MC, (FILE [SHARE;FFT MAIL]) at MIT-MC
  99 The FFT you get by doing LOAD("FFT") is the natural complex one:
 100
 101       N - 1
 102       ====                     - 1
 103       \          2 %I %PI I K N
 104  T  =  >    A  %E
 105   K   /      I
 106       ====
 107       I = 0
 108
 109 So if you do FFT(A_REAL,A_IMAG); what you get back is the
 110 coefficients of an exponential series which is the inverse
 111 transform:
 112
 113           N - 1
 114           ====                        - 1
 115       - 1 \           - 2 %I %PI I K N
 116 A  = N     >     T  %E
 117  I        /       K
 118           ====
 119           K = 0
 120
 121 Your question: "What do you do if all you have is a real function,
 122 and what you want is the coefficients of the SIN and COS terms?"
 123 A simple way is to put the function in A_REAL, and set A_IMAG to
 124 all zeros. Then FFT. Using the exponential formula for SIN and
 125 COS you get that (T[K]-T[-K])/2 gives the SIN term, and
 126 (T[K]+T[-K])/2 gives the COS term.
 127 A clever way starts out by putting the odd part of the function in A_REAL,
 128 and the even part in A_IMAG.
 129
 130 Say you do FFT(REAL,IMAG); then:
 131
 132 COSTERM[J]:(REAL[J]+REAL[N-J])/2;
 133 SINTERM[J]:(IMAG[J]-IMAG[N-J])/2;
 134
 135 Now, there is this question about the factor N^(-1), which
 136 some people write as SQRT(2*%PI*N)^(-1).
 137 All I can suggest is that you take the above formuli and use
 138 macsyma to carefully check that what you want to do is
 139 represented correctly. Using the PLOT2 package can also help in
 140 debugging.
 141
 142
 143 -gjc