/*

indent ltf_plan.c && \
gcc -g -O9 -ffast-math -Wconversion -Wall -W -o ltf_plan ltf_plan.c -lm && \
./ltf_plan > plan.tex

*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include <complex.h>

#define TEX 0

struct ltf_plan_one_t
{
  double f;			/* frequency at bin center [Hz] */
  double fres;			/* frequency resolution (bin width) [Hz] */
  double bin;			/* bin number (possibly non-integer) */
  unsigned dftlen;		/* length of one DFT segment */
  unsigned nseg;		/* number of DFT segments */
  unsigned *segstart;		/* starting position of each of the nseg segments */
};
typedef struct ltf_plan_one_t ltf_plan_one;

struct ltf_plan_t
{
/* copy of input data */
  unsigned ndata;		/* total number of input data */
  double fsamp;			/* sampling frequency [Hz] */
  double ovlap;			/* window overlap in % */
  double binmin;		/* smallest bin-number to use */
  unsigned mindftlen;		/* minimum length of one segment */
  unsigned freqdes;		/* desired number of frequencies */
  unsigned avgdes;		/* desired number of averages */

/* scheduler results */
  unsigned nfreq;		/* final number of frequencies */
  unsigned maxdftlen;		/* maximal length of one DFT segment */
  unsigned maxnseg;		/* maximal number of averages */
  ltf_plan_one *fplan;		/* detailed plan for each of nfreq frequencies */
};
typedef struct ltf_plan_t ltf_plan;

void
gerror (const char *err)
{
  fprintf (stderr, "\n*** Error: %s\n", err);
  exit (1);
}

int
myround (double x)
{
  return ((int) (floor (x + 0.5)));
}


void
ltf_plan_alloc (ltf_plan * p, unsigned size)
{
  char err[160];

  p->fplan = realloc (p->fplan, size * sizeof (ltf_plan_one));

  if (p->fplan == NULL)
    {
      sprintf (err,
	       "ltf_plan: error in allocating %u*%u=%u bytes for p.fplan",
	       size, sizeof (ltf_plan_one), size * sizeof (ltf_plan_one));
      gerror (err);
    }
}

