home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / PASCAL / MADTRB16.ZIP / FFT.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1985-02-19  |  4.7 KB  |  166 lines
  1. program fft;
  2. type
  3.    bigstr   =  string[12];
  4.    range    =  array[1..256] of real;
  5. var
  6.    filename,nstr,fileout                :bigstr;
  7.    filnam,fiout                         :string[8];
  8.    xform                                :string[25];
  9.    inverse,output                       :char;
  10.    npoints,invt,d                       :integer;
  11.    reel,imaginary                       :range;
  12.  
  13. procedure readdata(filename:bigstr;var x,y :range;var a:integer);
  14.  
  15. {range and bigstr are set up in the TYPE identifier
  16. range  = array[1..?] of real
  17. bigstr = string[12]
  18. a returns the transform length or file size}
  19.  
  20. type
  21.     data  =
  22.       record
  23.         w,z :real;
  24.       end;
  25. var
  26.    index    :integer;
  27.    xy       :array[1..256] of data;
  28.    datafile :file of data;
  29. begin
  30.      assign(datafile,filename);
  31.      reset(datafile);
  32.      a:=filesize(datafile);
  33.      for index:=1 to a do begin
  34.          read(datafile,xy[index]);
  35.          with xy[index] do begin
  36.               x[index]:=w;
  37.               y[index]:=z;
  38.          end;
  39.      end;
  40.      close(datafile);
  41. end;
  42.  
  43. procedure writedata(filename:bigstr;var x,y:range;var a:integer);
  44. type
  45.     data  =
  46.       record
  47.         w,z :real;
  48.       end;
  49. var
  50.    index    :integer;
  51.    xy       :array[1..256] of data;
  52.    datafile :file of data;
  53. begin
  54.      assign(datafile,filename);
  55.      rewrite(datafile);
  56.      for index:=1 to a do begin
  57.          with xy[index] do begin
  58.               w:=x[index];
  59.               z:=y[index];
  60.          end;
  61.      write(datafile,xy[index]);
  62.      end;
  63.      close(datafile);
  64. end;
  65.  
  66. procedure fft(inv,n:integer;var re,imag:range);
  67.  
  68. {range is set up in the TYPE identifier and is = array[1..?] of real
  69. inv  = 1 if direct transform and =-1 for inverse
  70. n    = transform length power of 2}
  71.  
  72. var
  73.    le,ip,m,i,l,nm1,j,t,le1              :integer;
  74.    ur,ui,tr,ti,tmpr,tmpi,nv2,k          :real;
  75.  
  76. begin
  77.      m:=round(ln(n)/ln(2));
  78.      if abs(m-int(m))=0 then begin
  79.         nv2:=n/2;
  80.         nm1:=n-1;
  81.         j:=1;
  82.         for i:=1 to nm1 do begin
  83.             if i<j then begin
  84.                tmpr:=re[j];tmpi:=imag[j];
  85.                re[j]:=re[i];imag[j]:=imag[i];
  86.                re[i]:=tmpr;imag[i]:=tmpi;
  87.             end;
  88.             k:=nv2;
  89.             while k<j do begin
  90.                   j:=round(j-k);
  91.                   k:=k/2;
  92.             end;
  93.             j:=round(j+k);
  94.         end;
  95.         for l:=1 to m do begin
  96.             le:=1;
  97.             for t:=1 to l do
  98.                 le:=le*2;
  99.             le1:=round(le/2);
  100.             ur:=1;ui:=0;
  101.             for j:=1 to le1 do begin
  102.                 i:=j;
  103.                 while i<=n do begin
  104.                       ip:=i+le1;
  105.                       tr:=re[ip]*ur-imag[ip]*ui;ti:=imag[ip]*ur+re[ip]*ui;
  106.                       re[ip]:=re[i]-tr;imag[ip]:=imag[i]-ti;
  107.                       re[i]:=re[i]+tr;imag[i]:=imag[i]+ti;
  108.                       i:=i+le;
  109.                 end;
  110.                 ur:=cos(pi*j/le1);ui:=-inv*sin(pi*j/le1);
  111.             end;
  112.         end;
  113.         if inv=-1 then begin
  114.            for i:=1 to n do begin
  115.                re[i]:=re[i]/n;
  116.                imag[i]:=imag[i]/n;
  117.            end;
  118.         end;
  119.      end;
  120. end;
  121.  
  122. begin
  123.      gotoxy(29,1);writeln('Fast Fourier Transform');
  124.      gotoxy(39,3);writeln('by');
  125.      gotoxy(31,4);writeln('Michael F. Griffin');
  126.      writeln;writeln;writeln;
  127.      gotoxy(19,6);writeln('*----------------------------------------*');
  128.      gotoxy(8,8);writeln('- All file names are assumed to have .DAT as the extension.');
  129.      gotoxy(8,9);writeln('- Output goes to either a file or a printer.');
  130.      writeln;writeln;
  131.      write('Enter data file name w/o extension : ');
  132.      readln(filnam);
  133.      filename:=filnam+'.dat';
  134.      writeln('Input File name ==> : ',filename);
  135.      write('Do you want a direct transform <Y/N>?');
  136.      read(kbd,inverse);writeln;
  137.      if (inverse='Y') or (inverse='y') then begin
  138.         xform:='Direct Fourier Transform';
  139.         invt:=1
  140.      end
  141.      else begin
  142.         xform:='Inverse Fourier Transform';
  143.         invt:=-1
  144.      end;
  145.      write('Do you want the output to a file <Y/N>?');
  146.      read(kbd,output);writeln;
  147.      if (output='Y') or (output='y') then begin
  148.         write('Enter output data file name w/o extension : ');
  149.         readln(fiout);
  150.         fileout:=fiout+'.dat';
  151.         writeln('Output File name ==> : ',fileout);
  152.      end;
  153.      readdata(filename,reel,imaginary,npoints);
  154.      fft(invt,npoints,reel,imaginary);
  155.      if (output='Y') or (output='y') then
  156.         writedata(fileout,reel,imaginary,npoints)
  157.      else begin
  158.         writeln(lst,xform);writeln(lst);
  159.         str(npoints,nstr);
  160.         for d:=1 to npoints do
  161.             writeln(lst,'(',d:length(nstr),')',reel[d],imaginary[d])
  162.      end;
  163. end.
  164.  
  165.  
  166.