Fast Fourier Transform

[_vince]

A fast algorithm for computing the DFT for N=2^k

To use in your own programs you just need SUB fft or SUB rfft if dealing with real only signals.
i.e.:

Code: QB64: [Select]

1. const n = 512
2.  
3. dim shared pi as double
4. pi = 4*atn(1)
5.  
6. dim x_r(n-1), x_i(n-1)
7. dim xx_r(n-1), xx_i(n-1)
8.  
9. 'create signal
10. for i=0 to n-1
11.         x_r(i) = 100*sin(2*pi*62.27*i/n) + 25*cos(2*pi*132.27*i/n)
12.         x_i(i) = 0
13. next
14.  
15. fft xx_r(), xx_i(), x_r(), x_i(), n
16.  

SUB fft
x_r, x_i is the input signal real and complex values, must be arrays of size n
xx_r, xx_i must also be arrays of size n and will contain output

Code: QB64: [Select]

1. sub fft(xx_r(), xx_i(), x_r(), x_i(), n)
2.         dim w_r as double, w_i as double, wm_r as double, wm_i as double
3.         dim u_r as double, u_i as double, v_r as double, v_i as double
4.  
5.         log2n = log(n)/log(2)
6.  
7.         'bit rev copy
8.         for i=0 to n - 1
9.                 rev = 0
10.                 for j=0 to log2n - 1
11.                         if i and (2^j) then rev = rev + (2^(log2n - 1 - j))
12.                 next
13.  
14.                 xx_r(i) = x_r(rev)
15.                 xx_i(i) = x_i(rev)
16.         next
17.  
18.  
19.         for i=1 to log2n
20.                 m = 2^i
21.                 wm_r = cos(-2*pi/m)
22.                 wm_i = sin(-2*pi/m)
23.  
24.                 for j=0 to n - 1 step m
25.                         w_r = 1
26.                         w_i = 0
27.  
28.                         for k=0 to m/2 - 1
29.                                 p = j + k
30.                                 q = p + (m \ 2)
31.  
32.                                 u_r = w_r*xx_r(q) - w_i*xx_i(q)
33.                                 u_i = w_r*xx_i(q) + w_i*xx_r(q)
34.                                 v_r = xx_r(p)
35.                                 v_i = xx_i(p)
36.  
37.                                 xx_r(p) = v_r + u_r
38.                                 xx_i(p) = v_i + u_i
39.                                 xx_r(q) = v_r - u_r
40.                                 xx_i(q) = v_i - u_i
41.  
42.                                 u_r = w_r
43.                                 u_i = w_i
44.                                 w_r = u_r*wm_r - u_i*wm_i
45.                                 w_i = u_r*wm_i + u_i*wm_r
46.                         next
47.                 next
48.         next
49. end sub
50.  

SUB rfft
This is about 2x faster than above fft if dealing with only real values
x_r is the input signal, must be an array of size n
xx_r, xx_i must also be arrays of size n and will contain output
Since in the real FFT the second half of the output are complex conjugates of the first half, you may wish to not compute the second half and remove the last 4 lines in this SUB. I left it in just in case.

Code: QB64: [Select]

