--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/test/new_ltf_code/ltf_plan.c	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,419 @@
+
+/*
+
+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;
+}