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 [ a₁ , .......... , a_n ] , "DFT" computes the complex components z_k = x_k + i y_k
of the Discrete Fourier Transform [ z₁ , .......... , z_n ]

-We have z_k = Sum_j=1,2,...,n a_j 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 = a₁ • R02 = a₂ ................ • Rnn = a_n
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	/	y_k
X	k	x_k
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" >>>>   x₁ = 12   X<>Y   y₁ = 0
      2       R/S          >>>>   x₂ = 4    X<>Y   y₂ = -2
      3       R/S          >>>>   x₃ = 4    X<>Y   y₃ = 0
      4       R/S          >>>>   x₄ = 4    X<>Y   y₄ = 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	/	y_k
X	k	x_k
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" >>>>   x₁ = 3    X<>Y   y₁ =   0
      2       R/S           >>>>   x₂ = 1    X<>Y   y₂ = 0.5
      3       R/S           >>>>   x₃ = 1    X<>Y   y₃ =   0
      4       R/S           >>>>   x₄ = 1    X<>Y   y₄ = -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) = a₀ + a₁.cos ( 360° x / n ) + b₁.sin ( 360° x / n ) + .... + a_k.cos ( 360° k.x / n ) + b_k.sin ( 360° k.x / n ) + ....
( n = period: y(x+n) = y(x) for all x )

