Discrete fourier transform on HP-41

Coert Vonk

The discrete Fourier transform transforms a time domain samples to the frequency domain.  It is used in digital signal processing applications such as audio compression (i.e. AAC), linear filtering and spectrum analysis.\(\) I wrote this in 1987.

Theory

Earlier, we introduced the Fourier Transform. The Discrete Fourier Transform (DFT) is used when all you have available are samples of the function, rather than the function itself. The Discrete Fourier Transform (DTF) is a numerical approximation to the Fourier Transform. [1]

$$ F(\omega_k)\ \triangleq \sum_{n=0}^{N-1}e^{-j\omega_kt_n}\ f(t_n)\label{eq:dft}\tag{4} $$

where

$$ \begin{aligned} k&=0,1,2,\ldots,N-1\\ \sum_{n=0}^{N-1}f(n)\ &\triangleq\ f(0)+f(1)+\ldots+f(N-1)\\ f(t_n)\ &\triangleq\ \text{input signal amplitude (}\mathbb{R},\mathbb{C}\text{) at time }t_n\\ t_n &\triangleq\ nT=n^{th}\text{ sampling instant [s], }n\in \mathbb{N_0}\\ T\ &\triangleq\ \text{sample interval [s]}\\ F(\omega_k)\ &\triangleq\ \text{spectrum of }f\text{ (complex valued) at frequency }\omega_k\\ \omega_k\ &\triangleq\ \frac{2k\pi}{NT} = k^{th}\text{ frequency sample [rad/s]}\\ f_s\ &\triangleq\ \frac{1}{T} = \text{sampling rate [Hz]}\\ N\ &\triangleq\ \text{number of samples} \end{aligned} $$

Examples

Example 1

Sampling a sinus function using 8 samples a period. 

01  LBL "FIE"
02  45
03	 *
04  SIN
05  VIEW X
06  RTN

The calculation steps are shown in the table below.

Sampling a sinus function using 8 samples a period
n An Bn Mn Φn
1 0.9003 -0.3729 0.975 -0.125*π
2 0.0000 0.0000 0.000 0.250*π
3 0.0000 0.0000 0.000 0.000*π
4 0.0000 0.0000 0.000 0.000*π
5 0.0000 0.0000 0.000 0.000*π
6 0.0000 0.0000 0.000 -0.125*π
7 0.1286 0.0533 0.139 0.125*π
8 0.0000 0.0000 0.0000 0.000*π

Example 2

Same example, but with a hamming correction. 

01  LBL "FIE"
02  45
03	 *
04  SIN
05  XEQ "HC"
06  VIEW X
07  RTN

The calculation steps are shown in the table below.

Same example, but with a hamming correction
n An Bn Mn Φn
1 0.4505  -0.1865  0.487  -0.125*π
2 -0.1592  0.1592  0.225  0.750*π
3 0.0000 0.0000  0.000  -0.687*π
4 0.0000 0.0000  0.000  0.000*π
5 0.0000 0.0000  0.000  -0.687*π
6  -0.0531  -0.531  0.075  -0.750*π
7  0.0643  0.266  0.070  0.125*π
8 0.0000 0.0000  0.000  0.000*π

Example 3

Sampling a sinus function using 10 samples with a sample time of 1/8 * Tsignal: 

01  LBL "FIE"
02  360
03  *
04  8
05  /
06  SIN
07  VIEW X
08  RTN

The calculation steps are shown in the table below.

sinus function using 10 samples with a sample time of 1/8 * Tsignal
n An Bn Harm MODn Φn
1 0.7201 0.3802 0.8 0.814 0.155*π
2 -0.3335 -0.0787 1.6 0.343 -0.926*π
3 -0.1650 -0.0206 2.4 0.166 -0.961*π
4 -0.1146 0.0064 3.2 0.115 -0.982*π
5 -0.0900 0.0000 4.0 0.090 1.000*π
6 -0.0764 0.0043 4.8 0.077 0.982*π
7 -0.0707 0.0088 5.6 0.071 0.961*π
8 -0..0834 0.0197 6.4 0.086 0.926*π
9 -0.0800 -0.0422 7.2 0.090 -0.155*π
10 0.0000 0.0000 8.0 0.000 0.000*π

