Text Only | Text with Attachments

QB64.org Forum

Active Forums => Programs => Topic started by: _vince on November 30, 2019, 02:34:17 pm

Title: Fast Fourier Transform
Post by: _vince on November 30, 2019, 02:34:17 pm
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.  
  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
  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
  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
  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
  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)
  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
  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.  
  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)
  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
  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)
  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)
  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
  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)
  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
  58.  
  59. t = timer
  60. for i=0 to 50
  61. rfft xx_r(), xx_i(), x_r(), sw
  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)
  70. line (0, 100+3*sh/4)-step(sw,0), _rgb(255,255,0),,&h5555
  71.  
  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
  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
  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
  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.  
  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)
  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
  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)
  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)
  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
  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)
  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)
  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
  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)
  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.  
  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
  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.  
  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)
  15.  
  16. for i=0 to sw-1
  17.         x_r(i) = 100*sin(2*pi*(sw*2000/44000)*i/sw) + (100*rnd - 50)
  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)
  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)
  41. line (0, 50+3*sh/4)-step(sw,0), _rgb(255,255,0),,&h5555
  42.  
  43.  
  44. 'find peak
  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
  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
  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.  
  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
  151.  
Title: Re: Fast Fourier Transform
Post by: Ashish on December 01, 2019, 11:04:03 am
I don't understand it but I think it is something good, something advance... :)
Title: Re: Fast Fourier Transform
Post by: bplus on December 01, 2019, 11:16:03 am
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
Title: Re: Fast Fourier Transform
Post by: anttiylikoski on December 25, 2019, 09:37:35 am
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/   

Title: Re: Fast Fourier Transform
Post by: _vince on December 29, 2019, 11:03:55 pm
free to use for anything, no credit required but a link to this forum or post appreciated
Title: Re: Fast Fourier Transform
Post by: Petr on April 30, 2020, 02:30:42 pm
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?
Title: Re: Fast Fourier Transform
Post by: _vince on April 30, 2020, 10:17:34 pm
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.
Title: Re: Fast Fourier Transform
Post by: Petr on May 01, 2020, 06:50:36 am
_Vince, thank you very much for this detailed answer, I will do some experiments to it for maximal possible understandt.
Title: Re: Fast Fourier Transform
Post by: Sanmayce on December 23, 2021, 02:17:43 am
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...
Title: Re: Fast Fourier Transform
Post by: MasterGy on December 23, 2021, 07:20:38 am
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 (http://tukorgravirozas.fw.hu/mastergy/4pohar/pg.avi)
Text Only | Text with Attachments