SPECTRAL Manual

 Overview
 

-This module put together several programs to solve spectral analysis problems:
-The Discrete Fourier Transform ( DFT ) & its inverse, in dimension 1, 2, 3.
-The Fast Fourier Transform ( FFT ) is also included ( dimension 1 only )

-The DFT is essentially equivalent to using the trapezoidal rule to compute the integrals that appear in the Fourier coefficients of a periodic function.
-The Filon's integration formulas are employed in "FILON" for a better accuracy.
-Other Filon-type formulae ( order 4 & 6 ) give an even better precision in "FIL4" & "FIL6"

-A direct method to evaluate these integrals is also available with FFOUR and FOURN

-Unlike the DFT, the Discrete Hartley Transform ( DHT ) has the advantage of returning real numbers only.
-Strictly symmetric DHTs are computed in dimensions 1, 2, 3.
-Thus, these DHTs are their own inverses.
-The Fast Hartley Transform ( FHT ) is also included ( dimension 1 only, not strictly symmetric )

-Finally, the Unilateral Z-Transform may be computed too.
 
 
 

XROM	Function	Desciption
08,00  08,01  08,02  08,03  08,04  08,05  08,06  08,07  08,08  08,09  08,10  08,11  08,12  08,13  08,14  08,15  08,16  08,17  08,18  08,19  08,20  08,21  08,22  08,23  08,24  08,25  08,26  08,27  08,28  08,29  08,30  08,31  08,32  08,33  08,34  08,35  08,36  08,37  08,38  08,39	-SPECTRAL C   AINT   E3/   E3/E+   EVEN?  "IN"   FFOUR   FOURN  "#FN"   ODD?  "OUT"  -FOURIER TF  "DFT1"  "DFT1Z"  "DFT2"  "DFT2Z"  "DFT3"  "DFT3Z"  "IDFT1"  "IDFT1Z"  "IDFT2"  "IDFT2Z"  "IDFT3"  "IDFT3Z"  "FFT"  "IFFT"  "FIL4"  "FIL6"  "FILON"  "SPA"  -HARTLEY TF   DHT  "DHT1"  "DHT2"  "DHT3"  "FHT"  "HRTL"  "ZTF"  "ZTFZ"  "360"	Section Header  Stores an integer into alpha  Divides X by 1000  Divides X by 1000 and add 1  Tests if X is even  Storing numbers in a block of registers  Caculates the Fourier Coefficients of a Periodic Function  Gives the value of a specific pair of coefficients an, bn  A subroutine called by FOURN  Tests if X is odd  Displays the elements of a block of registers  Section Header  Discrete Fourier Transform dim 1  Discrete Fourier Transform dim 1 ( complex )  Discrete Fourier Transform dim 2  Discreter Fourier Transform dim 2 ( complex )  Discrete Fourier Transform dim 3  Discrete Fourier Transform dim 3 ( complex )  Inverse Discrete Fourier Transform dim 1  Inverse Discrete Fourier Transform dim 1 ( complex )  Inverse Discrete Fourier Transform dim 2  Inverse Discrete Fourier Transform dim 2 ( complex )  Inverse Discrete Fourier Transform dim 3  Inverse Discrete Fourier Transform dim 3 ( complex )  Fast Fourier Transform  Inverse Fast Fourier Transform  Filon-type Integration ( order 4 )  Filon-type Integration ( order 6 )  Filon Integration  Spectral Analysis  Section Header  Discrete Hartley Transform Dim 1&2 ( M-Code )  Discrete Hartley Transform Dim 1  Discrete Hartley Transform Dim 2  Discrete Hartley Transform Dim 3  Fast Hartley Transform  Continuous Hartley Transform  Z-Transform  Complex Z-Transform  A Subroutine

 

Storing an Integer into Alpha
 

-The M-code function AINT actually appends the integer part of X to the alpha "register"
 

Dividing X by 1000
 

-The M-code function E3/ simply divides X by 1000 without disturbing registers P & Q
 
 

Dividing X by 1000 & Adding 1
 

  E3/E+  add one to the result of  E3/
 

Even ?
 

-The M-code function EVEN? tests if X-register contains an even number. The answer is YES or NO
-In a program, the following line is skipped if the answer is NO.
 
 

Odd ?
 

-The M-code function ODD? tests if X-register contains an odd number. The answer is YES or NO
-In a program, the following line is skipped if the answer is NO.
 
 

Storing Numbers
 

