hp41programs

Quat-DFT

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 )