# hp41programs

Spectral Analysis

# Spectral Analysis for the HP-41

Overview

1°) One Dimensional Problems - Real Arguments

1-a)  The Discrete Fourier Transform  ( program#1 )
1-b)  The Discrete Fourier Transform  ( program#2 )
1-c)  Filon's Formula
1-d)  More Accurate Filon-Type Formulae

1-d-1)  n = multiple of 4
1-d-2)  n = multiple of 6

2°) One Dimensional Problems - Complex Arguments

2-a)  The Discrete Fourier Transform
2-b)  The Fast Fourier Transform

3°) Two Dimensional Problems - Real Arguments

3-a)  The Discrete Fourier Transform  ( program#1 )
3-b)  The Discrete Fourier Transform  ( program#2 )

4°)  Two Dimensional Problems - Complex Arguments

5°)  Three Dimensional Problems - Real Arguments

6°)  Three Dimensional Problems - Complex Arguments

1°) One-dimensional Problems - Real Arguments

a) The Discrete Fourier Transform ( Program#1 )

-Given an array of n real numbers  [ a1 , .......... , an ] ,  "DFT" computes the complex components zk = xk + i yk
of the Discrete Fourier Transform   [ z1 , .......... , zn ]

-We have  zk = Sum j=1,2,...,n   aj exp( -2i(pi) (j-1).(k-1)/n )             i = sqrt(-1)

Data Registers:         •  R00 = n                                           ( Registers R00 thru Rnn are to be initialized before executing "DFT" )

•   R01 = a1   •   R02 = a2  ................  •   Rnn = an
Flags: /
Subroutines:  /

 01  LBL "DFT"  02  DEG  03  STO O  04  360  05  ST* Y  06  -  07  RCL 00  08  STO M  09  /  10  STO N  11  CLST  12  LBL 01  13  RCL M  14  RCL N  15  ST* Y  16  -  17  RCL IND M  18  P-R  19  ST+ T  20  RDN  21  -  22  DSE M  23  GTO 01  24  RCL O  25  SIGN  26  X<> Z  27  CLA  28  END

( 51 bytes / SIZE nnn+1 )

 STACK INPUTS OUTPUTS Y / yk X k xk L / k

Example:    Calculate  DFT  [ 6  3  2  1 ]

n = 4  STO 00            6  STO 01    3  STO 02    2  STO 03    1  STO 04

1  XEQ "DFT"  >>>>   x1 = 12   X<>Y   y1 =  0
2       R/S          >>>>   x2 =  4    X<>Y   y2 = -2
3       R/S          >>>>   x3 =  4    X<>Y   y3 =  0
4       R/S          >>>>   x4 =  4    X<>Y   y4 =  2

-So,   DFT [ 6  3  2  1 ]  =  [ 12  4-2.i  4  4+2.i ]

Notes:

-For n = 100 , execution time = 85 seconds / component
-The discrete inverse Fourier transform may be computed by an almost identical program:

 01  LBL "IDFT"  02  DEG 03  STO O 04  360 05  ST* Y 06  - 07  RCL 00 08  STO M        09  / 10  STO N 11  CLST 12  LBL 01 13  RCL M 14  RCL N 15  ST* Y 16  - 17  RCL IND M 18  P-R 19  ST+ T 20  RDN 21  + 22  DSE M 23  GTO 01        24  RCL 00 25  ST/ Z 26  / 27  RCL O 28  SIGN 29  X<> Z        30  CLA 31  END

( 56 bytes / SIZE nnn+1 )

 STACK INPUTS OUTPUTS Y / yk X k xk L / k

Example:    Calculate  IDFT  [ 6  3  2  1 ]

n = 4  STO 00            6  STO 01    3  STO 02    2  STO 03    1  STO 04

1  XEQ "IDFT"  >>>>   x1 =  3    X<>Y   y1 =   0
2       R/S           >>>>   x2 =  1    X<>Y   y2 =  0.5
3       R/S           >>>>   x3 =  1    X<>Y   y3 =   0
4       R/S           >>>>   x4 =  1    X<>Y   y4 = -0.5

-So,   IDFT [ 6  3  2  1 ]  =  [ 3  1+i/2  1  1-i/2 ]

Note:

-I've used the formulae that produce the same results as the FFT & IFFT of the HP-48
-Some authors use another sign for the exponential.

b) The Discrete Fourier Transform ( Program#2 )

-A periodic function y may be written     y(x) = a0 + a1.cos ( 360° x / n ) + b1.sin ( 360° x / n ) + .... +  ak.cos ( 360° k.x / n ) + bk.sin ( 360° k.x / n ) + ....
( n = period:  y(x+n) = y(x)  for all x )

-"SPA" computes the coefficients  a0 , a1 , b1 , ..... , ak , bk  ( 0 <= k <= n/2 )  using the Discrete Fourier Transform.

Data Registers:          •  R00 = n                                                 ( Registers R00 thru Rnn are to be initialized before executing "SPA" )

•   R01 = y(1)   •   R02 = y(2)  ................  •   Rnn = y(n) = y(0)
Flags: /
Subroutines:  /

 01  LBL "SPA"  02  DEG 03  STO O 04  360 05  * 06  RCL 00 07  STO M 08  / 09  STO N 10  CLST 11  LBL 01 12  RCL N 13  RCL M 14  * 15  RCL IND L 16  P-R 17  ST+ T 18  RDN 19  + 20  DSE M       21  GTO 01 22  RCL O 23  ST+ X 24  RCL 00 25  MOD 26  X#0? 27  SIGN 28  1 29  + 30  RCL 00       31  / 32  ST* Z 33  * 34  RCL O       35  SIGN 36  X<> Z 37  CLA 38  END

( 62 bytes / SIZE nnn+1 )

 STACK INPUTS OUTPUTS Y / bk X k ak L / k

( 0 <= k <= n/2 )

Example:      y  is a periodic function ( period = 4 ) and we have the samples:

 x 1 2 3 4 y(x) 6 3 2 1

-We store these data into registers R00 thru R04       4  STO 00    6 STO 01    3  STO 02   2  STO 03   1  STO 04
-Then,

0 XEQ "SPA"  >>>      3
1      R/S          >>>     -1   X<>Y  2
2      R/S          >>>     -1   X<>Y  0

-Thus    y(x)  ~  3 - cos(90°x) + 2 sin(90°x) - cos(180°x)

c) Filon's formula ( assuming n is even )

-The Discrete Fourier Transform provides a trigonometric sum which does collocate with the function y at the given data points.
-However, the true Fourier coefficients are equal to integrals of the form:     § y(x) cos(kx) dx   or   § y(x) sin(kx) dx   and the DFT is only an approximation:
it's almost equivalent to compute these integrals by the trapezoidal rule.
-One may think, for instance, to use Simpson's rule to obtain a better accurracy: it's often true for small k-values, but in the neighborhood of n/2 it doesn't work!
-Filon has developed special formulae for such integrals. The following program uses simplified versions of this method
and we assume that f is continuous over [ 0 n ] ( not only over [ 0 n [ )

Data Registers:                •  R00 = n                                                 ( Registers R00 thru Rnn are to be initialized before executing "FIL" )

•   R01 = y(1)   •   R02 = y(2)  ................  •   Rnn = y(n) = y(0)
Flag:  F10
Subroutines:  /

-Synthetic registers M N O P a   are used by this program and may be replaced by any unused data registers.
-Since synthetic register a   is used, "FIL" must not be called as more than a first level subroutine.

 01  LBL "FIL"    02  DEG 03  STO a 04  360 05  * 06  RCL 00 07  STO M 08  / 09  STO N 10  X=0? 11  GTO 00 12  COS 13  CHS 14  STO Y 15  LASTX 16  SIN 17  LASTX        18  D-R 19  / 20  ST+ T 21  ST+ X 22  + 23  * 24  1 25  + 26  RCL N 27  D-R 28  2 29  ST/ Z 30  / 31  X^2 32  ST/ Z 33  / 34  GTO 01        35  LBL 00 36  3 37  1/X 38  ENTER^ 39  ST+ Y 40  LBL 01 41  STO P 42  X<>Y 43  STO O 44  CLST 45  LBL 02 46  RCL O 47  X<> P 48  STO O 49  RCL IND M 50  * 51  SIGN 52  CLX 53  RCL M 54  RCL N 55  ST* Y 56  X<> L 57  P-R 58  ST+ T 59  RDN 60  + 61  DSE M 62  GTO 02        63  RCL 00 64  2 65  / 66  ST/ Z 67  / 68  RCL a 69  SIGN 70  X<> Z 71  CLA 72  END

( 110 bytes / SIZE nnn+1 )

 STACK INPUTS OUTPUTS Y / bk X k ak L / k

Example:      Suppose y is defined as   y(x) = ( 10-x ) ln ( 1+x )   over [ 0,10 ]  with a period 10

n = 10  STO 00   and we store the exact values  y(1) into R01 , y(2) into R02 , ........... , y(10) into R10

Filon's formula ( FIL )            Exact results ( rounded to 4D )        Discrete Fourier Transform ( SPA )
and,
0  XEQ "FIL"  >>>     6.5005                                              6.5073                                                6.4066
1     R/S          >>>    -3.2200   X<>Y   2.1912                 -3.1971     2.1834                               -3.4022       2.1780
2     R/S          >>>    -1.1338   X<>Y   0.5145                 -1.0982     0.5133                               -1.3152       0.5015
3     R/S          >>>    -0.5979   X<>Y   0.1855                 -0.5562     0.2012                               -0.7953       0.1809
4     R/S          >>>    -0.3706   X<>Y   0.0612                 -0.3353     0.0993                               -0.6120       0.0658
5     R/S          >>>    -0.2285   X<>Y   0                          -0.2238     0.0561                               -0.2818       0

-Thus    y(x)  ~  6.5005  - 3.2200 cos(36°x) + 2.1912  sin(36°x) - 1.1338 cos(72°x) +  0.5145 sin(72°x) - ......................

Note:

-"FIL" gives exact values when y is a polynomial of degree < 3.

d)  More Accurate Filon-Type Formulae

d-1)  n = multiple of 4

