// Define functions to do Fourier transform, 
// taking care of abscissas
// and normalization of ordinates
//
// invoke this file from scilab with:  exec "fourier-norm.sci";
// invoke this file from shell  with:  scilab -f fourier-norm.sci

pi = 4*atan(1);
i = sqrt(-1);

// First, some general utility functions ...

// Calculate length of a simple vector.
// Works for 1xN row vectors as well as Nx1 column vectors.
function nnn=veclen(vv)
  nnn = size(vv,1);
  if (nnn == 1) then	// must be row vector
    nnn = size(vv,2);
  end
endfunction

function rslt=expand(vv,fctr)
  nn = veclen(vv);
  rslt = zeros(nn*fctr,1);
  rslt(1:nn) = vv;	// assign to submatrix
endfunction

function [r1,r2]=myFT(xxx,yyy);
  istep = xxx(2) - xxx(1);	// input abscissa spacing
  if (istep == 0) then 
    r1 = "bad step";
    return; 
  end;
  nnn = veclen(yyy);
  r2 = istep * fft(yyy);
//..  twister = exp(2*pi*i*xxx(1)/istep/nnn);
//..  twist = 1;
  for kk = (1:nnn)
//..    r2(kk) = r2(kk) * twist;
//..    twist = twist * twister;
//.. The following is less efficient than the above
//.. due to more transcendental functions, 
//.. but it makes sure numerical errors don't accumulate.
    r2(kk) = r2(kk) * exp(-2*pi*i*xxx(1)*(kk-1)/istep/nnn);
  end
  r1 = (0:nnn-1)'/istep/nnn;
endfunction;

// Same code as above, except ifft instead of fft,
// and ifft needs explicit factor of nnn
function [r1,r2]=myIFT(xxx,yyy);
  istep = xxx(2) - xxx(1);	// abscissa spacing
  if (istep == 0) then 
    r1 = "bad step";
    return; 
  end;
  nnn = veclen(yyy);
  r2 = nnn * istep * ifft(yyy);
  for kk = (1:nnn)
    r2(kk) = r2(kk) * exp( 2*pi*i*xxx(1)*(kk-1)/istep/nnn);
  end
  r1 = (0:nnn-1)'/istep/nnn;
endfunction;

// should produce a spike at positive frequency
function posfreq;
   nnn = 128;
   xx = 0:nnn-1
   xx = xx/nnn;
   yy = cos(3*2*pi*xx) + i*sin(3*2*pi*xx);
//   xbasc;
//   plot(xx,yy);
   [r1,r2] = myFT(xx,yy);
   xbasc;
   r1 = r1(1:10)
   r2 = r2(1:10)
   plot(r1,r2);   
endfunction;

function play;
  shift = 10;
  nnn = 128;
  xx = (shift : shift+nnn-1);   // method 1: shifted abscissas
  yy = xx*0;
  yy(1) = 1;			// populate first slot
  [f1,f2] = myIFT(xx,yy);

  xx = (0 : nnn-1);		// method 2: unshifted abscissas
  yy = xx*0;
  yy(1+shift) = 1;		// populate later slot
  [g1,g2] = myIFT(xx,yy);
  xbasc;
  plot(f1,i*f2);
  plot(g1,i*g2,'r');
  d2 = (g2 - f2);		// both methods should give same answer
  plot(g1,abs(d2)*1e14);
endfunction;