1. sub rfft(xx_r(), xx_i(), x_r(), n)
2.         dim w_r as double, w_i as double, wm_r as double, wm_i as double
3.         dim u_r as double, u_i as double, v_r as double, v_i as double
4.  
5.         log2n = log(n/2)/log(2)
6.  
7.         for i=0 to n/2 - 1
8.                 rev = 0
9.                 for j=0 to log2n - 1
10.                         if i and (2^j) then rev = rev + (2^(log2n - 1 - j))
11.                 next
12.  
13.                 xx_r(i) = x_r(2*rev)
14.                 xx_i(i) = x_r(2*rev + 1)
15.         next
16.  
17.         for i=1 to log2n
18.                 m = 2^i
19.                 wm_r = cos(-2*pi/m)
20.                 wm_i = sin(-2*pi/m)
21.  
22.                 for j=0 to n/2 - 1 step m
23.                         w_r = 1
24.                         w_i = 0
25.  
26.                         for k=0 to m/2 - 1
27.                                 p = j + k
28.                                 q = p + (m \ 2)
29.  
30.                                 u_r = w_r*xx_r(q) - w_i*xx_i(q)
31.                                 u_i = w_r*xx_i(q) + w_i*xx_r(q)
32.                                 v_r = xx_r(p)
33.                                 v_i = xx_i(p)
34.  
35.                                 xx_r(p) = v_r + u_r
36.                                 xx_i(p) = v_i + u_i
37.                                 xx_r(q) = v_r - u_r
38.                                 xx_i(q) = v_i - u_i
39.  
40.                                 u_r = w_r
41.                                 u_i = w_i
42.                                 w_r = u_r*wm_r - u_i*wm_i
43.                                 w_i = u_r*wm_i + u_i*wm_r
44.                         next
45.                 next
46.         next
47.  
48.         xx_r(n/2) = xx_r(0)
49.         xx_i(n/2) = xx_i(0)
50.  
51.         for i=1 to n/2 - 1
52.                 xx_r(n/2 + i) = xx_r(n/2 - i)
53.                 xx_i(n/2 + i) = xx_i(n/2 - i)
54.         next
55.  
56.         dim xpr as double, xpi as double
57.         dim xmr as double, xmi as double
58.  
59.         for i=0 to n/2 - 1
60.                 xpr = (xx_r(i) + xx_r(n/2 + i)) / 2
61.                 xpi = (xx_i(i) + xx_i(n/2 + i)) / 2
62.  
63.                 xmr = (xx_r(i) - xx_r(n/2 + i)) / 2
64.                 xmi = (xx_i(i) - xx_i(n/2 + i)) / 2
65.  
66.                 xx_r(i) = xpr + xpi*cos(2*pi*i/n) - xmr*sin(2*pi*i/n)
67.                 xx_i(i) = xmi - xpi*sin(2*pi*i/n) - xmr*cos(2*pi*i/n)
68.         next
69.  
70.         'symmetry, complex conj
71.         for i=0 to n/2 - 1
72.                 xx_r(n/2 + i) = xx_r(n/2 - 1 - i)
73.                 xx_i(n/2 + i) =-xx_i(n/2 - 1 - i)
74.         next
75. end sub
76.  

As an example FYI, here is a direct unoptimized DFT for comparison. It is much slower but can be used with any n

Code: QB64: [Select]

1. sub dft(xx_r(), xx_i(), x_r(), x_i(), n)
2.         for i=0 to n-1
3.                 xx_r(i) = 0
4.                 xx_i(i) = 0
5.                 for j=0 to n-1
6.                         xx_r(i) = xx_r(i) + x_r(j)*cos(2*pi*i*j/n) + x_i(j)*sin(2*pi*i*j/n)
7.                         xx_i(i) = xx_i(i) - x_r(j)*sin(2*pi*i*j/n) + x_i(j)*cos(2*pi*i*j/n)
8.                 next
9.         next
10. end sub
11.  

if you wish to take the inverse FFT then do the following to get back x_r and x_i from xx_r and xx_i. I chose not to include a inverse flag in the sub for my own reasons, it should be easy enough to modify above SUBs to add an inverse option if you so wish.

Code: QB64: [Select]

1. 'inverse fft
2. for i=0 to n-1
3.         xx_i(i) = -xx_i(i)
4. next
5.  
6. fft x_r(), x_i(), xx_r(), xx_i(), n
7.  
8. for i=0 to n-1
9.         x_r(i) = x_r(i)/n
10.         x_i(i) = x_i(i)/n
11. next
12.  

Here is a basic example using and plotting all the SUBs:

Code: QB64: [Select]