-In order to compute integrals of the form    § y(x) cos(µ.x) dx   or   § y(x) sin(µ.x) dx  , we seek a formula that produces perfect accuracy
when y(x) is a polynomial p(x) with deg p < 5

-Successive integrations by parts give:

§ab y(x) cos(µ.x) dx = B cos µ.b - A cos µ.a + D sin µ.b - C sin µ.a

§ab y(x) sin(µ.x) dx  = -D cos µ.b + C cos µ.a + B sin µ.b - A sin µ.a

A = y'(a)/µ2 - y'''(a)/µ4 + ......              C = y(a)/µ - y''(a)/µ3 + y'''''(a)/µ5 - ..........
with
B = y'(b)/µ2 - y'''(b)/µ4 + ......              D = y(b)/µ - y''(b)/µ3 + y'''''(b)/µ5 - ..........

-So, we can use numerical differentiation to evaluate  A , B , C , D , namely:

y'(a) = [ -25 y(a) + 48 y(a+h) - 36 y(a+2h) + 16 y(a+3h) - 3 y(a+4h) ] / (12h)

y''(a) = [ 35 y(a) - 104 y(a+h) + 114 y(a+2h) -56 y(a+3h) + 11 y(a+4h) ] / (12h2)               and similar expressions for y'(b)  y''(b)  y'''(b)

