41Z Routine of the week: DFT

[Ángel Martin]

Discrete Fourier Transform with the 41Z Module.

Not much of a groundbreaking algorithm, the program below uses the straightforward definition to obtain the DFT and inverse DFT of a set of N complex data points, stored in complex data registers ZR01 to ZR0n. You can use the sub-function ZINPT from the module to enter all values.

Just enter the sample size in X, and execute ZDFT. The results are stored in complex registers ZRn+1 to ZR(2n). You can use the sub-function ZOUPT to see the values.

If you want to apply the inverse DFT to the results you can use REGSWAP to reset the locations: 2,00(2+2n)00(2n). i.e. 2,010|008 for n=4. Then run ZIDFT to obtain the original values (except rounding).

Example:

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

IDFT [ 16+20.i ; -8 ; -4-4.i ; -8.i ] = [ 1+2i ; 3+4i ; 5+6i ; 7+8i ]

Code:

01   LBL "ZDFT" 
02   CF 01
03   GTO 02
04   LBL "ZIDFT"
05   SF 01
06   LBL 02
07   STO N
08   E3/E+
09   STO M
10   LBL 01
11   VIEW M
12   CLX
13   STO 00
14   STO 01
15   RCL N
16   E3/E+
17   STO O
18   RCL M
19   INT 
20   E
21   -
22   PI
23   *
24   ST+ X
25   RCL N
26   /
27   STO P
28   LBL 00
29   RCL O
30   INT
31   E
32   -
33   RCL P
34   *
35   FC? 01
36   CHS
37   1
38   P-R
39   ZRC* IND O
40   ZST+ 00
41   ISG O
42   GTO 00
43   RCL M
44   INT
45   RCL N
46   +
47   ZRCL 00
48   FC? 01
40   GTO 02
41   RCL N
42   ST/ Z
43   /
44   LBL 02
45   ZSTO IND Z
46   ISG M
47   GTO 01
48   END

As always, the execution time can be rather long for large size samples. Using the 41-CL or an emulator in Turbo mode is recommended.
Feedback is welcome - would be great to get an FFT implementation as well... any volunteers?

[gjmcclure]

'Angel,

I am having a problem executing this program... specifically steps 27 thru 33. You use the P register for storage at step 27, then enter E at step 31. I think this blows part of the P register, so your recall at step 33 won't work properly, or do I understand your use of E as a step properly? I am thinking you are using E as a faster "1" entry.

Greg.

[Ángel Martin]

Good catch Greg, thanks for the warning. Indeed using "E" is not a good idea if register P is also used for storage, I had completely forgotten about that... since I moved to the all-MCODE version in the updated module attached ;-)

So just change the "E" to a "1" and the program should work fine.

Here's an alternative FOCAL version that should always work, just to close this issue. It does the intermediate multiplication in-situ so registers R00 and R01 are free to use now.

Code:

01    LBL "ZDFT"
02    CF 01
03    GTO 00
04    LBL "ZIDFT"
05    SF 01
06    LBL 00
07    STO 00
08    E3/E+
09    STO M(5)
10    LBL 01
11    VIEW M(5)
12    RCL 00
13    STO N(6)
14    E3/3+
15    STO O(7)
16    RCL 5(M)
17    INT
18    ST+  N(6)
19    E
20    -
21    PI
22    *
23    ST+ X(3)
24    RCL 00
25     /
26    STO 01
27    CLZ
28    ZSTO IND N(6)
29    LBL 02
30    RCL 0(7)
31    INT
32    E
33    -
34    RCL 01
35    *
36    FC? 01
37    CHS
38    E
39    P-R
40    ZRC* IND O(7)
41    ZST+  IND N(6)
42    ISG O(7)
43    GTO 02
44    FC? 01
45    GTO 00
46    ZRCL IND 01
47    RCL 00
48    ST/ Z
49     /
50    ZSTO IND 01
51    LBL 00
52    ISG M(5)
53    GTO 01
54    END

Latest revision of the 41Z deluxe is attached.


Cheers,
Ángel