1. 'defdbl a-z
2.  
3. const sw = 1024
4. const sh = 600
5.  
6. dim shared pi as double
7. 'pi = 2*asin(1)
8. pi = 4*atn(1)
9.  
10. declare sub rfft(xx_r(), xx_i(), x_r(), n)
11. declare sub fft(xx_r(), xx_i(), x_r(), x_i(), n)
12. declare sub dft(xx_r(), xx_i(), x_r(), x_i(), n)
13.  
14.  
15. dim x_r(sw-1), x_i(sw-1)
16. dim xx_r(sw-1), xx_i(sw-1)
17. dim t as double
18.  
19. for i=0 to sw-1
20.         x_r(i) = 100*sin(2*pi*62.27*i/sw) + 25*cos(2*pi*132.27*i/sw)
21.         x_i(i) = 0
22. next
23.  
24.  
25. 'screenres sw, sh, 32
26. screen _newimage(sw, sh, 32)
27.  
28. pset (0, sh/4 - x_r(0))
29. for i=0 to sw-1
30.         line -(i, sh/4 - x_r(i)), _rgb(100,100,100)
31. next
32.  
33. dft xx_r(), xx_i(), x_r(), x_i(), sw
34. pset (0, 3*sh/4 - 0.1*sqr(xx_r(0)*xx_r(0) + xx_i(0)*xx_i(0)) ), _rgb(0,255,0)
35. for i=0 to sw - 1
36.         line -(i, 3*sh/4 - 0.1*sqr(xx_r(i)*xx_r(i) + xx_i(i)*xx_i(i)) ), _rgb(0,255,0)
37. next
38. line (0, 3*sh/4)-step(sw,0), _rgb(0,255,0),,&h5555
39.  
40. t = timer
41. for i=0 to 50
42. fft xx_r(), xx_i(), x_r(), x_i(), sw
43. next
44. locate 1,1
45. print "50x fft ";timer - t
46.  
47. pset (0, 50+3*sh/4 - 0.1*sqr(xx_r(0)*xx_r(0) + xx_i(0)*xx_i(0)) ), _rgb(255,0,0)
48. for i=0 to sw - 1
49.         line -(i, 50+3*sh/4 - 0.1*sqr(xx_r(i)*xx_r(i) + xx_i(i)*xx_i(i)) ), _rgb(255,0,0)
50. next
51. line (0, 50+3*sh/4)-step(sw,0), _rgb(255,0,0),,&h5555
52.  
53.  
54. for i=0 to sw-1
55.         xx_r(i) = 0
56.         xx_i(i) = 0
57. next
58.  
59. t = timer
60. for i=0 to 50
61. rfft xx_r(), xx_i(), x_r(), sw
62. next
63. locate 2,1
64. print "50x rfft ";timer - t
65.  
66. pset (0, 100+3*sh/4 - 0.1*sqr(xx_r(0)*xx_r(0) + xx_i(0)*xx_i(0)) ), _rgb(255,255,0)
67. for i=0 to sw - 1
68.         line -(i, 100+3*sh/4 - 0.1*sqr(xx_r(i)*xx_r(i) + xx_i(i)*xx_i(i)) ), _rgb(255,255,0)
69. next
70. line (0, 100+3*sh/4)-step(sw,0), _rgb(255,255,0),,&h5555
71.  
72. sleep
73. system
74.  
75.  
76. sub rfft(xx_r(), xx_i(), x_r(), n)
77.         dim w_r as double, w_i as double, wm_r as double, wm_i as double
78.         dim u_r as double, u_i as double, v_r as double, v_i as double
79.  
80.         log2n = log(n/2)/log(2)
81.  
82.         for i=0 to n/2 - 1
83.                 rev = 0
84.                 for j=0 to log2n - 1
85.                         if i and (2^j) then rev = rev + (2^(log2n - 1 - j))
86.                 next
87.  
88.                 xx_r(i) = x_r(2*rev)
89.                 xx_i(i) = x_r(2*rev + 1)
90.         next
91.  
92.         for i=1 to log2n
93.                 m = 2^i
94.                 wm_r = cos(-2*pi/m)
95.                 wm_i = sin(-2*pi/m)
96.  
97.                 for j=0 to n/2 - 1 step m
98.                         w_r = 1
99.                         w_i = 0
100.  
101.                         for k=0 to m/2 - 1
102.                                 p = j + k
103.                                 q = p + (m \ 2)
104.  
105.                                 u_r = w_r*xx_r(q) - w_i*xx_i(q)
106.                                 u_i = w_r*xx_i(q) + w_i*xx_r(q)
107.                                 v_r = xx_r(p)
108.                                 v_i = xx_i(p)
109.  
110.                                 xx_r(p) = v_r + u_r
111.                                 xx_i(p) = v_i + u_i
112.                                 xx_r(q) = v_r - u_r
113.                                 xx_i(q) = v_i - u_i
114.  
115.                                 u_r = w_r
116.                                 u_i = w_i
117.                                 w_r = u_r*wm_r - u_i*wm_i
118.                                 w_i = u_r*wm_i + u_i*wm_r
119.                         next
120.                 next
121.         next
122.  
123.         xx_r(n/2) = xx_r(0)
124.         xx_i(n/2) = xx_i(0)
125.  
126.         for i=1 to n/2 - 1
127.                 xx_r(n/2 + i) = xx_r(n/2 - i)
128.                 xx_i(n/2 + i) = xx_i(n/2 - i)
129.         next
130.  
131.         dim xpr as double, xpi as double
132.         dim xmr as double, xmi as double
133.  
134.         for i=0 to n/2 - 1
135.                 xpr = (xx_r(i) + xx_r(n/2 + i)) / 2
136.                 xpi = (xx_i(i) + xx_i(n/2 + i)) / 2
137.  
138.                 xmr = (xx_r(i) - xx_r(n/2 + i)) / 2
139.                 xmi = (xx_i(i) - xx_i(n/2 + i)) / 2
140.  
141.                 xx_r(i) = xpr + xpi*cos(2*pi*i/n) - xmr*sin(2*pi*i/n)
142.                 xx_i(i) = xmi - xpi*sin(2*pi*i/n) - xmr*cos(2*pi*i/n)
143.         next
144.  
145.         'symmetry, complex conj
146.         for i=0 to n/2 - 1
147.                 xx_r(n/2 + i) = xx_r(n/2 - 1 - i)
148.                 xx_i(n/2 + i) =-xx_i(n/2 - 1 - i)
149.         next
150. end sub
151.  
152. sub fft(xx_r(), xx_i(), x_r(), x_i(), n)
153.         dim w_r as double, w_i as double, wm_r as double, wm_i as double
154.         dim u_r as double, u_i as double, v_r as double, v_i as double
155.  
156.         log2n = log(n)/log(2)
157.  
158.         'bit rev copy
159.         for i=0 to n - 1
160.                 rev = 0
161.                 for j=0 to log2n - 1
162.                         if i and (2^j) then rev = rev + (2^(log2n - 1 - j))
163.                 next
164.  
165.                 xx_r(i) = x_r(rev)
166.                 xx_i(i) = x_i(rev)
167.         next
168.  
169.  
170.         for i=1 to log2n
171.                 m = 2^i
172.                 wm_r = cos(-2*pi/m)
173.                 wm_i = sin(-2*pi/m)
174.  
175.                 for j=0 to n - 1 step m
176.                         w_r = 1
177.                         w_i = 0
178.  
179.                         for k=0 to m/2 - 1
180.                                 p = j + k
181.                                 q = p + (m \ 2)
182.  
183.                                 u_r = w_r*xx_r(q) - w_i*xx_i(q)
184.                                 u_i = w_r*xx_i(q) + w_i*xx_r(q)
185.                                 v_r = xx_r(p)
186.                                 v_i = xx_i(p)
187.  
188.                                 xx_r(p) = v_r + u_r
189.                                 xx_i(p) = v_i + u_i
190.                                 xx_r(q) = v_r - u_r
191.                                 xx_i(q) = v_i - u_i
192.  
193.                                 u_r = w_r
194.                                 u_i = w_i
195.                                 w_r = u_r*wm_r - u_i*wm_i
196.                                 w_i = u_r*wm_i + u_i*wm_r
197.                         next
198.                 next
199.         next
200. end sub
201.  
202. sub dft(xx_r(), xx_i(), x_r(), x_i(), n)
203.         for i=0 to n-1
204.                 xx_r(i) = 0
205.                 xx_i(i) = 0
206.                 for j=0 to n-1
207.                         xx_r(i) = xx_r(i) + x_r(j)*cos(2*pi*i*j/n) + x_i(j)*sin(2*pi*i*j/n)
208.                         xx_i(i) = xx_i(i) - x_r(j)*sin(2*pi*i*j/n) + x_i(j)*cos(2*pi*i*j/n)
209.                 next
210.         next
211. end sub
212.  

