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:

  Patrice Denis - Thèse - "Quaternions et Algèbres Géométriques, de nouveaux outils pour les images numériques couleur" ( in French )