-"SPA" computes the coefficients a₀ , a₁ , b₁ , ..... , a_k , b_k ( 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	/	b_k
X	k	a_k
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	/	b_k
X	k	a_k
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:

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

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

A = y'(a)/µ² - y'''(a)/µ⁴ + ...... C = y(a)/µ - y''(a)/µ³ + y'''''(a)/µ⁵ - ..........
with
B = y'(b)/µ² - y'''(b)/µ⁴ + ...... D = y(b)/µ - y''(b)/µ³ + y'''''(b)/µ⁵ - ..........

-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) ] / (12h²) 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) ] / (2h³)

y''''(a) = [ y(a) - 4 y(a+h) + 6 y(a+2h) - 4 y(a+3h) + y(a+4h) ] / h⁴ = 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 a₀ = 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	/	b_k
X	k	a_k

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" >>>>   a₀ =   8.455556                                                                  --- execution time = 6s ---
1        R/S         >>>>    a₁ = -1.049462   X<>Y   b₁ = -1.073589                         --- execution time = 37s ---
2        R/S         >>>>    a₂ = 1.838711   X<>Y   b₂ = -4.050732                                        idem
   .....................................................................................................

6        R/S         >>>>    a₆ = 0.157933   X<>Y   b₆ = -0.975730             Unlike "ASP" & "FIL",
7        R/S         >>>>    a₇ = 0.118962   X<>Y   b₇ = -0.861560             "FIL4" may produce a better approximation than b_n/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:

§₀ⁿ y(x) cos µ.x dx = Sum_{i=0,4,8,.....,n-4} y(i) A₀ cos µ.i + y(i+1) [ A₁ cos (µ.i+µ) + B₁ sin (µ.i+µ) ] + y(i+2) A₂ cos (µ.i+2µ)
+ y(i+3) [ A₁ cos (µ.i+3µ) - B₁ sin (µ.i+3µ) ]

     §₀ⁿ y(x) sin µ.x dx = Sum_{i=0,4,8,.....,n-4} y(i) A₀ sin µ.i + y(i+1) [ A₁ sin (µ.i+µ) - B₁ cos (µ.i+µ) ] + y(i+2) A₂ sin (µ.i+2µ)
                                                                                     + y(i+3) [ A₁ sin (µ.i+3µ) + B₁ cos (µ.i+3µ) ]
   with

     A₀ = (1/2µ²-3/µ⁴) cos 4µ + (25/6µ²-5/µ⁴) + (-11/6µ³+2/µ⁵) sin 4µ
     A₁ = (-4/3µ²+7/µ⁴) cos 3µ + (-4/µ²+9/µ⁴) cos µ + (14/3µ³-4/µ⁵) sin 3µ + (26/3µ³-4/µ⁵) sin µ
     A₂ = (6/µ²-24/µ⁴) cos 2µ + (-19/µ³+12/µ⁵) sin 2µ

and B₁ = (4/3µ²-7/µ⁴) sin 3µ + (-4/µ²+9/µ⁴) sin µ + (14/3µ³-4/µ⁵) cos 3µ + (-26/3µ³+4/µ⁵) 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:

            A₀ = 28/45 + 16µ²/63 - 128µ⁴/405 + 13568µ⁶/93555 - .............
            A₁ = 64/45 - 32µ²/63 + 712µ⁴/2835 - 1124µ⁶/18711 + .............
            A₂ = 8/15 + 16µ²/21 - 176µ⁴/945 + 544µ⁶/31185 - .............

B₁ = 64µ³/189 - 1888µ⁵/14175 + 2152µ⁷/93555 - .............

-In the following routine, the Taylor series are used for | µ | < 1/6 ( lines 30-31 ) and the terms of orders above µ⁵ 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 = a_k , R02 = b_k

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 X<Y?
33 GTO 00
34 RCL 05
35 RCL 05
36 128
37 *
38 CHS
39 405
40 /
41 16
42 63
43 /
44 STO 09
45 +
46 *
47 RCL 06
48 ST+ X
49 +

50 STO 06
51 CLX
52 712
53 *
54 2835
55 /
56 RCL 09
57 ST+ X
58 -
59 *
60 RCL 07
61 ST+ X
62 +
63 STO 07
64 CLX
65 176
66 *
67 CHS
68 945
69 /
70 16
71 21
72 /
73 +
74 *
75 RCL 08
76 ST+ X
77 +
78 STO 08
79 CLX
80 1888
81 *
82 14175
83 /
84 64
85 189
86 /
87 -
88 *
89 RCL 04
90 *
91 STO 09
92 GTO 01
93 LBL 00
94 5
95 RCL 05
96 /
97 CHS
98 .24

99 1/X
100 +
101 3
102 RCL 05
103 /
104 .5
105 -
106 RCL 03
107 4
108 *
109 STO 09
110 COS
111 *
112 -
113 2
114 RCL 05
115 /
116 11
117 6
118 /
119 -
120 RCL 04
121 /
122 RCL 09
123 SIN
124 *
125 +
126 STO 06
127 RCL 03
128 3
129 *
130 STO 08
131 7
132 RCL 05
133 /
134 .75
135 1/X
136 -
137 P-R
138 STO 07
139 X<>Y
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	/	b_k
X	k	a_k

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+" >>>>    a₀ =   4.770370                                                                  --- execution time = 12s ---
1        R/S            >>>>    a₁ = -5.802627   X<>Y   b₁ = 3.282884                          --- execution time = 23s ---
2        R/S            >>>>    a₂ = 1.060168   X<>Y   b₂ = 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+" >>>>    a₀ =   6.083282                                                                  --- execution time = 45s ---
      1        R/S            >>>>    a₁ = -2.250492   X<>Y   b₁ = -2.782105                         --- execution time = 64s --- ( 147s with "FIL4" )
      2        R/S            >>>>    a₂ = -0.940576   X<>Y   b₂ = -0.907591                                       ..........
        .....................................................................................................                                          ..........

-In this example, the coefficients for a₁ & b₁ are computed via the Taylor series ( 2.PI/48 < 1/6 )
whereas those for a₂ & b₂ 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:

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

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

A = y'(a)/µ² - y'''(a)/µ⁴ + y'''''(a)/µ⁶ - ...... C = y(a)/µ - y''(a)/µ³ + y'''''(a)/µ⁵ - y''''''(a)/µ⁷ + ..........
with
B = y'(b)/µ² - y'''(b)/µ⁴ + y'''''(b)/µ⁶ - ...... D = y(b)/µ - y''(b)/µ³ + y'''''(b)/µ⁵ - y''''''(b)/µ⁴ + ..........

-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) ] / h²

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) ] / h³

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) ] / h⁴

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) ] / h⁵

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) ] / h⁶ = 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 a₀ = 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	/	b_k
X	k	a_k

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" >>>>   a₀ =   8.475595                                                                  --- execution time = 6s ---
1        R/S         >>>>    a₁ = -1.059182   X<>Y   b₁ = -1.090833                         --- execution time = 46s ---
2        R/S         >>>>    a₂ = 1.776544   X<>Y   b₂ = -4.029470                                        ........
   .....................................................................................................

