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:

 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"  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