
 -- Fast Fourier Transform Examples / fft.lua
 -- 
 -- Copyright (C) 2009, 2010 Francesco Abbate
 -- 
 -- This program is free software; you can redistribute it and/or modify
 -- it under the terms of the GNU General Public License as published by
 -- the Free Software Foundation; either version 3 of the License, or (at
 -- your option) any later version.
 -- 
 -- This program is distributed in the hope that it will be useful, but
 -- WITHOUT ANY WARRANTY; without even the implied warranty of
 -- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 -- General Public License for more details.
 -- 
 -- You should have received a copy of the GNU General Public License
 -- along with this program; if not, write to the Free Software
 -- Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.

function demo1()
   local n, ncut = 256, 24

   local NS = 5
   local sqa = new(n, NS, |i,j| i < n/(j+1) and 0 or (i < 2*n/(j+1) and 1 or 0))

   for j= 1, NS do
      local sq = sqa:col(j)

      local pt = plot('Original signal / reconstructed')
      local pf = plot('FFT Power Spectrum')

      pt:addline(filine(|i| sq[i], n), 'black')

      fft(sq)

      pf:add(ibars(isample(|k| abs(sq:get(k)), 0, 60)), 'black')

      for k=ncut, n/2 do sq:set(k,0) end
      fft_inv(sq)

      pt:addline(filine(|i| sq[i], n), 'red')

      pf:show()
      pt:show()
   end
end

function demo2()
   local n, ncut, order = 512, 11, 8
   local x1 = besselJ_zero(order, 14)
   local xsmp = |k| x1*(k-1)/(n-1)

   local bess = new(n, 1, |i| besselJ(order, xsmp(i)))

   local p = plot('Original signal / reconstructed')
   p:addline(filine(|i| bess[i], n), 'black')

   fft(bess)

   fftplot = plot('FFT power spectrum')
   bars = ibars(isample(|k| abs(bess:get(k)), 0, 60))
   fftplot:add(bars, 'darkgreen')
   fftplot:addline(bars, 'black')
   fftplot:show()

   for k=ncut, n/2 do bess:set(k,0) end
   fft_inv(bess)

   p:addline(filine(|i| bess[i], n), 'red', {{'dash', 7, 3}})
   p:show()

   return p, fftplot
end

echo 'demo1() - GSL example with square function and frequency cutoff'
echo 'demo2() - frequency cutoff example on bessel function'
