The Museum of HP Calculators


The Discrete Hartley Transform for the HP-41

Updated 12/12/05. New functions and Improvements noted in red.

This program is Copyright © 2004 by Jean-Marc Baillard and is used here by permission.

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

Overview

-The Discrete Hartley Transform is closely related to the Discrete Fourier Transform,
  but unlike the DFT, the DHT has the advantage of producing real numbers.
-Furthermore, it is (quasi-)symmetrical.
 

   1°) One-dimensional data ( vector )
   2°) Two-dimensional data ( matrix )
   3°) Three-dimensional data ( hypermatrix )
   4°) The Fast Hartley Transform ( one-dimensional only ; X-Functions module required )
   5°) The Hartley Transform and the Filon's Integration formula  ( new )

-This last paragraph deals with the continuous Hartley Transform.
 

1°) One-dimensional Data
 

-The DHT of a vector  [ a1 , a2 , ........ , an ]  is a vector   [ b1 , b2 , ........ , bn ]

  where   bi =  SUMq=1,2,...,n   aq. cas  360°(q-1)(i-1)/n    with   cas µ = sin µ + cos µ
 

Data Registers:             R00 = n      R01 = a1 ;    R02 = a2 ; .......... ;    Rnn = an     ( these registers are to be initialized before executing "DHT" )
Flags: /
Subroutines: /
 

01  LBL "DHT"
02  DEG
03  STO N
04  360
05  ST* Y
06  -
07  RCL 00
08  STO M
09  /
10  0
11  LBL 01
12  RCL M
13  ENTER^
14  SIGN
15  -
16  R^
17  *
18  RCL IND M
19  P-R
20  +
21  +
22  DSE M
23  GTO 01
24  RCL N
25  SIGN
26  X<>Y
27  CLA
28  END

( 46 bytes / SIZE n+1 )
 
 
      STACK        INPUTS      OUTPUTS
           X             i             bi
           L             /             i

 ( 1 <= i <= n )
 

Example:    Find the DHT of  [ 2  4  7  6 ]

  n = 4   therefore  4  STO 00   and   2  STO 01   4 STO 02   7  STO 03   6  STO 04

  1  XEQ "DHT"  >>>   19
  2      R/S           >>>    -7
  3      R/S           >>>    -1
  4      R/S           >>>    -3     whence   DHT [ 2  4  7  6 ]  =  [ 19  -7  -1  -3 ]

Notes:   -If you call  a0 ; a1 ; ...... ; an-1  the components of these vectors  ( all subscripts decremented by 1 ) , replace lines 05-06 by the single line *
             -Execution time = 83 seconds per component if  n = 100.
             -The DFT may be obtained from the DHT as follows: keep the first element unchanged ( 19 ) and reverse the order of the others  ( -3 -1 -7 )

 then  perform   E =  (   [ 19  -7  -1  -3  ]  +  [ 19  -3  -1  -7 ]  ) / 2  =  [ 19 -5  -1  -5 ]  = the even part of the DHT
            and       O = (   [ 19  -7  -1  -3  ]  -  [ 19  -3  -1  -7 ]  ) / 2  =  [ 0 -2  0  2 ]        = the odd part of the DHT

-The DFT = E - i.O =  [ 19  -5+2.i  -1  -5-2.i ]

-We have the property   DHToDHT = n.Id
-Thus, applying the DHT twice yields the original vector multiplied by n.
 

2°) Two-dimensional Data

-The DHT may be generalized to  nxm matrices  [ ai,j ]  ( 1 <= i <= n ; 1 <= j <= m ). The result is also a  nxm matrix  [ bi,j ]

  where  bi,j =  SUMq=1,2,...,n  ; r=1,2,....,m   aq,r. cas 360°[(q-1)(i-1)/n+(r-1)(j-1)/m]            ( cas = sin + cos )
 

Data Registers:       R00 = n.mmm  =  n + m/1000     ( these  nm+1 registers are to be initialized before executing "DHT2" )
                                  R01 = a1,1        Rn+1 = a1,2   .......................................     Rnm-n+1 = a1,m
                                  R02 = a2,1        Rn+2 = a2,2   .......................................     Rnm-n+2 = a2,m
                                       ..........................................................................................................
                                  Rnn = an,1        R2n  =  an,2   ........................................     Rnm = an,m
Flag: /
Subroutine:  /
 