ltf_plan
plan (unsigned ndata, double fsamp, double ovlap, double binmin,
      unsigned mindftlen, unsigned freqdes, unsigned avgdes,
      unsigned getsegstart)
{
  ltf_plan p;
  ltf_plan_one *p1;
  double fmin, fmax, fres, f;
  unsigned j, k, dftlen, segstartpos;
  double bin;
  unsigned nseg;
  double xov;
  double freslim, fresmin;
  unsigned nalloc;
  char err[160];
  double logfact;
  unsigned nsegstart;		/* number of segstart entries */
  unsigned *segstartmemory;	/* begin of segstart entries */
  double segincr;

  if (ndata < 500)
    {
      sprintf (err, "ltf_plan: illegal value for ndata: %u<500", ndata);
      gerror (err);
    }

  if (ndata < 500)
    {
      sprintf (err, "ltf_plan: illegal value for ndata: %u<500", ndata);
      gerror (err);
    }

  if (fsamp <= 0)
    {
      sprintf (err, "ltf_plan: illegal value for fsamp: %g<=0", fsamp);
      gerror (err);
    }

  if (binmin < 0)
    {
      sprintf (err, "ltf_plan: illegal value for binmin: %g<0", binmin);
      gerror (err);
    }

  if (ovlap < 0)
    {
      sprintf (err, "ltf_plan: illegal value for ovlap: %g<=0", ovlap);
      gerror (err);
    }

  if (ovlap > 99)
    {
      sprintf (err, "ltf_plan: illegal value for ovlap: %g>99%%", ovlap);
      gerror (err);
    }

  if (freqdes < 3)
    {
      sprintf (err, "ltf_plan: illegal value for freqdes: %u<3", freqdes);
      gerror (err);
    }

  if (freqdes > ndata / 2)
    {
      sprintf (err, "ltf_plan: illegal value for freqdes: %u>%u", freqdes,
	       ndata / 2);
      gerror (err);
    }

  /* copy input data */
  p.ndata = ndata;
  p.fsamp = fsamp;
  p.ovlap = ovlap;
  if (binmin < 1)
    binmin = 1;
  p.binmin = binmin;
  p.mindftlen = mindftlen;
  p.freqdes = freqdes;
  p.avgdes = avgdes;
  p.maxdftlen = p.maxnseg = 0;
  nsegstart = 0;

  xov = (1. - ovlap / 100.);
  fmin = fsamp / ndata * binmin;
  fmax = fsamp / 2;
  fresmin = fsamp / (double) ndata;	/* absolute minimal allowed frequency resolution */
  freslim = fresmin * (1 + xov * (avgdes - 1));	/* limit above which fres leads to desired number of averages */
  logfact = (pow (ndata / 2., 1. / freqdes) - 1.);

  /* initial allocation */
  p.fplan = NULL;
  nalloc = freqdes;
  ltf_plan_alloc (&p, nalloc);

  for (j = 0, f = fmin; f < fmax; j++)
    {
      if (j >= nalloc)
	{
	  /* enlarge allocation */
	  nalloc *= 2;
	  ltf_plan_alloc (&p, nalloc);
	}

      p1 = &(p.fplan[j]);

      fres = f * logfact;
      fres = (fres > freslim) ? fres : sqrt (fres * freslim);
      if (fres < fresmin)
	fres = fresmin;

      bin = (f / fres);
      if (bin < binmin)
	{
	  bin = binmin;
	  fres = f / bin;
	}

      dftlen = myround (fsamp / fres);
      if (dftlen > ndata)
	dftlen = ndata;
      if (dftlen < mindftlen)
	dftlen = mindftlen;
      nseg = myround (((double) (ndata - dftlen) /
		       (xov * (double) dftlen)) + 1);
      if (nseg == 1)
	dftlen = ndata;

      fres = fsamp / dftlen;
      bin = f / fres;

      p1->f = f;
      p1->fres = fres;
      p1->bin = bin;
      p1->dftlen = dftlen;
      p1->nseg = nseg;
      if (nseg > p.maxnseg)
	p.maxnseg = nseg;
      nsegstart += nseg;
      if (dftlen > p.maxdftlen)
	p.maxdftlen = dftlen;
      f += fres;
    }
  p.nfreq = j;
  /* allocation shrink to fit  */
  nalloc = p.nfreq;
  ltf_plan_alloc (&p, nalloc);

  /* determine starting indices of DFT segments */
  if (!getsegstart)
    {
      for (j = 0; j < p.nfreq; j++)
	p.fplan[j].segstart = NULL;
    }
  else
    {
      segstartmemory = malloc (nsegstart * sizeof (unsigned));
      if (segstartmemory == NULL)
	{
	  sprintf (err,
		   "ltf_plan: error allocating %u bytes for segment start indices",
		   nsegstart * sizeof (unsigned));
	  gerror (err);
	}
      segstartpos = 0;
      for (j = 0; j < p.nfreq; j++)
	{
	  p1 = &(p.fplan[j]);
	  p1->segstart = &(segstartmemory[segstartpos]);
	  p1->segstart[0] = 0;
	  if (p1->nseg > 1)
	    {
	      segincr =
		(double) (ndata - p1->dftlen) / (double) (p1->nseg - 1);
	      for (k = 1; k < p1->nseg - 1; k++)
		p1->segstart[k] = myround ((double) k * segincr);
	      p1->segstart[p1->nseg - 1] = ndata - p1->dftlen;
	    }
	  segstartpos += p1->nseg;
	}
    }
  return p;
}

void
ltf_plan_free (ltf_plan p)
{
  if ((p.fplan[0]).segstart != NULL)
    free ((p.fplan[0]).segstart);
  free (p.fplan);
}

void
fprint (double f, char *s)
{
  if (f < 1e-5)
    sprintf (s, "%.3f$\\mu$", 1e6 * f);
  else if (f < 1e-4)
    sprintf (s, "%.2f$\\mu$", 1e6 * f);
  else if (f < 1e-3)
    sprintf (s, "%.1f$\\mu$", 1e6 * f);
  else if (f < 1e-2)
    sprintf (s, "%.3fm", 1e3 * f);
  else if (f < 1e-1)
    sprintf (s, "%.2fm", 1e3 * f);
  else if (f < 1)
    sprintf (s, "%.1fm", 1e3 * f);
  else
    sprintf (s, "%.3f", f);
}

