implementation
//PROGRAM D10R1
//Driver for routine FOUR1
uses
  unit2;
  {$R *.DFM}
procedure PRNTFT(DATA:array of real;  NN2:integer);
const
  s1 = '%13.6f'; s2 = '#0';
var
  NN,N,M,MM:integer;
begin
  Writeln(F,'   n    Real H(n)     Imag H(n)   Real H(N-n)  Imag H(N-n)');
  Writeln(F,'  ',FormatFloat('0',0),Format(s1,[DATA[1]]),
    Format(s1,[DATA[2]]),Format(s1,[DATA[1]]),Format(s1,[DATA[2]]));
  For NN:=2 To (((NN2 div 2) + 2) div 2) do
  begin
    N:= 2 * NN -1;
    M:= (N-1) div 2;
    MM:=NN2 + 2 - N;
    Writeln(F,'  ',FormatFloat(s2,M),Format(s1,[DATA[N]]),
     Format(s1,[DATA[N+1]]),Format(s1,[DATA[MM]]),Format(s1,[DATA[MM+1]]));
  end;
end;
procedure TForm1.Button1Click(Sender: TObject);
const
  s1='%13.6f'; s2 = '##';  NN = 32;  NN2 = 2 * NN;
var
  DATA,DCMP:array[0..64] of real;
  I,II,ISIGN,J:integer;  
begin
  //輸出計算結果到文件
  AssignFile(F, 'd:\delphi_shu\p10\d10r1.dat');
  Rewrite(F);
  Writeln(F, 'h(t)=real-valued even-function');
  Writeln(F, 'H(n)=H(N-n) and real?');
  For II:=1 To NN do
  begin
    I:= II * 2 - 1;
    DATA[I]:=1 / (Sqr((I - NN - 1) / NN) + 1);
    DATA[I + 1]:=0;
  end;
  ISIGN:=1;
  FOUR1(DATA, NN, ISIGN);
  PRNTFT(DATA, NN2);
  Writeln(F, 'h(t)=imagnary-valued even-function');
  Writeln(F, 'H(n)=H(N-n) and imaginary?');
  For II:=1 To  NN  do
  begin
    I:= II * 2 - 1;
    DATA[I + 1]:=1 /(Sqr((I - NN - 1) / NN) + 1);
    DATA[I]:=0;
  end;
  ISIGN:=1;
  FOUR1(DATA, NN, ISIGN);
  PRNTFT(DATA, NN2);
  Writeln(F, 'h(t)=imagnary-valued even-function');
  Writeln(F, 'H(n)=H(N-n) and imaginary?');
  For II:=1 To  NN  do
  begin
    I:= II * 2 - 1;
    DATA[I]:=(I - NN - 1) / NN /(Sqr((I - NN - 1) / NN)+ 1);
    DATA[I + 1]:=0;
  end;
  DATA[1]:=0;
  ISIGN:=1;
  FOUR1(DATA, NN, ISIGN);
  PRNTFT(DATA, NN2);
  Writeln(F, 'h(t)=imagnary-valued odd-function');
  Writeln(F, 'H(n)=-H(N-n) and real?');
  For II:=1 To  NN  do
  begin
    I:= II * 2 - 1;
    DATA[I + 1]:=(I - NN - 1) / NN / (Sqr((I - NN - 1) / NN) + 1);
    DATA[I]:=0;
  end;
  DATA[2]:=0;
  ISIGN:=1;
  FOUR1(DATA, NN, ISIGN);
  PRNTFT(DATA, NN2);
  //Transrorm, inverse-transform test
  For II:=1 To  NN  do
  begin
    I:= II * 2 - 1;
    DATA[I]:=1 / (Sqr(0.5 * (I - NN - 1) / NN) + 1);
    DCMP[I]:=DATA[I];
    DATA[I + 1]:=(0.25 * (I - NN - 1) / NN);
    DATA[I + 1]:=DATA[I + 1] * Exp(-Sqr(0.5 * (I - NN - 1) / NN));
    DCMP[I + 1]:=DATA[I + 1];
  end;
  ISIGN:=1;
  FOUR1(DATA, NN, ISIGN);
  ISIGN:=-1;
  FOUR1(DATA, NN, ISIGN);
  Writeln(F,'       Double Fourier Transform:           Original Data:');
  Writeln(F,'   k      Real h(k)   Imag h(k)    Real h(k)   Imag h(k)');
  For II:=1 To NN div 2 do
  begin
    I:= II * 2 - 1;
    J:=(I + 1) div 2;
    Writeln(F,'  ', FormatFloat(S2,J),Format(S1,[DCMP[I]]),
                   Format(S1,[DCMP[I + 1]]),Format(S1,[DATA[I] / NN]),
                   Format(S1,[DATA[I + 1] / NN]));
  end;
  CloseFile(F);
  //屏幕顯示計算結果
  memo1.Lines.LoadFromFile('d:\delphi_shu\p10\d10r1.dat');
end;