hartley.pas

The frequency domain plays an important role in image processing to smooth, enhance, and detect edg

{-----------------------------------------------------------------------------
     FAST HARTLEY TRANSFORM ROUTINE,
     VERSION 1.05,

     AUTHOR: M.A. O'NEILL,
             DEPARTMENT OF PHOTGRAMMETRY AND SURVEYING,
             UNIVERSITY COLLEGE LONDON,
             GOWER STREET,
             LONDON WC1E 6BT.

     DATED:    22ND NOVEMBER 1986.
     MODIFIED: 3RD MARCH 1987.
     MODIFIED FOR MAC (TML) PASCAL 5TH MARCH 1987.
     MODIFIED FOR TRANSFORM TIMINGS ON MAC 6TH MARCH 1987.
-----------------------------------------------------------------------------}

program fast_Hartley_test(input,output);
uses MacIntf;
label
   1000;
const
   datasize = 512;
type
   direction_type = (forwrd,reverse);
   data_array_type = array[1..datasize] of real;
var
   opt_set: Boolean;
   dir, test_option: char;
   i,j,syze,iter,demo: integer;
   data_array: data_array_type;
   transform_direction: direction_type;
   in_file,out_file,time_file: string;
   writefile: text;
   start_time: DateTimeRec;
   
{---------------------------------------------------------------------------
   External System bindings to Macintosh OS here ...
---------------------------------------------------------------------------

procedure SetTime(d: DateTimeRec);
EXTERNAL;}

{--------------------------------------------------------------------------

procedure GetTime(var d: DateTimeRec);
EXTERNAL;}

{--------------------------------------------------------------------------}

procedure initialise_stopwatch;

begin
   GetTime(start_time);
end;

{--------------------------------------------------------------------------}

procedure read_stopwatch;
var
   end_time: DateTimeRec;
   elapsed_time: DateTimeRec;
begin
   with elapsed_time do
   begin
      GetTime(end_time);
      Hour := end_time.Hour - start_time.Hour;
      Minute := end_time.Minute - start_time.Minute;
      Second := end_time.Second - start_time.Second;
      writeln(writefile,syze,'    ',Hour*3600+Minute*60+Second);                                     
   end;
end;

{-----------------------------------------------------------------------------
   Fast Hartley transform routine ...
   version 1.00, Dated 22nd November 1986,

   transform := foward; Forward transform from time frequency domain.
   transform := reverse; Reverse transform from frequency to time domain.

   power_index: Index to which 2 must be raised to generate a transform
                containing 'syze' elements:
                if syze = 8  then power_index = 3;
                if syze = 16 then power_index = 4; etc etc ...

   syze:  Number of element in the input data array.
-----------------------------------------------------------------------------}

procedure fht(var data_array: data_array_type;
                              power_index,
                              syze: integer;
                              transform: direction_type);

var
   i,
   j,
   k,
   trg_ind,
   trg_inc,
   power,
   t_a,
   f_a,
   i_temp,
   section,
   s_start,
   s_end: integer;
   sne,csn: array[1..datasize] of real;
   accu: array[1..2,1..datasize] of real;

{-----------------------------------------------------------------------------
   Permutation routine. This routine re-orders the data before the butterly 
   transform routine is called ...
----------------------------------------------------------------------------}

function permute(index: integer): integer;
var
  i,j,s: integer;
begin
   j := 0;
   index := index - 1;
   for i := 1 to power_index do
   begin
      s := index div 2;
      j := j + j + index - s - s;
      index := s;
   end;
   permute := j + 1;
end;

{-----------------------------------------------------------------------------
   Calculate the trignometric functions required by the FHT and store values.
   For a N point transform, the trignometric functions will be calculated  at
   intervals of Nths of a turn ...
-----------------------------------------------------------------------------}

procedure trig_table(npts: integer);
const
   pi = 3.14159265;
var
   i: integer;
   angle,omega: real;
begin
   angle := 0;
   omega := 2 * pi / npts;
   for i := 1 to npts do
   begin
      sne[i] := sin(angle);
      csn[i] := cos(angle);
      angle := angle + omega;
   end;
end;

{-----------------------------------------------------------------------------
   Calculate the address of the retrograde index for the sine term for the
   dual place alogrithm, if it is required ...
-----------------------------------------------------------------------------}

function modify(power,s_start,s_end,index: integer): integer;
begin
   if (s_start = index) or (power < 3) then
      modify := index
   else
      modify := s_start + s_end - index + 1;
end;

{------------------------------------------------------------------------------
   Butterfly transform an index pair ...
------------------------------------------------------------------------------}

procedure butterfly(trig_ind,i_1,i_2,i_3: integer);
begin
   accu[t_a,i_1] := accu[f_a,i_1] + accu[f_a,i_2] * csn[trig_ind] +
                                     accu[f_a,i_3] * sne[trig_ind];
   trig_ind := trig_ind + syze div 2;
   accu[t_a,i_2] := accu[f_a,i_1] + accu[f_a,i_2] * csn[trig_ind] +
                                     accu[f_a,i_3] * sne[trig_ind];
end;

{------------------------------------------------------------------------------
   Main program for the fast Hartley transform.
------------------------------------------------------------------------------}

begin

   if test_option = 't' then initialise_stopwatch;
   power := 1;
   f_a := 1;
   t_a := 2;
   trig_table(syze);
   for i := 1 to syze do accu[f_a,permute(i)] := data_array[i];

