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"
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 2.4154 1.2025 0.593 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