A common use for FFT is filtering unwanted frequencies. The following example demonstrates this, with screenshots attached

Code: QB64: [Select]

1. 'defdbl a-z
2.  
3. const sw = 512
4. const sh = 600
5.  
6. dim shared pi as double
7. 'pi = 2*asin(1)
8. pi = 4*atn(1)
9.  
10. declare sub fft(xx_r(), xx_i(), x_r(), x_i(), n)
11.  
12.  
13. dim x_r(sw-1), x_i(sw-1)
14. dim xx_r(sw-1), xx_i(sw-1)
15. dim t as double
16.  
17. for i=0 to sw-1
18.         'x_r(i) = 100*sin(2*pi*62.27*i/sw) + 25*cos(2*pi*132.27*i/sw)
19.         x_r(i) = 100*sin(0.08*i) + 25*cos(i)
20.         x_i(i) = 0
21. next
22.  
23.  
24. 'screenres sw, sh, 32
25. screen _newimage(sw*2, sh, 32)
26.  
27. 'plot input signal
28. pset (0, sh/4 - x_r(0))
29. for i=0 to sw-1
30.         line -(i, sh/4 - x_r(i)), _rgb(255,0,0)
31. next
32. line (0, sh/4)-step(sw,0), _rgb(255,0,0),,&h5555
33. color _rgb(255,0,0)
34. _printstring (0,0), "input signal"
35.  
36. fft xx_r(), xx_i(), x_r(), x_i(), sw
37.  
38. 'plot its fft
39. pset (0, 50+3*sh/4 - 0.01*sqr(xx_r(0)*xx_r(0) + xx_i(0)*xx_i(0)) ), _rgb(255,255,0)
40. for i=0 to sw/2
41.         line -(i*2, 50+3*sh/4 - 0.01*sqr(xx_r(i)*xx_r(i) + xx_i(i)*xx_i(i)) ), _rgb(255,255,0)
42. next
43. line (0, 50+3*sh/4)-step(sw,0), _rgb(255,255,0),,&h5555
44.  
45.  
46. 'set unwanted frequencies to zero
47. for i=50 to sw/2
48.         xx_r(i) = 0
49.         xx_i(i) = 0
50.         xx_r(sw - i) = 0
51.         xx_i(sw - i) = 0
52. next
53.  
54. 'plot fft of filtered signal
55. pset (sw, 50+3*sh/4 - 0.01*sqr(xx_r(0)*xx_r(0) + xx_i(0)*xx_i(0)) ), _rgb(255,255,0)
56. for i=0 to sw/2
57.         line -(sw + i*2, 50+3*sh/4 - 0.01*sqr(xx_r(i)*xx_r(i) + xx_i(i)*xx_i(i)) ), _rgb(0,155,255)
58. next
59. line (sw, 50+3*sh/4)-step(sw,0), _rgb(0,155,255),,&h5555
60.  
61. 'take inverse fft
62. for i=0 to sw-1
63.         xx_i(i) = -xx_i(i)
64. next
65.  
66. fft x_r(), x_i(), xx_r(), xx_i(), sw
67.  
68. for i=0 to sw-1
69.         x_r(i) = x_r(i)/sw
70.         x_i(i) = x_i(i)/sw
71. next
72.  
73.  
74. 'plot filtered signal
75. pset (sw, sh/4 - x_r(0))
76. for i=0 to sw-1
77.         line -(sw + i, sh/4 - x_r(i)), _rgb(0,255,0)
78. next
79. line (sw, sh/4)-step(sw,0), _rgb(0,255,0),,&h5555
80.  
81. color _rgb(0,255,0)
82. _printstring (sw,0), "filtered signal"
83.  
84. sleep
85. system
86.  
87. sub fft(xx_r(), xx_i(), x_r(), x_i(), n)
88.         dim w_r as double, w_i as double, wm_r as double, wm_i as double
89.         dim u_r as double, u_i as double, v_r as double, v_i as double
90.  
91.         log2n = log(n)/log(2)
92.  
93.         'bit rev copy
94.         for i=0 to n - 1
95.                 rev = 0
96.                 for j=0 to log2n - 1
97.                         if i and (2^j) then rev = rev + (2^(log2n - 1 - j))
98.                 next
99.  
100.                 xx_r(i) = x_r(rev)
101.                 xx_i(i) = x_i(rev)
102.         next
103.  
104.  
105.         for i=1 to log2n
106.                 m = 2^i
107.                 wm_r = cos(-2*pi/m)
108.                 wm_i = sin(-2*pi/m)
109.  
110.                 for j=0 to n - 1 step m
111.                         w_r = 1
112.                         w_i = 0
113.  
114.                         for k=0 to m/2 - 1
115.                                 p = j + k
116.                                 q = p + (m \ 2)
117.  
118.                                 u_r = w_r*xx_r(q) - w_i*xx_i(q)
119.                                 u_i = w_r*xx_i(q) + w_i*xx_r(q)
120.                                 v_r = xx_r(p)
121.                                 v_i = xx_i(p)
122.  
123.                                 xx_r(p) = v_r + u_r
124.                                 xx_i(p) = v_i + u_i
125.                                 xx_r(q) = v_r - u_r
126.                                 xx_i(q) = v_i - u_i
127.  
128.                                 u_r = w_r
129.                                 u_i = w_i
130.                                 w_r = u_r*wm_r - u_i*wm_i
131.                                 w_i = u_r*wm_i + u_i*wm_r
132.                         next
133.                 next
134.         next
135. end sub
136.  

