home *** CD-ROM | disk | FTP | other *** search
/ HAM Radio 1 / HamRadio.cdr / tech / design2 / fft.pas < prev    next >
Encoding:
Pascal/Delphi Source File  |  1986-07-22  |  11.0 KB  |  379 lines
  1. program fft;
  2. {        Forward & Inverse FFT Program     --     Jeff Crawford
  3.  
  4. Reference:  " Add the FFT to Your Box of Design Tools", Robert Cushman,
  5.             EDN, September 16, 1981.
  6.  
  7.          Program Calculates Both Forward and Inverse FFT's.  Data May be
  8.          Entered Manually, Point by Point, Read off Disk, or Calculated by
  9.          a Program FFT.In and Inserted in Input Arrays.
  10.  
  11.          Plotting of Both Input & Output Waveforms is Provided Through the
  12.          Use of Turbo Graphix Toolbox.  At Present, the Largest FFT Capable of
  13.          Using the Plotting Package is a 256 Point FFT.
  14.                                                                             }
  15.  
  16. {$ITypedef.sys}
  17. {$Igraphix.sys}
  18. {$Ikernel.sys}
  19. {$iwindows.sys}
  20. {$Iaxis.hgh}
  21. {$IPolygon.hgh}
  22. {$IModpoly.hgh}
  23.  
  24. type    vector =        array[0..256] of real;
  25.         XYdata = record
  26.                  XX,YY: real;
  27.                  end;
  28.  
  29. var      Ans,Anz,An     : string[1];
  30.          k              : real;
  31.          N              : integer;
  32.          D1,D2,X1,Y1    : vector;
  33.          N3,N4,N6,N7,N8 : integer;
  34.          N1,N2          : integer;
  35.          a,c,s,t1,t2    : real; { Used in FFt }
  36.          i              : integer;
  37.          file_name      : string[10];
  38.          write_file     : file of XYdata;
  39.          read_file      : file of XYdata;
  40.          file_rec       : XYdata;
  41.          OK             : Boolean;
  42.          FI             : integer; { Flagged for Forward or Inverse FFT }
  43.          PlotTr,PlotTi  : PlotArray; { Arrays for Plotting Input }
  44.          Graph_title    : string[25];
  45.          xloc,yloc      : integer;
  46.          Title          : string[15];
  47.          PlotFr,PlotFi  : PlotArray;  { Arrays for Plotting Output }
  48.          tinc           : real; { Used in FFT_Input Subroutine }
  49.          jm             : integer; { Provides Plotting for > 100 Points }
  50.  
  51. label 10,20,30,99;
  52.  
  53. {$IFFT.In} { Program Automatically Fills Input Arrays }
  54.  
  55. Procedure PlotG( A,B:PlotArray);   { PlotArray Defined as Type in Graphics.SYs}
  56.  
  57. var      ymin,ymax   : real;   { For Real Components }
  58.          ymin1,ymax1 : real; { For Imag Components }
  59.          i           : integer;
  60.          Colornum    : integer;
  61.  
  62.  
  63. begin
  64.     Colornum:= 16;          { Sets Color of Graphics Content }
  65.     N:= N div jm;
  66.     initgraphic;
  67.     ymin:= 1.0E20;  ymax:= -1.0E20;
  68.          for i:= 1 to N do
  69.               begin
  70.               if A[i,2] > ymax then ymax:= A[i,2];
  71.               if A[i,2] < ymin then ymin:= A[i,2];
  72.               end;
  73.     ymin1:= 1.0E20;  ymax1:= -1.0E20;
  74.          for i:= 1 to N do
  75.               begin
  76.               if B[i,2] > ymax1 then ymax1:= B[i,2];
  77.               if B[i,2] < ymin1 then ymin1:= B[i,2];
  78.               end;
  79.     if ymax < ymax1 then ymax:= ymax1;
  80.     if ymin > ymin1 then ymin:= ymin1;
  81.     Definewindow(1,2,0,Xmaxglb,Ymaxglb);
  82.     SetforegroundColor(Colornum);
  83.     Defineworld(1,A[1,1],Ymax,A[N,1],Ymin);
  84.     Defineheader(1,Graph_title);
  85.     setheaderon;
  86.     setbackground(0);
  87.     Selectworld(1);
  88.     Selectwindow(1);
  89.     Setlinestyle(0);
  90.     Drawborder;
  91.     DrawAxis(7,7,8,2,2,10,0,0,false);
  92.     DrawPolygon(A,1,N,4,1,5);
  93.     SetLinestyle(1);
  94.     Drawaxis(7,7,8,2,2,10,0,-1,false); { Must Recall as per Graphix Book p128}
  95.     DrawPolygon(B,1,-N,3,1,5); { Draws Imaginary Portion }
  96.     SetLinestyle(0);
  97.     SetClippingOff;
  98.     Drawtext(300,190,1,'  C  E  L  L    #  ' );
  99.     xloc:= 50;
  100.     yloc:= 25;
  101.     Title[1]:= 'M';
  102.     Title[2]:= 'A';
  103.     Title[3]:= 'G';
  104.     Title[4]:= 'N';
  105.     Title[5]:= 'I';
  106.     Title[6]:= 'T';
  107.     Title[7]:= 'U';
  108.     Title[8]:= 'D';
  109.     Title[9]:= 'E';
  110.     for i:= 1 to 9 do
  111.          begin
  112.          Drawtext(xloc,yloc+i*15,1,Title[i]);
  113.          end;
  114.     SetClippingOn;
  115.     repeat until keypressed;
  116.     leavegraphic;
  117.     N:= N * jm;
  118. end;
  119.  
  120. procedure Read_Data;
  121.     label 2;
  122.  
  123.     begin
  124.   2:clrscr;
  125.     gotoXY(10,10);
  126.     write('Enter File Name  ');
  127.     {$I-} readln(file_name)  {$I+};
  128.     OK:= (IOresult = 0);
  129.     if NOT OK then goto 2;
  130.     assign(read_file,file_name);
  131.     {$I-} Reset(read_file) {$I+};  { Flags if Disk is Non-Existent }
  132.     OK:= (IOResult = 0);
  133.     if not OK then
  134.          begin
  135.          clrscr;
  136.          gotoXY(10,10);
  137.          writeln('File "',file_name,'" Is Not Available  ');
  138.          writeln;
  139.          Delay(1000);
  140.          clrscr;
  141.          goto 2;
  142.          end;
  143.     i:= 0;
  144.     clrscr;
  145.     writeln('                           Filename:     ',file_name:10);
  146.     writeln;
  147.     writeln('Index            Real              Imaginary');
  148.     writeln;
  149.     while not EOF(read_file) do { If File is Available, Reads From Disk }
  150.          begin
  151.          read(read_file,file_rec);
  152.          with file_rec do
  153.               begin
  154.               D1[i]:= XX;
  155.               D2[i]:= YY;
  156.               writeln('   ',i+1,'       ',D1[i]:12:6,'       ',D2[i]:12:6);
  157.               i:= i + 1;
  158.               end;
  159.          end;
  160.     N:= i;
  161.     close(read_file);
  162. end;
  163.  
  164. procedure Write_Data( X,Y: vector);
  165.  
  166. var           Ans1,Ans2: string[1];
  167.  
  168. label 4;
  169. begin
  170. 4:  clrscr;
  171.     gotoXY(10,10);
  172.     write('Write Data to What File  ?');
  173.     readln(file_name);
  174.     assign(write_file,file_name);
  175.     {$I-} Reset(write_file) {$I+};
  176.     OK:= (Ioresult = 0); { Determines if Write-To File Already Exists }
  177.     if OK then
  178.          begin
  179.          clrscr;
  180.          gotoXY(10,10);
  181.          write('File Already Exists. . . Overwrite ? Y / N  ? ');
  182.          readln(Ans2);
  183.          Ans2:= upcase(Ans2);
  184.          if Ans2 <> 'Y' then goto 4;
  185.          end;
  186.     assign(write_file,file_name);
  187.     rewrite(write_file);
  188.     clrscr;
  189.     writeln('                     Filename:   ',file_name:12);
  190.     writeln;
  191.     writeln('Index            Real              Imaginary');
  192.     writeln;
  193.     with file_rec do
  194.          begin
  195.          for i:= 0 to N-1 do { Writes to Disk File "write_file"}
  196.               begin
  197.               writeln(' ',i+1,'         ',X[i]:12:6,'         ',Y[i]:12:6);
  198.               XX:= X[i];
  199.               YY:= Y[i];
  200.               write(write_file,file_rec);
  201.               end;
  202.          end;
  203.     close(write_file);
  204. end;
  205.  
  206. function XtoY(x,y: real): real;
  207.     begin
  208.     XtoY:= exp(y*Ln(x));
  209.     end;
  210.  
  211. begin
  212. 99: clrscr;
  213.     gotoxy(25,10);
  214.     write('     Fft     Ift    "I" or "F" ');
  215.     for i:= 0 to 256 do
  216.          begin
  217.          D1[i]:= 0.0;
  218.          D2[i]:= 0.0;
  219.          end;
  220.     readln(An);
  221.     An:= upcase(An);
  222.     if An = 'F' then FI:= 1 else FI:= -1;
  223.     gotoxy(5,10);
  224.     clrscr;
  225.     write(' Number of Steps  FFT/IFT  ');
  226.     readln(N);
  227.     k:= Trunc(Ln(N)/Ln(2));
  228.     clrscr;
  229.     writeln;
  230.     write(' Read Diskfile = F     Manual Input = M  ');
  231.     readln(Anz);
  232.     Anz:= upcase(Anz);
  233.     if Anz = 'M' then
  234.          begin
  235.          clrscr;
  236.          writeln;
  237.          write('        M - Manual Input     F - File Program Entry ');
  238.          readln(Ans);
  239.          Ans:= upcase(Ans);
  240.          if Ans = 'M' then
  241.               begin
  242.               writeln('Input Data ');
  243.               writeln;
  244.               writeln('   Data #     Real   Imag ');
  245.               tinc:= N;
  246.               tinc:= 1/tinc;
  247.               for i:= 0 to N-1 do
  248.                    begin
  249.                    write('     ',i+1,'     ');
  250.                    readln(D1[i], D2[i]);
  251.                    end;
  252.               end;
  253.          if Ans = 'F' then
  254.               begin
  255.               Input_Data;
  256.               end;
  257.          end
  258.          else Read_Data;
  259. writeln;
  260. write('                         Write Data to Disk  Y/N ?  ');
  261. readln(Anz);
  262. Anz:= upcase(Anz);
  263. if Anz = 'Y' then Write_Data(D1,D2);
  264. writeln;
  265. write('   Plot of Input Data  [real]  Y/N  ?  ');
  266. readln(Anz);
  267. Anz:= upcase (Anz);
  268. if Anz = 'Y' then
  269.     begin
  270.     clrscr;
  271.     gotoxy(10,10);
  272.     write('Input Graph Title  ');
  273.     readln(Graph_title);
  274.     if N > 100 then jm:= 4 else jm:= 1;
  275.     i:= 1;
  276.     while i <= N div jm do
  277.          begin
  278.          PlotTr[i,1]:= i*jm;
  279.          PlotTr[i,2]:= D1[(i-1)*jm];
  280.          PlotTi[i,1]:= i*jm;
  281.          PlotTi[i,2]:= D2[(i-1)*jm];
  282.          i:= i + 1;
  283.          end;
  284.     PlotG(PlotTr,PlotTi);
  285.     end;
  286.     N1:= 0;
  287.     N2:= N-1;
  288.     for N3:= 1 to N2 do
  289.          begin
  290.          N4:= N;
  291. 10:      N4:= N4 div 2;
  292.          if N1 + N4 > N2 then goto 10;
  293.          N1:= N1 - ( N1 div N4 ) * N4 + N4;
  294.          if N1 <= N3 then goto 20;
  295.          t1:= D1[N3];
  296.          D1[N3]:= D1[N1];
  297.          D1[N1]:= t1;
  298.          t2:= D2[N3];
  299.          D2[N3]:= D2[N1];
  300.          D2[N1]:= t2;
  301. 20:      end;
  302.          N4:= 1;
  303. 30:      N6:= 2*N4;
  304.          for N3:= 0 to N4 - 1 do
  305.               begin
  306.               a:= Fi*N3*Pi/N4;
  307.               c:= Cos(a);
  308.               s:= Sin(a);
  309.               N7:= N3;
  310.               while N7 <= N-1 do
  311.                    begin
  312.                    N8:= N7 + N4;
  313.                    t1:= c*D1[N8] - s*D2[N8];
  314.                    t2:= c*D2[N8] + s*D1[N8];
  315.                    D1[N8]:= D1[N7] - t1;
  316.                    D2[N8]:= D2[N7] - t2;
  317.                    D1[N7]:= D1[N7] + t1;
  318.                    D2[N7]:= D2[N7] + t2;
  319.                    N7:= N7 + N6;
  320.                    end;
  321.               end;
  322.               N4:= N6;
  323.               if N4 < N then goto 30;
  324.          clrscr;
  325.          write('                   Output Results ');
  326.          if Fi= 1 then writeln('--   Forward Transform');
  327.          if Fi= -1 then writeln('--  Inverse Transform');
  328.          writeln;
  329.          writeln('Index         Real           Imaginary');
  330.          writeln;
  331.          if An = 'F' then
  332.               begin
  333.               for i:= 0 to N-1 do
  334.                    begin
  335.                    X1[i]:= D1[i]/N; Y1[i]:= D2[i]/N;
  336.                    writeln('   ',i+1,'    ',X1[i]:12:6,'    ',Y1[i]:12:6);
  337.                    end;
  338.               end;
  339.          if An = 'I' then
  340.               begin
  341.               for i:= 0 to N-1 do
  342.                    begin
  343.                    X1[i]:= D1[i]; Y1[i]:= D2[i];
  344.                    writeln('   ',i+1,'    ',X1[i]:12:6,'     ',Y1[i]:12:6);
  345.                    end;
  346.               end;
  347.          gotoxy(60,24);
  348.          write('Any Key to Cont');
  349.          repeat until keypressed;
  350.          clrscr;
  351.          gotoxy(10,10);
  352.          write('                Write Data to Disk  Y/N ?  ');
  353.          readln(Anz);
  354.          Anz:= upcase(Anz);
  355.          if Anz = 'Y' then Write_Data(X1,Y1);
  356.          writeln;
  357.          write('                Plot of Input Data  [real]  Y/N  ?  ');
  358.          readln(Anz);
  359.          Anz:= upcase (Anz);
  360.          if Anz = 'Y' then
  361.               begin
  362.               clrscr;
  363.               gotoxy(10,10);
  364.               write('Input Graph Title  ');
  365.               readln(Graph_title);
  366.               if N > 100 then jm:= 4 else jm:= 1;
  367.               i:= 1;
  368.               while i <= N div jm do
  369.                    begin
  370.                    PlotFr[i,1]:= i*jm;
  371.                    PlotFr[i,2]:= X1[(i-1)*jm];
  372.                    PlotFi[i,1]:= i*jm;
  373.                    PlotFi[i,2]:= Y1[(i-1)*jm];
  374.                    i:= i + 1;
  375.               end;
  376.               PlotG(PlotFr,PlotFi);
  377.               end;
  378.          goto 99;
  379.    end.$