Mercurial > hg > ltpda
view m-toolbox/test/new_ltf_code/ltf_plan.c @ 14:6d43f39633b8 database-connection-manager
Remove unused functions from utils.jmysql
author | Daniele Nicolodi <nicolodi@science.unitn.it> |
---|---|
date | Mon, 05 Dec 2011 16:20:06 +0100 |
parents | f0afece42f48 |
children |
line wrap: on
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; }