Programmer 7500/MAX_PROGRAMMERS.iso/PASCAL/FFT_PROC.ZIP/FFT-PROC.PAS

Pascal/Delphi Source File | 1989-03-20 | 11.5 KB | 372 lines

Program FOURIER;
{ Real FFT with single sine look-up per pass }
{      This is an almost exact copy of the FFT given in "Introduction to Pascal
  for Scientists" by James W. Cooper, John Wiley & Sons, 1981.  The differences
  are that the Procedure POST has been changed  to fix an error in the
  calculation of the DC term and a factor of 2 error in the other terms.
       The results now agree with another FFT I have which uses a different
  algorithm and is written in FORTRAN (FILE-FFT).  The program is
  composed of an actual FFT (Procedure FFT) which finds the complex transform
  of data in the order of ascending real followed by ascending imaginary.
  Shuffl is used in front of FFT to convert data in the form of real(0),
  imag(0),real(1),imag(1),...,imag(N) to the order required.  The post routine
  needs explaination.  An array of all real data can be transformed using
  a half-sized FFT by treating the all-real data as alternate real,complex and
  doing a final fast and simple conversion with Post.  This technique is
  described in "The Fast Fourier Transform" by E. O. Brigman.  For complex
  input data, the procedure FFT is used alone.  It can be used as a forward
  or inverse transform.  For real input only, use as shown.  The example
  used here is an exponential which is also in the book.
    This FFT seems to work, but I encourage you to check it out since
    algorithms like these can be tricky.

                    Roger Coleman, Palm Bay, Fl.}

Const
  Asize = 4096;           { Size of array goes here }
  PI2   = 1.570796327;    { PI/2 }
  F     = 16;             { Format constants }
  D     = 6;              {   "      " }

Type
  Xform = (Fwd,Inverse);  { Transform types }
  Xary  = Array[1..Asize] of Real;

Var
  I,j,N,modulo   : Integer;
  F1    : Text;    {File of Real;}
  X     : Xary;           { Array to transform }
  Inv   : Xform;          { Transform type - Forward or Inverse }
  w     : Real;           { Frequency of wave }

{*****************************************************************************}

Procedure Debug;          { Used to print intermediate results }

Var I3 : Integer;

Begin
  For I3 := 1 to N Do
    Writeln((I3-1):3,X[I3]:F:D);
End;    { Debug }

{*****************************************************************************}

Function Ibitr(j,nu : Integer) : Integer;
  { Function to bit invert the number j by nu bits }

Var
  i,j2,ib    : Integer;

Begin
  ib := 0;   { Default return integer }
  For i := 1 to nu do
  Begin
    j2 := j Div 2;          { Divide by 2 and compare lowest bits }
    { ib is doubled and bit 0 set to 1 if j is odd }
    ib := ib*2 + (j - 2*j2);
    j := j2;                { For next pass }
  End;  {For}
  ibitr := ib;              { Return bit inverted value }
End;    { Ibitr }

{*****************************************************************************}

Procedure Expand;

Var
  i,nn2,nx2    : Integer;

Begin
  nn2 := n div 2;
  nx2 := n + n;
  For i := 1 to nn2 do x[i+n] := x[i+nn2]; { Copy IM to 1st half IM position }
  For i := 2 to nn2 do
  Begin
    x[nx2-i+2] := -x[i+n];    { Replicate 2nd half Imag as negative }
    x[n-i+2] := x[i];         { Replicate 2nd half Real as mirror image }
  End;
  n := nx2;       { We have doubled number of points }
End;

{*****************************************************************************}

Procedure Post(inv : Xform);
  { Post processing for forward real transforms and pre-processing for inverse
    real transforms, depending on state of the variable inv.
    Note: Lines marked by * are modified from Cooper's orgignal to agree with
          Brigham's algorithm in The Fast Fourier Transform, P.169. Lines
          marked by ** are added for the same reason. }

Var
  nn2,nn4,l,i,m,ipn2,mpn2   : Integer;
  arg,rmsin,rmcos,ipcos,ipsin,ic,is1,rp,rm,ip,im   : Real;

