The Museum of HP Calculators
The Museum of HP Calculators
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.
-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