6        R/S        >>>>    a₆ = 0.171103   X<>Y   b₆ = -0.979868
7        R/S         >>>>    a₇ = 0.048603   X<>Y   b₇ = -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:

§₀ⁿ y(x) cos µ.x dx = Sum_{i=0,6,12,.....,n-6} y(i) A₀ cos µ.i + y(i+1) [ A₁ cos (µ.i+µ) + B₁ sin (µ.i+µ) ] + y(i+2) [ A₂ cos (µ.i+2µ) + B₂ sin (µ.i+2µ) ]

+ y(i+3) A₃ cos (µ.i+3µ) + y(i+4) [ A₂ cos (µ.i+4µ) - B₂ sin (µ.i+4µ) ] + y(i+5) [ A₁ cos (µ.i+5µ) - B₁ sin (µ.i+5µ) ]

§₀ⁿ y(x) sin µ.x dx = Sum_{i=0,6,12,.....,n-6} y(i) A₀ sin µ.i + y(i+1) [ A₁ sin (µ.i+µ) - B₁ cos (µ.i+µ) ] + y(i+2) [ A₂ sin (µ.i+2µ) - B₂ cos (µ.i+2µ) ]

+ y(i+3) A₃ sin (µ.i+3µ) + y(i+4) [ A₂ sin (µ.i+4µ) + B₂ cos (µ.i+4µ) ] + y(i+5) [ A₁ sin (µ.i+5µ) + B₁ cos (µ.i+5µ) ]
with

      A₀ = (1/3µ²-15/4µ⁴+5/µ⁶) cos 6µ + (49/10µ²-49/4µ⁴+7/µ⁶) + (-137/90µ³+17/3µ⁵-2/µ⁷) sin 6µ
      A₁ = (-6/5µ²+13/µ⁴-16/µ⁶) cos 5µ - (6/µ²-29/µ⁴+20/µ⁶) cos µ + (27/5µ³-19/µ⁵+6/µ⁷) sin 5µ + (87/5µ³-31/µ⁵+6/µ⁷) sin µ
      A₂ = (15/4µ²-307/8µ⁴+85/2µ⁶) cos 4µ - (-15/2µ²+461/8µ⁴-95/2µ⁶) cos 2µ + (-33/2µ³+107/2µ⁵-15/µ⁷) sin 4µ + (-117/4µ³+137/2µ⁵-15/µ⁷) sin 2µ
      A₃ = 2(-20/3µ²+62/µ⁴-60/µ⁶) cos 3µ + 2(254/9µ³-242/3µ⁵+20/µ⁷) sin 3µ

and

B₁ = -(-6/5µ²+13/µ⁴-16/µ⁶) sin 5µ - (6/µ²-29/µ⁴+20/µ⁶) sin µ + (27/5µ³-19/µ⁵+6/µ⁷) cos 5µ - (87/5µ³-31/µ⁵+6/µ⁷) cos µ
B₂ = -(15/4µ²-307/8µ⁴+85/2µ⁶) sin 4µ - (-15/2µ²+461/8µ⁴-95/2µ⁶) sin 2µ + (-33/2µ³+107/2µ⁵-15/µ⁷) cos 4µ - (-117/4µ³+137/2µ⁵-15/µ⁷) cos 2µ

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

A₀ = 41/70 + 9µ²/25 - 567µ⁴/550 + 28053µ⁶/25025 - 16281µ⁸/25025 + 498636µ¹⁰/2127125 - 0.057560072125µ¹² + 0.010259175386µ¹⁴ - .............

A₁ = 54/35 - 27µ²/25 + 1917µ⁴/1100 - 256947µ⁶/200200 + 809283µ⁸/1601600 - 33654633µ¹⁰/272272000 + 0.020638591235µ¹²
- 0.0025064271605µ¹⁴ + .............