01  LBL "DHT2"
02  DEG
03  360
04  ST* Y
05  ST* Z
06  ST- Z
07  -
08  RCL 00
09  INT
10  STO M
11  /
12  STO N
13  CLX
14  RCL 00
15  FRC
16   E3
17  *
18  ST* M
19  /
20  STO O
21  CLX
22  LBL 01
23  RCL N
24  RCL M
25  ENTER^
26  SIGN
27  -
28  ST* Y
29  RCL 00
30  INT
31  /
32  INT
33  RCL O
34  *
35  +
36  RCL IND M
37  P-R
38  +
39  +
40  DSE M
41  GTO 01
42  CLA
43  END

( 69 bytes / SIZE nm+1 )
 
 
      STACK        INPUTS      OUTPUTS
           Y             j             /
           X             i           bi,j

 ( 1 <= i <= n ; 1 <= j <= m )
 

                                   [ [  1  3  4  10  ]
 Example:    Let  A =   [  4  5  7  14  ]       Calculate   b2,3
                                     [  2  9  6  11  ] ]

 We have 3 rows and 4 columns , therefore  3.004  STO 00  and we store the 12 elements in registers  R01 thru R12  in column order:
  1  STO 01    4  STO 02   2  STO 03  ..........   14  STO 11   11 STO 12

-Then    3  ENTER^
             2  XEQ "DHT2"  >>>  b2,3 =  5.464101608  ( in 14 seconds )

                                                [ [    76             -28         -28          8         ]
-Likewise, we find  DHT (A) =   [  -9.2679    5.9282    5.4641   -3.1962  ]        ( rounded to 4 decimals )
                                                  [ -12.7321  -7.9282   -1.4641    7.1962  ] ]

Notes:   -If all subscripts are decremented by 1 ( a0,0 to an-1,m-1) , replace lines 06-07 by the single line *  and delete line 04.
             -Execution time = 42 seconds per component if  n = 5 and m = 7.
-We have the property   DHToDHT = nm.Id
-Applying the DHT twice yields the original matrix multiplied by nm.
 

3°) Three-dimensional Data
 

-In this case, the DHT of a  nxmxp  hypermatrix [ ai,j,k ]  ( 1 <= i <= n ; 1 <= j <= m ; 1 <= k <= p ) is a  nxmxp hypermatrix  [ bi,j,k ]

  where  bi,j,k =  SUMq=1,2,...,n  ; r=1,2,....,m ; s=1,2,....,p   aq,r,s. cas 360°[(q-1)(i-1)/n+(r-1)(j-1)/m+(s-1)(k-1)/p]         ( cas = sin + cos )

Data Registers:       R00 = n.mmmppp  =  n + m/1,000 + p/1,000,000    ( these  nmp+1 registers are to be initialized before executing "DHT3" )
                                  R01 = a1,1,1        Rn+1 = a1,2,1   .......................................     Rnm-n+1 = a1,m,1
                                  R02 = a2,1,1        Rn+2 = a2,2,1   .......................................     Rnm-n+2 = a2,m,1
                                       ..................................................................................................................
                                  Rn = an,1,1          R2n  =  an,2,1   ........................................    Rnm = an,m,1

                                            Rnm+1 = a1,1,2        Rnm+n+1 = a1,2,2   .......................................     R2nm-n+1 = a1,m,2
                                            Rnm+2 = a2,1,2        Rnm+n+2 = a2,2,2   .......................................     R2nm-n+2 = a2,m,2
                                              ..................................................................................................................................
                                            Rnm+n = an,1,2        Rnm+2n  =  an,2,2   ........................................    R2nm = an,m,2

                                                ..........................................................................................................
                                                    ..........................................................................................................
                                                        ..........................................................................................................

                                                      Rnm(p-1)+1 = a1,1,p        Rnm(p-1)+n+1 = a1,2,p   .......................................     Rnmp-n+1 = a1,m,p
                                                      Rnm(p-1)+2 = a2,1,p        Rnm(p-1)+n+2 = a2,2,p   .......................................     Rnmp-n+2 = a2,m,p
                                                       ....................................................................................................................................................
                                                      Rnm(p-1)+n = an,1,p        Rnm(p-1)+2n  =  an,2,p   ........................................    Rnmp = an,m,p
Flag: /
Subroutine:  /
 

