line source
/*
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;
}