y'''(a) = [ -5 y(a) + 18 y(a+h) - 24 y(a+2h) +14 y(a+3h) - 3 y(a+4h) ] / (2h3)

y''''(a) = [ y(a) - 4 y(a+h) + 6 y(a+2h) - 4 y(a+3h) + y(a+4h) ] / h4 = y''''(b)

-These formulas are applied with  h = 1 and  b = a+4  over the intervals [0,4]  [4,8]  , ............. , [n-4,n]    So n must be a multiple of 4
-They are exact if y(x) is a polynomial p(x) with deg p < 5

-Bode's rule is used to compute the first coefficient a0 = the mean of y(x) over the interval [0,n]
-This integration formula also gives perfect accuracy if deg p = 5

Data Registers:            •  R00 = n = mulptiple of 4              ( Registers R00 & R10 thru Rnn+10 are to be initialized before executing "FIL4" )

•   R10 = y(0+)   •   R11 = y(1)  ................  •   Rnn+10 = y(n-)          R01 to R09: temp   R08 = 360 k/n
Flags:  /
Subroutines:  /

-Lines 40-117 are three-byte  GTOs

 01  LBL "FIL4"    02  DEG   03  360   04  *   05  RCL 00   06  STO 09   07  /   08  STO 08   09  D-R   10  X^2   11  STO 07   12  10.009   13  ST+ 09   14  CLA   15  GTO 01   16  LBL 00   17  RCL 01   18  RCL 05   19  +   20  7   21  *   22  RCL 02   23  RCL 04   24  +   25  32   26  *   27  +   28  RCL 03   29  12   30  *   31  +   32  45   33  /   34  ST+ M   35  LBL 01 36  RCL IND 09   37  STO 05   38  DSE 09   39  FS? 30   40  GTO 04    41  RCL IND 09   42  STO 04   43  DSE 09   44  RCL IND 09   45  STO 03   46  DSE 09   47  RCL IND 09   48  STO 02   49  DSE 09   50  RCL IND 09   51  STO 01   52  RCL 08   53  X=0?   54  GTO 00   55  XEQ 02   56  RCL 09   57  INT   58  10   59  -   60  RCL 08   61  *   62  STO 06   63  X<>Y   64  P-R   65  ST- M   66  X<>Y   67  ST- N   68  RCL 01   69  RCL 02   70  RCL 04 71  +   72  4   73  *   74  -   75  RCL 03          76  6   77  *   78  +   79  RCL 05   80  +   81  RCL 07   82  /   83  STO O   84  XEQ 03   85  RCL 06   86  X<>Y   87  P-R   88  ST+ N   89  X<>Y   90  ST- M   91  RCL 01   92  X<> 05   93  STO 01   94  RCL 02   95  X<> 04   96  STO 02   97  RCL O   98  XEQ 03   99  RCL 08 100  4 101  * 102  RCL 06 103  + 104  STO 06 105  X<>Y 106  P-R 107  ST- N 108  X<>Y 109  ST+ M 110  XEQ 02        111  RCL 06 112  X<>Y 113  P-R 114  ST- M 115  X<>Y 116  ST- N 117  GTO 01 118  LBL 02 119  RCL 01 120  5 121  * 122  RCL 02 123  18 124  * 125  - 126  RCL 03 127  24 128  * 129  + 130  RCL 04 131  14 132  * 133  - 134  RCL 05 135  3 136  * 137  + 138  RCL 07 139  ST+ X 140  / 141  RCL 01 142  25 143  * 144  RCL 02        145  48 146  * 147  - 148  RCL 03 149  36 150  * 151  + 152  RCL 04 153  16 154  * 155  - 156  RCL 05 157  3 158  * 159  + 160  12 161  / 162  - 163  RCL 07 164  / 165  RTN 166  LBL 03 167  RCL 01 168  35 169  * 170  RCL 02 171  104 172  * 173  - 174  RCL 03 175  114 176  * 177  + 178  RCL 04         179  56 180  * 181  - 182  RCL 05 183  11 184  * 185  + 186  12 187  / 188  - 189  RCL 07 190  / 191  RCL 01 192  + 193  RCL 08 194  D-R 195  / 196  RTN 197  LBL 04 198  RCL N 199  RCL M 200  RCL 00 201  2 202  / 203  ST/ Z 204  / 205  CLA 206  END

( 284 bytes / SIZE nn+11 )

 STACK INPUTS OUTPUTS Y / bk X k ak

Example:          y(x) is periodic ( period = 12 )

 x 0+ 1 2 3 4 5 6 7 8 9 10 11 12- y 1 2 4 7 10 12 12 8 6 7 10 14 20

n = 12   STO 00   and store these 13 y-values into  R10  R11  ..................  R22

0  XEQ "FIL4"  >>>>   a0 =   8.455556                                                                  --- execution time =  6s  ---
1        R/S         >>>>    a1 = -1.049462   X<>Y   b1 =  -1.073589                         --- execution time = 37s ---
2        R/S         >>>>    a2 =  1.838711   X<>Y   b2 =  -4.050732                                        idem
.....................................................................................................

6        R/S         >>>>    a6 =  0.157933   X<>Y   b6 =  -0.975730             Unlike "ASP" & "FIL",
7        R/S         >>>>    a7 =  0.118962   X<>Y   b7 =  -0.861560             "FIL4"  may produce a better approximation than  bn/2 = 0
......................................................................................................

Notes:

-The integrals are computed over the intervals [0,4] , [4,8] , ...etc...
so, this method works even if the derivatives have discontinuities for x = 0 , 4 , 8 , ...etc...
-However, the formulas involve numerical differentiation which may produce inaccurate results, especially if µ = 2(pi)k/n  is small.
-So, we can re-write these formulae in order to reduce roundoff errors, but we assume hereunder that the function y(x) is continuous for x = 0

-After some trigonometry, it yields:

§0n y(x) cos µ.x dx = Sum i=0,4,8,.....,n-4  y(i) A0 cos µ.i + y(i+1) [ A1 cos (µ.i+µ) + B1 sin (µ.i+µ) ] + y(i+2) A2 cos (µ.i+2µ)
+ y(i+3) [ A1 cos (µ.i+3µ) - B1 sin (µ.i+3µ) ]

§0n y(x) sin µ.x dx = Sum i=0,4,8,.....,n-4  y(i) A0 sin µ.i + y(i+1) [ A1 sin (µ.i+µ) - B1 cos (µ.i+µ) ] + y(i+2) A2 sin (µ.i+2µ)
+ y(i+3) [ A1 sin (µ.i+3µ) + B1 cos (µ.i+3µ) ]
with

A0 = (1/2µ2-3/µ4) cos 4µ + (25/6µ2-5/µ4) + (-11/6µ3+2/µ5) sin 4µ
A1 = (-4/3µ2+7/µ4) cos 3µ + (-4/µ2+9/µ4) cos µ + (14/3µ3-4/µ5) sin 3µ + (26/3µ3-4/µ5) sin µ
A2 = (6/µ2-24/µ4) cos 2µ + (-19/µ3+12/µ5) sin 2µ

and     B1 = (4/3µ2-7/µ4) sin 3µ + (-4/µ2+9/µ4) sin µ + (14/3µ3-4/µ5) cos 3µ + (-26/3µ3+4/µ5) cos µ

-Thus, these coefficients may be computed in the beginning and the program will be faster for large n-values.
-However, the formulas again produce important roundoff-errors when µ is small - by cancellation of leading digits -
but we can use Taylor series in this case, namely:

A0 = 28/45 + 16µ2/63 - 128µ4/405 + 13568µ6/93555 - .............
A1 = 64/45 - 32µ2/63 + 712µ4/2835 - 1124µ6/18711 + .............
A2 =  8/15 + 16µ2/21 - 176µ4/945 + 544µ6/31185 - .............

B1 = 64µ3/189 - 1888µ5/14175 + 2152µ7/93555 - .............

-In the following routine, the Taylor series are used for | µ | < 1/6  ( lines 30-31 ) and the terms of orders above µ5 are neglected
-If µ > 1/6 , the trigonometric expressions are employed.

-So, these coefficients are correct to 5D or 6D at least.
-Change the threshold and add more terms in the Taylor series if you need a better accuracy.

Data Registers:          •  R00 = n = mulptiple of 4              ( Registers R00 & R11 thru Rnn+10 are to be initialized before executing "FIL4+" )

•   R11 = y(1)  ................  •   Rnn+10 = y(n) = y(0)         R01 to R09: temp      R10 is unused

>>>> When the program stops,   R01 = ak ,  R02 = bk

Flags:  /
Subroutines:  /

-Lines 28-92 are three-byte  GTOs

 01  LBL "FIL4+"   02  DEG   03  360   04  *   05  RCL 00   06  /   07  STO 03   08  D-R   09  STO 04   10  X^2   11  STO 05   12  14   13  STO 06   14  32   15  STO 07   16  12   17  STO 08   18  CLX   19  STO 01   20  STO 02   21  STO 09   22  45   23  ST/ 06   24  ST/ 07   25  ST/ 08   26  RCL 04   27  X=0?   28  GTO 01   29  ABS   30  6   31  1/X   32  XY 140  STO 09 141  RCL 08 142  4 143  RCL 05 144  / 145  14 146  3 147  / 148  - 149  RCL 04        150  / 151  P-R 152  ST+ 09 153  X<>Y 154  ST- 07 155  RCL 03 156  9 157  RCL 05 158  / 159  4 160  - 161  P-R 162  ST+ 07 163  X<>Y 164  ST- 09 165  RCL 03 166  4 167  RCL 05 168  / 169  26 170  3 171  / 172  - 173  RCL 04 174  / 175  P-R 176  ST- 09 177  X<>Y 178  ST- 07 179  6 180  24 181  RCL 05 182  / 183  - 184  RCL 03 185  ST+ X 186  COS 187  * 188  12 189  RCL 05 190  / 191  19 192  - 193  RCL 04 194  / 195  RCL 03 196  ST+ X 197  SIN 198  * 199  + 200  RCL 05        201  ST/ 06 202  ST/ 07 203  ST/ 09 204  / 205  STO 08 206  LBL 01 207  10 208  RCL 00 209  STO 04 210  + 211  STO 05 212  LBL 02 213  RCL 03 214  RCL 04 215  * 216  RCL IND 05 217  RCL 06 218  * 219  P-R 220  ST+ 01 221  X<>Y 222  ST+ 02 223  DSE 04 224  DSE 05 225  RCL 03 226  RCL 04 227  * 228  RCL IND 05 229  P-R 230  RCL X 231  RCL 07 232  * 233  ST+ 01 234  CLX 235  RCL 09 236  * 237  ST- 02 238  X<> L 239  * 240  ST+ 01 241  CLX 242  RCL 07 243  * 244  ST+ 02 245  DSE 04 246  DSE 05 247  RCL 03 248  RCL 04 249  * 250  RCL IND 05 251  RCL 08 252  * 253  P-R 254  ST+ 01 255  X<>Y 256  ST+ 02 257  DSE 04 258  DSE 05 259  RCL 03 260  RCL 04 261  * 262  RCL IND 05 263  P-R 264  RCL X 265  RCL 07 266  * 267  ST+ 01 268  CLX 269  RCL 09 270  * 271  ST+ 02 272  X<> L 273  * 274  ST- 01 275  CLX 276  RCL 07 277  * 278  ST+ 02 279  DSE 05 280  DSE 04 281  GTO 02 282  RCL 00 283  2 284  / 285  ST/ 01 286  ST/ 02 287  RCL 02 288  RCL 01 289  END

( 405 bytes / SIZE nnn+11 )

 STACK INPUTS OUTPUTS Y / bk X k ak

Example1:         y(x) is periodic ( period = 12 )    y(0) = y(12) = 1

 x 1 2 3 4 5 6 7 8 9 10 11 12 y 2 4 7 10 12 12 8 5 1 -2 -2 1

n = 12   STO 00   and store these 12 y-values into   R11  ..................  R22

0  XEQ "FIL4+"  >>>>    a0 =   4.770370                                                                  --- execution time = 12s  ---
1        R/S            >>>>    a1 = -5.802627   X<>Y   b1 =  3.282884                          --- execution time = 23s ---
2        R/S            >>>>    a2 =  1.060168   X<>Y   b2 =  0.111899                                        idem
.....................................................................................................

Example2:      y(x)  periodic ( period = 48 )  and  y(x) = (x/10) Ln(49-x)  over  [0,48]

-To store y(n) in the proper registers, load for instance this short routine:

LBL 01     RCL Y      *              /                           DSE X                         and key in    48  XEQ 01
STO Y      -               10            STO IND Y         GTO 01
49             LN           ST+ Z       R^

-Then  n = 48  STO 00

0  XEQ "FIL4+"  >>>>    a0 =   6.083282                                                                  --- execution time = 45s  ---
1        R/S            >>>>    a1 = -2.250492   X<>Y   b1 =  -2.782105                         --- execution time = 64s  ---  ( 147s with "FIL4" )
2        R/S            >>>>    a2 = -0.940576   X<>Y   b2 =  -0.907591                                       ..........
.....................................................................................................                                          ..........

-In this example, the coefficients for  a1 & b1  are computed via the Taylor series ( 2.PI/48 < 1/6 )
whereas those for  a2 & b2  use the trigonometric expressions.

d-2)  n = multiple of 6

-To compute the integrals of the form    § y(x) cos(µ.x) dx   or   § y(x) sin(µ.x) dx  , we now seek a formula that produces perfect accuracy
when y(x) is a polynomial p(x) with deg p < 7

-More successive integrations by parts give:

§ab y(x) cos(µ.x) dx = B cos µ.b - A cos µ.a + D sin µ.b - C sin µ.a

§ab y(x) sin(µ.x) dx  = -D cos µ.b + C cos µ.a + B sin µ.b - A sin µ.a

A = y'(a)/µ2 - y'''(a)/µ4 + y'''''(a)/µ6 - ......              C = y(a)/µ - y''(a)/µ3 + y'''''(a)/µ5 - y''''''(a)/µ7  + ..........
with
B = y'(b)/µ2 - y'''(b)/µ4 + y'''''(b)/µ6  - ......             D = y(b)/µ - y''(b)/µ3 + y'''''(b)/µ5 - y''''''(b)/µ4 + ..........

-And we can use numerical differentiation to evaluate  A , B , C , D , namely:

y'(a) = [ -(49/20) y(a) + 6 y(a+h) - (15/2) y(a+2h) + (20/3) y(a+3h) - (15/4) y(a+4h) + (6/5) y(a+5h) - (1/6) y(a+6h) ] / h

y''(a) = [ (203/45) y(a) - (87/5) y(a+h) + (117/4) y(a+2h) - (254/9) y(a+3h) + (33/2) y(a+4h) - (27/5) y(a+5h) + (137/180) y(a+6h) ] / h2

y'''(a) = [ -(49/8) y(a) + 29 y(a+h) - (461/8) y(a+2h) +62 y(a+3h) - (307/8) y(a+4h) + 13 y(a+5h) - (15/8) y(a+6h) ] / h3

y''''(a) = [ (35/6) y(a) - 31 y(a+h) + (137/2) y(a+2h) - (242/3) y(a+3h) + (107/2) y(a+4h) - 19 y(a+5h) + (17/6) y(a+6h) ] / h4

y'''''(a) = [ -(7/2) y(a) + 20 y(a+h) - (95/2) y(a+2h) + 60 y(a+3h) - (85/2) y(a+4h) + 16 y(a+5h) - (5/2) y(a+6h) ] / h5

y''''''(a) = [  y(a) - 6 y(a+h) + 15 y(a+2h) - 20 y(a+3h) + 15 y(a+4h) - 6 y(a+5h) + y(a+6h) ] / h6 = y''''''(b)

and similar expressions for y'(b)  y''(b)  y'''(b)  y''''(b)  y'''''(b)

-Newton-Cote 7-point formula is used to compute the first coefficient a0 = the mean of y(x) over the interval [0,n]
-This integration formula also gives a perfect accuracy if deg p = 7