Another common thing seems to be wanting to find an exact frequency of a signal. The following demonstrates calculating the frequency of a sine corrupted in noise and applying bin to frequency correction for extra accuracy, screenshot attached.

Code: QB64: [Select]

1. 'defdbl a-z
2.  
3. const sw = 1024
4. const sh = 600
5.  
6. dim shared pi as double
7. 'pi = 2*asin(1)
8. pi = 4*atn(1)
9.  
10. declare sub rfft(xx_r(), xx_i(), x_r(), n)
11.  
12. dim x_r(sw-1), x_i(sw-1)
13. dim xx_r(sw-1), xx_i(sw-1)
14. dim t as double
15.  
16. for i=0 to sw-1
17.         x_r(i) = 100*sin(2*pi*(sw*2000/44000)*i/sw) + (100*rnd - 50)
18. next
19.  
20.  
21. 'screenres sw, sh, 32
22. screen _newimage(sw, sh, 32)
23.  
24. 'plot signal
25. pset (0, sh/4 - x_r(0))
26. for i=0 to sw-1
27.         line -(i, sh/4 - x_r(i)), _rgb(255,0,0)
28. next
29. line (0, sh/4)-step(sw,0), _rgb(255,0,0),,&h5555
30.  
31. _printstring (0, 0), "2000 kHz signal with RND noise sampled at 44 kHz in 1024 samples"
32.  
33.  
34. rfft xx_r(), xx_i(), x_r(), sw
35.  
36. 'plot its fft
37. pset (0, 50+3*sh/4 - 0.005*sqr(xx_r(0)*xx_r(0) + xx_i(0)*xx_i(0)) ), _rgb(255,255,0)
38. for i=0 to sw/2
39.         line -(i*2, 50+3*sh/4 - 0.005*sqr(xx_r(i)*xx_r(i) + xx_i(i)*xx_i(i)) ), _rgb(255,255,0)
40. next
41. line (0, 50+3*sh/4)-step(sw,0), _rgb(255,255,0),,&h5555
42.  
43.  
44. 'find peak
45. dim max as double, d as double
46. max = 0
47. m = 0
48. for i=0 to sw/2
49.         d = 0.01*sqr(xx_r(i)*xx_r(i) + xx_i(i)*xx_i(i))
50.         if d > max then
51.                 max = d
52.                 m = i
53.         end if
54. next
55.  
56. _printstring (0, sh/2), "m_peak ="+str$(m)
57. _printstring (0, sh/2 + 16), "f_peak = m_peak * 44 kHz / 1024 samples = "+str$(m*44000/1024)+" Hz"
58.  
59. 'apply frequency correction, only works for some signals
60. dim c as double
61. dim u_r as double, u_i as double
62. dim v_r as double, v_i as double
63.  
64. u_r = xx_r(m - 1) - xx_r(m + 1)
65. u_i = xx_i(m - 1) - xx_i(m + 1)
66. v_r = 2*xx_r(m) - xx_r(m - 1) - xx_r(m + 1)
67. v_i = 2*xx_i(m) - xx_i(m - 1) - xx_i(m + 1)
68. c = (u_r*v_r + u_i*v_i)/(v_r*v_r + v_i*v_i)
69.  
70. _printstring (0, sh/2 + 2*16), "f_corrected = "+str$((m+c)*44000/1024)+" Hz"
71.  
72. sleep
73. system
74.  
75.  
76. sub rfft(xx_r(), xx_i(), x_r(), n)
77.         dim w_r as double, w_i as double, wm_r as double, wm_i as double
78.         dim u_r as double, u_i as double, v_r as double, v_i as double
79.  
80.         log2n = log(n/2)/log(2)
81.  
82.         for i=0 to n/2 - 1
83.                 rev = 0
84.                 for j=0 to log2n - 1
85.                         if i and (2^j) then rev = rev + (2^(log2n - 1 - j))
86.                 next
87.  
88.                 xx_r(i) = x_r(2*rev)
89.                 xx_i(i) = x_r(2*rev + 1)
90.         next
91.  
92.         for i=1 to log2n
93.                 m = 2^i
94.                 wm_r = cos(-2*pi/m)
95.                 wm_i = sin(-2*pi/m)
96.  
97.                 for j=0 to n/2 - 1 step m
98.                         w_r = 1
99.                         w_i = 0
100.  
101.                         for k=0 to m/2 - 1
102.                                 p = j + k
103.                                 q = p + (m \ 2)
104.  
105.                                 u_r = w_r*xx_r(q) - w_i*xx_i(q)
106.                                 u_i = w_r*xx_i(q) + w_i*xx_r(q)
107.                                 v_r = xx_r(p)
108.                                 v_i = xx_i(p)
109.  
110.                                 xx_r(p) = v_r + u_r
111.                                 xx_i(p) = v_i + u_i
112.                                 xx_r(q) = v_r - u_r
113.                                 xx_i(q) = v_i - u_i
114.  
115.                                 u_r = w_r
116.                                 u_i = w_i
117.                                 w_r = u_r*wm_r - u_i*wm_i
118.                                 w_i = u_r*wm_i + u_i*wm_r
119.                         next
120.                 next
121.         next
122.  
123.         xx_r(n/2) = xx_r(0)
124.         xx_i(n/2) = xx_i(0)
125.  
126.         for i=1 to n/2 - 1
127.                 xx_r(n/2 + i) = xx_r(n/2 - i)
128.                 xx_i(n/2 + i) = xx_i(n/2 - i)
129.         next
130.  
131.         dim xpr as double, xpi as double
132.         dim xmr as double, xmi as double
133.  
134.         for i=0 to n/2 - 1
135.                 xpr = (xx_r(i) + xx_r(n/2 + i)) / 2
136.                 xpi = (xx_i(i) + xx_i(n/2 + i)) / 2
137.  
138.                 xmr = (xx_r(i) - xx_r(n/2 + i)) / 2
139.                 xmi = (xx_i(i) - xx_i(n/2 + i)) / 2
140.  
141.                 xx_r(i) = xpr + xpi*cos(2*pi*i/n) - xmr*sin(2*pi*i/n)
142.                 xx_i(i) = xmi - xpi*sin(2*pi*i/n) - xmr*cos(2*pi*i/n)
143.         next
144.  
145.         'symmetry, complex conj
146.         for i=0 to n/2 - 1
147.                 xx_r(n/2 + i) = xx_r(n/2 - 1 - i)
148.                 xx_i(n/2 + i) =-xx_i(n/2 - 1 - i)
149.         next
150. end sub
151.  