-"IN" stores real numbers in a block of registers.
-The input is the address of the first element.
-After each Prompt, place the required number in X and R/S.

-When all the elements are stored, respond R/S ( without any digit entry ) to the last PROMPT: the control number of the array is returned in X

--->  If you want to store a number like PI, set flag F22 before R/S ( otherwise, the routine would stop prematurely )
 

Displaying the elements of a Block of Registers
 

-"OUT"  displays successively the elements ( numbers or alpha strings ) in a block of registers.
-It requires an input of the form  bbb.eeeii where Rbb is the 1st register, Ree is the last one and ii = the increment ( or 0 if ii = 1 )

-Set flag F21 if you want that the program stops at each AVIEW
 

Fourier Coefficients of a Periodic Function
 

-The function FFOUR calculates the Fourier coefficients for a periodic function F(x), defined as:

    a_n = (1/L) §₀^2L  f(x')  cos (n.PI.x'/L) dx'
    b_n = (1/L) §₀^2L  f(x')  sin (n.PI.x'/L) dx'

  with the following characteristics:

-centered in x = x₀
-with period 2L on an interval [x₀, x₀+2L]
-with a given precision for calculations (significant decimal places)

The function must be programmed in main memory under its own global label.
The program prompts for the first index to calculate, and the number of desired coefficients.

The program also calculates the approximate value of the function at a given argument applying the summation of the terms, using the obtained coefficients:

   f(x') = a₀/2 + Sum_{n=1 to infinity}  a_n cos (n.PI.x'/L) + Sum_{n=1 to infinity}  b_n sin (n.PI.x'/L)

To use it simply enter the value of x and press “A” (XEQ A) in user mode on – this assumes that no function is assigned to the key.
The approximation will be more correct when a sufficient number of terms is included.
 

Example:  Calculate the first six coefficients for F(x) = x^2, assuming a period T=2pi, centered in x0 = 0.  As it’s known,

 X^2 = 4/3 Pi^2 + SUM{ 4 cos(nx) /n^2 - 4 Pi sin(nx) /n } | n = 0,1,…

Thus using an accuracy of 6 decimal places:
 
 

a0 = 13,1595	b0 = 0
a1 = 4	b1 = -12,5664
a2 = 1	b2 = -6,2832
a3 = 0,4444	b3 = -4,1888
a4 = 0,2500	b4 = -3,1416
a5= 0,1600	b5 = -2.5133

 

Calculating the value of a specific pair of coefficients an, bn
 

-The function FOURN  uses the Advantage Pac INTEG function
  to calculate by the direct method the value of a specific pair of coefficients a_n, b_n.
-It uses the subroutine #FN for the integrand.
 

Discrete Fourier Transform dim 1
 

-Given an array of n real numbers  [ a₁ , .......... , a_n ] ,  "DFT1" 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 "DFT1" )

                                    •   R01 = a₁   •   R02 = a₂  ................  •   Rnn = a_n
 
 

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

Inverse Discrete Fourier Transform dim 1
 

-It's almost the same program ( F01 is set ).
 

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

                                    •   R01 = a₁   •   R02 = a₂  ................  •   Rnn = a_n
 
 

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 "IDFT1"  >>>>   x₁ =  3    X<>Y   y₁ =   0
      2  XEQ "IDFT1"  >>>>   x₂ =  1    X<>Y   y₂ =  0.5
      3  XEQ "IDFT1"  >>>>   x₃ =  1    X<>Y   y₃ =   0
      4  XEQ "IDFT1"  >>>>   x₄ =  1    X<>Y   y₄ = -0.5

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

Discrete Fourier Transform dim 1 ( complex )
 

-Given an array of n complex numbers  [ a₁ , .......... , a_n ]  where  a_k = b_k + i c_k ,
 "DFT1Z" 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 "DFT1Z" )

                                    •   R01 = b₁   •   R03 = b₂  ................  •   R2n-1 = b_n
                                    •   R02 = c₁   •   R04 = c₂  ................   •   R2n = c_n
 
 

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

Inverse Discrete Fourier Transform dim 1 ( complex )
 

 -Almost the same calculations, F01 is set.
 

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

                                    •   R01 = b₁   •   R03 = b₂  ................  •   R2n-1 = b_n
                                    •   R02 = c₁   •   R04 = c₂  ................   •   R2n = c_n
 
 

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 "IDFT1Z"   >>>>   x₁ = 1   X<>Y   y₁ = 2
      2  XEQ "IDFT1Z"   >>>>   x₂ = 3    X<>Y   y₂ =  4
      3  XEQ "IDFT1Z"   >>>>   x₃ = 5    X<>Y   y₃ =  6
      4  XEQ "IDFT1Z"  >>>>   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 ]
 