Data Registers:         •  R00 = n = mulptiple of 6              ( Registers R00 & R12 thru Rnn+12 are to be initialized before executing "FIL6" )

•   R12 = y(0+)   •   R13 = y(1)  ................  •   Rnn+12 = y(n-)          R01 to R11: temp   R08 = 360 k/n
Flags:  /
Subroutines:  /

-Lines 46-138 are three-byte  GTOs

 01  LBL "FIL6"   02  DEG   03  360   04  *   05  RCL 00   06  STO 11   07  /   08  STO 08   09  D-R   10  X^2   11  STO 09   12  12.011   13  ST+ 11   14  CLA   15  GTO 01   16  LBL 00   17  RCL 01   18  RCL 07   19  +   20  41   21  *   22  RCL 02   23  RCL 06   24  +   25  216   26  *   27  +   28  RCL 03   29  RCL 05   30  +   31  27   32  *   33  +   34  RCL 04   35  272   36  *   37  +   38  280   39  /   40  ST+ M   41  LBL 01   42  RCL IND 11   43  STO 07   44  DSE 11   45  FS? 30   46  GTO 04     47  RCL IND 11   48  STO 06   49  DSE 11   50  RCL IND 11   51  STO 05   52  DSE 11   53  RCL IND 11 54  STO 04   55  DSE 11   56  RCL IND 11   57  STO 03   58  DSE 11   59  RCL IND 11   60  STO 02   61  DSE 11   62  RCL IND 11   63  STO 01   64  RCL 08   65  X=0?   66  GTO 00   67  XEQ 02   68  RCL 11   69  INT   70  12   71  -   72  RCL 08   73  *   74  STO 10   75  X<>Y   76  P-R   77  ST- M   78  X<>Y   79  ST- N   80  RCL 02   81  RCL 06   82  +   83  6   84  *   85  RCL 01   86  -   87  RCL 03   88  RCL 05   89  +   90  15   91  *   92  -   93  RCL 04   94  20   95  *   96  +   97  RCL 07   98  -   99  RCL 09 100  / 101  STO O 102  XEQ 03 103  RCL 10 104  X<>Y 105  P-R 106  ST+ N 107  X<>Y 108  ST- M 109  RCL 01 110  X<> 07 111  STO 01 112  RCL 02 113  X<> 06 114  STO 02 115  RCL 03 116  X<> 05 117  STO 03 118  RCL O 119  XEQ 03 120  RCL 08 121  6 122  * 123  RCL 10 124  + 125  STO 10 126  X<>Y 127  P-R 128  ST- N 129  X<>Y 130  ST+ M 131  XEQ 02 132  RCL 10 133  X<>Y 134  P-R 135  ST- M 136  X<>Y 137  ST- N 138  GTO 01 139  LBL 02 140  RCL 02 141  40 142  * 143  RCL 01 144  7 145  * 146  - 147  RCL 03 148  95 149  * 150  - 151  RCL 04 152  120 153  * 154  + 155  RCL 05 156  85 157  * 158  - 159  RCL 06 160  32 161  * 162  + 163  RCL 07 164  5 165  * 166  - 167  RCL 09 168  ST+ X 169  / 170  RCL 01 171  49 172  * 173  RCL 02 174  232 175  * 176  - 177  RCL 03 178  461 179  * 180  + 181  RCL 04 182  496 183  * 184  - 185  RCL 05 186  307 187  * 188  + 189  RCL 06 190  104 191  * 192  - 193  RCL 07 194  15 195  * 196  + 197  8 198  / 199  + 200  RCL 09 201  / 202  RCL 02 203  360 204  * 205  RCL 01 206  147 207  * 208  - 209  RCL 03        210  450 211  * 212  - 213  RCL 04 214  400 215  * 216  + 217  RCL 05 218  225 219  * 220  - 221  RCL 06 222  72 223  * 224  + 225  RCL 07 226  10 227  * 228  - 229  60 230  / 231  + 232  RCL 09 233  / 234  RTN 235  LBL 03 236  RCL 01 237  35 238  * 239  RCL 02 240  186 241  * 242  - 243  RCL 03 244  411 245  * 246  + 247  RCL 04 248  484 249  * 250  - 251  RCL 05 252  321 253  * 254  + 255  RCL 06 256  114 257  * 258  - 259  RCL 07        260  17 261  * 262  + 263  6 264  / 265  + 266  RCL 09 267  / 268  RCL 01 269  812 270  * 271  RCL 02 272  3132 273  * 274  - 275  RCL 03 276  5265 277  * 278  + 279  RCL 04 280  5080 281  * 282  - 283  RCL 05 284  2970 285  * 286  + 287  RCL 06 288  972 289  * 290  - 291  RCL 07 292  137 293  * 294  + 295  PI 296  R-D 297  / 298  - 299  RCL 09        300  / 301  RCL 01 302  + 303  RCL 08 304  D-R 305  / 306  RTN 307  LBL 04 308  RCL N 309  RCL M 310  RCL 00 311  2 312  / 313  ST/ Z 314  / 315  CLA 316  END

( 451 bytes / SIZE nnn+13 )

 STACK INPUTS OUTPUTS Y / bk X k ak

Example:             y(x) is periodic ( period = 12 )

 x 0+ 1 2 3 4 5 6 7 8 9 10 11 12- y 1 2 4 7 10 12 12 8 6 7 10 14 20

n = 12   STO 00   and store these 13 y-values into  R12  R13  ..................  R24

0  XEQ "FIL6"  >>>>   a0 =   8.475595                                                                  --- execution time =  6s  ---
1        R/S         >>>>    a1 = -1.059182   X<>Y   b1 =  -1.090833                         --- execution time = 46s ---
2        R/S         >>>>    a2 =  1.776544   X<>Y   b2 =  -4.029470                                        ........
.....................................................................................................

6        R/S        >>>>    a6 =  0.171103   X<>Y   b6 =  -0.979868
7        R/S         >>>>    a7 =  0.048603   X<>Y   b7 =  -0.834744
......................................................................................................

Notes:

-The results are accurate when the y-values are small integers, but otherwise, unfortunately, roundoff-errors may be huge if µ = 2.k.PI/n is small.
and one can find examples where the 3rd decimal is already wrong!
-So, another approach is usually preferable:

-We assume hereunder that y(x) is also continuous for x = 0 , then the formulae may be re-written:

§0n y(x) cos µ.x dx = Sum i=0,6,12,.....,n-6          y(i) A0 cos µ.i       + y(i+1) [ A1 cos (µ.i+µ) + B1 sin (µ.i+µ) ]    + y(i+2) [ A2 cos (µ.i+2µ) + B2 sin (µ.i+2µ) ]

+ y(i+3) A3 cos (µ.i+3µ) + y(i+4) [ A2 cos (µ.i+4µ) - B2 sin (µ.i+4µ) ] + y(i+5) [ A1 cos (µ.i+5µ) - B1 sin (µ.i+5µ) ]

§0n y(x) sin µ.x dx = Sum i=0,6,12,.....,n-6          y(i) A0 sin µ.i       + y(i+1) [ A1 sin (µ.i+µ) - B1 cos (µ.i+µ) ]      + y(i+2) [ A2 sin (µ.i+2µ) - B2 cos (µ.i+2µ) ]

+ y(i+3) A3 sin (µ.i+3µ) + y(i+4) [ A2 sin (µ.i+4µ) + B2 cos (µ.i+4µ) ] + y(i+5) [ A1 sin (µ.i+5µ) + B1 cos (µ.i+5µ) ]
with

A0 = (1/3µ2-15/4µ4+5/µ6) cos 6µ + (49/10µ2-49/4µ4+7/µ6) + (-137/90µ3+17/3µ5-2/µ7) sin 6µ
A1 = (-6/5µ2+13/µ4-16/µ6) cos 5µ - (6/µ2-29/µ4+20/µ6) cos µ + (27/5µ3-19/µ5+6/µ7) sin 5µ + (87/5µ3-31/µ5+6/µ7) sin µ
A2 = (15/4µ2-307/8µ4+85/2µ6) cos 4µ - (-15/2µ2+461/8µ4-95/2µ6) cos 2µ + (-33/2µ3+107/2µ5-15/µ7) sin 4µ + (-117/4µ3+137/2µ5-15/µ7) sin 2µ
A3 =  2(-20/3µ2+62/µ4-60/µ6) cos 3µ + 2(254/9µ3-242/3µ5+20/µ7) sin 3µ

and

B1 = -(-6/5µ2+13/µ4-16/µ6) sin 5µ - (6/µ2-29/µ4+20/µ6) sin µ + (27/5µ3-19/µ5+6/µ7) cos 5µ - (87/5µ3-31/µ5+6/µ7) cos µ
B2 = -(15/4µ2-307/8µ4+85/2µ6) sin 4µ - (-15/2µ2+461/8µ4-95/2µ6) sin 2µ + (-33/2µ3+107/2µ5-15/µ7) cos 4µ - (-117/4µ3+137/2µ5-15/µ7) cos 2µ

-When µ is small, we must use Taylor series to reduce roundoff-errors, namely:

