Mercurial > hg > ltpda

comparison m-toolbox/test/new_ltf_code/ltf_plan.c @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1
2 /*
3
4 indent ltf_plan.c && \
5 gcc -g -O9 -ffast-math -Wconversion -Wall -W -o ltf_plan ltf_plan.c -lm && \
6 ./ltf_plan > plan.tex
7
8 */
9
10 #include <stdio.h>
11 #include <stdlib.h>
12 #include <math.h>
13 #include <assert.h>
14 #include <complex.h>
15
16 #define TEX 0
17
18 struct ltf_plan_one_t
19 {
20 double f; /* frequency at bin center [Hz] */
21 double fres; /* frequency resolution (bin width) [Hz] */
22 double bin; /* bin number (possibly non-integer) */
23 unsigned dftlen; /* length of one DFT segment */
24 unsigned nseg; /* number of DFT segments */
25 unsigned *segstart; /* starting position of each of the nseg segments */
26 };
27 typedef struct ltf_plan_one_t ltf_plan_one;
28
29 struct ltf_plan_t
30 {
31 /* copy of input data */
32 unsigned ndata; /* total number of input data */
33 double fsamp; /* sampling frequency [Hz] */
34 double ovlap; /* window overlap in % */
35 double binmin; /* smallest bin-number to use */
36 unsigned mindftlen; /* minimum length of one segment */
37 unsigned freqdes; /* desired number of frequencies */
38 unsigned avgdes; /* desired number of averages */
39
40 /* scheduler results */
41 unsigned nfreq; /* final number of frequencies */
42 unsigned maxdftlen; /* maximal length of one DFT segment */
43 unsigned maxnseg; /* maximal number of averages */
44 ltf_plan_one *fplan; /* detailed plan for each of nfreq frequencies */
45 };
46 typedef struct ltf_plan_t ltf_plan;
47
48 void
49 gerror (const char *err)
50 {
51 fprintf (stderr, "\n*** Error: %s\n", err);
52 exit (1);
53 }
54
55 int
56 myround (double x)
57 {
58 return ((int) (floor (x + 0.5)));
59 }
60
61
62 void
63 ltf_plan_alloc (ltf_plan * p, unsigned size)
64 {
65 char err[160];
66
67 p->fplan = realloc (p->fplan, size * sizeof (ltf_plan_one));
68
69 if (p->fplan == NULL)
70 {
71 sprintf (err,
72 "ltf_plan: error in allocating %u*%u=%u bytes for p.fplan",
73 size, sizeof (ltf_plan_one), size * sizeof (ltf_plan_one));
74 gerror (err);
75 }
76 }
77
78 ltf_plan
79 plan (unsigned ndata, double fsamp, double ovlap, double binmin,
80 unsigned mindftlen, unsigned freqdes, unsigned avgdes,
81 unsigned getsegstart)
82 {
83 ltf_plan p;
84 ltf_plan_one *p1;
85 double fmin, fmax, fres, f;
86 unsigned j, k, dftlen, segstartpos;
87 double bin;
88 unsigned nseg;
89 double xov;
90 double freslim, fresmin;
91 unsigned nalloc;
92 char err[160];
93 double logfact;
94 unsigned nsegstart; /* number of segstart entries */
95 unsigned *segstartmemory; /* begin of segstart entries */
96 double segincr;
97
98 if (ndata < 500)
99 {
100 sprintf (err, "ltf_plan: illegal value for ndata: %u<500", ndata);
101 gerror (err);
102 }
103
104 if (ndata < 500)
105 {
106 sprintf (err, "ltf_plan: illegal value for ndata: %u<500", ndata);
107 gerror (err);
108 }
109
110 if (fsamp <= 0)
111 {
112 sprintf (err, "ltf_plan: illegal value for fsamp: %g<=0", fsamp);
113 gerror (err);
114 }
115
116 if (binmin < 0)
117 {
118 sprintf (err, "ltf_plan: illegal value for binmin: %g<0", binmin);
119 gerror (err);
120 }
121
122 if (ovlap < 0)
123 {
124 sprintf (err, "ltf_plan: illegal value for ovlap: %g<=0", ovlap);
125 gerror (err);
126 }
127
128 if (ovlap > 99)
129 {
130 sprintf (err, "ltf_plan: illegal value for ovlap: %g>99%%", ovlap);
131 gerror (err);
132 }
133
134 if (freqdes < 3)
135 {
136 sprintf (err, "ltf_plan: illegal value for freqdes: %u<3", freqdes);
137 gerror (err);
138 }
139
140 if (freqdes > ndata / 2)
141 {
142 sprintf (err, "ltf_plan: illegal value for freqdes: %u>%u", freqdes,
143 ndata / 2);
144 gerror (err);
145 }
146
147 /* copy input data */
148 p.ndata = ndata;
149 p.fsamp = fsamp;
150 p.ovlap = ovlap;
151 if (binmin < 1)
152 binmin = 1;
153 p.binmin = binmin;
154 p.mindftlen = mindftlen;
155 p.freqdes = freqdes;
156 p.avgdes = avgdes;
157 p.maxdftlen = p.maxnseg = 0;
158 nsegstart = 0;
159
160 xov = (1. - ovlap / 100.);
161 fmin = fsamp / ndata * binmin;
162 fmax = fsamp / 2;
163 fresmin = fsamp / (double) ndata; /* absolute minimal allowed frequency resolution */
164 freslim = fresmin * (1 + xov * (avgdes - 1)); /* limit above which fres leads to desired number of averages */
165 logfact = (pow (ndata / 2., 1. / freqdes) - 1.);
166
167 /* initial allocation */
168 p.fplan = NULL;
169 nalloc = freqdes;
170 ltf_plan_alloc (&p, nalloc);
171
172 for (j = 0, f = fmin; f < fmax; j++)
173 {
174 if (j >= nalloc)
175 {
176 /* enlarge allocation */
177 nalloc *= 2;
178 ltf_plan_alloc (&p, nalloc);
179 }
180
181 p1 = &(p.fplan[j]);
182
183 fres = f * logfact;
184 fres = (fres > freslim) ? fres : sqrt (fres * freslim);
185 if (fres < fresmin)
186 fres = fresmin;
187
188 bin = (f / fres);
189 if (bin < binmin)
190 {
191 bin = binmin;
192 fres = f / bin;
193 }
194
195 dftlen = myround (fsamp / fres);
196 if (dftlen > ndata)
197 dftlen = ndata;
198 if (dftlen < mindftlen)
199 dftlen = mindftlen;
200 nseg = myround (((double) (ndata - dftlen) /
201 (xov * (double) dftlen)) + 1);
202 if (nseg == 1)
203 dftlen = ndata;
204
205 fres = fsamp / dftlen;
206 bin = f / fres;
207
208 p1->f = f;
209 p1->fres = fres;
210 p1->bin = bin;
211 p1->dftlen = dftlen;
212 p1->nseg = nseg;
213 if (nseg > p.maxnseg)
214 p.maxnseg = nseg;
215 nsegstart += nseg;
216 if (dftlen > p.maxdftlen)
217 p.maxdftlen = dftlen;
218 f += fres;
219 }
220 p.nfreq = j;
221 /* allocation shrink to fit */
222 nalloc = p.nfreq;
223 ltf_plan_alloc (&p, nalloc);
224
225 /* determine starting indices of DFT segments */
226 if (!getsegstart)
227 {
228 for (j = 0; j < p.nfreq; j++)
229 p.fplan[j].segstart = NULL;
230 }
231 else
232 {
233 segstartmemory = malloc (nsegstart * sizeof (unsigned));
234 if (segstartmemory == NULL)
235 {
236 sprintf (err,
237 "ltf_plan: error allocating %u bytes for segment start indices",
238 nsegstart * sizeof (unsigned));
239 gerror (err);
240 }
241 segstartpos = 0;
242 for (j = 0; j < p.nfreq; j++)
243 {
244 p1 = &(p.fplan[j]);
245 p1->segstart = &(segstartmemory[segstartpos]);
246 p1->segstart[0] = 0;
247 if (p1->nseg > 1)
248 {
249 segincr =
250 (double) (ndata - p1->dftlen) / (double) (p1->nseg - 1);
251 for (k = 1; k < p1->nseg - 1; k++)
252 p1->segstart[k] = myround ((double) k * segincr);
253 p1->segstart[p1->nseg - 1] = ndata - p1->dftlen;
254 }
255 segstartpos += p1->nseg;
256 }
257 }
258 return p;
259 }
260
261 void
262 ltf_plan_free (ltf_plan p)
263 {
264 if ((p.fplan[0]).segstart != NULL)
265 free ((p.fplan[0]).segstart);
266 free (p.fplan);
267 }
268
269 void
270 fprint (double f, char *s)
271 {
272 if (f < 1e-5)
273 sprintf (s, "%.3f$\\mu$", 1e6 * f);
274 else if (f < 1e-4)
275 sprintf (s, "%.2f$\\mu$", 1e6 * f);
276 else if (f < 1e-3)
277 sprintf (s, "%.1f$\\mu$", 1e6 * f);
278 else if (f < 1e-2)
279 sprintf (s, "%.3fm", 1e3 * f);
280 else if (f < 1e-1)
281 sprintf (s, "%.2fm", 1e3 * f);
282 else if (f < 1)
283 sprintf (s, "%.1fm", 1e3 * f);
284 else
285 sprintf (s, "%.3f", f);
286 }
287
288 int
289 main (void)
290 {
291 ltf_plan p;
292 ltf_plan_one *p1;
293 unsigned ndata = 1000000;
294 double fsamp = 10.;
295 double ovlap = 75;
296 double binmin = 5.3;
297 unsigned mindftlen = 100;
298 unsigned freqdes = 200;
299 unsigned avgdes = 100;
300 unsigned j, k;
301 double eq1, eq2;
302 double dfdown, dfup;
303 double logfact, logres;
304
305 logfact = (pow (ndata / 2., 1. / freqdes) - 1.);
306
307 p = plan (ndata, fsamp, ovlap, binmin, mindftlen, freqdes, avgdes, 1);
308
309 if (TEX)
310 {
311
312 printf ("\\section{Sample output}\n");
313 printf
314 ("As reference for other implementations, here are the results for the case\\\\\n");
315 printf ("$N_{\\rm data}$ ({\\tt ndata}) $=%d$\\\\\n", ndata);
316 printf ("$f_{\\rm samp}$ ({\\tt fsamp}) $=%g$\\,Hz\\\\\n", fsamp);
317 printf ("$\\xi_{\\rm ovlap}$ ({\\tt ovlap}) $=%g$\\,\\%%\\\\\n", ovlap);
318 printf ("$b_{\\rm min}$ ({\\tt binmin}) $=%g$\\\\\n", binmin);
319 printf ("$L_{\\rm min}$ ({\\tt mindftlen}) $=%d$\\\\\n", mindftlen);
320 printf ("$J_{\\rm des}$ ({\\tt freqdes}) $=%d$\\\\\n", freqdes);
321 printf ("$K_{\\rm des}$ ({\\tt avgdes}) $=%d$\n", avgdes);
322 printf ("\n");
323 printf
324 ("Usually, $J_{\\rm des}$ ({\\tt freqdes}) will be chosen higher (1000\
325 is a typical value); %d has been chosen here to save space in the\
326 printed table.\n",
327 freqdes);
328 printf ("\n");
329 printf
330 ("Targets that are not met are printed {\\bfseries \\slshape like this}\n");
331 printf ("\n");
332
333 printf ("\\small{\n");
334 printf
335 ("\\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");
336 printf ("\\begin{supertabular}{|r|r|r|r|r|r|r|r|r|l|}\n");
337 }
338
339 for (j = 0; j < p.nfreq; j++)
340 {
341 p1 = &(p.fplan[j]);
342
343 /* check required relations */
344 eq1 = p1->f - p1->fres * p1->bin;
345 assert (fabs (eq1 / p1->f) < 1e-10);
346 eq2 = p.fsamp - p1->fres * p1->dftlen;
347 assert (fabs (eq2 / p.fsamp) < 1e-10);
348
349 if (j > 0)
350 dfdown = p.fplan[j].f - p.fplan[j - 1].f;
351 else
352 dfdown = 0;
353 if (j < p.nfreq - 1)
354 dfup = p.fplan[j + 1].f - p.fplan[j].f;
355 else
356 dfup = 0;
357 logres = p1->fres / (p1->f * logfact);
358
359 if (TEX)
360 {
361 char sf[80], sfres[80], snseg[80], sdown[80], sup[80], slogres[80];
362 fprint (p1->f, sf);
363 fprint (p1->fres, sfres);
364 if (p1->nseg < avgdes)
365 sprintf (snseg, "{\\bfseries \\slshape %7d}", p1->nseg);
366 else
367 sprintf (snseg, "%7d", p1->nseg);
368 if (j > 0)
369 sprintf (sdown, "%.3f", dfdown / p1->fres);
370 else
371 sprintf (sdown, "---");
372 if (j < p.nfreq - 1)
373 sprintf (sup, "%.3f", dfup / p1->fres);
374 else
375 sprintf (sup, "---");
376 if (logres < .995 || logres > 1.005)
377 sprintf (slogres, "{\\bfseries \\slshape %.3f}", logres);
378 else
379 sprintf (slogres, "%.3f", logres);
380 printf
381 ("%4d & %s & %s & %.2f & %s & %7d & %s & %s & %s & ",
382 j, sf, sfres, p1->bin, snseg, p1->dftlen, sdown, sup, slogres);
383 }
384 else
385 {
386 printf ("%4d %20.20f %20.20f %20.20f %20d %20d %20.20f %20.20f %20.20f ", j,
387 p1->f, p1->fres, p1->bin, p1->nseg, p1->dftlen,
388 dfdown / p1->fres, dfup / p1->fres, logres);
389 }
390 if (p1->segstart != NULL)
391 {
392 // if (p1->nseg < 5)
393 // {
394 // for (k = 0; k < p1->nseg; k++)
395 // printf ("%u ", p1->segstart[k]);
396 // printf ("\\\\\n");
397 // }
398 // else
399 // {
400 // printf ("%u %u %u ... %u\\\\\n",
401 // p1->segstart[0],
402 // p1->segstart[1], p1->segstart[2],
403 // p1->segstart[p1->nseg - 1]);
404 // }
405 printf("\n");
406 assert (p1->segstart[p1->nseg - 1] + p1->dftlen == ndata);
407 }
408 else
409 printf (" (NULL)\n");
410 }
411 if (TEX)
412 {
413 printf ("\\end{supertabular}}\n");
414 }
415 // else
416 // printf ("% max. %u segments %u dftlen\n", p.maxnseg, p.maxdftlen);
417 ltf_plan_free (p);
418 return 0;
419 }