The Discrete Hartley Transform for the HP-41

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.

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 )

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.

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 )

-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.

-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 )