program ALG081;
{  FAST FOURIER TRANSFORM

   To compute the coefficients in the discrete approximation
   for the data (x(J),y(J)), 0<=J<=2m-1 where m=2**p and
   x(J)=-pi+J*pi/m for 0<=J<=2m-1.

   INPUT:  m; y(0),y(1),...y(2m-1).

   OUTPUT: complex numbers c(0),...,c(2m-1); real numbers
           a(0),...,a(m); b(1),...,b(m-1).

   NOTE:   The multiplication by EXP(-K*PI*I) is done within the
           program.
}
const
   ZERO = 1.0e-20;
var
   CR, CI, WR, WI, Y : array[1..64] of real;
   WWR, WWI, T1R, T1I, T3R, T3I : real;
   Z, XR, XI, YR, YI, TW : real;
   NG,N, N2, NU1, I, K, L, M, NP, J, M1 : integer;
   FLAG : integer;
   OK : boolean ;
   A : char;
   INP,OUP : text;
   NAME : string [ 14 ];
{  Change F if program is to calculate y                               }
function F ( Z : real ) : real;
   begin
      F := exp(-1.0-Z/pi)
   end;
procedure INPUT;
   begin
      writeln('This is the Fast Fourier Transform.');
      writeln(' ');
      writeln('The user must make provisions if the');
      writeln('interval is not [-pi,pi].');
      writeln('The example illustrates the required');
      writeln('provisions under input method 3.');
      OK := false;
      while ( not OK ) do
         begin
            writeln ('Choice of input method: ');
            writeln ('1. Input entry by entry from keyboard ');
            writeln ('2. Input data from a text file ');
            writeln ('3. Generate data using a function F ');
            writeln ('Choose 1, 2, or 3 please. ');
            readln ( FLAG );
            if ( FLAG = 1 ) or ( FLAG = 2 ) or ( FLAG = 3 ) then OK := true
         end;
      case FLAG of
         1 : begin
                OK := false;
                while ( not OK ) do
                   begin
                      writeln ('Input m. ');
                      readln ( M );
                      if ( M > 0 ) then
                         begin
                            OK := true;
                            N := 2*M;
                            for I := 1 to N do
                               begin
                                  J := I - 1;
                                  writeln ('Input y(',J,').');
                                  readln(Y[I])
                               end
                         end
                      else writeln ('Number must be a positive integer. ')
                   end
             end;
         2 : begin
                write ('Has a text file been created with the ');
                writeln ('entries y(0),...,y(2m-1)');
                writeln ('separated by a blank? ');
                writeln ('Enter Y or N. ');
                readln ( A );
                if ( A = 'Y' ) or ( A = 'y' ) then
                   begin
                      write ('Input the file name in the form - ');
                      writeln ('drive:name.ext, ');
                      writeln ('for example:   A:DATA.DTA ');
                      readln ( NAME );
                      assign ( INP, NAME );
                      reset ( INP );
                      OK := false;
                      while ( not OK ) do
                         begin
                            writeln ('Input number m. ');
                            readln ( M );
                            N := 2*M;
                            if ( N > 0 ) then
                               begin
                                  for I := 1 to N do
                                     readln ( INP, Y[I]);
                                  close ( INP );
                                  OK := true
                               end
                            else writeln ('Number must be a positive integer.')
                         end
                   end
                else
                   begin
                      write ('The program will end so the input file can ');
                      writeln ('be created. ');
                      OK := false
                   end
             end;
         3 : begin
                write ('Has the function F been created in the program ');
                writeln ('immediately preceding ');
                writeln ('the INPUT procedure? ');
                writeln ('Enter Y or N. ');
                readln ( A );
                if (( A = 'Y' ) or ( A = 'y' )) then
                   begin
                      OK := false;
                      while ( not OK ) do
                         begin
                            writeln ('Input the number m.');
                            readln ( M );
                            N := 2*M;
                            if ( N > 0 ) then
                               begin
                                  for I := 1 to N do
                                     begin
                                        Z := -PI+(I-1)*PI/M;
                                        Y[I] := F(Z)
                                     end;
                                  OK := true
                               end
                            else writeln ('Number must be a positive integer.')
                         end
                   end
                else
                   begin
                      write ('The program will end so that the function F ');
                      writeln ('can be created. ');
                      OK := false
                   end
             end
      end
   end;
procedure OUTPUT;
   begin
      writeln ('Select output destination ');
      writeln ('1. Screen ');
      writeln ('2. Text file ');
      writeln ('Enter 1 or 2. ');
      readln ( FLAG );
      if ( FLAG = 2 ) then
         begin
            write ('Input the file name in the form - ');
            writeln ('drive:name.ext, ');
            writeln ('for example:   A:OUTPUT.DTA ');
            readln ( NAME );
            assign ( OUP, NAME )
         end
      else assign ( OUP, 'CON');
      rewrite ( OUP );
      writeln(OUP,'FAST FOURIER TRANSFORM');
      writeln(OUP)
   end;