[Ashish]

I don't understand it but I think it is something good, something advance... :)

[bplus]

What is DFT? https://en.wikipedia.org/wiki/Discrete_Fourier_transform

Ah, maybe you haven't heard, it's discreet.

Or it is Density Function Theory, nope not this code...

with an S it means Don't Forget To Smile :D

[anttiylikoski]

I have planned to author an article about the Fourier Transform for a Scandinavian radio amateur magazine.

It would be nice, if I could print your FFT programs in that article:

1.  giving the full copyrights credits, according to the immaterial rights laws

2.  only with the author's permission.

I would list the BASIC code, and, explain the code in the article.

See the Finland Radio Amateurs Association  https://sral.info/   

[_vince]

free to use for anything, no credit required but a link to this forum or post appreciated

[Petr]

Hi _Vince,

i try using your last source code, but it is not QB64 source. What is SHL statement? I'm playing with your program (with the sources above). The program obtains one signal from the signals. Could you please show me how to calculate the frequency for this acquired signal? I assume that a time window must be used to obtain the frequency, and the frequency will be calculated, for example, in a sample of 441 samples for a recording length of 10 milliseconds (for a sampling frequency of 44100). Or it will be enough to calculate the number of changes (ie if the previous sample increases and the next decreases, it is one change, if the previous increases and the next also, it is not a change), then the frequency would be calculated as 1/441 * number of changes - can this be count this way?