A0 = 41/70 + 9µ2/25 - 567µ4/550 + 28053µ6/25025 - 16281µ8/25025 + 498636µ10/2127125 - 0.057560072125µ12 + 0.010259175386µ14 - .............

A1 = 54/35 - 27µ2/25 + 1917µ4/1100 - 256947µ6/200200 + 809283µ8/1601600 - 33654633µ10/272272000 + 0.020638591235µ12
- 0.0025064271605µ14 +  .............

A2 =  27/140 + 27µ2/10 - 513µ4/220 + 9903µ6/10010 - 11859µ8/50050 + 76134µ10/2127125 - 0.0037150884891µ12
+ 0.000281765004581µ14 - .............

A3 = 68/35 - 18µ2/5 + 243µ4/110 - 2133µ6/4004 + 55161µ8/800800 - 763263µ10/136136000 + 0.000315589112894µ12
- 0.0000130645273645µ14 +  .............

B1 = 36µ3/25 - 18µ5/11 + 1508589µ7/1751750 - 5137µ9/19500 + 0.052739990401µ11 - 0.0074605626534µ13 +  .............

B2 = -9µ3/5 + 414µ5/275 - 88758µ7/175175 + 290µ9/3003 - 0.0120217037863µ11 + 0.00106051813438µ13 +  .............

-The coefficients in decimal form are only approximate values: the last decimals are not guaranteed...

-Taylor's series are used if µ < 0.43 ( line 33 ) so that Ai and Bj have at least 5 or 6 exact decimals.
-Choose another threshold and add more terms in the Taylor series if you need a better accuracy.

Data Registers:          •  R00 = n = mulptiple of 6              ( Registers R00 & R11 thru Rnn are to be initialized before executing "FIL6+" )

•   R11 = y(1)  ................  •   Rnn+10 = y(n) = y(0)         R01 to R09: temp      R10 is unused

Flags:  /
Subroutines:  /

-Lines 31-35-157-500 are three-byte  GTOs

 01  LBL "FIL6+"   02  DEG   03  CLA   04  360   05  *   06  RCL 00   07  /   08  STO 07   09  D-R   10  STO 08   11  X^2   12  STO 09   13  82   14  STO 01   15  216   16  STO 02   17  27   18  STO 03   19  272   20  STO 04   21  CLX   22  STO 05   23  STO 06   24  280   25  ST/ 01   26  ST/ 02   27  ST/ 03   28  ST/ 04   29  RCL 08   30  X=0?   31  GTO 01    32  ABS   33  .43   34  XY 223  STO 05 224  RCL 07 225  20 226  RCL 09 227  / 228  29 229  - 230  RCL 09 231  / 232  6 233  + 234  P-R 235  ST- 02 236  X<>Y 237  ST- 05 238  RCL 06 239  6 240  RCL 09 241  / 242  19 243  - 244  RCL 09 245  / 246  5.4 247  + 248  RCL 08 249  / 250  P-R 251  ST+ 05 252  X<>Y 253  ST+ 02 254  RCL 07        255  6 256  RCL 09 257  / 258  31 259  - 260  RCL 09 261  / 262  17.4 263  + 264  RCL 08 265  / 266  P-R 267  ST- 05 268  X<>Y 269  ST+ 02 270  RCL 07 271  4 272  * 273  STO 04 274  42.5 275  RCL 09 276  / 277  38.375 278  - 279  RCL 09 280  / 281  3.75 282  + 283  P-R 284  STO 03 285  X<>Y 286  CHS 287  STO 06 288  RCL 04 289  15 290  RCL 09 291  / 292  53.5 293  - 294  RCL 09 295  / 296  16.5 297  + 298  RCL 08 299  / 300  P-R 301  ST- 06 302  X<>Y 303  ST- 03 304  RCL 07 305  ST+ X 306  STO 04 307  47.5 308  RCL 09 309  / 310  57.625 311  - 312  RCL 09 313  / 314  7.5 315  + 316  P-R 317  ST+ 03 318  X<>Y 319  ST+ 06 320  RCL 04 321  15 322  RCL 09 323  / 324  68.5 325  - 326  RCL 09 327  / 328  29.25 329  + 330  RCL 08 331  / 332  P-R 333  ST+ 06 334  X<>Y 335  ST- 03 336  20 337  RCL 09        338  / 339  242 340  3 341  / 342  - 343  RCL 09 344  / 345  254 346  9 347  / 348  + 349  RCL 08 350  / 351  RCL 07 352  3 353  * 354  STO 04 355  SIN 356  * 357  60 358  RCL 09 359  / 360  62 361  - 362  RCL 09 363  / 364  20 365  3 366  / 367  + 368  RCL 04 369  COS 370  * 371  - 372  ST+ X 373  RCL 09 374  ST/ 01 375  ST/ 02 376  ST/ 03 377  ST/ 05 378  ST/ 06 379  / 380  STO 04 381  LBL 01 382  10 383  RCL 00 384  STO 08 385  + 386  STO 09 387  LBL 02 388  RCL 07 389  RCL 08 390  * 391  RCL IND 09 392  RCL 01 393  * 394  P-R 395  ST+ M 396  X<>Y 397  ST+ N 398  DSE 08 399  DSE 09 400  RCL 07 401  RCL 08 402  * 403  RCL IND 09 404  P-R 405  RCL X 406  RCL 02 407  * 408  ST+ M 409  CLX 410  RCL 05 411  * 412  ST+ N 413  X<> L 414  * 415  ST- M 416  CLX 417  RCL 02 418  * 419  ST+ N 420  DSE 08 421  DSE 09 422  RCL 07 423  RCL 08 424  * 425  RCL IND 09 426  P-R 427  RCL X 428  RCL 03 429  * 430  ST+ M 431  CLX 432  RCL 06 433  * 434  ST+ N 435  X<> L 436  * 437  ST- M 438  CLX 439  RCL 03 440  * 441  ST+ N 442  DSE 08 443  DSE 09 444  RCL 07 445  RCL 08 446  * 447  RCL IND 09 448  RCL 04 449  * 450  P-R 451  ST+ M 452  X<>Y 453  ST+ N 454  DSE 08 455  DSE 09 456  RCL 07 457  RCL 08 458  * 459  RCL IND 09 460  P-R 461  RCL X 462  RCL 03 463  * 464  ST+ M 465  CLX 466  RCL 06 467  * 468  ST- N 469  X<> L 470  * 471  ST+ M 472  CLX 473  RCL 03 474  * 475  ST+ N 476  DSE 08 477  DSE 09 478  RCL 07 479  RCL 08 480  * 481  RCL IND 09 482  P-R 483  RCL X 484  RCL 02 485  * 486  ST+ M 487  CLX 488  RCL 05 489  * 490  ST- N 491  X<> L 492  * 493  ST+ M 494  CLX 495  RCL 02 496  * 497  ST+ N 498  DSE 09 499  DSE 08 500  GTO 02 501  RCL N 502  RCL M 503  RCL 00 504  2 505  / 506  ST/ Z 507  / 508  CLA 509  END

( 796 bytes / SIZE nnn+11 )

 STACK INPUTS OUTPUTS Y / bk X k ak

Example:      y(x)  periodic ( period = 48 )  and  y(x) = (x/10) Ln(49-x)  over  [0,48]

-To store y(n) in the proper registers, load the same routine:

LBL 01     RCL Y      *              /                           DSE X                         and key in    48  XEQ 01
STO Y      -               10            STO IND Y         GTO 01
49             LN           ST+ Z       R^

-Then  n = 48  STO 00

0  XEQ "FIL6+"  >>>>    a0 =   6.083399                                                                  --- execution time = 48s  ---
1        R/S            >>>>    a1 = -2.250254   X<>Y   b1 =  -2.782055                         --- execution time = 71s  ---
2        R/S            >>>>    a2 = -0.940317   X<>Y   b2 =  -0.907481                                       ..........
.....................................................................................................                                          ..........

7         R/S           >>>>    a7 = -0.158130   X<>Y   b7 =  -0.086882                                       ..........

Notes:

-If we use "FIL6" instead of "FIL6+", a0 is correctly computed but "FIL6" gives  a1 = -2.243751  &  b1 = -2.784534 ,  a2 = -0.940469  &  b2 =  -0.907335
-However, a7 & b7 are accurately calculated ... but in 3mn07s.
-If you don't want to use synthetic registers, replace registers M & N by R10 & R11
delete line 508, replace line 382 by 11 ( instead of 10 ), add  STO 10   STO 11  after line 23, and delete line 03.
-And store  y(1) , ..... , y(n)  into  R12 , ........... , Rnn+11

A few Comparisons:

 ASP FIL FIL4+ FIL6+ exact(6D) a0 6.074830 6.083005 6.083282 6.083399 6.083605 a1 -2.267371 -2.251065 -2.250492 -2.250254 -2.249808 b1 -2.781891 -2.782237 -2.782105 -2.782055 -2.781991 a2 -0.957391 -0.941185 -0.940576 -0.940317 -0.939786 b2 -0.907200 -0.907842 -0.907590 -0.907481 -0.907391 a7 -0.175917 -0.160233 -0.158926 -0.158130 -0.157669 b7 -0.086422 -0.087694 -0.086800 -0.086882 -0.087134