A₂ = 27/140 + 27µ²/10 - 513µ⁴/220 + 9903µ⁶/10010 - 11859µ⁸/50050 + 76134µ¹⁰/2127125 - 0.0037150884891µ¹²
+ 0.000281765004581µ¹⁴ - .............

A₃ = 68/35 - 18µ²/5 + 243µ⁴/110 - 2133µ⁶/4004 + 55161µ⁸/800800 - 763263µ¹⁰/136136000 + 0.000315589112894µ¹²
- 0.0000130645273645µ¹⁴ + .............

B₁ = 36µ³/25 - 18µ⁵/11 + 1508589µ⁷/1751750 - 5137µ⁹/19500 + 0.052739990401µ¹¹ - 0.0074605626534µ¹³ + .............

B₂ = -9µ³/5 + 414µ⁵/275 - 88758µ⁷/175175 + 290µ⁹/3003 - 0.0120217037863µ¹¹ + 0.00106051813438µ¹³ + .............

-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 A_i and B_j 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 X<Y?
35 GTO 00
36 RCL 09
37 RCL 09
38 RCL 09
39 17
40 /
41 CHS
42 .2344
43 +
44 *
45 .6506
46 -
47 *
48 1.121
49 +
50 *
51 567
52 550
53 /
54 -
55 *
56 .36
57 +
58 *
59 RCL 01
60 +
61 ST+ 01
62 CLX
63 .1236
64 *
65 CHS
66 .5053
67 +
68 *
69 1.28345
70 -
71 *
72 19.17
73 11
74 /
75 +
76 *
77 1.08
78 -
79 *
80 RCL 02
81 +
82 ST+ 02
83 CLX
84 27.94
85 /

86 .2369
87 -
88 *
89 .9893
90 +
91 *
92 513
93 220
94 /
95 -
96 *
97 2.7
98 +
99 *
100 RCL 03
101 +
102 ST+ 03
103 CLX
104 14.52
105 /
106 .5327
107 -
108 *
109 243
110 110
111 /
112 +
113 *
114 3.6
115 -
116 *
117 RCL 04
118 +
119 ST+ 04
120 CLX
121 19
122 /
123 .2634
124 -
125 *
126 .8612
127 +
128 *
129 18
130 11
131 /
132 -
133 *
134 1.44
135 +
136 *
137 RCL 08
138 *
139 STO 05
140 CLX
141 10.36
142 /
143 .5067
144 -
145 *
146 414
147 275
148 /
149 +
150 *
151 1.8
152 -
153 *
154 RCL 08
155 *
156 STO 06
157 GTO 01
158 LBL 00
159 RCL 09
160 1/X
161 1.75
162 -
163 RCL 09
164 /
165 .7
166 +
167 7
168 *
169 5
170 RCL 09

171 /
172 3.75
173 -
174 RCL 09
175 /
176 3
177 1/X
178 +
179 RCL 07
180 6
181 *
182 STO 01
183 COS
184 *
185 +
186 2
187 RCL 09
188 /
189 17
190 3
191 /
192 -
193 RCL 09
194 /
195 137
196 90
197 /
198 +
199 RCL 08
200 /
201 RCL 01
202 SIN
203 *
204 -
205 STO 01
206 RCL 07
207 5
208 *
209 STO 06
210 16
211 RCL 09
212 /
213 13
214 -
215 RCL 09
216 /
217 1.2
218 +
219 P-R
220 CHS
221 STO 02
222 X<>Y
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	/	b_k
X	k	a_k

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:

-Then n = 48 STO 00

      0 XEQ "FIL6+" >>>>    a₀ =   6.083399                                                                  --- execution time = 48s ---
      1        R/S            >>>>    a₁ = -2.250254   X<>Y   b₁ = -2.782055                         --- execution time = 71s ---
      2        R/S            >>>>    a₂ = -0.940317   X<>Y   b₂ = -0.907481                                       ..........
        .....................................................................................................                                          ..........

7 R/S >>>> a₇ = -0.158130 X<>Y b₇ = -0.086882 ..........

Notes:

