The Museum of HP Calculators
Spectral Analysis for the HP-41
This program is Copyright © 2004 by
Jean-Marc Baillard
and is used here by permission.
This program is supplied without representation or warranty of any kind.
Jean-Marc Baillard and The Museum of HP Calculators therefore assume no
responsibility and shall have no liability, consequential or otherwise, of
any kind arising from the use of this program material or any part thereof.
Overview
-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 = the period of y : 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.
-"FIL" calculates these coefficients by Filon's formula.
-"SPA2" generalizes "SPA" to 2-dimensional problems
-All these programs use data points given at equally spaced arguments.
1°) One-dimensional Problem
a) The Discrete Fourier Transform
Data Registers: Registers R00 thru Rnn are to
be initialized before executing "SPA":
• R00 = n •
R01 = y(1) • R02 = y(2) ................
• Rnn = y(n) = y(0)
Flags: /
Subroutine: /
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:
-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" yields 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)
b) Filon's formula (
assuming n is even )
-The Discrete Fourier Transform (DFT) 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 f continuous
over [ 0;n ]
( not only over [ 0;n [ )
Data Registers: Registers R00 thru Rnn are to
be initialized before executing "FIL":
• R00 = n •
R01 = y(1) • R02 = y(2) ................
• Rnn = y(n) = y(0)
Flag: F10
Subroutine: /
N.B: 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 )
Then 0 XEQ "FIL" yields
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) - ......................
Notes: -"FLT" gives the exact values when y is
a polynomial of degree < 4
-Unlike "SPA" , "FIL"
does not collocate the samples.
2°) Two-dimensional Problem
-We assume z(x,y) is a doubly periodic function:
z(x+n;y) = z(x;y) for all x;y
and 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: Registers R00 thru Rnm
are to be initialized before executing "SPA2":
• R00 = n.mmm (
I mean n + m/1000 )
• 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
Subroutine: /
N.B: 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 ( synthetic )
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
and 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
1 3 7 10
R01 R04 R07 R10
3.004 STO 00 and we store the 12
numbers: 4 5 8
14 into 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)
-Thus 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 )
References: -Abramowitz & Stegun
"Handbook of Mathematical Functions" Dover Publications
ISBN 0-486-61272-4
-Ronald N. Bracewell "The Fourier Transform and its Applications"
McGraw-Hill ISBN 0-07-116043-4
-F.Scheid - "Numerical Analysis" - Mc Graw Hill - ISBN 07-055197-9
Go back to the HP-41 software library
Go back to the general software library
Go
back to the main exhibit hall