1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | |
8 | |
9 | |
10 | |
11 | |
12 | |
13 | |
14 | |
15 | |
16 | |
17 | |
18 | |
19 | |
20 | |
21 | |
22 | |
23 | |
24 | |
25 | |
26 | |
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 | |
44 | |
45 | |
46 | |
47 | |
48 | |
49 | |
50 | |
51 | |
52 | |
53 | |
54 | |
55 | |
56 | |
57 | |
58 | |
59 | |
60 | |
61 | |
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 | |
68 | static void make_bands(struct meteor_working_data_s *s, int i) |
69 | { |
70 | int j; |
71 | int kmax; |
72 | |
73 | |
74 | |
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 | |
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 | |
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 | |
91 | } |
92 | |
93 | s->spec->spec[i].last_col = s->spec->spec[i].first_col + kmax; |
94 | } |
95 | |
96 | |
97 | static double trig0(struct meteor_working_data_s *s, int i, double freq) |
98 | { |
99 | double res; |
100 | |
101 | |
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 | |
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 | |
117 | } |
118 | |
119 | return res; |
120 | } |
121 | |
122 | |
123 | static double trig2(struct meteor_working_data_s *s, int i, double freq) |
124 | { |
125 | double res; |
126 | |
127 | |
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 | |
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 | |
143 | } |
144 | |
145 | return res; |
146 | } |
147 | |
148 | |
149 | static void convex(struct meteor_working_data_s *s, int i) |
150 | { |
151 | int row; |
152 | int col; |
153 | |
154 | |
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 | |
159 | s->c[col] = 0.0; |
160 | for (row = 0; row < s->m; row++) |
161 | { |
162 | |
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 | |
168 | } |
169 | |
170 | s->tab[s->m][col] = 0.0; |
171 | } |
172 | |
173 | } |
174 | |
175 | |
176 | static void limit(struct meteor_working_data_s *s, int i) |
177 | { |
178 | int row; |
179 | int col; |
180 | |
181 | |
182 | |
183 | |
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 | |
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 | |
206 | } |
207 | |
208 | if (s->spec->spec[i].sense == sense_lower) |
209 | s->c[col] = -s->c[col]; |
210 | |
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 | |
216 | s->tab[s->m][col] = (s->spec->spec[i].hug) ? 0.0 : 1.0; |
217 | } |
218 | |
219 | } |
220 | |
221 | |
222 | static void setup(struct meteor_working_data_s *s) |
223 | { |
224 | int i; |
225 | |
226 | |
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 | |
239 | } |
240 | |
241 | s->num_cols = s->spec->spec[s->spec->num_specs - 1].last_col; |
242 | } |
243 | |
244 | |
245 | static void column_search(struct meteor_working_data_s *s) |
246 | { |
247 | int i; |
248 | int col; |
249 | double cost; |
250 | |
251 | |
252 | |
253 | |
254 | for (i = 0; i <= s->m; i++) |
255 | s->price[i] = -s->carry[0][i + 1]; |
256 | |
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 | |
266 | if (s->cbar > cost) |
267 | { |
268 | s->cbar = cost; |
269 | s->pivot_col = col + 1; |
270 | } |
271 | |
272 | } |
273 | |
274 | if (s->cbar > -EPS1.0e-8) |
275 | s->optimal = true1; |
276 | |
277 | } |
278 | |
279 | |
280 | static 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 | |
288 | |
289 | for (i = 1; i <= (s->m + 1); i++) |
290 | { |
291 | |
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 | |
296 | } |
297 | |
298 | |
299 | s->cur_col[0] = s->cbar; |
300 | s->pivot_row = -1; |
301 | min_ratio = LARGE1.0e+31; |
302 | |
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 | |
311 | min_ratio = ratio; |
312 | s->pivot_row = i; |
313 | s->pivot_element = s->cur_col[i + 1]; |
314 | } |
315 | else |
316 | { |
317 | |
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 | |
324 | } |
325 | |
326 | } |
327 | |
328 | } |
329 | |
330 | s->unbounded = (s->pivot_row == -1); |
331 | } |
332 | |
333 | |
334 | static 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 | |
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 | |
350 | } |
351 | |
352 | } |
353 | |
354 | return -s->carry[0][0]; |
355 | } |
356 | |
357 | |
358 | static double change_phase(struct meteor_working_data_s *s) |
359 | { |
360 | int i; |
361 | int j; |
362 | int b; |
363 | |
364 | |
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 | |
371 | } |
372 | |
373 | |
374 | for (i = 0; i < s->num_cols; i++) |
375 | s->d[i] = s->c[i]; |
376 | |
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 | |
383 | b = s->basis[i]; |
384 | if (b >= 1) |
385 | s->carry[0][j] -= s->c[b - 1]*s->carry[i + 1][j]; |
386 | |
387 | } |
388 | |
389 | } |
390 | |
391 | return -s->carry[0][0]; |
392 | } |
393 | |
394 | |
395 | static double magnitude_response(struct meteor_working_data_s *s, double freq) |
396 | { |
397 | int i; |
398 | double temp; |
399 | |
400 | |
401 | temp = 0.0; |
402 | for (i = 0; i < s->m; i++) |
403 | temp += s->coeff[i]*trig0(s, i, freq); |
404 | |
405 | return temp; |
406 | } |
407 | |
408 | |
409 | static double half_magnitude_response(struct meteor_working_data_s *s, double freq) |
410 | { |
411 | int i; |
412 | double temp; |
413 | |
414 | |
415 | temp = 0.0; |
416 | for (i = 0; i < (s->m + 1)/2; i++) |
417 | temp += s->coeff[i]*trig0(s, i, freq); |
418 | |
419 | return temp; |
420 | } |
421 | |
422 | |
423 | static 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 | |
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 | |
439 | } |
440 | |
441 | |
442 | for (i = 1; i <= (s->m + 1); i++) |
443 | s->carry[i][i] = 1.0; |
444 | |
445 | |
446 | s->carry[0][0] = -1.0; |
447 | s->cur_cost = -s->carry[0][0]; |
448 | |
449 | s->carry[s->m + 1][0] = 1.0; |
450 | |
451 | for (i = 0; i <= s->m; i++) |
452 | s->basis[i] = -i; |
453 | |
454 | |
455 | if (s->num_cols <= NCOL_MAX6000) |
456 | { |
457 | |
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 | |
464 | } |
465 | |
466 | } |
467 | else |
468 | { |
469 | printf("...termination: too many columns for storage\n"); |
470 | done = true1; |
471 | s->result = too_many_columns; |
472 | } |
473 | |
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 | |
485 | |
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 | |
494 | printf("Phase 1 successfully completed\n"); |
495 | s->cur_cost = change_phase(s); |
496 | } |
497 | |
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 | |
504 | printf("Phase 2 successfully completed\n"); |
505 | done = true1; |
506 | s->result = optimum_obtained; |
507 | } |
508 | |
509 | } |
510 | else |
511 | { |
512 | row_search(s); |
513 | if (s->unbounded) |
514 | { |
515 | |
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 | |
526 | } |
527 | |
528 | } |
529 | |
530 | } |
531 | |
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 | |
537 | s->result = infeasible_primal; |
538 | } |
539 | |
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 | |
546 | |
547 | |
548 | return s->result; |
549 | } |
550 | |
551 | |
552 | static void print_result(int result) |
553 | { |
554 | |
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 | |
586 | } |
587 | |
588 | |
589 | static 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 | |
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 | |
612 | s->m = left_m + (right_m - left_m)/2; |
613 | } |
614 | |
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 | |
624 | } |
625 | else |
626 | { |
627 | length = s->m*2; |
628 | } |
629 | |
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 | |
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 | |
649 | } |
650 | else |
651 | { |
652 | length = s->best_m*2; |
653 | } |
654 | |
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 | |
660 | } |
661 | |
662 | |
663 | if (result != optimum_obtained) |
664 | { |
665 | left_m = s->m; |
666 | |
667 | checked_left = true1; |
668 | } |
669 | |
670 | |
671 | if (right_m > left_m + 1) |
672 | s->m = left_m + (right_m - left_m)/2; |
673 | |
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 | |
692 | } |
693 | |
694 | |
695 | if (right_m == left_m) |
696 | found_m = true1; |
697 | |
698 | } |
699 | |
700 | |
701 | if (!s->found_feasible_solution) |
702 | return no_feasible_solution_found; |
703 | |
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 | |
714 | } |
715 | |
716 | |
717 | if (!s->odd_length) |
718 | printf("Best length L=%d\n", s->best_m*2); |
719 | |
720 | return 0; |
721 | } |
722 | |
723 | |
724 | static 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 | |
735 | one_space = 0.5/s->spec->grid_points; |
736 | stop_space = one_space/10.0; |
737 | |
738 | if (s->which_way == rr) |
| 2 | | Assuming the condition is false | |
|
| |
739 | { |
740 | |
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 | |
747 | } |
748 | |
749 | right_edge = 0.5; |
750 | } |
751 | else |
752 | { |
753 | |
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 | |
761 | } |
762 | |
763 | } |
764 | |
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 | |
778 | } |
779 | |
780 | setup(s); |
781 | s->result = simplex(s); |
782 | print_result(s->result); |
783 | if (s->result == optimum_obtained) |
| |
784 | { |
785 | if (s->which_way == rr) |
786 | left_edge = newe; |
787 | else |
788 | right_edge = newe; |
789 | |
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 | |
795 | } |
796 | else |
797 | { |
798 | if (s->which_way == rr) |
| 10 | | Assuming the condition is false | |
|
| |
799 | right_edge = newe; |
800 | else |
801 | left_edge = newe; |
802 | |
803 | } |
804 | |
805 | } |
806 | |
807 | putchar('\n'); |
808 | if (!s->found_feasible_solution) |
| 13 | | Assuming the condition is false | |
|
| |
809 | return no_feasible_band_edge_found; |
810 | |
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 | |
819 | } |
820 | |
821 | for (i = 0; i < s->spec->num_specs; i++) |
822 | make_bands(s, i); |
823 | |
824 | return 0; |
825 | } |
826 | |
827 | |
828 | static int get_max_dist(struct meteor_working_data_s *s) |
829 | { |
830 | int i; |
831 | |
832 | |
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 | |
840 | printf("Final cost = distance from constraints = %.5f\n", s->cur_cost); |
841 | |
842 | for (i = 0; i < s->m; i++) |
843 | s->coeff[i] = -s->carry[0][i + 1]; |
844 | |
845 | return 0; |
846 | } |
847 | |
848 | |
849 | static 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 | |
860 | coeffs[j++] = s->coeff[0]; |
861 | for (i = 1; i < s->m; i++) |
862 | coeffs[j++] = s->coeff[i]/2.0; |
863 | |
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 | |
870 | for (i = 0; i < s->m; i++) |
871 | coeffs[j++] = s->coeff[i]/2.0; |
872 | |
873 | } |
874 | else if (s->odd_length && s->spec->symmetry_type == symmetry_sine) |
875 | { |
876 | |
877 | |
878 | for (i = s->m - 1; i >= 0; i--) |
879 | coeffs[j++] = -s->coeff[i]/2.0; |
880 | |
881 | |
882 | coeffs[j++] = 0.0; |
883 | for (i = 0; i < s->m; i++) |
884 | coeffs[j++] = s->coeff[i]/2.0; |
885 | |
886 | } |
887 | else if (!s->odd_length && s->spec->symmetry_type == symmetry_sine) |
888 | { |
889 | |
890 | for (i = s->m - 1; i >= 0; i--) |
891 | coeffs[j++] = -s->coeff[i]/2.0; |
892 | |
893 | for (i = 0; i < s->m; i++) |
894 | coeffs[j++] = s->coeff[i]/2.0; |
895 | |
896 | } |
897 | |
898 | return j; |
899 | } |
900 | |
901 | |
902 | static 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 | |
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 | |
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 | |
938 | } |
939 | else |
940 | { |
941 | s->smallest_m = s->spec->shortest/2; |
942 | s->largest_m = s->spec->longest/2; |
943 | } |
944 | |
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 | |
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 | |
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 | |
985 | putchar('\n'); |
986 | } |
987 | |
988 | } |
989 | |
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 | |
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 | |
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 | |
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 | |
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 | |
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 | |
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 | |
1068 | } |
1069 | |
1070 | |
1071 | if (all_hugged) |
1072 | { |
1073 | printf("All constraints are hugged: ill-posed problem\n"); |
1074 | return badly_formed_requirements; |
1075 | } |
1076 | |
1077 | return 0; |
1078 | } |
1079 | |
1080 | |
1081 | void 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 | |
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 | |
1101 | } |
1102 | |
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 | |
1109 | fprintf(magnitude_file, "Frequency, Gain (dB), Gain (linear), Half gain (linear)\n"); |
1110 | |
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 | |
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 | |
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 | |
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 | |
1136 | } |
1137 | |
1138 | if (s->log_fd == NULL((void*)0)) |
1139 | fclose(magnitude_file); |
1140 | |
1141 | } |
1142 | |
1143 | |
1144 | int 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 | |
1154 | |
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 | |
1169 | if (res < 0) |
1170 | return res; |
1171 | |
1172 | return meteor_get_coefficients(s, coeffs); |
1173 | } |
1174 | |
1175 | |