-If we use "FIL6" instead of "FIL6+", a₀ is correctly computed but "FIL6" gives a₁ = -2.243751 & b₁ = -2.784534 , a₂ = -0.940469 & b₂ = -0.907335
-However, a₇ & b₇ 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)
a₀	6.074830	6.083005	6.083282	6.083399	6.083605
a₁	-2.267371	-2.251065	-2.250492	-2.250254	-2.249808
b₁	-2.781891	-2.782237	-2.782105	-2.782055	-2.781991
a₂	-0.957391	-0.941185	-0.940576	-0.940317	-0.939786
b₂	-0.907200	-0.907842	-0.907590	-0.907481	-0.907391
a₇	-0.175917	-0.160233	-0.158926	-0.158130	-0.157669
b₇	-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 [ a₁ , .......... , a_n ] where a_k = b_k + i c_k ,
"DFTZ" computes the complex components z_k = x_k + i y_k of the Discrete Fourier Transform [ z₁ , .......... , z_n ]

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

• R01 = b₁ • R03 = b₂ ................ • R2n-1 = b_n
• R02 = c₁ • R04 = c₂ ................ • R2n = c_n
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	/	y_k
X	k	x_k
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" >>>>   x₁ = 16   X<>Y   y₁ = 20
      2       R/S            >>>>   x₂ = -8    X<>Y   y₂ =   0
      3       R/S            >>>>   x₃ = -4    X<>Y   y₃ = -4
      4       R/S            >>>>   x₄ = 0    X<>Y   y₄ = -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	/	y_k
X	k	x_k
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"   >>>>   x₁ = 1   X<>Y   y₁ = 2
      2       R/S              >>>>   x₂ = 3    X<>Y   y₂ = 4
      3       R/S              >>>>   x₃ = 5    X<>Y   y₃ = 6
      4       R/S              >>>>   x₄ = 7    X<>Y   y₄ = 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 [ z₁ , .......... , z_n ] z_k = x_k + i y_k
of an array of n complex numbers [ a₁ , .......... , a_n ] a_k = b_k + i c_k

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

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

Inputs: • R01 = b₁ • R03 = b₂ ................ • R2n-1 = b_n
• R02 = c₁ • R04 = c₂ ................ • R2n = c_n R2n+1 to R4n: temp

Outputs: • R01 = x₁ • R03 = x₂ ................ • R2n-1 = x_n
• R02 = y₁ • R04 = y₂ ................ • R2n = y_n
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, n_max = 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 [ a_i,j ] ( 1 <= i <= n , 1 <= j <= m ). The result is a complex nxm matrix [ z_i,j ]

where z_i,j = SUM_{q=1,2,...,n
; r=1,2,....,m} a_q,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 = a_1,1      • Rn+1 = a_1,2   .......................................   • Rnm-n+1 = a_1,m
                                • R02 = a_2,1      • Rn+2 = a_2,2   .......................................   • Rnm-n+2 = a_2,m
                                       ..........................................................................................................
                                • Rnn = a_n,1      • R2n = a_n,2   ........................................   • Rnm = a_n,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	y_i,j
X	i	x_i,j

with z_i,j = x_i,j + i y_i,j

Example:

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

-Calculate z_2,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 z_2,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	y_i,j
X	i	x_i,j