2°) One-dimensional Problems - Complex Arguments

a) The Discrete Fourier Transform

-Given an array of n complex numbers  [ a1 , .......... , an ]  where  ak = bk + i ck ,
"DFTZ" computes the complex components zk = xk + i yk  of the Discrete Fourier Transform   [ z1 , .......... , zn ]

Data Registers:         •  R00 = n                                            ( Registers R00 thru R2n are to be initialized before executing "DFTZ" )

•   R01 = b1   •   R03 = b2  ................  •   R2n-1 = bn
•   R02 = c1   •   R04 = c2  ................   •   R2n = cn
Flags: /
Subroutines:  /

 01  LBL "DFTZ"  02  DEG 03  CLA 04  STO O 05  PI 06  R-D 07  ST* Y 08  - 09  RCL 00        10  STO M 11  ST+ M 12  / 13  STO N 14  CLX 15  LBL 01 16  RCL M 17  RCL N 18  ST* Y 19  ST+ X 20  - 21  RCL IND M 22  DSE M 23  RCL IND M 24  R-P 25  RDN 26  - 27  R^ 28  P-R 29  ST+ P 30  RDN 31  - 32  DSE M 33  GTO 01 34  RCL O         35  SIGN 36  X<> P 37  CLA 38  END

( 64 bytes / SIZE 2n+1 )

 STACK INPUTS OUTPUTS Y / yk X k xk L / k

Example:    Calculate  DFT  [ 1+2i  3+4i  5+6i  7+8i ]

n = 4  STO 00

1  STO 01    3  STO 03    5  STO 05    7  STO 07
2  STO 02    4  STO 04    6  STO 06    8  STO 08

1  XEQ "DFTZ"  >>>>   x1 = 16   X<>Y   y1 =  20
2       R/S            >>>>   x2 = -8    X<>Y   y2 =   0
3       R/S            >>>>   x3 = -4    X<>Y   y3 =  -4
4       R/S            >>>>   x4 =  0    X<>Y   y4 =  -8

-So,   DFT [ 1+2i  3+4i  5+6i  7+8i ]  =  [ 16+20.i  -8  -4-4.i   -8.i ]

-The discrete inverse Fourier transform may be obtained by the following routine:

 01  LBL "IDFTZ" 02  DEG 03  CLA 04  STO O 05  PI 06  R-D 07  ST* Y 08  - 09  RCL 00         10  STO M 11  ST+ M 12  / 13  STO N 14  CLX 15  LBL 01 16  RCL M         17  RCL N 18  ST* Y 19  ST+ X 20  - 21  RCL IND M 22  DSE M 23  RCL IND M 24  R-P 25  RDN 26  + 27  R^ 28  P-R 29  ST+ P 30  RDN 31  + 32  DSE M         33  GTO 01 34  RCL 00 35  ST/ P 36  / 37  RCL O         38  SIGN 39  X<> P 40  CLA 41  END

( 69 bytes / SIZE 2n+1 )

 STACK INPUTS OUTPUTS Y / yk X k xk L / k

Example:   Calculate  IDFT [ 16+20.i  -8  -4-4.i   -8.i ]

n = 4  STO 00

16  STO 01    -8  STO 03    -4  STO 05     0  STO 07
20  STO 02     0  STO 04    -4  STO 06    -8  STO 08

1  XEQ "IDFTZ"   >>>>   x1 = 1   X<>Y   y1 =  2
2       R/S              >>>>   x2 = 3    X<>Y   y2 =  4
3       R/S              >>>>   x3 = 5    X<>Y   y3 =  6
4       R/S              >>>>   x4 = 7    X<>Y   y4 =  8

-So,   IDFT  [ 16+20.i  -8  -4-4.i   -8.i ] = [ 1+2i  3+4i  5+6i  7+8i ]    which is quite logical!

b) The Fast Fourier Transform  ( X-Functions Module required )

-"FFT" computes the Discrete Fourier Transform   [ z1 , .......... , zn ]       zk = xk + i yk
of  an array of n complex numbers  [ a1 , .......... , an ]      ak = bk + i ck

>>>   n must be an integer power of 2 , n > 1

Data Registers:         •  R00 = n = 2j > 1                              ( Registers R00 thru R2n are to be initialized before executing "FFT" )

Inputs:                  •   R01 = b1   •   R03 = b2  ................  •   R2n-1 = bn
•   R02 = c1   •   R04 = c2  ................   •   R2n = cn                   R2n+1 to R4n:   temp

Outputs:               •   R01 = x1   •   R03 = x2  ................   •   R2n-1 = xn
•   R02 = y1   •   R04 = y2  ................   •   R2n = yn
Flags: /
Subroutines:  /

 01  LBL "FFT"    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  ST+ X   35  RCL M   36  ST+ X   37  RCL IND X   38  STO IND Z   39  CLX   40  SIGN   41  ST- Z   42  -   43  RCL IND X   44  STO IND Z   45  DSE M   46  GTO 01   47  360   48  STO O   49  SIGN   50  STO M 51  GTO 04   52  LBL 03   53  RCL 00   54  ST+ X   55  .1   56  STO Z   57  %   58  ISG X   59  +   60  %   61  1   62  +   63  REGMOVE   64  LBL 04   65  2   66  ST/ O   67  RCL 00   68  STO P   69  LBL 05   70  RCL P   71  RCL X   72  RCL O   73  STO Q   74  SIGN   75  - 76  ENTER^   77  ST* Q   78  RCL M   79  ST+ X   80  MOD   81  ST- Z   82  CLX   83  RCL M   84  MOD   85  +   86  RCL 00   87  +   88  ST+ X   89  RCL X   90  RCL M   91  ST+ X   92  ST+ Z   93  X<> Q   94  RCL IND Z   95  DSE T   96  RCL IND T   97  R-P   98  X<> Z   99  - 100  X<>Y 101  P-R 102  DSE T 103  RCL IND T 104  + 105  X<>Y 106  RCL IND Z 107  + 108  RCL P 109  ST+ X 110  RDN 111  STO IND T 112  DSE T 113  X<>Y 114  STO IND T 115  DSE P 116  GTO 05 117  RCL M 118  ST+ M 119  DSE N 120  GTO 03 121  CLA 122  END

( 197 bytes / SIZE 4n+1 )

 STACK INPUT OUTPUT X / /

Example:    Calculate  DFT [ 1+2i  3+4i  5+6i  7+8i ]

n = 4  STO 00

1  STO 01    3  STO 03    5  STO 05    7  STO 07
2  STO 02    4  STO 04    6  STO 06    8  STO 08

XEQ "FFT"   and 27 seconds later, we have in registers R01 thru R08

16   20    -8   0    -4   -4    0   -8       thus,   DFT [ 1+2i  3+4i  5+6i  7+8i ]  =  [ 16+20.i  -8  -4-4.i   -8.i ]

-In fact, I've rounded these values, for example  R04 = 0.000000006 instead of exactly 0. Roundoff errors are unavoidable...

Execution time:

 n 2 4 8 16 32 64 t 8s 28s 1m14s 3m11s 7m57s 18m54s

-It's of course not so fast as the built-in  FFT on the HP-48, but it's faster than  n x "DFTZ"
- t is approximately proportional to  n Log n

Notes:

-This program doesn't work if n = 1, but in this case,  DFT [z] = [z]
-n must be an integer power of 2.
-Since there are at most 319 registers, nmax = 64

-If you want to calculate the inverse Fourier transform IFFT,

replace line 99 by  +  ( instead of - )   and divide the registers R01 thru R2n by n:  for instance,  add

RCL 00  ENTER^  ST+ Y  LBL 06  ST/ IND Y  DSE Y  GTO 06   after line 120

3°) Two-dimensional Problems - Real Arguments

