Program Fourier;
{$path "ram:include/","pascal:include/" }
{$incl "graphics.lib", "intuition/intuition.h" }
Const
  N = 511;
Type
  MAS = Array [0..512] of Real;
  RPort = ^RastPort;
Var
  Win,Win1,Win2,Win3: ^Window;
  Rp: RPort; { dieser Typ wird in "graphics/rastport.h" definiert }
  err: boolean;
  Msg: Ptr;
  i,nn: integer;
  xr,xi,a,f:MAS;
  m:real;

Procedure InputSignal;
 var i:integer;
 begin
   for i:=0 to N do
    begin
      xr[i]:=round(60*sin(i/8)+random*50);
      xi[i]:=0;
    end;
 end;

Procedure DrawSignal(Rp:RPort, xr:MAS, m:real);
 const
  x0 = 10;
  y0 = 45;
 var i:integer;
 begin
  SetAPen(Rp,3);
  Move(Rp,x0,0);Draw(Rp,x0,70);
  Move(Rp,x0,75-y0);Draw(Rp,570,75-y0);
  SetAPen(Rp,1);
  move(Rp,x0,75-y0-round(xr[0]));
  for i:=0 to N do Draw(Rp,x0+i,75-y0-round(xr[i]*m));
 end;

Procedure DPF(n:real; var  r,i:MAS);
 var a,b:real;p,j:integer;
 gr,gi,hr,hi:MAS;
 begin
   if n=2 then begin a:=r[0]+r[1]; b:=r[0]-r[1];
                     r[0]:=a; r[1]:=b;
                     a:=i[0]+i[1]; b:=i[0]-i[1];
                     i[0]:=a; i[1]:=b;
               end
          else begin
                 for j:=0 to round((n/2)-1) do begin
                                             gr[j]:=r[2*j];gi[j]:=i[2*j];
                                             hr[j]:=r[2*j+1];hi[j]:=i[2*j+1];
                                            end;
                 DPF(n/2,gr,gi);
                 DPF(n/2,hr,hi);
                 for j:=0 to round((n/2)-1) do begin r[j]:=gr[j]+hr[j]*cos(-2*pi*j/n)-hi[j]*sin(-2*pi*j/n);
                                                     i[j]:=gi[j]+hi[j]*cos(-2*pi*j/n)+hr[j]*sin(-2*pi*j/n);
                                        end;
                 p:=round(n/2);
                 for j:=p to round(n-1) do begin r[j]:=gr[j-p]+hr[j-p]*cos(-2*pi*j/n)-hi[j-p]*sin(-2*pi*j/n);
                                                 i[j]:=gi[j-p]+hi[j-p]*cos(-2*pi*j/n)+hr[j-p]*sin(-2*pi*j/n);
                                      end;
               end;
 end;
Begin
  OpenLib(GfxBase, 'graphics.library', 0);

  Win:= Open_Window(0, 0, 640, 256, 1, _CLOSEWINDOW,
        GIMMEZEROZERO+ACTIVATE+WINDOWCLOSE+WINDOWDRAG+WINDOWDEPTH,
        'DPF', Nil, 640, 200, 640, 200);
  Win1:= Open_Window(20,16, 600, 75, 1,0,GIMMEZEROZERO,'Input Signal:', Nil, 600,70, 600,70);
  Win2:= Open_Window(20,96, 600, 75, 1,0,GIMMEZEROZERO,'Amplitude:', Nil, 600,70, 600,70);
  Win3:= Open_Window(20,176, 600, 75, 1,0,GIMMEZEROZERO,'Phase:', Nil, 600,70, 600,70);
  Rp:= Win1^.RPort;
  InputSignal;
  DrawSignal(Rp,xr,6/25);
  Rp:= Win2^.RPort;
  DPF(N+1,xr,xi);
  for i:=0 to n do a[i]:=sqrt(sqr(xr[i])+sqr(xi[i]));
  m:=0;
  for i:=0 to n do if abs(a[i])>m then m:=abs(a[i]);
  m:=25/m;
  DrawSignal(Rp,a,m);
  for i:=0 to n do f[i]:=arctan(xi[i]/xr[i]);
  m:=0;
  for i:=0 to n do if abs(f[i])>m then m:=abs(f[i]);
  m:=25/m;
  DrawSignal(Win3^.RPort,f,m);
  Msg:= Wait_Port(Win^.UserPort);  { aufs Schließen warten }
  Msg:= Get_Msg(Win^.Userport);
  Reply_Msg(Msg);

  Close_Window(Win1);
  Close_Window(Win2);
  Close_Window(Win3);
  Close_Window(Win);
End.