01  LBL "DHT3"
02  DEG
03  360
04  ST* Y
05  ST* Z
06  ST* T
07  ST- T
08  ST- Z
09  -
10  RCL 00
11  INT
12  STO M
13  STO a
14  /
15  STO N
16  CLX
17  RCL 00
18  FRC
19   E3
20  *
21  INT
22  ST* M
23  ST* a
24  /
25  STO O
26  CLX
27  RCL 00
28   E6
29  ST* Y
30  SQRT
31  MOD
32  ST* M
33  /
34  STO P       ( not   STOP  but synthetic  STO P )
35  CLX
36  LBL 01
37  RCL M
38  ENTER^
39  SIGN
40  -
41  STO Z
42  RCL a
43  /
44  INT
45  RCL P
46  *
47  RCL N
48  R^
49  *
50  ST+ Y
51  X<> L
52  RCL 00
53  INT
54  /
55  INT
56  RCL O
57  *
58  +
59  RCL IND M
60  P-R
61  +
62  +
63  DSE M
64  GTO 01
65  CLA
66  END

( 104 bytes / SIZE nmp+1 )
 
 
 
      STACK        INPUTS      OUTPUTS
           Z             k             /
           Y             j             /
           X             i           bi,j,k

 ( 1 <= i <= n ; 1 <= j <= m ; 1 <= k <= p )
 

Example:    Let A =  [ [ [   1  25  40   5    2  ]                        To store these 60 coefficients into the proper registers, you can key in    60   XEQ 10
                                      [   4  36  18  32  37 ]                         after loading this short routine:
                                      [   9   8   39  20  33 ]
                                      [ 16  23  21  10  31 ] ]                       LBL 10   ENTER^   X^2   41   MOD   STO IND Y   X<>Y   DSE X   GTO 10
                                           [ [  31  10  21  23  16 ]
                                             [  33  20  39   8    9  ]
                                             [  37  32  18  36   4  ]
                                             [   2    5   40  25   1  ] ]
                                                    [ [   0  16  23  21  10 ]
                                                      [   1  25  40   5    2  ]
                                                      [   4  36  18  32  37 ]
                                                      [   9   8   39  20  33 ] ] ]                         -Evaluate   b2,4,3

-First,  we have a   4x5x3  hypermatrix  therefore:   4.005003   STO 00
-Then,  3  ENTER^
            4  ENTER^
            2  XEQ "DHT3"  >>>    b2,4,3 =  27.04299874

Notes:   -If all subscripts are decremented by 1 (  a0,0,0 to an-1,m-1,p-1 ) , replace lines 07-08-09 by the single line *  and delete line 04.
             -Execution time = 79 seconds per component if  n = 4 ; m = 5 ; p = 3.
-We have the property   DHToDHT = nmp.Id
-Applying the DHT twice yields the original hypermatrix multiplied by nmp.
 

4°) The Fast Hartley Transform ( n must be an integer power of 2 ; n > 1 )

 -If  n = 2N is an integer power of 2 , the DHT of  [ a1 , a2 , ........ , an ]  =  [ b1 , b2 , ........ , bn ]  may be perform more quickly by the following method:

 a)  First, the  ai  are placed into the jth position where j is calculated as follows ( lines 01 to 37 ):         i-1 is written in binary ,
      and the digits are reversed which yields  j-1   ( for instance, if  i = 4 , 4-1 = 3 = (011)2 >>> (110)2 = 6 = 7-1  thus, a4 is placed into the 7th position )

 b)  Then, for L = 1 to N,   registers  Rii  are replaced by  Rj1 + Rj2 .cos 360°(i-1)/2L + Rj3 .sin 360°(i-1)/2L    ( i = 1 ,2 , ...... , n )

      with  j1 = i - (i-1) mod 2L + (i-1) mod 2L-1  ;
              j2 = j1 + 2L-1  ;
              j3 =  j2           if   ( i-1) mod 2L-1 = 0  and   j3 =  2L + i - (i-1) mod 2L - (i-1) mod 2L-1  otherwise.
 

Data Registers:         R00 = n      R01 = a1 ;    R02 = a2 ; .......... ;    Rnn = an     ( these n+1 registers are to be initialized before executing "FHT" )
           and when the program stops,   R01 = b1 ;     R02 = b2 ; .......... ;     Rnn = bn
                                                       ( Rn+1 thru R2n are used for temporary data storage )