[_vince]

Hah, didn't think anyone would ever use this.  I've revised the original post with easier to use code and better examples.

To answer your question, when you take an FFT of, say, 1024 samples sampled at 44 kHz you only have 1024 bins giving you about 44 Hz of resolution. That's the accuracy with which you can pick out individual frequencies, the more samples the better accuracy.  If you want to know the exact frequency of a particular bin, say bin m, then

f = m * f_s / 1024

where f_s is the sampling frequency, ie 44100 Hz.  That is, the magnitude of that frequency f is sqr(xx_r(m)^2 + xx_i(m)^2).  I've included this concept as an example in above post.

[Petr]

_Vince, thank you very much for this detailed answer, I will do some experiments to it for maximal possible understandt.

[Sanmayce]

Many thanks for sharing!
If I am successful, in next days will share my BAR-JUMPER - to enable visual feedback to the MP3/WAV files...

[MasterGy]

Great !! Congratulations! It is a very useful thing to deal with. Rarely, but then very! About 15 years ago, I controlled an engraving machine with qb from dos. I experimented with using a microphone to tell the computer when the milling head would reach the object. This was necessary to be able to engrave on an irregular surface. Unfortunately, the solution was inaccurate due to the noises, and I remember looking for this fourier transform thing at the time, because I wanted to take out the frequency where the milling motor was buzzing.
That’s what you’re dealing with now, I should have really, really used to! Congratulations!

http://tukorgravirozas.fw.hu/mastergy/4pohar/pg.avi