Bug Summary

File:src/../tools/meteor-engine.c
Warning:line 811, column 5
Function call argument is an uninitialized value

Annotated Source Code

1/*
2 * SpanDSP - a series of DSP components for telephony
3 *
4 * meteor-engine.c - The meteor FIR design algorithm
5 *
6 * Written by Steve Underwood <steveu@coppice.org>
7 *
8 * Copyright (C) 2013 Steve Underwood
9 *
10 * All rights reserved.
11 *
12 * This program is free software; you can redistribute it and/or modify
13 * it under the terms of the GNU General Public License version 2, as
14 * published by the Free Software Foundation.
15 *
16 * This program is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 * GNU General Public License for more details.
20 *
21 * You should have received a copy of the GNU General Public License
22 * along with this program; if not, write to the Free Software
23 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
24 */
25
26/*! \file */
27
28#include <inttypes.h>
29#include <stdbool.h>
30#include <time.h>
31#include <stdio.h>
32#include <string.h>
33#include <ctype.h>
34#include <math.h>
35#include <setjmp.h>
36#include <assert.h>
37#include <stddef.h>
38#include <stdlib.h>
39
40#include "meteor-engine.h"
41#include "ae.h"
42
43/* version I: Wed Jun 27 13:36:06 EDT 1990 */
44
45/* copyright Prof. K. Steiglitz
46 Dept. of Computer Science
47 Princeton University
48 Princeton, NJ 08544 */
49
50/* Constraint-based design of linear-phase fir filters with
51 upper and lower bounds, and convexity constraints.
52 Finds minimum length, or optimizes fixed length, or pushes band-edges.
53 If L is the filter length, the models are
54
55 odd-length
56 cosine: sum(i from 0 to (L-1)/2) coeff[i]*cos(i*omega)
57 sine: sum(i from 0 to (L-3)/2) coeff[i]*sin((i + 1)*omega)
58
59 even-length
60 cosine: sum(i from 0 to L/2 - 1) coeff[i]*cos((i + 0.5)*omega)
61 sine: sum(i from 0 to L/2 - 1) coeff[i]*sin((i + 0.5)*omega) */
62
63#define MAX_PIVOTS1000 1000 /* Maximum no. of pivots */
64#define SMALL1.0e-8 1.0e-8 /* Small number used in defining band-edges */
65#define LARGE1.0e+31 1.0e+31 /* Large number used in search for minimum cost column */
66#define EPS1.0e-8 1.0e-8 /* For testing for zero */
67
68static void make_bands(struct meteor_working_data_s *s, int i)
69{
70 int j;
71 int kmax;
72
73 /* Fill in frequencies to make grid - frequencies are kept as reals
74 in radians, and each band has equally spaced grid points */
75 if (i == 0)
76 s->spec->spec[i].first_col = 1;
77 else
78 s->spec->spec[i].first_col = s->spec->spec[i - 1].last_col + 1;
79 /*endif*/
80 kmax = (int) ((s->spec->spec[i].right_freq - s->spec->spec[i].left_freq)*s->spec->grid_points/0.5 + SMALL1.0e-8);
81 /* kmax + 1 cols. in this band */
82 if (kmax == 0)
83 {
84 s->freq[s->spec->spec[i].first_col - 1] = 2.0*M_PI3.14159265358979323846*s->spec->spec[i].left_freq;
85 }
86 else
87 {
88 for (j = 0; j <= kmax; j++)
89 s->freq[s->spec->spec[i].first_col + j - 1] = 2.0*M_PI3.14159265358979323846*(s->spec->spec[i].left_freq + (s->spec->spec[i].right_freq - s->spec->spec[i].left_freq)*j/kmax);
90 /*endfor*/
91 }
92 /*endif*/
93 s->spec->spec[i].last_col = s->spec->spec[i].first_col + kmax;
94}
95/*- End of function --------------------------------------------------------*/
96
97static double trig0(struct meteor_working_data_s *s, int i, double freq)
98{
99 double res;
100
101 /* Trig function in filter transfer function */
102 if (s->odd_length)
103 {
104 if (s->spec->symmetry_type == symmetry_cosine)
105 res = cos(i*freq);
106 else
107 res = sin((i + 1.0)*freq);
108 /*endif*/
109 }
110 else
111 {
112 if (s->spec->symmetry_type == symmetry_cosine)
113 res = cos((i + 0.5)*freq);
114 else
115 res = sin((i + 0.5)*freq);
116 /*endif*/
117 }
118 /*endif*/
119 return res;
120}
121/*- End of function --------------------------------------------------------*/
122
123static double trig2(struct meteor_working_data_s *s, int i, double freq)
124{
125 double res;
126
127 /* Second derivative of trig function in filter transfer function */
128 if (s->odd_length)
129 {
130 if (s->spec->symmetry_type == symmetry_cosine)
131 res = -i*i*cos(i*freq);
132 else
133 res = -(i + 1.0)*(i + 1.0)*sin(i*freq);
134 /*endif*/
135 }
136 else
137 {
138 if (s->spec->symmetry_type == symmetry_cosine)
139 res = -(i + 0.5)*(i + 0.5)*cos((i + 0.5)*freq);
140 else
141 res = -(i + 0.5)*(i + 0.5)*sin((i + 0.5)*freq);
142 /*endif*/
143 }
144 /*endif*/
145 return res;
146}
147/*- End of function --------------------------------------------------------*/
148
149static void convex(struct meteor_working_data_s *s, int i)
150{
151 int row;
152 int col;
153
154 /* Set up tableau columns for convexity constraints on magnitude */
155 make_bands(s, i);
156 for (col = s->spec->spec[i].first_col - 1; col < s->spec->spec[i].last_col; col++)
157 {
158 /* For all frequencies in band */
159 s->c[col] = 0.0;
160 for (row = 0; row < s->m; row++)
161 {
162 /* Normal constraint is <= */
163 if (s->spec->spec[i].sense == sense_convex)
164 s->tab[row][col] = -s->tab[row][col];
165 else
166 s->tab[row][col] = trig2(s, row, s->freq[col]);
167 /*endif*/
168 }
169 /*endfor*/
170 s->tab[s->m][col] = 0.0;
171 }
172 /*endfor*/
173}
174/*- End of function --------------------------------------------------------*/
175
176static void limit(struct meteor_working_data_s *s, int i)
177{
178 int row;
179 int col;
180
181 /* Sets up tableau columns for upper or lower bounds on transfer
182 function for specification i; the bound is linearly interpolated
183 between the start and end of the band */
184 make_bands(s, i);
185 for (col = s->spec->spec[i].first_col - 1; col < s->spec->spec[i].last_col; col++)
186 {
187 /* For all frequencies in band */
188 if (s->spec->spec[i].first_col == s->spec->spec[i].last_col)
189 {
190 s->c[col] = s->spec->spec[i].left_bound;
191 }
192 else
193 {
194 switch (s->spec->spec[i].interpolation)
195 {
196 case interpolation_geometric:
197 s->c[col] = s->spec->spec[i].left_bound*exp((double) (col - s->spec->spec[i].first_col + 1) /
198 (s->spec->spec[i].last_col - s->spec->spec[i].first_col)*log(fabs(s->spec->spec[i].right_bound/s->spec->spec[i].left_bound)));
199 break;
200 case interpolation_arithmetic:
201 s->c[col] = s->spec->spec[i].left_bound + (double) (col - s->spec->spec[i].first_col + 1) /
202 (s->spec->spec[i].last_col - s->spec->spec[i].first_col)*(s->spec->spec[i].right_bound - s->spec->spec[i].left_bound);
203 break;
204 }
205 /*endswitch*/
206 }
207 /*endif*/
208 if (s->spec->spec[i].sense == sense_lower)
209 s->c[col] = -s->c[col];
210 /*endif*/
211 for (row = 0; row < s->m; row++)
212 {
213 s->tab[row][col] = (s->spec->spec[i].sense == sense_lower) ? -trig0(s, row, s->freq[col]) : trig0(s, row, s->freq[col]);
214 }
215 /*endfor*/
216 s->tab[s->m][col] = (s->spec->spec[i].hug) ? 0.0 : 1.0;
217 }
218 /*endfor*/
219}
220/*- End of function --------------------------------------------------------*/
221
222static void setup(struct meteor_working_data_s *s)
223{
224 int i;
225
226 /* Initialize constraints */
227 for (i = 0; i < s->spec->num_specs; i++)
228 {
229 switch (s->spec->spec[i].type)
230 {
231 case constraint_type_convexity:
232 convex(s, i);
233 break;
234 case constraint_type_limit:
235 limit(s, i);
236 break;
237 }
238 /*endswitch*/
239 }
240 /*endfor*/
241 s->num_cols = s->spec->spec[s->spec->num_specs - 1].last_col;
242}
243/*- End of function --------------------------------------------------------*/
244
245static void column_search(struct meteor_working_data_s *s)
246{
247 int i;
248 int col;
249 double cost;
250
251 /* Look for favorable column to enter basis.
252 returns lowest cost and its column number, or turns on the flag optimal */
253 /* Set up price vector */
254 for (i = 0; i <= s->m; i++)
255 s->price[i] = -s->carry[0][i + 1];
256 /*endfor*/
257 s->optimal = false0;
258 s->cbar = LARGE1.0e+31;
259 s->pivot_col = 0;
260 for (col = 0; col < s->num_cols; col++)
261 {
262 cost = s->d[col];
263 for (i = 0; i <= s->m; i++)
264 cost -= s->price[i]*s->tab[i][col];
265 /*endfor*/
266 if (s->cbar > cost)
267 {
268 s->cbar = cost;
269 s->pivot_col = col + 1;
270 }
271 /*endif*/
272 }
273 /*endfor*/
274 if (s->cbar > -EPS1.0e-8)
275 s->optimal = true1;
276 /*endif*/
277}
278/*- End of function --------------------------------------------------------*/
279
280static void row_search(struct meteor_working_data_s *s)
281{
282 int i;
283 int j;
284 double ratio;
285 double min_ratio;
286
287 /* Look for pivot row. returns pivot row number, or turns on the flag unbounded */
288 /* Generate column */
289 for (i = 1; i <= (s->m + 1); i++)
290 {
291 /* Current column = B inverse * original col. */
292 s->cur_col[i] = 0.0;
293 for (j = 0; j <= s->m; j++)
294 s->cur_col[i] += s->carry[i][j + 1]*s->tab[j][s->pivot_col - 1];
295 /*endfor*/
296 }
297 /*endfor*/
298 /* First element in current column */
299 s->cur_col[0] = s->cbar;
300 s->pivot_row = -1;
301 min_ratio = LARGE1.0e+31;
302 /* Ratio test */
303 for (i = 0; i <= s->m; i++)
304 {
305 if (s->cur_col[i + 1] > EPS1.0e-8)
306 {
307 ratio = s->carry[i + 1][0]/s->cur_col[i + 1];
308 if (min_ratio > ratio)
309 {
310 /* Favorable row */
311 min_ratio = ratio;
312 s->pivot_row = i;
313 s->pivot_element = s->cur_col[i + 1];
314 }
315 else
316 {
317 /* Break tie with max pivot */
318 if (min_ratio == ratio && s->pivot_element < s->cur_col[i + 1])
319 {
320 s->pivot_row = i;
321 s->pivot_element = s->cur_col[i + 1];
322 }
323 /*endif*/
324 }
325 /*endif*/
326 }
327 /*endif*/
328 }
329 /*endfor*/
330 s->unbounded = (s->pivot_row == -1);
331}
332/*- End of function --------------------------------------------------------*/
333
334static double pivot(struct meteor_working_data_s *s)
335{
336 int i;
337 int j;
338
339 s->basis[s->pivot_row] = s->pivot_col;
340 for (j = 0; j <= (s->m + 1); j++)
341 s->carry[s->pivot_row + 1][j] /= s->pivot_element;
342 /*endfor*/
343 for (i = 0; i <= (s->m + 1); i++)
344 {
345 if ((i - 1) != s->pivot_row)
346 {
347 for (j = 0; j <= (s->m + 1); j++)
348 s->carry[i][j] -= s->carry[s->pivot_row + 1][j]*s->cur_col[i];
349 /*endfor*/
350 }
351 /*endif*/
352 }
353 /*endfor*/
354 return -s->carry[0][0];
355}
356/*- End of function --------------------------------------------------------*/
357
358static double change_phase(struct meteor_working_data_s *s)
359{
360 int i;
361 int j;
362 int b;
363
364 /* Change phase from 1 to 2, by switching to the original cost vector */
365 s->phase = 2;
366 for (i = 0; i <= s->m; i++)
367 {
368 if (s->basis[i] <= 0)
369 printf("...artificial basis element %5d remains in basis after phase 1\n", s->basis[i]);
370 /*endif*/
371 }
372 /*endfor*/
373 /* Switch to original cost vector */
374 for (i = 0; i < s->num_cols; i++)
375 s->d[i] = s->c[i];
376 /*endfor*/
377 for (j = 0; j <= (s->m + 1); j++)
378 {
379 s->carry[0][j] = 0.0;
380 for (i = 0; i <= s->m; i++)
381 {
382 /* Ignore artificial basis elements that are still in basis */
383 b = s->basis[i];
384 if (b >= 1)
385 s->carry[0][j] -= s->c[b - 1]*s->carry[i + 1][j];
386 /*endif*/
387 }
388 /*endfor*/
389 }
390 /*endfor*/
391 return -s->carry[0][0];
392}
393/*- End of function --------------------------------------------------------*/
394
395static double magnitude_response(struct meteor_working_data_s *s, double freq)
396{
397 int i;
398 double temp;
399
400 /* Compute magnitude function, given radian frequency f */
401 temp = 0.0;
402 for (i = 0; i < s->m; i++)
403 temp += s->coeff[i]*trig0(s, i, freq);
404 /*endfor*/
405 return temp;
406}
407/*- End of function --------------------------------------------------------*/
408
409static double half_magnitude_response(struct meteor_working_data_s *s, double freq)
410{
411 int i;
412 double temp;
413
414 /* Compute magnitude function, given radian frequency f */
415 temp = 0.0;
416 for (i = 0; i < (s->m + 1)/2; i++)
417 temp += s->coeff[i]*trig0(s, i, freq);
418 /*endfor*/
419 return temp;
420}
421/*- End of function --------------------------------------------------------*/
422
423static int simplex(struct meteor_working_data_s *s)
424{
425 int i;
426 int j;
427 int col;
428 int row;
429 bool_Bool done;
430
431 /* Simplex for linear programming */
432 done = false0;
433 s->phase = 1;
434 for (i = 0; i <= (s->m + 1); i++)
435 {
436 for (j = 0; j <= (MAX_COEFFS64 + 1); j++)
437 s->carry[i][j] = 0.0;
438 /*endfor*/
439 }
440 /*endfor*/
441 /* Artificial basis */
442 for (i = 1; i <= (s->m + 1); i++)
443 s->carry[i][i] = 1.0;
444 /*endfor*/
445 /* - initial cost */
446 s->carry[0][0] = -1.0;
447 s->cur_cost = -s->carry[0][0];
448 /* Variable minimized in primal */
449 s->carry[s->m + 1][0] = 1.0;
450 /* Initial, artificial basis */
451 for (i = 0; i <= s->m; i++)
452 s->basis[i] = -i;
453 /*endfor*/
454 /* Check number of columns */
455 if (s->num_cols <= NCOL_MAX6000)
456 {
457 /* Initialize cost for phase 1 */
458 for (col = 0; col < s->num_cols; col++)
459 {
460 s->d[col] = 0.0;
461 for (row = 0; row <= s->m; row++)
462 s->d[col] -= s->tab[row][col];
463 /*endfor*/
464 }
465 /*endfor*/
466 }
467 else
468 {
469 printf("...termination: too many columns for storage\n");
470 done = true1;
471 s->result = too_many_columns;
472 }
473 /*endif*/
474 s->num_pivots = 0;
475 while (s->num_pivots < MAX_PIVOTS1000 && !done && (s->cur_cost > s->low_limit || s->phase == 1))
476 {
477 column_search(s);
478 if (s->optimal)
479 {
480 if (s->phase == 1)
481 {
482 if (s->cur_cost > EPS1.0e-8)
483 {
484 /* Dual of problem is infeasible */
485 /* This happens if all specs are hugged */
486 done = true1;
487 s->result = infeasible_dual;
488 }
489 else
490 {
491 if (s->num_pivots != 1 && s->num_pivots%10 != 0)
492 printf("Pivot %d cost = %.5f\n", s->num_pivots, s->cur_cost);
493 /*endif*/
494 printf("Phase 1 successfully completed\n");
495 s->cur_cost = change_phase(s);
496 }
497 /*endif*/
498 }
499 else
500 {
501 if (s->num_pivots != 1 && s->num_pivots%10 != 0)
502 printf("Pivot %d cost = %.5f\n", s->num_pivots, s->cur_cost);
503 /*endif*/
504 printf("Phase 2 successfully completed\n");
505 done = true1;
506 s->result = optimum_obtained;
507 }
508 /*endif*/
509 }
510 else
511 {
512 row_search(s);
513 if (s->unbounded)
514 {
515 /* Dual of problem is unbounded */
516 done = true1;
517 s->result = unbounded_dual;
518 }
519 else
520 {
521 s->cur_cost = pivot(s);
522 s->num_pivots++;
523 if (s->num_pivots == 1 || s->num_pivots%10 == 0)
524 printf("Pivot %d cost = %.5f\n", s->num_pivots, s->cur_cost);
525 /*endif*/
526 }
527 /*endif*/
528 }
529 /*endif*/
530 }
531 /*endwhile*/
532 if (s->cur_cost <= s->low_limit && s->phase == 2)
533 {
534 if (s->num_pivots != 1 && s->num_pivots%10 != 0)
535 printf("Pivot %d cost = %.5f\n", s->num_pivots, s->cur_cost);
536 /*endif*/
537 s->result = infeasible_primal;
538 }
539 /*endif*/
540 if (s->num_pivots >= MAX_PIVOTS1000)
541 {
542 printf("...termination: maximum number of pivots exceeded\n");
543 s->result = too_many_pivots;
544 }
545 /*endif*/
546
547 /* Optimal */
548 return s->result;
549}
550/*- End of function --------------------------------------------------------*/
551
552static void print_result(int result)
553{
554 /* Print enumerated result type */
555 switch (result)
556 {
557 case badly_formed_requirements:
558 printf("badly formed requirements\n");
559 break;
560 case optimum_obtained:
561 printf("optimum obtained\n");
562 break;
563 case too_many_columns:
564 printf("too many columns in specifications\n");
565 break;
566 case too_many_pivots:
567 printf("too many pivots\n");
568 break;
569 case unbounded_dual:
570 printf("infeasible (unbounded dual)\n");
571 break;
572 case infeasible_dual:
573 printf("infeasible or unbounded\n");
574 break;
575 case infeasible_primal:
576 printf("infeasible\n");
577 break;
578 case no_feasible_solution_found:
579 printf("no infeasible solution found\n");
580 break;
581 case no_feasible_band_edge_found:
582 printf("no infeasible bend edge found\n");
583 break;
584 }
585 /*endswitch*/
586}
587/*- End of function --------------------------------------------------------*/
588
589static int get_m(struct meteor_working_data_s *s)
590{
591 int left_m;
592 int right_m;
593 bool_Bool found_m;
594 bool_Bool checked_left;
595 bool_Bool checked_right;
596 int result;
597 int length;
598 int i;
599
600 /* Find best order (and hence length) */
601 s->found_feasible_solution = false0;
602 left_m = s->smallest_m;
603 right_m = s->largest_m;
604 found_m = false0;
605 checked_left = false0;
606 checked_right = false0;
607 for (s->iteration = 0; !found_m; s->iteration++)
608 {
609 if (s->iteration == 0)
610 {
611 /* First time through */
612 s->m = left_m + (right_m - left_m)/2;
613 }
614 /*endif*/
615 printf("\nIteration %d\n", s->iteration);
616
617 if (s->odd_length)
618 {
619 if (s->spec->symmetry_type == symmetry_cosine)
620 length = s->m*2 - 1;
621 else
622 length = s->m*2 + 1;
623 /*endif*/
624 }
625 else
626 {
627 length = s->m*2;
628 }
629 /*endif*/
630 printf("L=%d\n", length);
631
632 setup(s);
633 result = simplex(s);
634 print_result(result);
635 if (result == optimum_obtained)
636 {
637 s->found_feasible_solution = true1;
638 right_m = s->m;
639 s->best_m = s->m;
640 /* Right side of bracket has been checked */
641 checked_right = true1;
642 if (s->odd_length)
643 {
644 if (s->spec->symmetry_type == symmetry_cosine)
645 length = s->best_m*2 - 1;
646 else
647 length = s->best_m*2 + 1;
648 /*endif*/
649 }
650 else
651 {
652 length = s->best_m*2;
653 }
654 /*endif*/
655 printf("New best length L=%d\n", length);
656
657 for (i = 0; i < s->m; i++)
658 s->coeff[i] = -s->carry[0][i + 1];
659 /*endfor*/
660 }
661 /*endif*/
662
663 if (result != optimum_obtained)
664 {
665 left_m = s->m;
666 /* Left side of bracket has been checked */
667 checked_left = true1;
668 }
669 /*endif*/
670
671 if (right_m > left_m + 1)
672 s->m = left_m + (right_m - left_m)/2;
673 /*endif*/
674
675 if (right_m == left_m + 1)
676 {
677 if (!checked_left)
678 {
679 s->m = left_m;
680 checked_left = true1;
681 }
682 else if (!checked_right)
683 {
684 s->m = right_m;
685 checked_right = true1;
686 }
687 else
688 {
689 found_m = true1;
690 }
691 /*endif*/
692 }
693 /*endif*/
694
695 if (right_m == left_m)
696 found_m = true1;
697 /*endif*/
698 }
699 /*endfor*/
700
701 if (!s->found_feasible_solution)
702 return no_feasible_solution_found;
703 /*endif*/
704 s->m = s->best_m;
705
706 putchar('\n');
707 if (s->odd_length)
708 {
709 if (s->spec->symmetry_type == symmetry_cosine)
710 printf("Best length L=%d\n", s->best_m*2 - 1);
711 else
712 printf("Best length L=%d\n", s->best_m*2 + 1);
713 /*endif*/
714 }
715 /*endif*/
716
717 if (!s->odd_length)
718 printf("Best length L=%d\n", s->best_m*2);
719 /*endif*/
720 return 0;
721}
722/*- End of function --------------------------------------------------------*/
723
724static int get_edge(struct meteor_working_data_s *s)
725{
726 double left_edge;
727 double right_edge;
728 double newe;
729 double beste;
1
'beste' declared without an initial value
→
730 double one_space;
731 double stop_space;
732 int i;
733
734 /* Optimize band edge */
735 one_space = 0.5/s->spec->grid_points; /* Space between grid points */
736 stop_space = one_space/10.0;
737 /* Stop criterion is 1/10 nominal grid spacing */
738 if (s->which_way == rr)
2
←
Assuming the condition is false
→
3
←
Taking false branch
→
739 {
740 /* Start with rightmost left edge */
741 left_edge = s->spec->spec[s->spec->spec[0].band_pushed - 1].left_freq;
742 for (i = 1; i < s->num_pushed; i++)
743 {
744 if (s->spec->spec[s->spec->spec[i].band_pushed - 1].left_freq > left_edge)
745 left_edge = s->spec->spec[s->spec->spec[i].band_pushed - 1].left_freq;
746 /*endif*/
747 }
748 /*endfor*/
749 right_edge = 0.5;
750 }
751 else
752 {
753 /* Start with leftmost right edge */
754 left_edge = 0.0;
755 right_edge = s->spec->spec[s->spec->spec[0].band_pushed - 1].right_freq;
756 for (i = 1; i < s->num_pushed; i++)
4
←
Assuming the condition is false
→
5
←
Loop condition is false. Execution continues on line 765
→
757 {
758 if (s->spec->spec[s->spec->spec[i].band_pushed - 1].right_freq < right_edge)
759 right_edge = s->spec->spec[s->spec->spec[i].band_pushed - 1].right_freq;
760 /*endif*/
761 }
762 /*endfor*/
763 }
764 /*endif*/
765 s->found_feasible_solution = false0;
766 for (s->iteration = 0; (right_edge - left_edge) > stop_space; s->iteration++)
6
←
Loop condition is true. Entering loop body
→
12
←
Loop condition is false. Execution continues on line 807
→
767 {
768 newe = (right_edge + left_edge)/2.0;
769 printf("\nIteration %d\n", s->iteration);
770 printf("Trying new edge = %10.4f\n", newe);
771 for (i = 0; i < s->num_pushed; i++)
7
←
Assuming the condition is false
→
8
←
Loop condition is false. Execution continues on line 780
→
772 {
773 if (s->which_way == rr)
774 s->spec->spec[s->spec->spec[i].band_pushed - 1].right_freq = newe;
775 else
776 s->spec->spec[s->spec->spec[i].band_pushed - 1].left_freq = newe;
777 /*endif*/
778 }
779 /*endif*/
780 setup(s);
781 s->result = simplex(s);
782 print_result(s->result);
783 if (s->result == optimum_obtained)
9
←
Taking false branch
→
784 {
785 if (s->which_way == rr)
786 left_edge = newe;
787 else
788 right_edge = newe;
789 /*endif*/
790 s->found_feasible_solution = true1;
791 beste = newe;
792 for (i = 0; i < s->m; i++)
793 s->coeff[i] = -s->carry[0][i + 1];
794 /*endfor*/
795 }
796 else
797 {
798 if (s->which_way == rr)
10
←
Assuming the condition is false
→
11
←
Taking false branch
→
799 right_edge = newe;
800 else
801 left_edge = newe;
802 /*endif*/
803 }
804 /*endif*/
805 }
806 /*endfor*/
807 putchar('\n');
808 if (!s->found_feasible_solution)
13
←
Assuming the condition is false
→
14
←
Taking false branch
→
809 return no_feasible_band_edge_found;
810 /*endif*/
811 printf("Found edge = %10.4f\n", beste);
15
←
Function call argument is an uninitialized value
812 for (i = 0; i < s->num_pushed; i++)
813 {
814 if (s->which_way == rr)
815 s->spec->spec[s->spec->spec[i].band_pushed - 1].right_freq = beste;
816 else
817 s->spec->spec[s->spec->spec[i].band_pushed - 1].left_freq = beste;
818 /*endif*/
819 }
820 /*endfor*/
821 for (i = 0; i < s->spec->num_specs; i++)
822 make_bands(s, i);
823 /*endfor*/
824 return 0;
825}
826/*- End of function --------------------------------------------------------*/
827
828static int get_max_dist(struct meteor_working_data_s *s)
829{
830 int i;
831
832 /* Maximize distance from constraints */
833 printf("Optimization: maximize distance from constraints\n");
834 setup(s);
835 s->result = simplex(s);
836 print_result(s->result);
837 if (s->result != optimum_obtained)
838 return s->result;
839 /*endif*/
840 printf("Final cost = distance from constraints = %.5f\n", s->cur_cost);
841 /* Record coefficients */
842 for (i = 0; i < s->m; i++)
843 s->coeff[i] = -s->carry[0][i + 1];
844 /*endfor*/
845 return 0;
846}
847/*- End of function --------------------------------------------------------*/
848
849static int meteor_get_coefficients(struct meteor_working_data_s *s, double coeffs[])
850{
851 int i;
852 int j;
853
854 j = 0;
855 if (s->odd_length && s->spec->symmetry_type == symmetry_cosine)
856 {
857 for (i = s->m - 1; i >= 1; i--)
858 coeffs[j++] = s->coeff[i]/2.0;
859 /*endfor*/
860 coeffs[j++] = s->coeff[0];
861 for (i = 1; i < s->m; i++)
862 coeffs[j++] = s->coeff[i]/2.0;
863 /*endfor*/
864 }
865 else if (!s->odd_length && s->spec->symmetry_type == symmetry_cosine)
866 {
867 for (i = s->m - 1; i >= 0; i--)
868 coeffs[j++] = s->coeff[i]/2.0;
869 /*endfor*/
870 for (i = 0; i < s->m; i++)
871 coeffs[j++] = s->coeff[i]/2.0;
872 /*endfor*/
873 }
874 else if (s->odd_length && s->spec->symmetry_type == symmetry_sine)
875 {
876 /* L = length, odd */
877 /* Negative of the first m coefs. */
878 for (i = s->m - 1; i >= 0; i--)
879 coeffs[j++] = -s->coeff[i]/2.0;
880 /*endfor*/
881 /* Middle coefficient is always 0 */
882 coeffs[j++] = 0.0;
883 for (i = 0; i < s->m; i++)
884 coeffs[j++] = s->coeff[i]/2.0;
885 /*endfor*/
886 }
887 else if (!s->odd_length && s->spec->symmetry_type == symmetry_sine)
888 {
889 /* Negative of the first m coefs. */
890 for (i = s->m - 1; i >= 0; i--)
891 coeffs[j++] = -s->coeff[i]/2.0;
892 /*endfor*/
893 for (i = 0; i < s->m; i++)
894 coeffs[j++] = s->coeff[i]/2.0;
895 /*endfor*/
896 }
897 /*endif*/
898 return j;
899}
900/*- End of function --------------------------------------------------------*/
901
902static int vet_data(struct meteor_working_data_s *s)
903{
904 int i;
905 bool_Bool all_hugged;
906 char ch;
907
908 printf("Filter name: '%s'\n", s->spec->filter_name);
909
910 if (s->spec->shortest < 1 || s->spec->longest > MAX_TAPS129)
911 {
912 printf("Shortest or longest out of range\n");
913 return badly_formed_requirements;
914 }
915 /*endif*/
916
917 if ((s->spec->shortest & 1) != (s->spec->longest & 1))
918 {
919 printf("Parity of smallest andlongest unequal\n");
920 return badly_formed_requirements;
921 }
922 /*endif*/
923
924 s->odd_length = s->spec->shortest & 1;
925 if (s->odd_length)
926 {
927 if (s->spec->symmetry_type == symmetry_cosine)
928 {
929 s->smallest_m = (s->spec->shortest + 1)/2;
930 s->largest_m = (s->spec->longest + 1)/2;
931 }
932 else
933 {
934 s->smallest_m = (s->spec->shortest - 1)/2;
935 s->largest_m = (s->spec->longest - 1)/2;
936 }
937 /*endif*/
938 }
939 else
940 {
941 s->smallest_m = s->spec->shortest/2;
942 s->largest_m = s->spec->longest/2;
943 }
944 /*endif*/
945
946 if (s->spec->shortest != s->spec->longest)
947 {
948 s->what_to_do = find_len;
949 printf("Finding minimum length: range %d to %d\n", s->spec->shortest, s->spec->longest);
950 }
951 else
952 {
953 s->m = s->smallest_m;
954 s->length = s->spec->shortest;
955
956 printf("Fixed length of %4d\n", s->length);
957 scanf("%c%*[^\n]", &ch);
958 /* Right, left, or neither: edges to be pushed? */
959 getchar();
960
961 if (ch == 'n')
962 {
963 s->what_to_do = max_dist;
964 }
965 else
966 {
967 s->what_to_do = push_edge;
968
969 s->which_way = (ch == 'r') ? rr : ll;
970
971 scanf("%d%*[^\n]", &s->num_pushed);
972 getchar();
973 for (i = 0; i < s->num_pushed; i++)
974 scanf("%d", &s->spec->spec[i].band_pushed);
975 /*endfor*/
976 scanf("%*[^\n]");
977 getchar();
978
979 printf("Pushing band edges right\n", (s->which_way == rr) ? "right" : "left");
980
981 printf("Constraint numbers: ");
982 for (i = 0; i < s->num_pushed; i++)
983 printf("%3d ", s->spec->spec[i].band_pushed);
984 /*endfor*/
985 putchar('\n');
986 }
987 /*endif*/
988 }
989 /*endif*/
990
991 for (i = 0; i < s->spec->num_specs; i++)
992 {
993 printf("Constraint name '%s'\n", s->spec->spec[i].name);
994 switch (s->spec->spec[i].type)
995 {
996 case constraint_type_convexity:
997 switch (s->spec->spec[i].sense)
998 {
999 case sense_convex:
1000 printf("Constraint %2d: convexity, sense convex\n", i);
1001 break;
1002 case sense_concave:
1003 printf("Constraint %2d: convexity, sense concave\n", i);
1004 break;
1005 }
1006 /*endswitch*/
1007
1008 printf(" Band edges: %10.4f %10.4f\n", s->spec->spec[i].left_freq, s->spec->spec[i].right_freq);
1009 break;
1010 case constraint_type_limit:
1011 if (s->spec->spec[i].interpolation == interpolation_geometric && s->spec->spec[i].left_bound*s->spec->spec[i].right_bound == 0.0)
1012 {
1013 printf("Geometrically interpolated band edge in constraint %5d is zero\n", i);
1014 return badly_formed_requirements;
1015 }
1016 /*endif*/
1017
1018 switch (s->spec->spec[i].sense)
1019 {
1020 case sense_lower:
1021 printf(" Constraint %2d: lower limit\n", i);
1022 break;
1023 case sense_upper:
1024 printf(" Constraint %2d: upper limit\n", i);
1025 break;
1026 case sense_envelope:
1027 printf(" Constraint %2d: envelope limit\n", i);
1028 break;
1029 }
1030 /*endswitch*/
1031
1032 switch (s->spec->spec[i].interpolation)
1033 {
1034 case interpolation_geometric:
1035 printf(" Geometric interpolation\n");
1036 break;
1037 case interpolation_arithmetic:
1038 printf(" Arithmetic interpolation\n");
1039 break;
1040 }
1041 /*endswitch*/
1042
1043 if (s->spec->spec[i].hug)
1044 printf(" This constraint will be hugged\n");
1045 else
1046 printf(" This constraint will be optimized\n");
1047 /*endif*/
1048
1049 printf(" Band edges: %10.4f %10.4f\n", s->spec->spec[i].left_freq, s->spec->spec[i].right_freq);
1050 printf(" Bounds: %10.4f %10.4f\n", s->spec->spec[i].left_bound, s->spec->spec[i].right_bound);
1051 break;
1052 }
1053 /*endswitch*/
1054 make_bands(s, i);
1055 printf(" Initial columns: %10d %10d\n", s->spec->spec[i].first_col, s->spec->spec[i].last_col);
1056 }
1057 s->num_cols = s->spec->spec[s->spec->num_specs - 1].last_col;
1058
1059 printf("Number of specs = %5d\n", s->spec->num_specs);
1060 printf("Initial number of columns = %5d\n", s->num_cols);
1061
1062 all_hugged = true1;
1063 for (i = 0; i < s->spec->num_specs; i++)
1064 {
1065 if (s->spec->spec[i].type == constraint_type_limit && !s->spec->spec[i].hug)
1066 all_hugged = false0;
1067 /*endif*/
1068 }
1069 /*endfor*/
1070
1071 if (all_hugged)
1072 {
1073 printf("All constraints are hugged: ill-posed problem\n");
1074 return badly_formed_requirements;
1075 }
1076 /*endif*/
1077 return 0;
1078}
1079/*- End of function --------------------------------------------------------*/
1080
1081void output_filter_performance_as_csv_file(struct meteor_working_data_s *s, const char *file_name)
1082{
1083 int i;
1084 double mg;
1085 double mg2;
1086 FILE *magnitude_file;
1087
1088 if (s->log_fd)
1089 {
1090 magnitude_file = s->log_fd;
1091 }
1092 else
1093 {
1094 /* Print final frequency response */
1095 if ((magnitude_file = fopen(file_name, "wb")) == NULL((void*)0))
1096 {
1097 fprintf(stderrstderr, "Cannot open file '%s'\n", file_name);
1098 exit(2);
1099 }
1100 /*endif*/
1101 }
1102 /*endif*/
1103
1104 if (s->spec->filter_name && s->spec->filter_name[0])
1105 {
1106 fprintf(magnitude_file, "%s\n", s->spec->filter_name);
1107 }
1108 /*endif*/
1109 fprintf(magnitude_file, "Frequency, Gain (dB), Gain (linear), Half gain (linear)\n");
1110 /* Magnitude on regular grid */
1111 for (i = 0; i <= s->spec->grid_points; i++)
1112 {
1113 mg = fabs(magnitude_response(s, i*M_PI3.14159265358979323846/s->spec->grid_points));
1114 if (mg == 0.0)
1115 mg = SMALL1.0e-8;
1116 /*endif*/
1117 mg2 = fabs(half_magnitude_response(s, i*M_PI3.14159265358979323846/s->spec->grid_points));
1118 if (mg2 == 0.0)
1119 mg2 = SMALL1.0e-8;
1120 /*endif*/
1121 fprintf(magnitude_file, "%10.4lf, %.10lf, %.5lf, %.5lf\n", 0.5*s->spec->sample_rate*i/s->spec->grid_points, 20.0*log10(mg), mg, mg2);
1122 }
1123 /*endfor*/
1124 fprintf(magnitude_file, "\nMagnitude at band edges\n\n");
1125 for (i = 0; i < s->spec->num_specs; i++)
1126 {
1127 if (s->spec->spec[i].type == constraint_type_limit)
1128 {
1129 fprintf(magnitude_file, "%10.4f %.5E\n",
1130 s->freq[s->spec->spec[i].first_col - 1]*0.5/M_PI3.14159265358979323846, magnitude_response(s, s->freq[s->spec->spec[i].first_col - 1]));
1131 fprintf(magnitude_file, "%10.4f %.5E\n",
1132 s->freq[s->spec->spec[i].last_col - 1]*0.5/M_PI3.14159265358979323846, magnitude_response(s, s->freq[s->spec->spec[i].last_col - 1]));
1133 putchar('\n');
1134 }
1135 /*endif*/
1136 }
1137 /*endfor*/
1138 if (s->log_fd == NULL((void*)0))
1139 fclose(magnitude_file);
1140 /*endif*/
1141}
1142/*- End of function --------------------------------------------------------*/
1143
1144int meteor_design_filter(struct meteor_working_data_s *s, struct meteor_spec_s *t, double coeffs[])
1145{
1146 int res;
1147
1148 memset(s, 0, sizeof(*s));
1149
1150 s->spec = t;
1151 if ((res = vet_data(s)) < 0)
1152 return res;
1153 /*endif*/
1154 /* Dual cost negative => primal infeasible */
1155 s->low_limit = -EPS1.0e-8;
1156 switch (s->what_to_do)
1157 {
1158 case find_len:
1159 res = get_m(s);
1160 break;
1161 case push_edge:
1162 res = get_edge(s);
1163 break;
1164 case max_dist:
1165 res = get_max_dist(s);
1166 break;
1167 }
1168 /*endswitch*/
1169 if (res < 0)
1170 return res;
1171 /*endif*/
1172 return meteor_get_coefficients(s, coeffs);
1173}
1174/*- End of function --------------------------------------------------------*/
1175/*- End of file ------------------------------------------------------------*/