#! /usr/bin/gawk -f
# Last edited on 2013-01-07 23:20:54 by stolfilocal

BEGIN{

  # Generates plot data to illustrate the Lanczos filtering article
  # Client must define {rad} with "-v".

  if (rad == "") { arg_error(("must define {rad}")); }

  oprefix = "out/lanczos";
  pi = 3.1415926;
  
  xmin = 0;    # First sample to plot (int).
  xmax = 15    # Last sample to plot (int).
  ms = 5;      # Extra hiden samples at each end (int).
  kmin = xmin - ms; # First data sample (int). 
  kmax = xmax + ms; # Last data sample (int). 
  
  nh = 20;  # Subsampling ratio (int).

  pick_samples();
  write_samples();
  
  write_scaled_kernel(4,rad);
  write_scaled_kernel(11,rad)
  write_interpolant(rad);1
  
  printf "done.\n" > "/dev/stderr";
}

function pick_samples(   k)
{
  # Defines {smp[kmin..kmax]}
  srand(41700417);
  for (k = kmin; k <= kmax; k++) { smp[k] = 0.25 + rand() + rand(); }
}

function write_samples(   fname,k)
{
  # Writes to file "{oprefix}-samples.txt" the samples {smp[kmin..kmax]}.
  fname = (oprefix "-samples.txt");
  for (k = kmin; k <= kmax; k++)
    { printf "%+10.6f %+10.6f\n", k, smp[k] > fname;
    }
  close(fname);
}

function write_scaled_kernel(xs,rad,   fname,d,x,sincx,sincax,Lx)
{
  # Writes to file "{oprefix}-kern-t{xs}.txt" the Lanczos kernel
  # with size parameter {rad} scaled by smp[xs] and shifted by {xs}.
    fname = (oprefix "-r" sprintf("%02d",rad) "-kern-t" sprintf("%03d",xs) ".txt");
  for (d=-rad*nh; d <= +rad*nh; d++)
    { x = d/nh;
      Lx = lanczos(x,rad);
      printf "%+10.6f %+10.6f\n", xs + x, smp[xs]*Lx > fname;
    }
  close(fname);
}

function write_interpolant(rad,   fname,d,x,Sx)
{
  # Writes to file "{oprefix}-interp.txt" the Lanczos interpolant
  # of {smp[kmin..kmax]} using a kernel with size parameter {rad}.
  fname = (oprefix "-r" sprintf("%02d",rad) "-interp.txt");
  for (d=kmin*nh; d <= kmax*nh; d++)
    { x = d/nh;
      Sx = interpolant(x,rad);
      printf "%+10.6f %+10.6f\n", x, Sx > fname;
    }
  close(fname);
}

function lanczos(x,rad,  sx,sax)
{
  # Evaluates the Lanczos kernel with radius {rad} at {x}.
  if (x <= -rad) { return 0; }
  if (x >= +rad) { return 0; }
  sx = (x == 0 ? 1 : sin(pi*x)/(pi*x));
  sax = (x == 0 ? 1 : sin(pi*x/rad)/(pi*x/rad));
  return sax*sx;
}

function interpolant(x,rad,  xfloor,klo,khi,sum,k)
{
  # Evaluates the Lanczos interpolant
  # of {smp[kmin..kmax]} at {x}.
  xfloor = int(x + 100000) - 100000;
  klo = xfloor - rad + 1;
  khi = xfloor + rad;
  if (klo < kmin) { klo = kmin; }
  if (khi > kmax) { khi = kmax; }
  sum = 0;
  for (k = klo; k <= khi; k++)
    { sum = sum + smp[k]*lanczos(x - k,rad); }
  return sum;
}

function arg_error(msg)
{ 
  printf "** error - %s\n", msg > "/dev/stderr";
  exit 1;
}