a) The Discrete Fourier Transform ( Program#1 )

-The DFT may be generalized to  nxm matrices  [ ai,j ]  ( 1 <= i <= n , 1 <= j <= m ). The result is a complex  nxm matrix  [ zi,j ]

where  zi,j =  SUMq=1,2,...,n  ; r=1,2,....,m   aq,r  exp { -2i(pi).[ (q-1)(i-1)/n+(r-1)(j-1)/m ] }

Data Registers:     •  R00 = n.mmm  =  n + m/1000                  ( Registers R00 thru Rnm are to be initialized before executing "DFT2" )

•  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
Flags: /
Subroutines: /

 01  LBL "DFT2" 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  STO P 23  LBL 01        24  RCL N 25  RCL M 26  ENTER^ 27  SIGN 28  - 29  ST* Y 30  RCL 00 31  INT 32  / 33  INT 34  RCL O 35  * 36  + 37  RCL IND M 38  P-R 39  ST+ P 40  RDN 41  - 42  DSE M 43  GTO 01        44  RCL P 45  CLA 46  END

( 75 bytes / SIZE n.m+1 )

 STACK INPUTS OUTPUTS Y j yi,j X i xi,j

with   zi,j = xi,j + i yi,j

Example:

[ [  1  3  4  10  ]
A =   [  4  5  7  14  ]
[  2  9  6  11  ] ]

-Calculate   z2,3

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 "DFT2"  >>>>      2                                                                                       ( in 14 seconds )
X<>Y    -3.464101617       So  z2,3 =  2 - 3.464101617 i
-Likewise,

[ [         76                -10+18 i                       -28                     -10-18 i         ]
DFT (A) =     [   -11-1.7321 i     6.5622+0.6340 i        2-3.4641 i       -5.5622-2.3660 i   ]        ( rounded to 4 decimals )
[  -11+1.7321 i    -5.5622+2.3660 i       2+3.4641 i       6.5622-0.6340 i   ] ]

Note:

-To compute the Inverse Fourier Transform:

add   RCL a    ST/ Z   /     after line 44
replace line 41 by +
add   RCL M   STO a      after line 20     ( 51 lines / 85 bytes )

-However, since status register a is used, this routine cannot be called as more than a first-level subroutine.
-If you want to avoid register a, use the following variant:

 01  LBL "IDFT2" 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  RCL M 22  STO P 23  CLST 24  LBL 01 25  RCL M 26  ENTER^ 27  SIGN 28  - 29  STO Q        30  RCL N 31  * 32  X<> Q 33  RCL 00 34  INT 35  / 36  INT 37  RCL O 38  * 39  RCL Q 40  + 41  RCL IND M 42  P-R 43  ST+ T 44  RDN 45  + 46  DSE M 47  GTO 01        48  X<>Y 49  RCL P 50  ST/ Z 51  / 52  CLA 53  END

( 87 bytes / SIZE n.m+1 )

 STACK INPUTS OUTPUTS Y j yi,j X i xi,j

with   zi,j = xi,j + i yi,j

-The instructions are identical.

b) The Discrete Fourier Transform ( Program#2 )

-We assume  z(x,y)  is a doubly periodic function:

z(x+n;y) = z(x;y)   for all  x;y
z(x;y+m) = z(x;y)  for all  x;y           n and m  being integers.

-We have   z(x;y) = SUMi;j  [ ai;j cos ( 360 ( ix/n+jy/m ) )  + bi;j sin ( 360 ( ix/n+jy/m ) )  ]

( 0 <= i <= n/2 ; 0 <= j <= m/2 and  0 < i < n/2  ;  -m/2 < j < 0 )

Data Registers:          •  R00 = n.mmm  =  n + m/1000                ( Registers  R00 thru Rnm   are to be initialized before executing "SPA2" )

•  R01 = z(1;1)     •  Rn+1 = z(1;2)    ..................................................     •  Rnm-n+1 = z(1;m)
•  R02 = z(2;1)     •  Rn+2 = z(2;2)    ..................................................     •  Rnm-n+2 = z(2;m)

...............................................................................................................................................

•  Rnn = z(n;1)      •  R2n  = z(n;2)     ..................................................     •     Rnm      = z(n;m)

Flag:  F10
Subroutines:  /

-Synthetic registers M N O P a   are used by this program and may be replaced by any unused data registers.
-Since synthetic register a   is used, "SPA2" must not be called as more than a first level subroutine.

 01  LBL "SPA2" 02  DEG 03  SF 10 04  PI 05  R-D 06  STO N 07  STO O 08  CLX 09  2 10  ST* Z 11  * 12  ST* N 13  RCL 00 14  INT 15  ST/ N 16  MOD 17  X<>Y 18  ST* O 19  RCL 00        20  FRC 21   E3 22  * 23  ST/ O 24  MOD 25  LASTX 26  RCL 00 27  INT 28  * 29  STO M 30  STO a 31  RDN 32  + 33  X=0? 34  CF 10 35  CLX 36  STO P  37  LBL 01 38  RCL N        39  RCL M 40  ST* Y 41  ENTER^ 42  SIGN 43  - 44  RCL 00 45  INT 46  / 47  INT 48  RCL O 49  ST* Y 50  + 51  + 52  RCL IND M 53  P-R 54  ST+ P 55  RDN 56  + 57  DSE M 58  GTO 01 59  RCL P        60  R-P 61  RCL a 62  / 63  FS?C 10 64  ST+ X 65  P-R 66  CLA 67  END

( 102 bytes / SIZE n.m+1 )

 STACK INPUTS OUTPUTS Y j bi,j X i ai,j

( 0 <= i <= n/2 ;  0 <= j <= m/2
0 < i < n/2    ;   -m/2 < j < 0 )

Example:        z(x;y)  is doubly periodic  with    n = 3   ;   m = 4    and takes the following values:

z(1;1) = 1   z(1;2) = 3  z(1;3) = 7  z(1;4) = 10
z(2;1) = 4   z(2;2) = 5  z(2;3) = 8  z(2;4) = 14
z(3;1) = 2   z(3;2) = 9  z(3;3) = 6  z(3;4) = 11

3.004   STO 00    and we store the 12 numbers:

1   3   7   10            R01   R04   R07   R10
4   5   8   14    in     R02   R05   R08   R11      respectively
2   9   6   11           R03   R06   R09   R12

-There are seven (i;j) pairs satisfying the required properties  ( 0 <= i <= 1.5 ; 0 <= j <= 2   and  0 < i < 1.5 ; -2 < j < 0 )
namely:     (0;0) (0;1) (0;2) ,  (1;0) (1;1) (1;2)  and  (1;-1)

0  ENTER^    XEQ "SPA2"    >>>   a0,0 =   6.6667
1  ENTER^       0    R/S           >>>   a0;1 =      3         X<>Y   b0;1 = -2.3333
2  ENTER^       0    R/S           >>>   a0;2 =      2         X<>Y   b0;2 =      0
0  ENTER^       1    R/S           >>>   a1;0 =   0.3333   X<>Y    b1;0 = -1.4434
1  ENTER^             R/S           >>>   a1;1 =  -0.7113   X<>Y    b1;1 = -0.1220
2  ENTER^       1    R/S           >>>   a1;2 =      1          X<>Y   b1;2 = -0.2887
-1  ENTER^      1    R/S           >>>   a1;-1 = -1.2887    X<>Y   b1;-1 = -0.4553

-Actually,

z(x;y) =    20/3   +  3 cos ( 90y ) - 7/3 sin ( 90y )    +  2 cos ( 180y )
+ 1/3 cos ( 120x ) - 5/r sin ( 120x )   +  ( -1 + 1/r ) cos ( 120x + 90y ) + ( 1/6 - 1/r ) sin ( 120x + 90y )  +   cos ( 120x + 180y ) - 1/r sin ( 120x + 180y )
+ ( -1 - 1/r ) cos ( 120x - 90y ) + ( -1/6 - 1/r ) sin ( 120x - 90y )

where  r is the square root of 12

4°) Two-dimensional Problems - Complex Arguments

-"DFTZ2"  calculates the components  zi,j = xi,j + i yi,j   of  DFT [ ai,j ]   with ai,j = bi,j + i ci,j   ( 1 <= i <= n , 1 <= j <= m ).

Data Registers:     •  R00 = n.mmm  =  n + m/1000             ( Registers R00 thru R2nm are to be initialized before executing "DFTZ2" )

•  R01 = b1,1      •  R2n+1 = b1,2   .......................................   •  R2nm-2n+1 = b1,m
•  R02 = c1,1      •  R2n+2 = c1,2   .......................................   •  R2nm-2n+2 = c1,m

•  R03 = b2,1      •  R2n+3 = b2,2   .......................................   •  R2nm-2n+3 = b2,m
•  R04 = c2,1      •  R2n+4 = c2,2   .......................................   •  R2nm-2n+4 = c2,m
..........................................................................................................
..........................................................................................................

•  R2n-1 = bn,1  •  R4n-1  =  bn,2   .......................................  •  R2nm-1 = bn,m
•  R2n  = cn,1     •  R4n    =  cn,2   .......................................    •  R2nm   = cn,m
Flags: /
Subroutines: /

 01  LBL "DFTZ2" 02  DEG 03  360 04  ST* Y 05  ST* Z 06  ST- Z 07  - 08  RCL 00 09  INT 10  ST+ X 11  STO M 12  / 13  STO N         14  CLX 15  RCL 00 16  FRC 17   E3 18  * 19  ST* M 20  / 21  STO O 22  CLX 23  STO P 24  LBL 01         25  RCL N 26  RCL M 27  ENTER^ 28  SIGN 29  ST+ X 30  - 31  ST* Y 32  RCL 00 33  INT 34  ST+ X 35  / 36  INT 37  RCL O 38  * 39  + 40  RCL IND M 41  DSE M 42  RCL IND M 43  R-P 44  RDN 45  - 46  R^ 47  P-R 48  ST+ P 49  RDN 50  - 51  DSE M 52  GTO 01         53  RCL P 54  CLA 55  END

( 90 bytes / SIZE 2n.m+1 )

 STACK INPUTS OUTPUTS Y j yi,j X i xi,j

with   zi,j = xi,j + i yi,j

Example:

[ [  1+2.i  3+4.i  4+5.i  6+7.i  ]
A =   [  4+6.i  5+7.i  7+8.i  3+4.i  ]
[  2+3.i  9+7.i  6+5.i  1+4.i  ] ]

-Calculate   z2,3

We have 3 rows and 4 columns , therefore  3.004  STO 00

