Quaternionic Discrete Fourier Transform for the HP-41
Overview
-"QDFT" calculates the discrete Fourier Transform of a vector or a matrix of quaternions.
-Due to the non-commutativity of the multiplication, many Fourier-Transforms
may be defined.
-Here we use the right-sided formulae:
1-Dimension: F (p) = ( 1/sqrt(M) ) SUMm=0,1,....,M-1 f(m) exp ( -2PI µ m p / M )
2-Dimensions: F (p,q) = ( 1/sqrt(M.N) ) SUMm=0,1,....,M-1 SUMn=0,1,....,N-1 f(m,n) exp [ -2PI µ { ( m p / M ) + ( n q / N ) } ]
where µ = µ' i + µ" j + µ''' k is unitary and purely imaginary: | µ | = 1 & Re(µ) = 0 ( this program does not check )
-The inverse transform is simply obtained by replacing -2PI
by +2PI in the above formulae.
Program Listing
Data Registers: • R00 = M.NNN ( Registers R00 & R18 thru R20+4M.N are to be initialized before executing "QDFT" )
R01 thru R17: temp R09 to R12: partial sums
• R18 = µ'
• R21 to R24 = f(0,0)
....................
• R21+4M(N-1) to R24+4M(N-1) = f(0,N-1)
• R19 = µ"
..............................
....................
...............................................................
• R20 = µ''' •
R17+4M to R20+4M = f(M-1,0) ....................
• R17+4M.N to R20+4M.N = f(M-1,N-1)
Flags: F00 & F01
CF 00 All the coefficients of the DFT are computed
and stored into R21+4M.N ................... R20+8M.N
SF 00 Only 1 coefficient F(p,q)
is computed and returned in the stack.
CF 01 = Fourier Transform
SF 01 = Inverse Fourier Transform
-Line 141 is a three-byte GTO 12
| 01 LBL "QDFT"
02 STO 13 03 X<>Y 04 STO 14 05 RCL 00 06 INT 07 ST- L 08 LASTX 09 E3 10 * 11 X=0? 12 SIGN 13 * 14 STO 16 15 FS? 00 16 GTO 00 17 STO 17 18 LBL 12 19 RCL 17 20 1 21 - 22 STO Y 23 RCL 00 24 INT 25 MOD 26 STO 13 27 X<>Y 28 RCL 00 29 INT 30 / 31 INT 32 STO 14 |
33 LBL 00
34 RCL 00 35 INT 36 ST/ 13 37 ST- L 38 LASTX 39 E3 40 * 41 X=0? 42 SIGN 43 ST/ 14 44 * 45 4 46 * 47 20.02 48 + 49 STO 15 50 CLX 51 STO 09 52 STO 10 53 STO 11 54 STO 12 55 LBL 01 56 RCL 15 57 INT 58 4 59 / 60 6 61 - 62 STO Y 63 RCL 00 64 INT |
65 MOD
66 RCL 13 67 * 68 X<>Y 69 RCL 00 70 INT 71 / 72 INT 73 RCL 14 74 * 75 + 76 ST+ X 77 PI 78 * 79 RCL 20 80 RCL 19 81 RCL 18 82 R^ 83 FC? 01 84 CHS 85 ST* T 86 ST* Z 87 * 88 0 89 XEQ "E^Q" 90 STO 05 91 RDN 92 STO 06 93 RDN 94 STO 07 95 X<>Y 96 STO 08 |
97 RCL IND 15
98 STO 04 99 DSE 15 100 RCL IND 15 101 STO 03 102 DSE 15 103 RCL IND 15 104 STO 02 105 DSE 15 106 RCL IND 15 107 STO 01 108 XEQ "Q*Q" 109 ST+ 09 110 RDN 111 ST+ 10 112 RDN 113 ST+ 11 114 X<>Y 115 ST+ 12 116 DSE 15 117 GTO 01 118 RCL 16 119 SQRT 120 ST/ 09 121 ST/ 10 122 ST/ 11 123 ST/ 12 124 FS? 00 125 GTO 00 126 RCL 16 127 RCL 17 128 + |
129 4
130 * 131 17 132 + 133 .004 134 + 135 E3 136 / 137 9 138 + 139 REGMOVE 140 DSE 17 141 GTO 12 142 RCL 16 143 4 144 * 145 ENTER^ 146 ST+ X 147 E3 148 / 149 + 150 21.02 151 + 152 GTO 09 153 LBL 00 154 RCL 12 155 RCL 11 156 RCL 10 157 RCL 09 158 LBL 09 159 END |
( 239 bytes / SIZE 21+4 M.N or 21+8
M.N )
| SF 00 STACK | INPUTS | OUTPUTS |
| T | / | t |
| Z | / | z |
| Y | q | y |
| X | p | x |
with F(p,q) = x + y i + z j + t k 0 <= p < M , 0 <= q < N
Example: µ = ( i + 2 j + 4 k ) / sqrt(21) and Matrix A =
( 1 + 2 i + 3 j + 4 k ) ( 2 + 3 i + j + 2 k
) ( 1 + 3 i + 4 j + 2 k )
( 3 + 4 i + 2 j + k )
( 2 + i + 4 j + 3 k ) ( 1 + 2
i + j + 3 k ) ( 2 + i + 3 j +
2 k ) ( 2 + 3 i + 4 j + 5 k )
( 3 + 4 i + 2 j + k ) ( 3 + 4
i + 5 j + 6 k ) ( 4 + 2 i + 3 j + 4 k )
( 1 + 2 i + 3 j + 4 k )
-M = 3 , N = 4 so 3.004 STO 00
-Key in 1 STO 18 2 STO 19 4 STO 20 21 SQRT ST/ 18 ST/ 19 ST/ 20 ( the program does not check that µ is unitary )
-Store the 12 quaternions ( 48 real numbers ) into R21 thru R68 as shown below:
R21 R22 R23 R24 R33
R34 R35 R36 R45 R46 R47
R48 R57 R58 R59 R60
R25 R26 R27 R28 R37
R38 R39 R40 R49 R50 R51
R52 R61 R62 R63 R64
R29 R30 R31 R32 R41
R42 R43 R44 R53 R54 R55
R56 R65 R66 R67 R68
-To compute, for instance, F(2;3)
SF 00 CF 01
3 ENTER^
2 XEQ "QDFT" >>>> -0.182134053
RDN -0.119263223
RDN 1.977063795
RDN 1.221466790
Whence F(2;3) = -0.182134053 - 0.119263223
i + 1.977063795 j + 1.221466790 k
Notes:
-If N = 1 ( vector ) you may also store M in R00
instead of M.001 and q is arbitrary
provided it is not an alpha string.
-Register R17 is unused if flag F00 is set.
| CF 00 STACK | INPUTS | OUTPUTS |
| X | / | bbb.eee |
where bbb.eee is the control number of the DFT
Example: With the same µ and A, which are still in the same registers ( SIZE 117 or greater )
CF 00 CF 01 XEQ "QDFT" >>>> 69.116 the DFT is in registers R69 to R116 at the end
-The first column ( in registers R69 to R80 ) is
7.217 + 8.949 i + 10.104 j + 10.681
k
-1.396 + 0.940 i - 1.267 j - 0.080
k
( rounded to 3D )
0.241 + 0.503 i - 0.176 j - 2.807
k
-The last element in registers R113 to R116 - which has been computed
earlier -
is F(2;3) = -0.182134053 - 0.119263223
i + 1.977063795 j + 1.221466790 k
Notes:
-If N = 1 ( vector ), you may also store M in R00 instead of M.001
-May be you will want to calculate the Inverse DFT from these results:
-Key in 69.021048 REGMOVE SF 01 R/S
>>>> 69.116
-Compare the numbers in R69 to R116 to the initial matrix to get an
idea of the roundoff-errors...
-A good emulator in turbo mode is not superfluous...
Reference:
[1] Patrice Denis - Thèse - "Quaternions
et Algèbres Géométriques, de nouveaux outils pour
les images numériques couleur" ( in French )