share/numeric/fft.dem

   1 if not ?fboundp (fft) then load (fft);
   2
   3 "The fft of a list:"$
   4 (N:32,
   5   g:makelist(if i<N/2 then 1.0 else 0.0,i,0,N-1),
   6   plot2d([discrete,makelist([i,g[i]],i,1,N)],['y,-0.1,1.1]),
   7   gt:fft(g));
   8
   9 "The inverse_fft of a list:"$
  10 gti:inverse_fft(gt);
  11 "The norm-difference between list G and output of inverse_fft"$
  12 apply(max,map(cabs,gti-g));
  13
  14 "The fft of an array:"$
  15 (N:8,
  16   g:g+%i*g,
  17   ga:make_array(any,N), fillarray(ga,g),
  18   gta:fft(ga));
  19
  20 "The inverse_fft of an array:"$
  21 gtia:inverse_fft(gta);
  22 "The norm-difference between array GA and output of inverse_fft"$
  23 (mx:0, for i:0 thru N-1 do mx:max(mx,cabs(gtia[i]-ga[i])), mx);
  24
  25 "Example from info pages:"$
  26 L : [1, 2, 3, 4, 5, 6, 7, 8];
  27 (n : length (L),
  28   x : make_array (any, n),
  29   fillarray (x, L) );
  30 y : fft (x) ;
  31 (a : make_array (any, n/2 + 1) ,
  32   b : make_array (any, n/2 + 1) ,
  33   a[0] : realpart (y[0]) ,
  34   b[0] : 0 ,
  35   for k : 1 thru n/2 - 1 do
  36   (a[k] : realpart (y[k] + y[n - k]),
  37     b[k] : imagpart (y[n - k] - y[k])),
  38   a[n/2] : y[n/2] ,
  39   b[n/2] : 0 )$
  40 listarray (a);
  41 listarray (b);
  42 f(j) := sum (a[k]*cos(-2*%pi*j*k/n) + b[k]*sin(-2*%pi*j*k/n), k, 0, n/2) $
  43 "The following should equal the list L"$
  44 makelist (float(f (j)), j, 0, n - 1);