Fast Fourier Transform
 

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

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 we get 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 ]

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
 

Inverse Fast Fourier Transform
 

-Almost the same calculations, F01 is set.
 

Discrete Fourier Transform dim 2
 

-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
 
 

STACK	INPUTS	OUTPUTS
Y	j	yi,j
X	i	xi,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
                             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   ] ]
 

Inverse Discrete Fourier Transform dim 2
 

-Almost the same calculations, F01 is set.
 

Discrete Fourier Transform dim 2 ( complex )
 

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

                                •  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
 
 

STACK	INPUTS	OUTPUTS
Y	j	yi,j
X	i	xi,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 "DFT2Z"  >>>>>    0.696152398
                               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  ] ]
 

Notes:

-Synthetic register a is used.
-So do not call this program as more than a first level subroutine.
-Same remark for "IDFT2Z" below.
 

Inverse Discrete Fourier Transform dim 2 ( complex )
 

-Almost the same calculations, F01 is set.
 

Discrete Fourier Transform dim 3
 

-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
 
 

STACK	INPUTS	OUTPUTS
Z	k	/
Y	j	yi,j,k
X	i	xi,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
                               X<>Y    10.96762920

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

Notes:

-Synthetic register a is used.
-So do not call this program as more than a first level subroutine.
-Same remark for "IDFT3" below.
 

Inverse Discrete Fourier Transform dim 3
 

-Almost the same calculations, F01 is set.
 

Discrete Fourier Transform dim 3 ( complex )
 

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

                                •  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
 
 

STACK	INPUTS	OUTPUTS
Z	k	/
Y	j	yi,j,k
X	i	xi,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 "DFT3Z"  >>>>    20.47040998
                                X<>Y   -159.3203501

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

Notes:

-Synthetic register a is used.
-So do not call this program as more than a first level subroutine.
-Same remark for "IDFT3Z" below.
 

Inverse Discrete Fourier Transform dim 3 ( complex )
 

-Almost the same calculations, F01 is set.
 

Spectral Analysis
 

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

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)
 

Filon Integration
 

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

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

STACK	INPUTS	OUTPUTS
Y	/	bk
X	k	ak

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 )

       0  XEQ "FILON"  >>>     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:

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

Filon-type Integration ( order 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)
 

-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 this 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.
 

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:  /
 
 

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"  >>>>    a₀ =   4.770370
  1        R/S         >>>>    a₁ = -5.802627   X<>Y   b₁ =  3.282884
  2        R/S         >>>>    a₂ =  1.060168   X<>Y   b₂ =  0.111899
   .....................................................................................................
 

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
      1        R/S         >>>>    a₁ = -2.250492   X<>Y   b₁ =  -2.782105
      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.
 