1  2    3  4    4  5    6  7              R01 R02   R07 R08   R13 R14   R19 R20
Store:    4  6    5  7    7  8    3  4     in      R03 R04   R09 R10   R15 R16   R21 R22      respectively
2  3    9  7    6  5    1  4              R05 R06   R11 R12   R17 R18   R23 R24

3  ENTER^
2  XEQ "DFTZ2"  >>>>>    0.696152398            --- execution time = 24s ---
X<>Y    -8.330126900

-So    z2,3 = 0.696152398 - 8.330126900 i

-Likewise,

[ [         51+62 i                  -7-14 i                   -3-4.i                      -13              ]
DFT [A] =   [  0.6962-4.8660 i   -0.3038+6.1340 i    0.6962-8.3301 i     1.3038-9.8660 i  ]       rounded to 4 decimals
[-9.6962-3.1340 i  -10.6962+7.8660 i  -9.6962+0.3301 i  11.6962-8.1340 i  ] ]

Note:

-To calculate the inverse Fourier transform:

add   RCL a    ST/ Z   /     after line 53
replace lines 45 and 50 by  +  ( instead of - )
add   ST* a   after line 19
and add   STO a   after line 09

5°) Three-dimensional Problems - Real Arguments

-Similar programs can be used to compute the DFT of a real hypermatrix  [ ai,j,k ]    ( 1 <= i <= n , 1 <= j <= m , 1 <= k <= p ).
-The result is a complex hypermatrix  [ zi,j,k ]   with  zi,j,k = xi,j,k + i yi,j,k

zi,j,k  = SUMq=1,2,...,n  ; r=1,2,....,m ; s=1,2,.....,p   aq,r,s  exp { -2i(pi).[ (q-1)(i-1)/n+(r-1)(j-1)/m+(s-1)(k-1)/p ] }

Data Registers:     •  R00 = n.mmmppp  =  n + m/1000 + p/1000000       ( Registers R00 thru Rnmp are to be initialized before executing "DFT3" )

•  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
Flags: /
Subroutines:  /

 01  LBL "DFT3" 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 35  CLST 36  LBL 01        37  RCL M 38  ENTER^ 39  SIGN 40  - 41  STO Q 42  RCL a 43  / 44  INT 45  RCL P 46  * 47  RCL N 48  X<> Q 49  ST* Q 50  X<> Q 51  + 52  X<> Q 53  RCL 00        54  INT 55  / 56  INT 57  RCL O 58  * 59  RCL Q 60  + 61  RCL IND M 62  P-R 63  ST+ T 64  RDN 65  - 66  DSE M 67  GTO 01 68  X<>Y 69  CLA 70  END

( 112 bytes / SIZE nmp+1 )

 STACK INPUTS OUTPUTS Z k / Y j yi,j,k X i xi,j,k

with   zi,j,k = xi,j,k + i yi,j,k

( 1 <= i <= n , 1 <= j <= m , 1 <= k <= p )

Example:

A =  [ [ [   1  25  40   5    2  ]                        To store these 60 coefficients into the proper registers, you can key in    60   XEQ 01
[   4  36  18  32  37 ]                         after loading this short routine:
[   9   8   39  20  33 ]
[ 16  23  21  10  31 ] ]                       LBL 01   ENTER^   X^2   41   MOD   STO IND Y   X<>Y   DSE X   GTO 01
[ [  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   z2,4,3

-We have a   4x5x3  hypermatrix  therefore:   4.005003   STO 00

-Then,

3  ENTER^
4  ENTER^
2  XEQ "DFT3"  >>>>>   38.01062796                              --- Execution time = 84 seconds ---
X<>Y    10.96762920

-Thus,     z2,4,3 =  38.01062796 + 10.96762920 i

Note:

-To compute the inverse Fourier transform:

replace line 65  by   +
and add   RCL 00    E6    ST* Y    SQRT    MOD    RCL a     *     ST/ Z     /      after line 68

6°) Three-dimensional Problems - Complex Arguments

-"DFTZ3"  computes the DFT of a complex hypermatrix  [ ai,j,k ]  with ai,j,k = bi,j,k + i ci,j,k    ( 1 <= i <= n , 1 <= j <= m , 1 <= k <= p ).
-The result is also a complex hypermatrix  [ zi,j,k ]   with  zi,j,k = xi,j,k + i yi,j,k

Data Registers:     •  R00 = n.mmmppp  =  n+m/1000+p/1000000       ( Registers R00 thru R2nmp are to be initialized before executing "DFTZ3" )

•  R01 = b1,1,1      ...................................................................   •  R2nm-2n+1 = b1,m,1
•  R02 = c1,1,1      ...................................................................   •  R2nm-2n+2 = c1,m,1
..........................................................................................................
..........................................................................................................

•  R2n-1 = bn,1     ....................................................................  •  R2nm-1 = bn,m
•  R2n  = cn,1        ..................................................................    •  R2nm   = cn,m

.......................................................................................................................
.......................................................................................................................

...................................................................................................................
...................................................................................................................

•  R2nm(p-1)+1 = b1,1,p    .....................................................................   •  R2nmp-2n+1 = b1,m,p
•  R2nm(p-1)+2 = c1,1,p    .....................................................................   •  R2nmp-2n+2 = c1,m,p
..........................................................................................................
..........................................................................................................

•  R2nm(p-1)+2n-1 = bn,1,p  ....................................................................  •  R2nmp-1 = bn,m,p
•  R2nm(p-1)+2n  = cn,1,p     ....................................................................  •  R2nmp   = cn,m,p
Flags: /
Subroutines: /

 01  LBL "DFTZ3" 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  ST+ X 13  STO M 14  / 15  STO N 16  CLX 17  RCL 00 18  FRC 19   E3 20  * 21  INT 22  ST* M 23  / 24  STO O 25  CLX 26  RCL 00        27   E6 28  ST* Y 29  SQRT 30  MOD 31  ST* M 32  / 33  STO P 34  CLX 35  STO a 36  LBL 01 37  RCL M 38  ENTER^ 39  SIGN 40  ST+ X 41  - 42  STO Q 43  RCL 00        44  FRC 45  PI 46  INT 47  10^X 48  * 49  INT 50  RCL 00 51  INT 52  * 53  ST+ X 54  / 55  INT 56  RCL P 57  * 58  RCL Q 59  RCL N 60  * 61  + 62  RCL Q 63  RCL 00 64  INT 65  ST+ X 66  / 67  INT 68  RCL O 69  * 70  + 71  RCL IND M 72  DSE M 73  RCL IND M 74  R-P 75  RDN 76  - 77  R^ 78  P-R 79  ST+ a 80  RDN 81  - 82  DSE M 83  GTO 01 84  RCL a 85  CLA 86  END

( 131 bytes / SIZE 2nmp+1 )

 STACK INPUTS OUTPUTS Z k / Y j yi,j,k X i xi,j,k

with   zi,j,k = xi,j,k + i yi,j,k

( 1 <= i <= n , 1 <= j <= m , 1 <= k <= p )

Example:

A =  [ [ [   1+4.i   40+18.i   2+37.i   10+20.i   23+8.i  ]
[  9+16.i  39+21.i  33+31.i   32+5.i   36+25.i ]
[ 25+36.i  5+32.i   31+33.i   21+39.i  16+9.i  ]
[  8+23.i  20+10.i   37+2.i    18+40.i    4+i     ] ]
[ [      i       23+40.i    10+2.i     2+10.i   40+23.i  ]
[   4+9.i   18+39.i   37+33.i   20+32.i   8+36.i   ]
[ 16+25.i  21+5.i    31+31.i    5+21.i   25+16.i  ]
[  36+8.i  32+20.i   33+37.i   39+18.i    9+4.i    ] ]
[ [      1        8+23.i    20+10.i    37+2.i   18+40.i  ]
[   1+4.i    40+18.i    2+37.i    10+20.i   23+8.i   ]
[  9+16.i   39+21.i   33+31.i    32+5.i   36+25.i  ]
[ 25+36.i   5+32.i    31+33.i   21+39.i   16+9.i   ]

-To store the 60 complex coefficients of this 4x5x3 hypermatrix into the proper registers, you can use the routine ( SIZE 121 ):

LBL 01       41                     X<>Y                     Then key in   120  XEQ 01   in RUN mode
ENTER^    MOD                DSE X
X^2           STO IND Y      GTO 01

-Evaluate    z2,4,3

n.mmmppp = 4.005003   STO 00

3   ENTER^
4   ENTER^
2   XEQ "DFTZ3"  >>>>    20.47040998                    --- Execution time = 2mn31s ---
X<>Y   -159.3203501

-Thus,     z2,4,3 =  20.47040998 - 159.3203501 i

Note:

-To compute the inverse Fourier transform:

replace lines  76 and 81  by   +
and add   RCL 00    FRC    E3    *    INT   ST- L    E3    ST* L    X<> L    *    RCL 00    INT     *     ST/ Z     /      after line 84

References:

[1]  Abramowitz & Stegun - "Handbook of Mathematical Functions" - Dover Publications - ISBN 0-486-61272-4
[2]  Ronald N. Bracewell - "The Fourier Transform and its Applications" - McGraw-Hill - ISBN 0-07-116043-4
[3]  F.Scheid - "Numerical Analysis" - Mc Graw Hill - ISBN  07-055197-9