with z_i,j = x_i,j + i y_i,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) = SUM_i;j [ a_i;j cos ( 360 ( ix/n+jy/m ) ) + b_i;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	b_i,j
X	i	a_i,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"    >>>   a_0,0 =   6.6667
      1 ENTER^       0    R/S           >>>   a_0;1 =      3         X<>Y   b_0;1 = -2.3333
      2 ENTER^       0    R/S           >>>   a_0;2 =      2         X<>Y   b_0;2 =      0
      0 ENTER^       1    R/S           >>>   a_1;0 =   0.3333   X<>Y    b_1;0 = -1.4434
      1 ENTER^             R/S           >>>   a_1;1 = -0.7113   X<>Y    b_1;1 = -0.1220
      2 ENTER^       1    R/S           >>>   a_1;2 =      1          X<>Y   b_1;2 = -0.2887
     -1 ENTER^      1    R/S           >>>   a_1;-1 = -1.2887    X<>Y   b_1;-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 z_i,j = x_i,j + i y_i,j of DFT [ a_i,j ] with a_i,j = b_i,j + i c_i,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 = b_1,1 • R2n+1 = b_1,2 ....................................... • R2nm-2n+1 = b_1,m
• R02 = c_1,1 • R2n+2 = c_1,2 ....................................... • R2nm-2n+2 = c_1,m

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

• R2n-1 = b_n,1 • R4n-1 = b_n,2 ....................................... • R2nm-1 = b_n,m
• R2n = c_n,1 • R4n = c_n,2 ....................................... • R2nm = c_n,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	y_i,j
X	i	x_i,j

with z_i,j = x_i,j + i y_i,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 z_2,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 z_2,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 [ a_i,j,k ] ( 1 <= i <= n , 1 <= j <= m , 1 <= k <= p ).
-The result is a complex hypermatrix [ z_i,j,k ] with z_i,j,k = x_i,j,k + i y_i,j,k

z_i,j,k = SUM_{q=1,2,...,n
; r=1,2,....,m ; s=1,2,.....,p} a_q,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 = a_1,1,1      • Rn+1 = a_1,2,1   .......................................   • Rnm-n+1 = a_1,m,1
                                • R02 = a_2,1,1      • Rn+2 = a_2,2,1   .......................................   • Rnm-n+2 = a_2,m,1
                                       ..................................................................................................................
                                • Rn = a_n,1,1        • R2n = a_n,2,1   ........................................ • Rnm = a_n,m,1

                                          • Rnm+1 = a_1,1,2      • Rnm+n+1 = a_1,2,2   .......................................   • R2nm-n+1 = a_1,m,2
                                          • Rnm+2 = a_2,1,2      • Rnm+n+2 = a_2,2,2   .......................................   • R2nm-n+2 = a_2,m,2
                                              ..................................................................................................................................
                                          • Rnm+n = a_n,1,2      • Rnm+2n = a_n,2,2   ........................................ • R2nm = a_n,m,2

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

                                                    • Rnm(p-1)+1 = a_1,1,p      • Rnm(p-1)+n+1 = a_1,2,p   .......................................   • Rnmp-n+1 = a_1,m,p
                                                    • Rnm(p-1)+2 = a_2,1,p      • Rnm(p-1)+n+2 = a_2,2,p   .......................................   • Rnmp-n+2 = a_2,m,p
                                                       ....................................................................................................................................................
                                                    • Rnm(p-1)+n = a_n,1,p      • Rnm(p-1)+2n = a_n,2,p   ........................................ • Rnmp = a_n,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	y_i,j,k
X	i	x_i,j,k

with z_i,j,k = x_i,j,k + i y_i,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 z_2,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, z_2,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 [ a_i,j,k ] with a_i,j,k = b_i,j,k + i c_i,j,k ( 1 <= i <= n , 1 <= j <= m , 1 <= k <= p ).
-The result is also a complex hypermatrix [ z_i,j,k ] with z_i,j,k = x_i,j,k + i y_i,j,k

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

                                • R01 = b_1,1,1      ...................................................................   • R2nm-2n+1 = b_1,m,1
                                • R02 = c_1,1,1      ...................................................................   • R2nm-2n+2 = c_1,m,1
                                       ..........................................................................................................
                                       ..........................................................................................................

• R2n-1 = b_n,1 .................................................................... • R2nm-1 = b_n,m
• R2n = c_n,1 .................................................................. • R2nm = c_n,m

                                                         • R2nm(p-1)+1 = b_1,1,p    .....................................................................   • R2nmp-2n+1 = b_1,m,p
                                                         • R2nm(p-1)+2 = c_1,1,p    .....................................................................   • R2nmp-2n+2 = c_1,m,p
                                                                          ..........................................................................................................
                                                                          ..........................................................................................................

• R2nm(p-1)+2n-1 = b_n,1,p .................................................................... • R2nmp-1 = b_n,m,p
• R2nm(p-1)+2n = c_n,1,p .................................................................... • R2nmp = c_n,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	y_i,j,k
X	i	x_i,j,k

with z_i,j,k = x_i,j,k + i y_i,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 z_2,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, z_2,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

hp41programs

Spectral Analysis for the HP-41