function IBR( J, NU : integer) : integer;
var
   K, I, J2, J1 : integer;
begin
   J1 := J;
   K := 0;
   for I := 1 to NU do
      begin
         J2 := J1 div 2;
         K := 2*K+(J1-2*J2);
         J1 := J2
      end;
   IBR := K
end;
procedure CMULT ( A,B,C,D : real; var E,F : real );
   var
      A1, B1, C1, D1 : real;
   begin
      if ( abs(A) <= ZERO)  then A1 := 0.0
                            else A1 := A;
      if ( abs(B) <= ZERO)  then B1 := 0.0
                            else B1 := B;
      if ( abs(C) <= ZERO)  then C1 := 0.0
                            else C1 := C;
      if ( abs(D) <= ZERO)  then D1 := 0.0
                            else D1 := D;
      E := ( A1*C1 ) - (B1*D1 );
      F := A1*D1 + B1*C1
   end;
procedure CEXP ( A,B : real; var C,D : real );
   var
      E : real;
   begin
      E := exp( A );
      C := E * cos( B );
      D := E * sin( B )
   end;
begin
   INPUT;
   if (OK) then
      begin
         OUTPUT;
         TW := ln(2.0);
{        STEP 1                                                        }
{        use N2 for m, NG for p, NU1 for q, WW for zeta                }
         N2 := N div 2;
{        STEP 2                                                        }
         for I := 1 to N do
            begin
               CR[I] := Y[I];
               CI[I] := 0.0
            end;
         Z := N;
         NG := round(ln(Z)/TW);
         NU1 := NG - 1;
         YR := 0.0;
         YI := 2.0*PI/N;
         CEXP(YR,YI,WWR,WWI);
{        STEP 3                                                        }
         for I := 1 to N2 do
            begin
               XR := 1.0;
               XI := 0.0;
               YR := 1.0;
               YI := 0.0;
               for J := 1 to I do
                  begin
                     CMULT(XR,XI,WWR,WWI,YR,YI);
                     XR := YR;
                     XI := YI
                  end;
               WR[I] := XR;
               WI[I] := XI;
               WR[N2+I] := -XR;
               WI[N2+I] := -XI
            end;
{        STEP 4                                                        }
         K := 0;
{        STEP 5                                                        }
         for L := 1 to NG do
            begin
{              STEP 6                                                  }
               while (K < N-1) do
                  begin
{                    STEP 7                                            }
                     for I := 1 to N2 do
                        begin
{                          STEP 8                                      }
                           Z := exp(NU1*TW);
                           M1 := round(Z);
                           M1 := K div M1;
{                          IBR does the bit reversal                   }
                           NP := IBR(M1,NG);
{                          T1 = T1R + T1i is eta                       }
                           T1R := CR[K+N2+1];
                           T1I := CI[K+N2+1];
{                          STEP 9                                      }
                           if (NP <> 0) then
                              begin
                                 XR := T1R;
                                 XI := T1I;
                                 CMULT(XR,XI,WR[NP],WI[NP],T1R,T1I)
                              end;
                           CR[K+N2+1] := CR[K+1] - T1R;
                           CI[K+N2+1] := CI[K+1] - T1I;
                           CR[K+1] := CR[K+1] + T1R;
                           CI[K+1] := CI[K+1] + T1I;
{                          STEP 10                                     }
                           K := K+1
                        end;
{                    STEP 11                                           }
                     K := K+N2
                  end;
{             STEP 12                                                  }
              K := 0;
              N2 := N2 div 2;
              NU1 := NU1 - 1;
            end;
{        STEP 13                                                       }
         while (K < N-1) do
            begin
{              STEP 14                                                 }
               I := IBR(K,NG);
{              STEP 15                                                 }
               if (I > K) then
                  begin
                     T3R := CR[K+1];
                     CR[K+1] := CR[I+1];
                     CR[I+1] := T3R;
                     T3I := CI[K+1];
                     CI[K+1] := CI[I+1];
                     CI[I+1] := T3I
                  end;
{              STEP 16                                                 }
               K := K+1
            end;
{        STEPS 17 and 18                                               }
         writeln(OUP,'Coefficients c(0), ... , c(2m-1)');
         writeln(OUP);
         for I := 1 to N do
            begin
               YR := 0.0;
               YI := -(I-1)*PI;
               CEXP(YR,YI,XR,XI);
               CMULT(XR,XI,CR[I],CI[I],YR,YI);
               CR[I] := YR/(0.5*N);
               CI[I] := YI/(0.5*N);
               K := I - 1;
               writeln(OUP,K:3,' ',CR[I]:14:8,CI[I]:14:8)
            end;
         writeln(OUP);
         writeln(OUP,'Coefficients a(0), ..., a(m)');
         writeln(OUP);
         for I := 1 to M+1 do writeln(OUP, CR[I]);
         writeln(OUP);
         writeln(OUP,'Coefficients b(1), ..., b(m-1)');
         writeln(OUP);
         for I := 2 to M do writeln(OUP, CI[I]);
         close (OUP)
      end
   end.