Begin
  nn2 := n div 2;   { n is global }
  nn4 := n div 4;
  { imax represents PI/2 }
  For l := 1 to nn4 Do
  { Start at ends of array and work towards middle }
  Begin
    i := l+1;       { Point near beginning }
    m := nn2-i+2;   { Point near end }
    ipn2 := i+nn2;  { Avoids recalculation each time }
    mpn2 := m+nn2;  { Calcs ptrs to imaginary part }
    rp := x[i]+x[m];
    rm := x[i]-x[m];
    ip := x[ipn2]+x[mpn2];
    im := x[ipn2]-x[mpn2];
    { Take cosine of PI/2 }
    arg := (Pi2/nn4)*(i-1);
    ic := Cos(arg);
    { Cosine term is minus if inverse }
    If inv = Inverse Then ic := -ic;
    Is1 := Sin(arg);
    ipcos := ip*ic;  { Avoid remultiplication below }
    ipsin := ip*is1;
    rmsin := rm*is1;
    rmcos := rm*ic;
    x[i] := (rp + ipcos - rmsin)/2.0;    {* Combine for real-1st pt }
    x[ipn2] := (im - ipsin - rmcos)/2.0; {* Imag-1st point }
    x[m] := (rp - ipcos + rmsin)/2.0;    {* Real - last pt }
    x[mpn2] := -(im +ipsin +rmcos)/2.0;  {* Imag, last pt }
  End;  { For }
  x[1] := x[1]+x[nn2+1];    {** For first complex pair}
  x[nn2+1] := 0.0;          {**}
End;    { Post }

{*****************************************************************************}

Procedure Shuffl(inv : Xform);
  { This procedure shuffels points from alternate real-imaginary to
    1st-half real, 2nd-half imaginary order if inv=Fwd, and reverses the
    process if inv=Inverse.  Algorithm is much like Cooley-Tukey:
      Starts with large cells and works to smaller ones for Fwd
      Starts with small cells and increases if inverse }

Var
  i,j,k,ipcm1,celdis,celnum,parnum  : Integer;
  xtemp                             : Real;

Begin
  { Choose whether to start with large cells and go down or start with small
    cells and increase }

  Case inv Of

  Fwd: Begin
    celdis := n div 2;     { Distance between cells }
    celnum := 1;           { One cell in first pass }
    parnum := n div 4;     { n/4 pairs per cell in 1st pass }
    End;   { Fwd case }

  Inverse: Begin
    celdis := 2;           { Cells are adjacent }
    celnum := n div 4;     { n/4 cells in first pass }
    parnum := 1;
    End;  { Inverse case }

  End;  { Case }

  Repeat      { Until cells large if Fwd or small if Inverse }
    i := 2;
    For j:= 1 to celnum do
    Begin
      For k := 1 to parnum do  { Do all pairs in each cell }
      Begin
        Xtemp := x[i];
        ipcm1 := i+celdis-1;   { Index of other pts }
        x[i] := x[ipcm1];
        x[ipcm1] := xtemp;
        i := i+2;
      End;  { For k }

      { End of cell, advance to next one }
      i := i+celdis;
    End;    { For j }

    { Change values for next pass }

    Case inv Of
    Fwd:    Begin
      celdis := celdis div 2;         { Decrease cell distance }
      celnum := celnum * 2;           { More cells }
      parnum := parnum div 2;         { Less pairs per cell }
      End;   { Case Fwd }

    Inverse: Begin
      celdis := celdis * 2;           { More distance between cells }
      celnum := celnum div 2;         { Less cells }
      parnum := parnum * 2;           { More pairs per cell }
      End;   { Case Inverse }
    End;  { Case }

  Until  (((inv = Fwd) And (Celdis < 2)) Or ((inv=Inverse) And (celnum = 0)));

End;  { Shuffl }

{*****************************************************************************}

Procedure FFT(inv : Xform);
  { Fast Fourier transform procedure operating on data in 1st half real,
    2nd half imaginary order and produces a complex result }

Var
  n1,n2,nu,celnum,celdis,parnum,ipn2,kpn2,jpn2,
  i,j,k,l,i2,imax,index      : Integer;
  arg,cosy,siny,r2cosy,r2siny,i2cosy,i2siny,picons,
  y,deltay,k1,k2,tr,ti,xtemp  : Real;