function example;
  nn = 32;
  resolution=1;
  xx = (0:nn-1)';
  yy = xx;
  for ii=1:nn
    yy(ii)  = sin(2*pi*xx(ii)/nn) + sin(9.6*2*pi*xx(ii)/nn);
  end
  
  xbasc;  xset("font size",3);
  [f1,f2] = myFT(xx,expand(yy,resolution));
  f1 = f1(1:1+resolution*nn/2);
  f2 = f2(1:1+resolution*nn/2);

  resolution=20;
  [fx1,fx2] = myFT(xx,expand(yy,resolution));
  fx1 = fx1(1:1+resolution*nn/2);
  fx2 = fx2(1:1+resolution*nn/2);

  plot(f1,abs(f2),'b-s');
  xs2ppm(0,"fourier-lores.ppm");

  xbasc;  xset("font size",3);
  plot(f1,abs(f2),'b-s');
  plot(fx1,abs(fx2),'r');
  xs2ppm(0,"fourier-2res.ppm");

  xbasc;  xset("font size",3);
  xxr = (0:resolution*nn-1)';
  xxr = xxr/resolution;
  yyr = xxr;
  for ii=1:nn*resolution
    yyr(ii)  = sin(2*pi*xxr(ii)/nn) + sin(9.6*2*pi*xxr(ii)/nn);
  end

  plot(xx,yy,'b.s');
  plot(xxr,yyr);
  xs2ppm(0,"fourier-xy.ppm");

  xbasc;  xset("font size",3);
//  plot(f1,abs(f2),'b.s');
  plot(fx1,abs(fx2),'r');
  xs2ppm(0,"fourier-hires.ppm");

  xbasc;  xset("font size",3);
  f2p=f2;
  for ii=1:veclen(f2);
    f2p(ii) = f2(ii) * conj(f2(ii));
  end;
  plot(f1,f2p,'b-s');
  fx2p = fx2;
  for ii=1:veclen(fx2);
    fx2p(ii) = fx2(ii) * conj(fx2(ii));
  end;
  plot(fx1,fx2p,'r');
  xs2ppm(0,"fourier-2res-power.ppm");

endfunction;
 
function rslt = maybe_abs(xx) ;
//  rslt = xx;
  rslt = abs(xx);
  return;
endfunction;

function radial_symmetry;
  nn = 64;
  phase = .2;
  frequency = 6;
  resolution=1;
  xx = (-nn/2:nn/2-1)';
  yy = xx;
  for ii=1:nn
    yy(ii)  = sin(phase + 2*pi*frequency*maybe_abs(xx(ii))/nn);
  end
  
  xbasc;  xset("font size",3);
  [f1,f2] = myFT(xx,expand(yy,resolution));
  f1 = f1(1:1+resolution*nn/2);
  f2 = f2(1:1+resolution*nn/2);

  resolution=20;
  [fx1,fx2] = myFT(xx,expand(yy,resolution));
  fx1 = fx1(1:1+resolution*nn/2);
  fx2 = fx2(1:1+resolution*nn/2);

  plot(f1,abs(f2),'b-s');
  xs2ppm(0,"radial_symmetry-lores.ppm");

  xbasc;  xset("font size",3);
  plot(f1,abs(f2),'b-s');
  plot(fx1,abs(fx2),'r');
  xs2ppm(0,"radial_symmetry-2res.ppm");

  xbasc;  xset("font size",3);
  xxr = (-nn*resolution/2:nn*resolution/2-1)';
  xxr = xxr/resolution;
  yyr = xxr;
  for ii=1:nn*resolution
    yyr(ii)  = sin(phase + 2*pi*frequency*maybe_abs(xxr(ii))/nn);
  end

//xx  plot(xx,yy,'b.s');        // blocks, showing lowres points
  plot(xxr,yyr);
  xs2ppm(0,"radial_symmetry-xy.ppm");

  xbasc;  xset("font size",3);
//  plot(f1,abs(f2),'b.s');
  plot(fx1,abs(fx2),'r');
  xs2ppm(0,"radial_symmetry-hires.ppm");

  xbasc;  xset("font size",3);
  f2p=f2;
  for ii=1:veclen(f2);
    f2p(ii) = f2(ii) * conj(f2(ii));
  end;
  plot(f1,f2p,'b-s');
  fx2p = fx2;
  for ii=1:veclen(fx2);
    fx2p(ii) = fx2(ii) * conj(fx2(ii));
  end;
  plot(fx1,fx2p,'r');
  xs2ppm(0,"radial_symmetry-2res-power.ppm");

endfunction;