Example 4

Same example, but with the hamming correction 

01	LBL "FIE"
02	360
03	*
04	8
05	/
06	SIN
07	VIEW X
08	XEQ "HC"
09	RTN

The calculation steps are shown in the table below.

Same example, but with the hamming correction
n An Bn Harm MODn Φn
0 0.4263 0.2038 0.8 0.472 0.142*π
1 -0.316 -0.0532 1.6 0.321 -0.947*π
2 -0.0262 -0.0050 2.4 0.027 -0.060*π
3 0.0041 -0.0019 3.2 0.005 -0.134*π
4 0.0017 0.0000 4.0 0.002 0.000*π
5 0.0028 0.0012 4.8 0.003 0.134*π
6 0.0112 0.0021 5.6 0.011 0.060*π
7 -0.0791 0.0133 6.4 0.080 0.947*π
8 -0.0474 -0.0226 7.2 0.052 -0.0142*π
9 0.0000 0.000 8.0 0.000 0.000*π

Example 4

Manual input 

01  LBL "FIE"
02  "N"
03  FIX 0
04  ARCL X
05  >"?"
06  PROMPT
07  RTN

Source code

  • Requires X-Functions module on the HP-41cv
  • Available as source code, raw for the V41 emulator and barcode for the HP Wand (HP82153A)

Listing

Available through the repository 

01	LBL "DFT"
02	DEG
03	sREG 11		; ΣREG
04	CLS  	  	; CLΣ
05	"HARM=?"
06	PROMPT
07	STO 00
08	"M=?"
09	PROMPT
10	"MONSTER:"
11	AVIEW
12	STO 03
13	STO 01
14	/
15	STO 02
16	180
17	ST* 02
18	1
19	ST- 01
20	RCL 01
21	.01
22	+
23	 E3
24	/
25	STO 01
26	LBL 00
27	RCL 01
28	INT
29	ST+ X
30	1
31	+
32	RCL 02
33	*
34	ENTER^
35	SIN
36	X<>Y
37	COS
38	STO 04
39	X<>Y
40	STO 05
41	RCL 01
42	INT
43	XEQ "FIE"
44	RCL 05
45	RCL 04
46	RCL Z
47	*
48	X<>Y
49	LASTX
50	*
51	S+		; Σ+
52	ISG 01
53	GTO 00
54	RCL 02
55	SIN
56	RCL 00
57	/
58	PI
59	/
60	ST+ X
61	ST* 11
62	ST* 13
63	FIX 0
64	"A"
65	ARCL 00
66	>"="
67	FIX 4
68	ARCL 11
69	PROMPT
70	"B"
71	FIX 0
72	ARCL 00
73	>"="
74	FIX 4
75	ARCL 13
76	PROMPT
77	RCL 13
78	RCL 11
79	R-P
80	"AMP="
81	FIX 3
82	ARCL X
83	PROMPT
84	X<>Y
85	180
86	/
87	FIX 3
88	"ARG="
89	ARCL X
90	>"*PI"
91	AVIEW
92	RTN

; Hamming Correction

93	LBL "HC"
94	1
95	RCL 01
96	INT
97	RCL 03
98	/
99	360
100	*
101	COS
102	-
103	2
104	/
105	*
106	END

References

[1] Mathematics of the Discrete Fourier Transform (DFT)
Julius O. Smith III, April 13, 2007
W3K Publishing; 2nd edition, ISBN 097456074X

Leave a Reply

Your email address will not be published. Required fields are marked *

 

This site uses Akismet to reduce spam. Learn how your comment data is processed.