{------------------------------------------------------------------------------
   Start of the Hartley butterfly transform ...
------------------------------------------------------------------------------}

   for i := 1 to power_index do
   begin
      j := 1;
      section := 1;
      trg_inc := syze div (power + power);
      repeat
         trg_ind := 1;
         s_start := section * power + 1;
         s_end   := (section + 1) * power;
         for k := 1 to power do
         begin
            butterfly(trg_ind,j,j + power,
                      modify(power,s_start,s_end,j + power));
            trg_ind := trg_ind + trg_inc;
            j := j + 1;
         end;
         j := j + power;
         section := section + 2;
      until j > syze;
      power := power + power;
      i_temp := t_a;
      t_a    := f_a;
      f_a    := i_temp;
   end;

{------------------------------------------------------------------------------
   End of Hartley butterfly. The results are scaled in necessary, and then
   placed in back into the array data ...
------------------------------------------------------------------------------}

   case transform_direction of
        forwrd: for i := 1 to syze do data_array[i] := accu[f_a,i] / syze;
        reverse: for i := 1 to syze do data_array[i] := accu[f_a,i];
   end;
   if test_option = 't' then
      read_stopwatch
   else
   begin
      open(writefile,out_file);
      rewrite(writefile);
      for i := 1 to syze do
          if test_option = 's' then
             writeln(writefile,i:1,'   ',data_array[i]);
      close(writefile);
   end;
end;

{------------------------------------------------------------------------------
   Main test program here ...
------------------------------------------------------------------------------}

begin
   writeln('[Fast Hartley Transform Demonstrator (MAC) Ver. 1.05]');
   writeln('(C) M.A. O''Neill 1987');
   repeat
      opt_set := true;
      write('Foward ''f'' or  reverse ''r'' transform > '); readln(dir);
      if dir = 'q' then
      begin
         writeln('Quitting FHT demonstrator');
         goto 1000;
      end;
      if (dir <> 'f') and (dir <> 'r') then
      begin
         opt_set := false;
         writeln('You must type ''f'', ''r'' or ''q''');
      end
      else
      begin
         if dir = 'f' then transform_direction := forwrd;
         if dir = 'r' then transform_direction := reverse;
      end;
   until opt_set;

   repeat
      opt_set := true;
      write('Demonstration function (1 or 2) > '); readln(demo);
      if (demo < 1) or (demo > 2) then
      begin
         opt_set := false;
         writeln('Options 1 or 2 only allowed');
      end;
   until opt_set;

{--------------------------------------------------------------------------
   If test_option = 's' then the Hartley Transform is computed for a
   specific number of points (which must be a power of two) ...

   Alternatively, if 't' is specified, the user is prompted for an
   iteration count. In this case, the FHT is iteratively computed for
   an increasing number of data points, and the time required to transform
   an N point dataset is written to file.
 -------------------------------------------------------------------------}

   repeat
      opt_set := true;
      write('Specimen transform (s) or timing test (t) > ');
      readln(test_option);
      if test_option = 's' then   
      begin
         writeln('Generating specimen Fast Hartley Transform ...');
         writeln;
         write('Power index > '); readln(iter);
         write('Input data file > '); readln(in_file);
         write('Output data file > '); readln(out_file);
      end;

      if test_option = 't' then
      begin
         writeln('Fast Hartley Transform timing test ...'); writeln;
         write('How many iterations > '); readln(iter);
         write('Timings file > '); readln(time_file);
         open(writefile,time_file);
         rewrite(writefile);
      end;

      if (test_option <> 's') and (test_option <> 't') then 
      begin
         opt_set := false;
         writeln('You must type ''s'' or ''t''');
      end;
   until opt_set;

{------------------------------------------------------------------------------
   Enter the function to be transformed here ...
   The examples given here are top hat function (1) or Y = X (2) ...
------------------------------------------------------------------------------}

   case test_option of
        's': writeln('Generating transform function ...');
        't': writeln('Running timing test ...');
   end;
   writeln;
   for j := 1 to iter do
   begin
      syze := round(exp(j * ln(2)));
      if test_option = 't' then writeln('Iteration: ',j:1,' Size: ',syze:1);
      if demo = 1 then
      begin
         if test_option = 's' then
         begin
            open(writefile,in_file);
            rewrite(writefile);
         end;
         for i := 1 to syze div 2 do
         begin
            data_array[i] := 1;
            if test_option = 's' then writeln(writefile,i,'  ',data_array[i]);
         end;
         for i := syze div 2 + 1 to syze do
         begin
            data_array[i] := -1;
            if test_option = 's' then writeln(writefile,i,'  ',data_array[i]);
         end;
         if test_option = 's' then close(writefile);
      end;

      if demo = 2 then
      begin
         if test_option = 's' then
         begin
            open(writefile,in_file);
            rewrite(writefile);
         end;
         for i := 1 to syze do
         begin
            data_array[i] := i;
            if test_option = 's' then writeln(writefile,i,'  ',data_array[i]);
         end;
        if test_option = 's' then close(writefile);
     end;

{------------------------------------------------------------------------------
   Perform transform ...
------------------------------------------------------------------------------}

      fht(data_array,j,syze,transform_direction);
   end;
   if test_option = 't' then close(writefile);
   write('*** Finished. Press any key to return to desktop ...');
   readln;
   writeln('Quitting to desktop ...');

1000: end.