Begin
  { Calculate nu = log2(n) }
  nu := 0;
  n1 := n div 2;
  n2 := n1;
  While n1 >= 2 Do
  Begin
    nu := nu + 1;            { Increment power-of-2 counter }
    n1 := n1 div 2;          { divide by 2 until zero }
  End;
  { Shuffel the data into bit-inverted order }
  For i := 1 to n2 Do
  Begin
    k := ibitr(i-1,nu)+1;  { Calc bit-inverted position in array }
    If i>k Then            { Prevent swapping twice }
    Begin
      ipn2 := i+n2;
      kpn2 := k+n2;
      tr := x[k];         { Temp storage of real }
      ti := x[kpn2];      { Temp imag }
      x[k] := x[i];
      x[kpn2] := x[ipn2];
      x[i] := tr;
      x[ipn2] := ti;
    End;  { If }
  End;  { For }

  { Do first pass specially, since it has no multiplications }
  i := 1;
  While i <= n2 Do
  Begin
    k := i+1;
    kpn2 := k+n2;
    ipn2 := i+n2;
    k1 := x[i]+x[k];        { Save this sum }
    x[k] := x[i]-x[k];      { Diff to k's }
    x[i] := k1;             { Sum to I's }
    k1 := x[ipn2]+x[kpn2];  { Sum of imag }
    x[kpn2] := x[ipn2]-x[kpn2];
    x[ipn2] := k1;
    i := i+2;
  End;  { While }

  { Set up deltay for 2nd pass, deltay=PI/2 }
    deltay := PI2;   { PI/2 }
    celnum := n2 div 4;
    parnum := 2;     { Number of pairs between cell }
    celdis := 2;     { Distance between cells }


  { Each pass after 1st starts here }
  Repeat             { Until number of cells becomes zero }

{ Writeln(Lst,'After Nth Pass:');  ### }
{ Debug; }

    { Each new cell starts here }
    index := 1;
    y := 0;      { Exponent of w }
    { Do the number of pairs in each cell }
    For i2 := 1 To parnum Do
    Begin
      If y <> 0 Then
      Begin                { Use sines and cosines if y is not zero }
        cosy := cos(y);    { Calc sine and cosine }
        siny := sin(y);
        { Negate sine terms if transform is inverse }
        If inv = Inverse then siny := -siny;
      End;   { If }
      { These are the fundamental equations of the FFT }
      For l := 0 to celnum-1 Do
      Begin
        i := (celdis*2)*l + index;
        j := i+celdis;
        ipn2 := i + n2;
        jpn2 := j + n2;
        If y = 0 Then   { Special case for y=0 -- No sine or cosine terms }
        Begin
          k1 := x[i]+x[j];
          k2 := x[ipn2]+x[jpn2];
          x[j] := x[i]-x[j];
          x[jpn2] := x[ipn2]-x[jpn2];
        End   { If-Then }
        Else
        Begin
          r2cosy := x[j]*cosy;   { Calc intermediate constants }
          r2siny := x[j]*siny;
          i2cosy := x[jpn2]*cosy;
          i2siny := x[jpn2]*siny;
          { Butterfly }
          k1 := x[i] + r2cosy + i2siny;
          k2 := x[ipn2] - r2siny + i2cosy;
          x[j] := x[i] - r2cosy - i2siny;
          x[jpn2] := x[ipn2] + r2siny - i2cosy;
        End;   { Else }
        { Replace the i terms }
        x[i] := k1;
        x[ipn2] := k2;

      { Advance angle for next pair }
      End;  { For l }

      Y := y + deltay;
      index := index + 1;
    End;  { For i2 }

    { Pass done - change cell distance and number of cells }
    celnum := celnum div 2;
    parnum := parnum * 2;
    celdis := celdis * 2;
    deltay := deltay/2;

  Until celnum = 0;

End;  { FFT }

{*****************************************************************************}


Begin    { * Main program * }
  For i := 0 To Asize-1 Do
  Begin
    x[i] := 0.0;
  End;

{  Write('Enter number of points: ');
  Readln(n);}
  n := 32;
  If (n > Asize) Then
  Begin
    Writeln('Too large, will use maximum');
    n := Round(asize/2.0);
  End;
  For i := 2 to n Do x[i] := Exp((1-i)*0.25); { Create Real array }
  x[1] := 0.5;

  Writeln('Input Array:');
  Debug;

  Shuffl(Fwd);

  FFT(Fwd);

  Post(Fwd);

  For i:= 1 to n Do x[i] := x[i]*8/n;
  Writeln('Forward FFT with real array first:');
  Debug;


End.