Filon-type Integration ( order 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

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

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:  /
 
 

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+"  >>>>    a₀ =   6.083399
      1        R/S            >>>>    a₁ = -2.250254   X<>Y   b₁ =  -2.782055
      2        R/S            >>>>    a₂ = -0.940317   X<>Y   b₂ =  -0.907481
        .....................................................................................................

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

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

 

Discrete Hartley Transform Dim 1
 

-The DHT of a vector  [ a₁ , a₂ , ........ , a_n ]  is a vector   [ b₁ , b₂ , ........ , b_n ]

  where   b_i = (1/sqrt(n))  SUM_q=1,2,...,n   a_q. cas  360°(q-1)(i-1)/n    with   cas µ = sin µ + cos µ
 

Data Registers:                                         ( Registers R00 thru Rnn are to be initialized before executing "DHT1" )

                                •  R00 = n    •  R01 = a₁ ;  •  R02 = a₂ ; .......... ;  •  Rnn = a_n
Flags:  /
Subroutines: /
 

STACK	INPUTS	OUTPUTS
X	i	bi
L	/	i

    ( 1 <= i <= n )

Example:    Find the DHT of  [ 2  4  7  6 ]

  n = 4   therefore  4  STO 00   and   2  STO 01   4 STO 02   7  STO 03   6  STO 04

  1  XEQ "DHT1"  >>>    9.5
  2      R/S             >>>   -3.5
  3      R/S             >>>   -0.5
  4      R/S             >>>   -1.5     whence   DHT [ 2  4  7  6 ]  =  [ 9.5  -3.5  -0.5  -1.5 ]

Note:

-We have the property   DHT_oDHT = Id
 

Fast Hartley Transform
 

-If  n = 2^N is an integer power of 2 , the DHT of  [ a₁ , a₂ , ........ , a_n ]  =  [ b₁ , b₂ , ........ , b_n ]  may be perform more quickly by the following method:

 a)  First, the  a_i  are placed into the jth position where j is calculated as follows ( lines 01 to 37 ):         i-1 is written in binary ,
      and the digits are reversed which yields  j-1   ( for instance, if  i = 4 , 4-1 = 3 = (011)₂ >>> (110)₂ = 6 = 7-1  thus, a₄ is placed into the 7th position )

 b)  Then, for L = 1 to N,   registers  Rii  are replaced by  Rj₁ + Rj₂ .cos 360°(i-1)/2^L+ Rj₃ .sin 360°(i-1)/2^L    ( i = 1 ,2 , ...... , n )   with

        j₁ = i - (i-1) mod 2^L + (i-1) mod 2^L-1  ;
        j₂ = j₁ + 2^L-1  ;
        j₃ =  j₂           if   ( i-1) mod 2^L-1 = 0  and   j₃ =  2^L + i - (i-1) mod 2^L - (i-1) mod 2^L-1  otherwise.
 

Data Registers:                                            ( Registers R00 thru Rnn are to be initialized before executing "FHT" )

                              •  R00 = n    •  R01 = a₁   •  R02 = a₂  ..........   •  Rnn = a_n
    and when the program stops,      R01 = b₁      R02 = b₂  ..........      Rnn = b_n

                                                       ( Rn+1 thru R2n are used for temporary data storage )
Flags: /
Subroutines:  /
 

Example:    Find   DHT [ 2  4  7  6 ]

n = 4   therefore  4  STO 00   and   2  STO 01   4 STO 02   7  STO 03   6  STO 04   XEQ "FHT"

 and we get  R01 = 19   R02 = -7   R03 = -1   R04 = -3  whence   DHT [ 2  4  7  6 ]  =  [ 19  -7  -1  -3 ]

Notes:

-This program doesn't work if  n = 1    ( but DHT [ a ] = [ a ] )
-Here, the DHT is not quite symmetric: divide all the coefficients by sqrt(n) if you want a strictly symmetric FHT.
-I've measured the following execution times:
 

n	2	4	8	16	32	64	128
nxDHT	4s	14s	53s	3m32s	14m10s	56m40s	3h47m
FHT	6s	21s	60s	2m38s	6m36s	15m59s	37m34s

-FHT execution time is approximately proportional to  n.Log n  ( instead of n² for n.DHT )
-Obviously, this "fast" algorithm remains relatively slow with an HP-41 and it's worthwhile only if  n > 8
-However, the advantage is substantial if  n = 128  ( about 6 times as fast )
-Furthermore, calculations are less complicated and roundoff errors are smaller.
 

Discrete Hartley Transform Dim 2
 

-The DHT may be generalized to  nxm matrices  [ a_i,j ]  ( 1 <= i <= n ; 1 <= j <= m ). The result is also a  nxm matrix  [ b_i,j ]

  where  b_i,j =  (1/sqrt(n.m)) SUM_{q=1,2,...,n  ; r=1,2,....,m}   a_q,r. cas 360°[(q-1)(i-1)/n+(r-1)(j-1)/m]            ( cas = sin + cos )
 

Data Registers:     •  R00 = n.mmm  =  n + m/1000     ( these  nm+1 registers are to be initialized before executing "DHT2" )

                                •  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
Flag: /
Subroutine:  /
 

STACK	INPUTS	OUTPUTS
Y	j	/
X	i	bi,j

   ( 1 <= i <= n ; 1 <= j <= m )

 Example:

               [ [  1  3  4  10  ]
 Let  A =   [  4  5  7  14  ]       Calculate   b_2,3
                 [  2  9  6  11  ] ]

 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 "DHT2"  >>>  b_2,3 =  1.577350267

-Likewise, we find

                       [ [  21.9393  -8.0829   -8.0829   -2.3094 ]
     DHT (A) =   [  -2.6754    1.7113    1.5774   -0.9226  ]        ( rounded to 4 decimals )
                         [  -3.6754  -2.2887   -0.4226     2.0774  ] ]
 

