1 /* Conditional Median Filter - Better Version
2 * Using fft and eliminating high frequencies
3 * to fill over the ticks.
4 *
5 * Copyright (C) 1998 J.A. Bezemer
6 * 2001 S.J. Tappin
7 *
8 * Licensed under the terms of the GNU General Public License.
9 * ABSOLUTELY NO WARRANTY.
10 * See the file `COPYING' in this directory.
11 */
12
13 /* Remove the `dont' to get the gate instead of the normal output - useful
14 for verifying properties. */
15 #define dontVIEW_INTERNALS
16
17 /* Choose the highpass filter: */
18 #define noSECOND_ORDER
19 #define FOURTH_ORDER
20 #define noSIXTH_ORDER
21
22 #define noDEBUGFILE
23
24 #include "signpr_cmf3.h"
25 #include "signpr_general.h"
26 #include "signpr_l1fit.h"
27 #include "errorwindow.h"
28 #include "stringinput.h"
29 #include "buttons.h"
30 #include "clrscr.h"
31 #include "boxes.h"
32 #include "helpline.h"
33 #include "yesnowindow.h"
34 #include <stdlib.h>
35 #include <stdio.h>
36 #include <math.h>
37 #ifndef OLD_CURSES
38 #include <ncurses.h>
39 #else
40 #include <curses.h>
41 #endif
42
43 #ifdef DEBUGFILE
44 static FILE *debugf=NULL;
45 #endif
46
47 #ifndef M_PIl
48 # define M_PIl 3.1415926535897932384626433832795029L /* pi */
49 #endif
50
51 void
cond_median3_param_defaults(parampointer_t parampointer)52 cond_median3_param_defaults (parampointer_t parampointer)
53 {
54 parampointer->postlength2 = 4;
55 parampointer->prelength2 = 4;
56 parampointer->postlength3 = 5;
57 parampointer->prelength3 = 5;
58 parampointer->postlength4 = 256;
59 parampointer->prelength4 = 255;
60
61 parampointer->int1 = 12;
62 parampointer->int2 = 9; /* This could be derived from
63 postlength4 and prelength4,
64 but it's messy */
65 #if defined (SECOND_ORDER)
66 parampointer->long1 = 1000; /* Threshold to detect precise tick length. */
67 parampointer->long2 = 2500; /* Must be above this to be a tick. */
68 #elif defined (FOURTH_ORDER)
69 parampointer->long1 = 2000;
70 parampointer->long2 = 8500;
71 #elif defined (SIXTH_ORDER)
72 parampointer->long1 = 1500;
73 parampointer->long2 = 7500;
74 #else
75 #error A Highpass version must be defined (signpr_cmf3.c)
76 #endif
77
78 }
79
80
81 #ifdef FOURTH_ORDER
82 #undef SIGNPR_CMF3_PARAMSCR_HEADERTEXT
83 #define SIGNPR_CMF3_PARAMSCR_HEADERTEXT "CMF IIF [FOURTH ORDER] - Parameters"
84 #endif
85
86 #ifdef SIXTH_ORDER
87 #undef SIGNPR_CMF3_PARAMSCR_HEADERTEXT
88 #define SIGNPR_CMF3_PARAMSCR_HEADERTEXT "CMF IIF [SIXTH ORDER] - Parameters"
89 #endif
90
91 void
cond_median3_param_screen(parampointer_t parampointer)92 cond_median3_param_screen (parampointer_t parampointer)
93 {
94 stringinput_t rmslengthstr, rmflengthstr, decimatestr, threshold1str;
95 stringinput_t threshold2str, fftstr;
96 button_t ok_button, cancel_button, defaults_button;
97 int dont_stop = TRUE;
98 int focus = 0;
99 int in_ch;
100 int i;
101 long helplong;
102
103 char *helplines[9] =
104 {
105 " ^: less ticks detected. v: not all of tick interpolated. ",
106 " ^: bad following of dynamics. v: less ticks detected. ",
107 " ^: bad following of dynamics. v: less ticks detected. ",
108 " ^: detected tick length too short v: detected tick length longer. ",
109 " ^: only strong ticks detected. v: music-ticks also filtered out. ",
110 " ^: Slower interpolation v: possible problems with long ticks ",
111 " Discard changes. ",
112 " Reset default values. ",
113 " Accept changes. "};
114
115 rmslengthstr.maxlen = 500;
116 rmslengthstr.string = (char *) malloc (rmslengthstr.maxlen *
117 sizeof (char));
118 sprintf (rmslengthstr.string, "%ld", parampointer->prelength2 +
119 parampointer->postlength2 + 1);
120 rmslengthstr.y = 6;
121 rmslengthstr.x = 59;
122 rmslengthstr.w = 19;
123 rmslengthstr.cursorpos = strlen (rmslengthstr.string);
124 rmslengthstr.firstcharonscreen = 0;
125
126 rmflengthstr.maxlen = 500;
127 rmflengthstr.string = (char *) malloc (rmflengthstr.maxlen *
128 sizeof (char));
129 sprintf (rmflengthstr.string, "%ld", parampointer->prelength3 +
130 parampointer->postlength3 + 1);
131 rmflengthstr.y = 8;
132 rmflengthstr.x = 59;
133 rmflengthstr.w = 19;
134 rmflengthstr.cursorpos = strlen (rmflengthstr.string);
135 rmflengthstr.firstcharonscreen = 0;
136
137 decimatestr.maxlen = 500;
138 decimatestr.string = (char *) malloc (decimatestr.maxlen *
139 sizeof (char));
140 sprintf (decimatestr.string, "%d", parampointer->int1);
141 decimatestr.y = 10;
142 decimatestr.x = 59;
143 decimatestr.w = 19;
144 decimatestr.cursorpos = strlen (decimatestr.string);
145 decimatestr.firstcharonscreen = 0;
146
147 threshold1str.maxlen = 500;
148 threshold1str.string = (char *) malloc (threshold1str.maxlen *
149 sizeof (char));
150 sprintf (threshold1str.string, "%ld", parampointer->long1);
151 threshold1str.y = 12;
152 threshold1str.x = 59;
153 threshold1str.w = 19;
154 threshold1str.cursorpos = strlen (threshold1str.string);
155 threshold1str.firstcharonscreen = 0;
156
157 threshold2str.maxlen = 500;
158 threshold2str.string = (char *) malloc (threshold2str.maxlen *
159 sizeof (char));
160 sprintf (threshold2str.string, "%ld", parampointer->long2);
161 threshold2str.y = 14;
162 threshold2str.x = 59;
163 threshold2str.w = 19;
164 threshold2str.cursorpos = strlen (threshold2str.string);
165 threshold2str.firstcharonscreen = 0;
166
167 fftstr.maxlen = 500;
168 fftstr.string = (char *) malloc (fftstr.maxlen *
169 sizeof (char));
170 sprintf (fftstr.string, "%d", parampointer->int2);
171 fftstr.y = 16;
172 fftstr.x = 59;
173 fftstr.w = 19;
174 fftstr.cursorpos = strlen (fftstr.string);
175 fftstr.firstcharonscreen = 0;
176
177
178 ok_button.text = " OK ";
179 ok_button.y = 20;
180 ok_button.x = 71;
181 ok_button.selected = FALSE;
182
183 cancel_button.text = " Cancel ";
184 cancel_button.y = 20;
185 cancel_button.x = 5;
186 cancel_button.selected = FALSE;
187
188 defaults_button.text = " Defaults ";
189 defaults_button.y = 20;
190 defaults_button.x = 36;
191 defaults_button.selected = FALSE;
192
193 clearscreen (SIGNPR_CMF3_PARAMSCR_HEADERTEXT);
194
195 do
196 {
197 header (SIGNPR_CMF3_PARAMSCR_HEADERTEXT);
198
199 if (focus == 6)
200 cancel_button.selected = TRUE;
201 else
202 cancel_button.selected = FALSE;
203
204 if (focus == 7)
205 defaults_button.selected = TRUE;
206 else
207 defaults_button.selected = FALSE;
208
209 if (focus == 8)
210 ok_button.selected = TRUE;
211 else
212 ok_button.selected = FALSE;
213
214 mvprintw (3, 2,
215 "See also the Signproc.txt file for the meaning of the parameters.");
216
217 stringinput_display (&rmslengthstr);
218 mvprintw (rmslengthstr.y, 2,
219 "Length of the RMS operation (samples):");
220
221 stringinput_display (&rmflengthstr);
222 mvprintw (rmflengthstr.y, 2,
223 "Length of the recursive median operation (samples):");
224
225 stringinput_display (&decimatestr);
226 mvprintw (decimatestr.y, 2,
227 "Decimation factor for the recursive median:");
228
229 stringinput_display (&threshold1str);
230 mvprintw (threshold1str.y, 2,
231 "Fine threshold for tick start/end (thousandths):");
232
233 stringinput_display (&threshold2str);
234 mvprintw (threshold2str.y, 2,
235 "Threshold for detection of tick presence (thousandths):");
236
237 stringinput_display (&fftstr);
238 mvprintw (fftstr.y, 2,
239 "Length for fft to interpolate (2^n):");
240
241 button_display (&cancel_button);
242 mybox (cancel_button.y - 1, cancel_button.x - 1,
243 3, strlen (cancel_button.text) + 2);
244 button_display (&defaults_button);
245 mybox (defaults_button.y - 1, defaults_button.x - 1,
246 3, strlen (defaults_button.text) + 2);
247 button_display (&ok_button);
248 mybox (ok_button.y - 1, ok_button.x - 1,
249 3, strlen (ok_button.text) + 2);
250
251 helpline (helplines[focus]);
252
253 switch (focus)
254 {
255 case 0:
256 stringinput_display (&rmslengthstr);
257 break;
258 case 1:
259 stringinput_display (&rmflengthstr);
260 break;
261 case 2:
262 stringinput_display (&decimatestr);
263 break;
264 case 3:
265 stringinput_display (&threshold1str);
266 break;
267 case 4:
268 stringinput_display (&threshold2str);
269 break;
270 case 5:
271 stringinput_display (&fftstr);
272 break;
273 default:
274 move (0, 79);
275 }
276
277 refresh ();
278
279 in_ch = getch ();
280
281 switch (focus)
282 {
283 case 0: /* rmslengthstr */
284 stringinput_stdkeys (in_ch, &rmslengthstr);
285 switch (in_ch)
286 {
287 case KEY_ENTER:
288 case 13:
289 i = sscanf (rmslengthstr.string, "%li", &helplong);
290 if (i < 1 || helplong < 1 || helplong % 2 == 0)
291 error_window ("A whole, odd number, greater than 0, must \
292 be specified.");
293 else
294 focus++;
295 break;
296
297 case KEY_UP:
298 focus--;
299 break;
300 case KEY_DOWN:
301 focus++;
302 break;
303 }
304 break;
305
306 case 1: /* rmflengthstr */
307 stringinput_stdkeys (in_ch, &rmflengthstr);
308 switch (in_ch)
309 {
310 case KEY_ENTER:
311 case 13:
312 i = sscanf (rmflengthstr.string, "%li", &helplong);
313 if (i < 1 || helplong < 1 || helplong % 2 == 0)
314 error_window ("A whole, odd number, greater than 0, must \
315 be specified.");
316 else
317 focus++;
318 break;
319
320 case KEY_UP:
321 focus--;
322 break;
323 case KEY_DOWN:
324 focus++;
325 break;
326 }
327 break;
328
329 case 2: /* decimatestr */
330 stringinput_stdkeys (in_ch, &decimatestr);
331 switch (in_ch)
332 {
333 case KEY_ENTER:
334 case 13:
335 i = sscanf (decimatestr.string, "%li", &helplong);
336 if (i < 1 || helplong < 1)
337 error_window ("A whole number, greater than 0, must \
338 be specified.");
339 else
340 focus++;
341 break;
342
343 case KEY_UP:
344 focus--;
345 break;
346 case KEY_DOWN:
347 focus++;
348 break;
349 }
350 break;
351
352 case 3: /* threshold1str */
353 stringinput_stdkeys (in_ch, &threshold1str);
354 switch (in_ch)
355 {
356 case KEY_ENTER:
357 case 13:
358 i = sscanf (threshold1str.string, "%li", &helplong);
359 if (i < 1 || helplong < 1)
360 error_window ("A whole, positive number must be specified.");
361 else
362 focus++;
363 break;
364
365 case KEY_UP:
366 focus--;
367 break;
368 case KEY_DOWN:
369 focus++;
370 break;
371 }
372 break;
373
374 case 4: /* threshold2str */
375 stringinput_stdkeys (in_ch, &threshold2str);
376 switch (in_ch)
377 {
378 case KEY_ENTER:
379 case 13:
380 i = sscanf (threshold2str.string, "%li", &helplong);
381 if (i < 1 || helplong < 1)
382 error_window ("A whole, positive number must be specified.");
383 else
384 focus++;
385 break;
386
387 case KEY_UP:
388 focus--;
389 break;
390 case KEY_DOWN:
391 focus++;
392 break;
393 }
394 break;
395
396 case 5: /* fft length */
397 stringinput_stdkeys (in_ch, &fftstr);
398 switch (in_ch)
399 {
400 case KEY_ENTER:
401 case 13:
402 i = sscanf (fftstr.string, "%li", &helplong);
403 if (i < 1 || helplong > 12 || helplong < 6)
404 error_window ("A number between 6 and 12 must be specified.");
405 else
406 focus = 8;
407 break;
408
409 case KEY_UP:
410 focus--;
411 break;
412 case KEY_DOWN:
413 focus++;
414 break;
415 }
416 break;
417
418 case 6: /* Cancel */
419 switch (in_ch)
420 {
421 case KEY_ENTER:
422 case 13:
423 dont_stop = FALSE;
424 break;
425
426 case KEY_LEFT:
427 case KEY_UP:
428 focus--;
429 break;
430 case KEY_RIGHT:
431 case KEY_DOWN:
432 focus++;
433 break;
434 }
435 break;
436
437 case 7: /* Defaults */
438 switch (in_ch)
439 {
440 case KEY_ENTER:
441 case 13:
442 if (yesno_window ("Restore default parameters?", " Yes ",
443 " No ", 0))
444 {
445 cond_median3_param_defaults (parampointer);
446 dont_stop = FALSE;
447 }
448 break;
449
450 case KEY_LEFT:
451 case KEY_UP:
452 focus--;
453 break;
454 case KEY_RIGHT:
455 case KEY_DOWN:
456 focus++;
457 break;
458 }
459 break;
460
461 case 8: /* OK */
462 switch (in_ch)
463 {
464 case KEY_ENTER:
465 case 13:
466
467 i = sscanf (rmslengthstr.string, "%li", &helplong);
468 if (i < 1 || helplong < 1 || helplong % 2 == 0)
469 {
470 error_window ("A whole, odd number, greater than 0, must \
471 be specified as RMS length.");
472 rmslengthstr.cursorpos =
473 strlen (rmslengthstr.string);
474 focus = 0;
475 break;
476 }
477
478 parampointer->prelength2 = (helplong - 1) / 2;
479 parampointer->postlength2 = (helplong - 1) / 2;
480
481 i = sscanf (rmflengthstr.string, "%li", &helplong);
482 if (i < 1 || helplong < 1 || helplong % 2 == 0)
483 {
484 error_window ("A whole, odd number, greater than 0, must \
485 be specified as length of the recursive median.");
486 rmflengthstr.cursorpos =
487 strlen (rmflengthstr.string);
488 focus = 1;
489 break;
490 }
491
492 parampointer->prelength3 = (helplong - 1) / 2;
493 parampointer->postlength3 = (helplong - 1) / 2;
494
495 i = sscanf (decimatestr.string, "%li", &helplong);
496 if (i < 1 || helplong < 1)
497 {
498 error_window ("A whole number, greater than 0, must \
499 be specified as decimation factor.");
500 decimatestr.cursorpos =
501 strlen (decimatestr.string);
502 focus = 2;
503 break;
504 }
505
506 parampointer->int1 = helplong;
507
508 i = sscanf (threshold1str.string, "%li", &helplong);
509 if (i < 1 || helplong < 1)
510 {
511 error_window ("A whole, positive number must be \
512 specified as threshold.");
513 threshold1str.cursorpos =
514 strlen (threshold1str.string);
515 focus = 3;
516 break;
517 }
518
519 parampointer->long1 = helplong;
520
521 i = sscanf (threshold2str.string, "%li", &helplong);
522 if (i < 1 || helplong < 1000)
523 {
524 error_window ("A whole, positive number must be \
525 specified as threshold.");
526 threshold2str.cursorpos =
527 strlen (threshold2str.string);
528 focus = 4;
529 break;
530 }
531
532 parampointer->long2 = helplong;
533
534 i = sscanf (fftstr.string, "%li", &helplong);
535 if (i < 1 || helplong > 12 || helplong < 6) {
536 error_window ("A number between 6 and 12 must be specified \
537 as fft length");
538 fftstr.cursorpos = strlen(fftstr.string);
539 focus = 5;
540 break;
541 }
542
543 parampointer->int2 = helplong;
544 parampointer->prelength4 = 1;
545 for (i=1; i < parampointer->int2; i++) parampointer->prelength4 *= 2;
546 parampointer->postlength4 = parampointer->prelength4-1;
547
548 dont_stop = FALSE;
549 break;
550
551 case KEY_LEFT:
552 case KEY_UP:
553 focus--;
554 break;
555 case KEY_RIGHT:
556 case KEY_DOWN:
557 focus++;
558 break;
559 }
560 break;
561 }
562
563 if (in_ch == 9) /* TAB */
564 focus++;
565
566 if (in_ch == 27)
567 dont_stop = FALSE;
568
569 if (focus > 8)
570 focus = 0;
571 if (focus < 0)
572 focus = 8;
573 }
574 while (dont_stop);
575
576 free (rmslengthstr.string);
577 free (rmflengthstr.string);
578 free (decimatestr.string);
579 free (threshold1str.string);
580 free (threshold2str.string);
581 }
582
583 void
init_cond_median3_filter(int filterno,parampointer_t parampointer)584 init_cond_median3_filter (int filterno, parampointer_t parampointer)
585 {
586 long total_post;
587 long total_pre;
588 long l;
589
590 total_post = parampointer->postlength4 + parampointer->prelength4 + 1 + 4;
591 /* +1+4=+5 : for highpass, max. 11th order */
592
593 total_pre = parampointer->postlength4 + parampointer->prelength4 + 1;
594 l = parampointer->prelength4 + parampointer->prelength3 *
595 parampointer->int1 + parampointer->prelength2 + 5;
596 /* + 5 : for highpass, max. 11th order */
597 if (l > total_pre)
598 total_pre = l;
599
600 parampointer->buffer = init_buffer (total_post, total_pre);
601 parampointer->buffer2 = init_buffer (parampointer->postlength2,
602 parampointer->prelength2);
603 parampointer->buffer3 = init_buffer (parampointer->postlength3,
604 parampointer->prelength3 * parampointer->int1);
605 parampointer->buffer4 = init_buffer (parampointer->postlength4,
606 parampointer->prelength4);
607
608 parampointer->filterno = filterno;
609
610 /* Set up the FFT plans here. Since we expect to do lots of FFTs,
611 we take the time to MEASURE the best way to do them [SJT] */
612
613 parampointer->planf = rfftw_create_plan(parampointer->postlength4 +
614 parampointer->prelength4 + 1,
615 FFTW_REAL_TO_COMPLEX,
616 FFTW_MEASURE);
617 parampointer->planr = rfftw_create_plan(parampointer->postlength4 +
618 parampointer->prelength4 + 1,
619 FFTW_COMPLEX_TO_REAL,
620 FFTW_MEASURE);
621
622 #ifdef DEBUGFILE
623 debugf = fopen("./gram.txt","w");
624 #endif
625
626
627 }
628
629
630 void
delete_cond_median3_filter(parampointer_t parampointer)631 delete_cond_median3_filter (parampointer_t parampointer)
632 {
633 delete_buffer (¶mpointer->buffer);
634 delete_buffer (¶mpointer->buffer2);
635 delete_buffer (¶mpointer->buffer3);
636 delete_buffer (¶mpointer->buffer4);
637
638 rfftw_destroy_plan(parampointer->planf);
639 rfftw_destroy_plan(parampointer->planr);
640
641 #ifdef DEBUGFILE
642 fclose(debugf);
643 #endif
644 }
645
646
647 sample_t
cond_median3_highpass(long offset,long offset_zero,parampointer_t parampointer)648 cond_median3_highpass (long offset, long offset_zero,
649 parampointer_t parampointer)
650 {
651 sample_t sample;
652 longsample_t sum;
653
654 offset += offset_zero; /* middle for highpass filter in
655 'big buffer' */
656 sum.left = 0;
657 sum.right = 0;
658
659 #if defined (SECOND_ORDER)
660 #define notTEST_DAVE_PLATT
661 #ifndef TEST_DAVE_PLATT
662 /* Original: */
663 sample = get_from_buffer (¶mpointer->buffer, offset - 1);
664 sum.left += sample.left;
665 sum.right += sample.right;
666 sample = get_from_buffer (¶mpointer->buffer, offset);
667 sum.left -= 2 * (long) sample.left;
668 sum.right -= 2 * (long) sample.right;
669 sample = get_from_buffer (¶mpointer->buffer, offset + 1);
670 sum.left += sample.left;
671 sum.right += sample.right;
672
673 sum.left /= 4;
674 sum.right /= 4;
675 #else /* TEST_DAVE_PLATT */
676 /* Testing, suggested by Dave Platt. Invert phase of one channel, then
677 do tick detection using the sum signal. This is because most ticks
678 are out-of-phase signals. I've not really tested this - it might
679 require other settings for thresholds etc.
680 Note: implemented for second_order only! */
681 sample = get_from_buffer (¶mpointer->buffer, offset - 1);
682 sum.left += sample.left;
683 sum.left -= sample.right;
684 sample = get_from_buffer (¶mpointer->buffer, offset);
685 sum.left -= 2 * (long) sample.left;
686 sum.left += 2 * (long) sample.right;
687 sample = get_from_buffer (¶mpointer->buffer, offset + 1);
688 sum.left += sample.left;
689 sum.left -= sample.right;
690
691 /* just in case L/R: 32000/-32000 -32000/32000 32000/-32000 : */
692 sum.left /= 8;
693 sum.right = sum.left;
694 #endif /* TEST_DAVE_PLATT */
695
696 #elif defined (FOURTH_ORDER)
697 sample = get_from_buffer (¶mpointer->buffer, offset - 2);
698 sum.left += sample.left;
699 sum.right += sample.right;
700 sample = get_from_buffer (¶mpointer->buffer, offset - 1);
701 sum.left -= 4 * (long) sample.left;
702 sum.right -= 4 * (long) sample.right;
703 sample = get_from_buffer (¶mpointer->buffer, offset);
704 sum.left += 6 * (long) sample.left;
705 sum.right += 6 * (long) sample.right;
706 sample = get_from_buffer (¶mpointer->buffer, offset + 1);
707 sum.left -= 4 * (long) sample.left;
708 sum.right -= 4 * (long) sample.right;
709 sample = get_from_buffer (¶mpointer->buffer, offset + 2);
710 sum.left += sample.left;
711 sum.right += sample.right;
712
713 sum.left /= 4;
714 sum.right /= 4;
715
716 #elif defined (SIXTH_ORDER)
717 sample = get_from_buffer (¶mpointer->buffer, offset - 3);
718 sum.left -= sample.left;
719 sum.right -= sample.right;
720 sample = get_from_buffer (¶mpointer->buffer, offset - 2);
721 sum.left += 6 * (long) sample.left;
722 sum.right += 6 * (long) sample.right;
723 sample = get_from_buffer (¶mpointer->buffer, offset - 1);
724 sum.left -= 15 * (long) sample.left;
725 sum.right -= 15 * (long) sample.right;
726 sample = get_from_buffer (¶mpointer->buffer, offset);
727 sum.left += 20 * (long) sample.left;
728 sum.right += 20 * (long) sample.right;
729 sample = get_from_buffer (¶mpointer->buffer, offset + 1);
730 sum.left -= 15 * (long) sample.left;
731 sum.right -= 15 * (long) sample.right;
732 sample = get_from_buffer (¶mpointer->buffer, offset + 2);
733 sum.left += 6 * (long) sample.left;
734 sum.right += 6 * (long) sample.right;
735 sample = get_from_buffer (¶mpointer->buffer, offset + 3);
736 sum.left -= sample.left;
737 sum.right -= sample.right;
738
739 /* Should be /64, but the signal is extremely soft, so divide by less to
740 get more quantization levels (more accurate) */
741 sum.left /= 4;
742 sum.right /= 4;
743 #endif
744
745 if (sum.left > 32767)
746 sample.left = 32767;
747 else if (sum.left < -32768)
748 sample.left = -32768;
749 else
750 sample.left = sum.left;
751
752
753 if (sum.right > 32767)
754 sample.right = 32767;
755 else if (sum.right < -32768)
756 sample.right = -32768;
757 else
758 sample.right = sum.right;
759
760
761 return sample;
762 }
763
764 fillfuncpointer_t cond_median3_highpass_pointer = cond_median3_highpass;
765
766 sample_t
cond_median3_rms(long offset,long offset_zero,parampointer_t parampointer)767 cond_median3_rms (long offset, long offset_zero,
768 parampointer_t parampointer)
769 {
770 sample_t sample;
771 doublesample_t doublesample;
772 doublesample_t sum;
773 long i;
774
775 advance_current_pos_custom (¶mpointer->buffer2,
776 cond_median3_highpass_pointer,
777 offset + offset_zero,
778 parampointer);
779
780 sum.left = 0;
781 sum.right = 0;
782
783 for (i = -parampointer->postlength2; i <= parampointer->prelength2;
784 i++)
785 {
786 sample = get_from_buffer (¶mpointer->buffer2, i);
787 doublesample.left = sample.left;
788 doublesample.right = sample.right;
789 sum.left += doublesample.left * doublesample.left;
790 sum.right += doublesample.right * doublesample.right;
791 }
792
793 sum.left /= (parampointer->postlength2 +
794 parampointer->prelength2 + 1);
795 sum.right /= (parampointer->postlength2 +
796 parampointer->prelength2 + 1);
797
798 sample.left = sqrt (sum.left + 1);
799 sample.right = sqrt (sum.right + 1);
800
801 return sample;
802 }
803
804 fillfuncpointer_t cond_median3_rms_pointer = cond_median3_rms;
805
806 sample_t
cond_median3_gate(long offset,long offset_zero,parampointer_t parampointer)807 cond_median3_gate (long offset, long offset_zero,
808 parampointer_t parampointer)
809 /* Well, not a `gate' any more - just (w[t]-b[t])/b[t], comparision to
810 make the real gate is done later. */
811 {
812 sample_t sample;
813 sample_t w_t;
814 sample_t b_t;
815 sample_t returnval;
816 signed short list1[parampointer->postlength3 +
817 parampointer->prelength3 * parampointer->int1 + 1];
818 signed short list2[parampointer->postlength3 +
819 parampointer->prelength3 * parampointer->int1 + 1];
820 long i, j;
821
822 advance_current_pos_custom (¶mpointer->buffer3,
823 cond_median3_rms_pointer,
824 offset + offset_zero,
825 parampointer);
826
827 w_t = get_from_buffer (¶mpointer->buffer3, 0);
828
829 /* The RMF Filter */
830
831 for (i = 0; i < parampointer->postlength3; i++)
832 {
833 sample = get_from_buffer (¶mpointer->buffer3,
834 i - parampointer->postlength3);
835 list1[i] = sample.left;
836 list2[i] = sample.right;
837 }
838
839 j = i;
840
841 for (; i <= parampointer->postlength3 +
842 parampointer->prelength3 * parampointer->int1;
843 i += parampointer->int1)
844 {
845 sample = get_from_buffer (¶mpointer->buffer3,
846 i - parampointer->postlength3);
847 list1[j] = sample.left;
848 list2[j] = sample.right;
849 j++;
850 }
851
852 b_t.left = median (list1, j);
853 b_t.right = median (list2, j);
854
855 put_in_buffer (¶mpointer->buffer3, 0, b_t);
856
857
858 i = (labs (w_t.left - b_t.left) * 1000)
859 /
860 b_t.left;
861 if (i > 32767)
862 i = 32767;
863 else if (i < -32768)
864 i = -32768;
865
866 returnval.left = i;
867
868 i = (labs (w_t.right - b_t.right) * 1000)
869 /
870 b_t.right;
871 if (i > 32767)
872 i = 32767;
873 else if (i < -32768)
874 i = -32768;
875 returnval.right = i;
876
877 return returnval;
878 }
879
880 fillfuncpointer_t cond_median3_gate_pointer = cond_median3_gate;
881
882 sample_t
cond_median3_filter(parampointer_t parampointer)883 cond_median3_filter (parampointer_t parampointer)
884 {
885 sample_t sample, gate, returnval;
886 /* Length of the fft we'll do to get the smoothed interpolate */
887
888 fftw_real list3[parampointer->postlength4 +
889 parampointer->prelength4 + 1];
890 fftw_real list4[parampointer->postlength4 +
891 parampointer->prelength4 + 1];
892
893 long fft_l = parampointer->postlength4 + parampointer->prelength4 + 1;
894 long i;
895 int toleft, toright, nfreq;
896 signed short maxval;
897
898 advance_current_pos (¶mpointer->buffer, parampointer->filterno);
899
900 advance_current_pos_custom (¶mpointer->buffer4,
901 cond_median3_gate_pointer,
902 0,
903 parampointer);
904
905 gate = get_from_buffer (¶mpointer->buffer4, 0);
906
907 #ifdef VIEW_INTERNALS
908 returnval.left = 0;
909 returnval.right = 0;
910 #else
911 /* 'Default' value - unchanged if there is no tick */
912 returnval = get_from_buffer (¶mpointer->buffer, 0);
913 #endif
914
915 if (gate.left > parampointer->long1)
916 {
917 maxval = gate.left;
918
919 toleft = -1;
920 sample.left = 0;
921 do
922 {
923 toleft++;
924 if (toleft < parampointer->postlength4)
925 {
926 sample = get_from_buffer (¶mpointer->buffer4, -toleft - 1);
927 if (sample.left > maxval)
928 /* if so, the `while' will continue anyway, so maxval may
929 be adjusted here already (if necessary) */
930 maxval = sample.left;
931 }
932 }
933 while (toleft < parampointer->postlength4 &&
934 sample.left > parampointer->long1);
935
936 toright = -1;
937 sample.left = 0;
938 do
939 {
940 toright++;
941 if (toright < parampointer->prelength4)
942 {
943 sample = get_from_buffer (¶mpointer->buffer4, toright + 1);
944 if (sample.left > maxval)
945 /* if so, the `while' will continue anyway, so maxval may
946 be adjusted here already (if necessary) */
947 maxval = sample.left;
948 }
949 }
950 while (toright < parampointer->prelength4 &&
951 sample.left > parampointer->long1);
952
953 /* only interpolate if there really is a tick */
954 if (maxval > parampointer->long2)
955 {
956
957 #ifdef VIEW_INTERNALS
958 returnval.left = (toright + toleft + 1) * 500;
959 #else
960
961 /* Use a HANNING window here for the time being; note that
962 the FFT is centred at the middle of the tick not at the
963 point we are interpolating. */
964 for (i = 0; i < fft_l; i++)
965 {
966 list3[i] = get_from_buffer(¶mpointer->buffer,
967 i - parampointer->prelength4 +
968 (toright - toleft + 1)/2).left *
969 (2.-cos(2.*M_PIl*(double) i/ ((double) fft_l - 1.)))/2.;
970 }
971 rfftw_one(parampointer->planf, list3, list4);
972 nfreq=floor((double) fft_l/(double) (2*(toleft+toright+1)));
973 for (i = 2*nfreq; i <= fft_l - 2*nfreq; i++) list4[i] = 0.;
974 for (i = nfreq; i<2*nfreq && i< fft_l/2; i++) {
975 list4[i] *= (1. - (double) (i-nfreq)/ (double) nfreq);
976 list4[fft_l-i] *= (1. - (double) (i-nfreq)/ (double) nfreq);
977 }
978 rfftw_one(parampointer->planr, list4, list3);
979 returnval.left = (signed short) (list3[parampointer->prelength4 -
980 (toright - toleft + 1)/2]/
981 (double) fft_l);
982
983 /* DON'T ASK !!! -- I have NO idea why I have to MULTIPLY by
984 the hanning window here. Everything sensible says DIVIDE, but
985 multiply works, divide doesn't */
986
987 returnval.left *= (2.-cos(2.*M_PIl*(double) ((toright - toleft + 1)/2) /
988 ((double) fft_l - 1.)))/2.;
989 #ifdef DEBUGFILE
990 fprintf(debugf, "L: %ld %d %ld %ld\n",
991 fft_l, (toleft+toright+1), nfreq,
992 fft_l-nfreq);
993 #endif
994
995 #endif
996 }
997 }
998
999 if (gate.right > parampointer->long1)
1000 {
1001 maxval = gate.right;
1002
1003 toleft = -1;
1004 sample.right = 0;
1005 do
1006 {
1007 toleft++;
1008 if (toleft < parampointer->postlength4)
1009 {
1010 sample = get_from_buffer (¶mpointer->buffer4, -toleft - 1);
1011 if (sample.right > maxval)
1012 /* if so, the `while' will continue anyway, so maxval may
1013 be adjusted here already (if necessary) */
1014 maxval = sample.right;
1015 }
1016 }
1017 while (toleft < parampointer->postlength4 &&
1018 sample.right > parampointer->long1);
1019
1020 toright = -1;
1021 sample.right = 0;
1022 do
1023 {
1024 toright++;
1025 if (toright < parampointer->prelength4)
1026 {
1027 sample = get_from_buffer (¶mpointer->buffer4, toright + 1);
1028 if (sample.right > maxval)
1029 /* if so, the `while' will continue anyway, so maxval may
1030 be adjusted here already (if necessary) */
1031 maxval = sample.right;
1032 }
1033 }
1034 while (toright < parampointer->prelength4 &&
1035 sample.right > parampointer->long1);
1036
1037 /* only interpolate if there really is a tick */
1038 if (maxval > parampointer->long2)
1039 {
1040
1041 #ifdef VIEW_INTERNALS
1042 returnval.right = (toright + toleft + 1) * 500;
1043 #else
1044 /* Use a HANNING window here for the time being; note that
1045 the FFT is centred at the middle of the tick not at the
1046 point we are interpolating. */
1047 for (i = 0; i < fft_l; i++)
1048 {
1049 list3[i] = get_from_buffer(¶mpointer->buffer,
1050 i - parampointer->prelength4 +
1051 (toright - toleft + 1)/2).right *
1052 (2.-cos(2.*M_PIl*(double) i/ ((double) fft_l - 1.)))/2.;
1053 }
1054 rfftw_one(parampointer->planf, list3, list4);
1055
1056 nfreq=floor((double) fft_l/(double) (2*(toleft+toright+1)));
1057 for (i = 2*nfreq; i <= fft_l - 2*nfreq; i++) list4[i] = 0.;
1058 for (i = nfreq; i<2*nfreq && i<fft_l/2; i++) {
1059 list4[i] *= (1. - (double) (i-nfreq)/ (double) nfreq);
1060 list4[fft_l-i] *= (1. - (double) (i-nfreq)/ (double) nfreq);
1061 }
1062
1063 rfftw_one(parampointer->planr, list4, list3);
1064 returnval.right = (signed short) (list3[parampointer->prelength4 -
1065 (toright - toleft + 1)/2]/
1066 (double) fft_l) ;
1067
1068 /* DON'T ASK !!! -- I have NO idea why I have to MULTIPLY by
1069 the hanning window here. Everything sensible says DIVIDE, but
1070 multiply works, divide doesn't */
1071
1072 returnval.right *= (2.-cos(2.*M_PIl*(double) ((toright - toleft + 1)/2) /
1073 ((double) fft_l - 1.)))/2.;
1074 #ifdef DEBUGFILE
1075 fprintf(debugf, "R: %ld %d %ld %ld\n",
1076 fft_l, (toleft+toright+1), nfreq,
1077 fft_l - nfreq);
1078 #endif
1079
1080 #endif
1081 }
1082 }
1083
1084 return returnval;
1085 }
1086