int
main (void)
{
  ltf_plan p;
  ltf_plan_one *p1;
  unsigned ndata = 1000000;
  double fsamp = 10.;
  double ovlap = 75;
  double binmin = 5.3;
  unsigned mindftlen = 100;
  unsigned freqdes = 200;
  unsigned avgdes = 100;
  unsigned j, k;
  double eq1, eq2;
  double dfdown, dfup;
  double logfact, logres;

  logfact = (pow (ndata / 2., 1. / freqdes) - 1.);

  p = plan (ndata, fsamp, ovlap, binmin, mindftlen, freqdes, avgdes, 1);

  if (TEX)
    {

      printf ("\\section{Sample output}\n");
      printf
	("As reference for other implementations, here are the results for the case\\\\\n");
      printf ("$N_{\\rm data}$ ({\\tt ndata}) $=%d$\\\\\n", ndata);
      printf ("$f_{\\rm samp}$ ({\\tt fsamp}) $=%g$\\,Hz\\\\\n", fsamp);
      printf ("$\\xi_{\\rm ovlap}$ ({\\tt ovlap}) $=%g$\\,\\%%\\\\\n", ovlap);
      printf ("$b_{\\rm min}$ ({\\tt binmin}) $=%g$\\\\\n", binmin);
      printf ("$L_{\\rm min}$ ({\\tt mindftlen}) $=%d$\\\\\n", mindftlen);
      printf ("$J_{\\rm des}$ ({\\tt freqdes}) $=%d$\\\\\n", freqdes);
      printf ("$K_{\\rm des}$ ({\\tt avgdes}) $=%d$\n", avgdes);
      printf ("\n");
      printf
	("Usually, $J_{\\rm des}$ ({\\tt freqdes}) will be chosen higher (1000\
	is a typical value); %d has been chosen here to save space in the\
	printed	table.\n",
	 freqdes);
      printf ("\n");
      printf
	("Targets that are not met are printed {\\bfseries \\slshape like this}\n");
      printf ("\n");

      printf ("\\small{\n");
      printf
	("\\tablehead{$j$ & $f_j$ & $r_j$ & $b_j$ & $K_j$ & $L_j$ & $\\frac{f_j-f_{j-1}}{r_j}$ & $\\frac{f_{j+1}-f_j}{r_j}$ & $\\frac{r_j}{c\\cdot f_j}$ & $q_k^{(j)}$\\\\\\hline}\n");
      printf ("\\begin{supertabular}{|r|r|r|r|r|r|r|r|r|l|}\n");
    }

  for (j = 0; j < p.nfreq; j++)
    {
      p1 = &(p.fplan[j]);

      /* check required relations */
      eq1 = p1->f - p1->fres * p1->bin;
      assert (fabs (eq1 / p1->f) < 1e-10);
      eq2 = p.fsamp - p1->fres * p1->dftlen;
      assert (fabs (eq2 / p.fsamp) < 1e-10);

      if (j > 0)
	dfdown = p.fplan[j].f - p.fplan[j - 1].f;
      else
	dfdown = 0;
      if (j < p.nfreq - 1)
	dfup = p.fplan[j + 1].f - p.fplan[j].f;
      else
	dfup = 0;
      logres = p1->fres / (p1->f * logfact);

      if (TEX)
	{
	  char sf[80], sfres[80], snseg[80], sdown[80], sup[80], slogres[80];
	  fprint (p1->f, sf);
	  fprint (p1->fres, sfres);
	  if (p1->nseg < avgdes)
	    sprintf (snseg, "{\\bfseries \\slshape %7d}", p1->nseg);
	  else
	    sprintf (snseg, "%7d", p1->nseg);
	  if (j > 0)
	    sprintf (sdown, "%.3f", dfdown / p1->fres);
	  else
	    sprintf (sdown, "---");
	  if (j < p.nfreq - 1)
	    sprintf (sup, "%.3f", dfup / p1->fres);
	  else
	    sprintf (sup, "---");
	  if (logres < .995 || logres > 1.005)
	    sprintf (slogres, "{\\bfseries \\slshape %.3f}", logres);
	  else
	    sprintf (slogres, "%.3f", logres);
	  printf
	    ("%4d & %s & %s & %.2f & %s & %7d & %s & %s & %s & ",
	     j, sf, sfres, p1->bin, snseg, p1->dftlen, sdown, sup, slogres);
	}
      else
	{
	  printf ("%4d %20.20f %20.20f %20.20f %20d %20d %20.20f %20.20f %20.20f ", j,
		  p1->f, p1->fres, p1->bin, p1->nseg, p1->dftlen,
		  dfdown / p1->fres, dfup / p1->fres, logres);
	}
      if (p1->segstart != NULL)
	{
// 	  if (p1->nseg < 5)
// 	    {
// 	      for (k = 0; k < p1->nseg; k++)
// 		printf ("%u ", p1->segstart[k]);
// 	      printf ("\\\\\n");
// 	    }
// 	  else
// 	    {
// 	      printf ("%u %u %u ... %u\\\\\n",
// 		      p1->segstart[0],
// 		      p1->segstart[1], p1->segstart[2],
// 		      p1->segstart[p1->nseg - 1]);
// 	    }
    printf("\n");
	  assert (p1->segstart[p1->nseg - 1] + p1->dftlen == ndata);
	}
      else
	printf (" (NULL)\n");
    }
  if (TEX)
    {
      printf ("\\end{supertabular}}\n");
    }
//   else
//     printf ("% max. %u segments %u dftlen\n", p.maxnseg, p.maxdftlen);
  ltf_plan_free (p);
  return 0;
}