Discrete Hartley Transform Dim 1&2 ( M-Code )
 

 >>>   DHT works for vectors and matrices.
 

Data Registers:     •  R00 = n.mmm  =  n + m/1000    or  n  if  m = 1        ( These  nm+1 registers are to be initialized before executing "DHT" )

                                •  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
Flag: /
Subroutine:  /
 
 

STACK	INPUT1	OUTPUT1	INPUT2	OUTPUT2
Y	j	/	/	/
X	i	bi,j	0	b1,1

  If  m = 1 ,  j  may take any real value.

  If  X = 0 , the whole DHT is computed and sored into register  Rmn+1 thru R2mn

Example1:           A = [ 1 2 4 7 ]     here, n = 4 & m = 1

   Store  1  2  4  7  into  R01  R02  R03  R04   and   n = 4   STO 00

   •  0  XEQ "DHT"  returns  b₁ = 7  and registers  [ R05  R06  R07  R08 ] = [ 7  -4  -2   1 ]

   •  3  XEQ "DHT"  gives  b₃ = -2

Example2:        Let  M  the 2x3 matrix defined by

    1  2  4      =      R01   R03   R05          Store  n.mmm = 2.003  STO 00
    3  5  6      =      R02   R04   R06

   •  0  XEQ "DHT"  returns  b₁ = 8.573214097

   R07   R09   R11    =    8.5732   -2.8978   -0.7765       ( rounded to 4D )
   R08   R10   R12    =   -2.8577  -0.1494    0.5577

 If you copy R07 ....  R12 to R01 .... R06  and press  0  XEQ "DHT"  again,
 you'll get the elements of the original matrix M with a mean error of about 3 E-9

   •  3  ENTER^
       2  XEQ "DHT"  gives  b_2,3 = 0.557677535
 

Discrete Hartley Transform Dim 3
 

-In this case, the DHT of a  nxmxp  hypermatrix [ a_i,j,k ]  ( 1 <= i <= n ; 1 <= j <= m ; 1 <= k <= p ) is a  nxmxp hypermatrix  [ b_i,j,k ]

  where  b_i,j,k = (1/sqrt(n.m.p)) SUM_{q=1,2,...,n  ; r=1,2,....,m ; s=1,2,....,p}   a_q,r,s. cas 360°[(q-1)(i-1)/n+(r-1)(j-1)/m+(s-1)(k-1)/p]         ( cas = sin + cos )

Data Registers:     •  R00 = n.mmmppp  =  n + m/1,000 + p/1,000,000  ( these  nmp+1 registers are to be initialized before executing "DHT3" )

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

STACK	INPUTS	OUTPUTS
Z	k	/
Y	j	/
X	i	bi,j,k

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