Flag: /
Subroutine:  /
 

  01  LBL "FHT"
  02  DEG
  03  RCL 00
  04  STO M
  05  1
  06  +
  07  LOG
  08  2
  09  LOG
  10  /
  11  INT
  12  STO N
  13  LBL 01
  14  RCL N
  15  STO O
  16  0
  17  RCL M
  18  DSE X
  19  LBL 02
  20  RCL X
  21  2
  22  ST* T
  23  MOD
  24  ST+ Z
  25  X<> L
  26  /
  27  INT
  28  DSE O
  29  GTO 02
  30  SIGN
  31  +
  32  RCL 00
  33  +
  34  RCL IND M
  35  STO IND Y
  36  DSE M
  37  GTO 01
  38  360
  39  STO O
  40  SIGN
  41  STO M
  42  GTO 04
  43  LBL 03
  44  RCL 00
  45  .1
  46  STO Z
  47  %
  48  ISG X
  49  +
  50  %
  51  1
  52  +
  53  REGMOVE
  54  LBL 04
  55  2
  56  ST/ O
  57  RCL 00
  58  STO P              ( synthetic  STO P )
  59  LBL 05
  60  RCL 00
  61  RCL P
  62  RCL X
  63  RCL O
  64  STO Q
  65  SIGN
  66  -
  67  ST* Q
  68  RCL M
  69  ST+ X
  70  ST+ T
  71  MOD
  72  -
  73  ST+ Y
  74  RCL 00
  75  +
  76  RCL P
  77  ENTER^
  78  SIGN
  79  -
  80  RCL M
  81  MOD
  82  ST- Z
  83  +
  84  RCL X
  85  RCL M
  86  ST+ Z
  87  X<> L
  88  X=0?
  89  ENTER^
  90  X=0?
  91  +
  92  CLX
  93  RCL IND T
  94  RCL Q
  95  SIN
  96  ST* Y
  97  X<> L
  98  COS
  99  RCL IND T
100  *
101  +
102  RCL IND Y
103  +
104  STO IND P
105  DSE P
106  GTO 05
107  RCL M
108  ST+ M
109  DSE N
110  GTO 03
111  CLA
112  END

( 176 bytes / SIZE 2n+1 )
 

Example:    Find   DHT [ 2  4  7  6 ]

n = 4   therefore  4  STO 00   and   2  STO 01   4 STO 02   7  STO 03   6  STO 04   XEQ "FHT"

-21 seconds later, we have  R01 = 19   R02 = -7   R03 = -1   R04 = -3  whence   DHT [ 2  4  7  6 ]  =  [ 19  -7  -1  -3 ]

Notes:  -This program doesn't work if  n = 1    ( but DHT [ a ] = [ a ] )
            - I've measured the following execution times:
 
       n      2      4      8     16     32     64   128
  nxDHT     4s    14s    53s  3m32s 14m10s 56m40s  3h47m
    FHT     6s    21s    60s  2m38s  6m36s 15m59s 37m34s

-FHT execution time is approximately proportional to  n.Log n  ( instead of n2 for n.DHT )
-Obviously, this "fast" algorithm remains relatively slow with an HP-41 and it's worthwhile only if  n > 8
-However, the advantage is substantial if  n = 128  ( about 6 times faster )
 and good emulators can reduce execution times to even more "reasonable" values...
-Furthermore, calculations are less complicated and roundoff errors are smaller.
 

5°) The Hartley Transform & the Filon's Integration Formula
 

-Actually, Hartley has defined the real transform  H(x) =  (2.pi) -1/2  §-infinity+infinity  f(t) ( cos x.t + sin x.t ) dt      ( § = integral )
-This transform is strictly reciprocal.

-In the following program, we assume  f(t)  is defined over an interval  [ a ; b ]    ( f = 0 elsewhere )  where n = b - a  is an even integer
 and we omit the factor (2.pi) -1/2      ( Add  PI  ST+ X  SQRT  /   after line 112 if you want to take it into account )

-The integral is estimated by the Filon's formula.
 

Data Registers:             R00 = f(a)      R01 = f(a+1)     R02 = f(a+2)  , .......... ,    Rnn = f(b) = f(a+n)

                                      ( these n+1 registers are to be initialized before executing "HRTL" )
Flags: /
Subroutines: /

