From news.kth.se!eru.mt.luth.se!bloom-beacon.mit.edu!gatech!udel!news.mathworks.com!news.ultranet.com!news.sprintlink.net!news.dgsys.com!news Fri Mar 31 08:22:51 1995
Path: news.kth.se!eru.mt.luth.se!bloom-beacon.mit.edu!gatech!udel!news.mathworks.com!news.ultranet.com!news.sprintlink.net!news.dgsys.com!news
From: foo@bar (Happy Hacker)
Newsgroups: comp.dsp
Subject: Re: Fast Hartley Transform - fhtfoo.txt [1/1]
Date: 28 Mar 1995 03:46:14 GMT
Organization: Digital Gateway Systems
Lines: 200
Distribution: world
Message-ID: <3l80q6$560@news.dgsys.com>
References: <3kpog7$frn@ninurta.fer.uni-lj.si>
NNTP-Posting-Host: rhp.com
Mime-Version: 1.0
Content-Type: multipart/mixed;
     Boundary="*-*-*- Next Section -*-*-*"
X-Newsreader: WinVN 0.93.11

--*-*-*- Next Section -*-*-*

In article <3kpog7$frn@ninurta.fer.uni-lj.si>, matic@gea says...
>
>
>Hi.
>
>A am looking for info about Fast Hartley Transform. Any response
>is highly appreciated.
>
>MatiC.
>
>
--*-*-*- Next Section -*-*-*

{...
Below is a Pascal program of a fast Hartley transform... If my
memory serves me it was originally derived from a BYTE magazine
article on the Hartley transform about... 10 years ago, sorry
that I do not have the exact reference.  As I recall the article
gave quite a good explanation of the theory.

Hope this is of some use to you.

Cheers!  Paul Williams
..}

PROGRAM FTHTEST(INPUT,OUTPUT);

CONST
  BUFFER_MAX = 256;
  DISKBLOCKSIZE = 256;
  NPTS = 16;
  ORDER = 4;

TYPE
  INTEGER_BUFFER = ARRAY[1..DISKBLOCKSIZE] OF INTEGER;
  REAL_BUFFER = ARRAY[1..BUFFER_MAX] OF REAL;

VAR
  DAFILE : FILE OF INTEGER_BUFFER;
  IBUFFER : INTEGER_BUFFER;
  ACCU: ARRAY[1..2,1..BUFFER_MAX] OF REAL;

  PI : REAL;
  INFILE : TEXT;
  I : INTEGER;
  J : INTEGER;
  SINE_TABLE : REAL_BUFFER;
  COSINE_TABLE : REAL_BUFFER;


PROCEDURE TRIG_TABLE;
VAR
   I: INTEGER;
   ANGLE,OMEGA: REAL;
BEGIN
   ANGLE := 0;
   OMEGA := 2 * PI / NPTS;
   FOR I := 1 TO NPTS DO
   BEGIN
      SINE_TABLE[I] := SIN(ANGLE);
      COSINE_TABLE[I] := COS(ANGLE);
      ANGLE := ANGLE + OMEGA;
   END;
END;


{-----------------------------------------------------------------------------
   Fast Hartley transform routine ...

-----------------------------------------------------------------------------}

PROCEDURE FHT;

VAR
   I : INTEGER;
   J : INTEGER;
   K : INTEGER;
   TRG_IND : INTEGER;
   TRIG_IND : INTEGER;
   TRG_INC : INTEGER;
   POWER : INTEGER;
   T_A : INTEGER;
   F_A : INTEGER;
   I_TEMP : INTEGER;
   SECTION : INTEGER;
   S_START : INTEGER;
   S_END : INTEGER;
   MODIFY : INTEGER;
   INDEX : INTEGER;

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

FUNCTION PERMUTE(INDEX: INTEGER): INTEGER;
VAR
  I : INTEGER;
  J : INTEGER;
  S : INTEGER;
BEGIN
   J := 0;
   INDEX := INDEX - 1;
   FOR I := 1 TO ORDER DO
   BEGIN
      S := INDEX DIV 2;
      J := J + J + INDEX - S - S;
      INDEX := S;
   END;
   PERMUTE := J + 1;
END;


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

BEGIN
   FOR I := 1 TO NPTS DO
      ACCU[1,PERMUTE(I)] := IBUFFER[I];

   POWER := 1;
   F_A := 1;
   T_A := 2;
{------------------------------------------------------------------------------
   Start of the Hartley butterfly transform ...
------------------------------------------------------------------------------}

   FOR I := 1 TO ORDER DO
   BEGIN
      J := 1;
      SECTION := 1;
      TRG_INC := NPTS DIV (POWER + POWER);
      REPEAT
         TRG_IND := 1;
         S_START := SECTION * POWER + 1;
         S_END   := (SECTION + 1) * POWER;
         FOR K := 1 TO POWER DO
         BEGIN
            INDEX := J+POWER;
            IF (S_START = INDEX) OR (POWER < 3) THEN
               MODIFY := INDEX
            ELSE
               MODIFY := S_START + S_END - INDEX + 1;
            ACCU[T_A,J] := ACCU[F_A,J] +
                           ACCU[F_A,INDEX] * COSINE_TABLE[TRG_IND] +
                           ACCU[F_A,MODIFY] * SINE_TABLE[TRG_IND];
            TRIG_IND := TRG_IND + NPTS DIV 2;
            ACCU[T_A,INDEX] := ACCU[F_A,J] +
                               ACCU[F_A,INDEX] * COSINE_TABLE[TRIG_IND] +
                               ACCU[F_A,MODIFY] * SINE_TABLE[TRIG_IND];
            TRG_IND := TRG_IND + TRG_INC;
            J := J + 1;
         END;
         J := J + POWER;
         SECTION := SECTION + 2;
      UNTIL J > NPTS;
      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 ...
------------------------------------------------------------------------------}
   FOR I := 1 TO NPTS DO
     BEGIN
       ACCU[F_A,I] := ABS(ACCU[F_A,I] / NPTS);
       WRITELN(ACCU[F_A,I]);
     END;
END;


BEGIN
  PI := 4*ARCTAN(1);
  TRIG_TABLE;
  {...
  ASSIGN(DAFILE,'F701.DA');
  RESET(DAFILE);
  READ(DAFILE,IBUFFER);
  ...}
      FOR I:=1 TO NPTS DO
          IBUFFER[I] := 0;
      IBUFFER[3] := 1;
      IBUFFER[5] := -1;

      FHT;
  {...
  CLOSE(DAFILE);
  ...}

END.


--*-*-*- Next Section -*-*-*--