Example:

    Let A =  [ [ [   1  25  40   5    2  ]                        To store these 60 coefficients into the proper registers, you can key in    60   XEQ 10
                       [   4  36  18  32  37 ]                         after loading this short routine:
                       [   9   8   39  20  33 ]
                       [ 16  23  21  10  31 ] ]                       LBL 10   ENTER^   X^2   41   MOD   STO IND Y   X<>Y   DSE X   GTO 10
                             [ [  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   b_2,4,3
-First,   we have a   4x5x3  hypermatrix  therefore:   4.005003   STO 00
-Then,  3  ENTER^
            4  ENTER^
            2  XEQ "DHT3"  >>>    b_2,4,3 =  3.491236125

Note:

-Synthetic register a is used.
-So do not call this program as more than a first level subroutine.
 

Continuous Hartley Transform
 

-Actually, Hartley has defined the real transform  H(x) =  (2.pi) ^-1/2  §_-infinity^+infinity  f(t) ( cos x.t + sin x.t ) dt      ( § = integral )
-This transform is strictly reciprocal.

-In the following program, we assume  f(t)  is defined over an interval  [ a ; b ]    ( f = 0 elsewhere )  where n = b - a  is an even integer
-The integral is estimated by the Filon's formula.
 

Data Registers:               ( Registers R00 thru Rnn are to be initialized before executing "HRTL" )

                                   •  R00 = f(a)    •  R01 = f(a+1)   •  R02 = f(a+2)  , .......... ,  •  Rnn = f(b) = f(a+n)

Flags: /
Subroutines: /

-This program uses synthetic register a
-So, it must not be called as more than a first level subroutine.
 
 

STACK	INPUTS	OUTPUTS
Z	a	/
Y	b	/
X	x	H(x)
L	/	x

  where  b-a =  n  must be  even  and  x  is expressed in  radians

Example1:      f  is only known by the following samples, and  f = 0  if  t < 3  or  t > 7
 
 

t	3	4	5	6	7
f(t)	1	2	1	-2	-7

-We store the f-values into registers  R00 thru R04:   1  STO 00   2  STO 01   1  STO 02   -2  STO 03   -7  STO 04
 and then, for instance:

  3  ENTER^
  7  ENTER^
  2  XEQ "HRTL"   >>>>   H(2) = -1.5466  and similarly:
 
 

x	-3	-2	-1	0	1	2	3
H(x)	0.3062	-1.4017	-1.1278	-0.5319	-3.8603	-1.5466	-1.4954

-Actually, we had   f(t) = -14 +8.t - t²   and these results are exact because Filon's formulae produce perfect accuracy if  f is a polynomial of degree 2

-If you compute more H(x)-values , you'll see the following graph ( approximately ):

                                                 H(x)
                                 *                |                           *
                               *   *             |        *               *   *
                    *       *      *        *  |     *   *           *       *        *
      --------*--*--*-----*---*--*|--*----*----- *-------*--*---*------------------------- x
                          *          *  *       |*          *      *               *
                                        *         |              *  *
                                                   |                *
 

Example2:     We take  f(t) =  e ^-t/2  for  t >= 0  and  f(t) = 0  if  t < 0

- H(x) may be expressed analytically and  H(x) = 2(1+2x)/(1+4x²) (2.pi) ^-1/2

-If we use f(0) , f(1) , .............. , f(16) to represent the function  ( these 17 numbers are to be stored in registers R00 thru R16 ), we obtain the following results:
 
 

x	-1	-0.5	-0.25	0	0.2071	1	2
HRTL	-0.1584	0.0013	0.3197	0.7979	0.9636	0.4797	0.2366
exact	-0.1596	0	0.3192	0.7979	0.9631	0.4787	0.2347

-The accuracy is quite satisfactory.
-The graph of the Hartley transform looks like this:

                                                              H(x)
                                                                 |
                                                                 |  *                               H(x) is maximum for x = 0.2071
                                                                 |*   *
                                                                 |          *
                                                                 |                       *
                                                                 |                                                          *
  ---------------------------------------* -|----------------------------------------------- x
            *                                                   |
                                     *                 *       |
                                                                 |
 

Z-Transform
 

-The unilateral  Z-transform of a sequence a₀ a₁ ............ a_n ........  is defined by  Z(x) = SUM _{k=0 to infinity}  a_k / x^k
-"ZTF" calculates the sum from k=0 to k=n , n must be specified by the user.

-The sequence (a_n) only contains real number.
 

Data Registers:                ( Registers R00 thru Rnn are to be initialized before executing "ZTF" )

                    •  R00 = a₀        •  R01 = a₁        •  R02 = a₂       ........................              •  Rnn = a_n

Flag:   F02
 
 

STACK	INPUTS	OUTPUTS
Y	n	n
X	x	Z(x)
L	/	x

Example:   Assuming  { a₀ a₁ ............ a_n } = { 1 2 3 2 1 1/2 1/3 1/4 1/5 1/6 1/7 } , store these 11 numbers into R00 thru R10   ( n = 10 )

-To compute  Z(PI)

       10    ENTER^
       PI    XEQ "ZTF"  >>>>   Z(PI) = 2.017443943
 

Complex Z-Transform
 

-The same method is used for complex arguments.
 

Data Registers:                ( Registers R00 thru Rnn are to be initialized before executing "ZTFZ" )

                    •  R00 = a₀        •  R01 = a₁        •  R02 = a₂       ........................              •  Rnn = a_n

Flag:   F02
 
 

STACK	INPUTS	OUTPUTS
Z	n	/
Y	y	V
X	x	U

  Where    Z(x+i.y) = U + i.V

Example:   With the same sequence, compute  Z(2+3.i)

     10  ENTER^
      3   ENTER^
      2   XEQ "ZTFZ"    >>>>     1.173222278
                                    RDN    -0.677712005

 So,  Z(2+3.i) = 1.173222278 - 0.677712005 i
 

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
[4]  http://mathworld.wolfram.com/Z-Transform.html