-This program uses synthetic register a
-So, it must not be called as more than a first level subroutine.
 

  01  LBL "HRTL"
  02  RAD
  03  RDN
  04  X<>Y
  05  STO a
  06  -
  07  STO M
  08  1.5
  09  1/X
  10  STO O
  11  R^
  12  STO N
  13  ENTER^
  14  ENTER^
  15  X=0?
  16  GTO 01
  17  COS
  18  X^2
  19  1
  20  +
  21  X<>Y
  22  ST+ X
  23  SIN
  24  R^
  25  /
  26  STO O
  27  -
  28  ST+ X
  29  X<>Y
  30  X^2
  31  /
  32  X<> O
  33  2
  34  /
  35  X<>Y
  36  SIN
  37  LASTX
  38  /
  39  X^2
  40  ST+ X
  41  -
  42  1
  43  +
  44  X<>Y
  45  /
  46  LBL 01
  47  RCL M
  48  RCL a
  49  +
  50  R^
  51  *
  52  RCL IND M
  53  P-R
  54  STO P                    ( synthetic )
  55  X<>Y
  56  ST+ P
  57  X<>Y
  58  -
  59  RCL a
  60  R^
  61  *
  62  RCL 00
  63  P-R
  64  ST- P
  65  X<>Y
  66  ST- P
  67  -
  68  +
  69  *
  70  RCL O
  71  RCL P
  72  *
  73  2
  74  /
  75  -
  76  RCL N
  77  X=0?
  78  GTO 00
  79  SIN
  80  LASTX
  81  ST/ Y
  82  COS
  83  -
  84  4
  85  *
  86  RCL N
  87  X^2
  88  /
  89  GTO 02
  90  LBL 00
  91  CLX
  92  RCL O
  93  ST+ X
  94  LBL 02
  95  STO P                    ( synthetic )
  96  X<>Y
  97  LBL 03
  98  RCL M
  99  RCL a
100  +
101  RCL N
102  *
103  RCL IND M
104  P-R
105  +
106  RCL P
107  X<> O
108  STO P                    ( synthetic )
109  *
110  +
111  DSE M
112  GTO 03
113  RCL N
114  SIGN
115  RDN
116  CLA
117  DEG
118  END

( 167 bytes / SIZE nnn+1 )
 
 
 
      STACK        INPUTS      OUTPUTS
           Z             a             /
           Y             b             /
           X             x          H(x)
           L             /            x

  where  b-a =  n  must be  even  and  x  is expressed in  radians

Example1:      f  is only known by the following samples, and  f = 0  if  t < 3  or  t > 7
 
 
     t     3     4     5     6     7
   f(t)     1     2     1    -2    -7

-We store the f-values into registers  R00 thru R04:   1  STO 00   2  STO 01   1  STO 02   -2  STO 03   -7  STO 04
 and then, for instance:

  3  ENTER^
  7  ENTER^
  2  XEQ "HRTL"   >>>>   H(2) = -3.8768  and similarly:
 
      x      -3      -2      -1       0       1       2       3
   H(x)   0.7675  -3.5134  -2.8271  -1.3333  -9.6763  -3.8768  -3.7485

-Actually, we had   f(t) = -14 +8.t - t2   and these results are exact because Filon's formulae produce perfect accuracy if  f is a polynomial of degree 2

-If you compute more H(x)-values , you'll see the following graph ( approximately ):

                                                 H(x)
                                 *                |                           *
                               *   *             |        *               *   *
                    *       *      *        *  |     *   *           *       *        *
      --------*--*--*-----*---*--*|--*----*----- *-------*--*---*------------------------- x
                          *          *  *       |*          *      *               *
                                        *         |              *  *
                                                   |                *
 

Example2:     We take  f(t) =  e -t/2  for  t >= 0  and  f(t) = 0  if  t < 0

- H(x) may be expressed analytically and  H(x) = 2(1+2x)/(1+4x2)

-If we use f(0) , f(1) , .............. , f(16) to represent the function  ( these 17 numbers are to be stored in registers R00 thru R16 ),
 we obtain the following results:
 
 
      x       -1     -0.5    -0.25       0   0.2071       1       2
  HRTL  -0.3971   0.0034   0.8015   2.0000   2.4154   1.2025   0.5930
   exact    -0.4       0      0.8       2   2.4142      1.2   0.5882

-The accuracy is quite satisfactory.
-The graph of the Hartley transform looks like this:

                                                              H(x)
                                                                 |
                                                                 |  *                               H(x) is maximum for x = 0.2071
                                                                 |*   *
                                                                 |          *
                                                                 |                       *
                                                                 |                                                          *
  ---------------------------------------* -|----------------------------------------------- x
            *                                                   |
                                     *                 *       |
                                                                 |
 

Reference:    Ronald N. Bracewell   "The Fourier Transform and its Applications"   McGraw-Hill        ISBN 0-07-116043-4
 

Go back to the HP-41 software library
Go back to the general software library
Go back to the main exhibit hall