home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / filter1d.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-03-13  |  27.3 KB  |  878 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)filter1d.c    2.19  3/13/95
  3.  *
  4.  *    Copyright (c) 1991-1995 by P. Wessel and W. H. F. Smith
  5.  *    See README file for copying and redistribution conditions.
  6.  *--------------------------------------------------------------------*/
  7. /*
  8.  * filter1d will read N colums of data from file/stdin and return
  9.  * filtered output at user-selected positions.  The time variable
  10.  * can be in any specified column of the input.  Several filters
  11.  * are available, with robustness as an option.
  12.  *
  13.  * Filters:    Boxcar, Cosine Arch, Gaussian, Median, or Mode.
  14.  * Robust:    Option replaces outliers with medians in filter.
  15.  * Output:    At input times, or from t_start to t_stop by t_int.
  16.  * Lack:    Option checks for data gaps in input series.
  17.  * Symmetry:    Option checks for asymmetry in filter window.
  18.  * Quality:    Option checks for low mean weight in window.
  19.  *
  20.  * Author:    Walter H. F. Smith
  21.  * Date:    13 January, 1989
  22.  * Version:    2.0, Developmental
  23.  */
  24.  
  25. #include "gmt.h"
  26.  
  27. #define BOXCAR 0
  28. #define COS_ARCH 1
  29. #define GAUSSIAN 2
  30. #define MEDIAN 3
  31. #define MODE 4
  32. #define N_FILTERS 5
  33. #define CONVOLVE 2        /* If filter_type > CONVOLVE then a MEDIAN or MODE is selected  */
  34.  
  35. PFD get_weight[3];        /* Selects desired weight function.  */
  36.  
  37. double    boxcar_weight ();    /* Returns the weight for a given delta_time  */
  38. double    cosine_weight ();
  39. double    gaussian_weight ();
  40.  
  41. int    get_median_scale();    /* Function that sets the median and L1 scale in a data window */
  42. int    median();        /* New routine that performs iterative search to update last median  */
  43. int    mode();            /* Searches for the mode of the array  */
  44. int    load_data_and_check_extrema ();        /* Does just what it says  */
  45. int    set_up_filter ();    /* Creates filter weights or reads them from file and sets start & stop  */
  46. int    do_the_filter ();    /* Does the job  */
  47. int    allocate_space (), free_space ();    /* Does just what it says  */
  48. BOOLEAN    lack_check ();        /* Tests for lack of data condition  (optional)  */
  49.  
  50. double    *f_wt;            /* Pointer for array of filter coefficients  */
  51. double    *wt_sum;        /* Pointer for array of weight sums [each column]  */
  52. double    *data_sum;        /* Pointer for array of data * weight sums [columns]  */
  53. double    *min_loc;        /* Pointer for array of values, one per [column]  */
  54. double    *max_loc;
  55. double    *last_loc;
  56. double    *this_loc;
  57. double    *min_scl;
  58. double    *max_scl;
  59. double    *last_scl;
  60. double    *this_scl;
  61. double    **work;            /* Pointer to array of pointers to doubles for work  */    
  62. double    **data;            /* Pointer to array of pointers to doubles for data  */    
  63. double    dt;            /* Delta time resolution for filter computation  */
  64. double    q_factor;        /* Quality level for mean weights or n in median  */
  65. double filter_width;        /* Full width of filter in user's units */
  66. double half_width;        
  67. double ignore_val;        /* Value to skip over in input stream  */    
  68. double t_start;            /* x-value of first output point */
  69. double t_stop;            /* x-value of last output point */
  70. double t_int;            /* Output interval */
  71. double sym_coeff = 0.0;        /* Symmetry coefficient  */
  72. double lack_width = 0.0;    /* Lack of data width  */
  73.  
  74. int    *n_this_col;        /* Pointer to array of counters [one per column]  */
  75. int    *n_left;        /* Pointer to array of counters [one per column]  */
  76. int    *n_right;        /* Pointer to array of counters [one per column]  */
  77.  
  78. int    n_rows;            /* Number of rows of input  */
  79. int    n_cols = 2;        /* Number of columns of input  */
  80. int    t_col = 0;        /* Column of time abscissae (independent variable)  */
  81. int    n_f_wts;        /* Number of filter weights  */
  82. int    half_n_f_wts;        /* Half the number of filter weights  */
  83. int    n_row_alloc = GMT_CHUNK;        /* Number of rows of data to allocate  */
  84. int    n_work_alloc = GMT_CHUNK;    /* Number of rows of workspace to allocate  */
  85. int    filter_type = 0;    /* Flag indicating desired filter type  */
  86.  
  87. BOOLEAN    *good_one;        /* Pointer to array of logicals [one per column]  */
  88.  
  89. BOOLEAN    use_ends = FALSE;    /* True to start/stop at ends of series instead of 1/2 width inside  */
  90. BOOLEAN check_asym = FALSE;    /* TRUE to test whether the data are asymmetric about the output time  */
  91. BOOLEAN check_lack = FALSE;    /* TRUE to test for lack of data (gap) in the filter window */
  92. BOOLEAN check_q = FALSE;    /* TRUE to test average weight or N in median */
  93. BOOLEAN robust = FALSE;        /* Look for outliers in data when TRUE */
  94. BOOLEAN equidist = TRUE;    /* Data is evenly sampled in t */
  95. BOOLEAN out_at_time = FALSE;    /* TRUE when output is required at evenly spaced intervals */
  96. BOOLEAN ignore_test = FALSE;    /* TRUE when input values that eqaul ignore_val should be discarded */
  97. BOOLEAN error = FALSE;
  98. BOOLEAN single_precision = FALSE;
  99. BOOLEAN b_in = FALSE;
  100. BOOLEAN b_out = FALSE;
  101.  
  102. char    buffer[512];        /* Used to scan input lines  */
  103.  
  104. FILE    *fp_in = NULL;        /* File pointer to data file or stdin */
  105. FILE    *fp_wt = NULL;        /* File pointer to custom weight coefficients (optional) */
  106.  
  107. main (argc, argv)
  108. int argc;
  109. char **argv;
  110. {
  111.  
  112.     int    i;
  113.     
  114.     argc = gmt_begin (argc, argv);
  115.  
  116.     for (i = 1; i < argc; i++) {
  117.         if (argv[i][0] == '-') {
  118.             switch (argv[i][1]) {
  119.         
  120.                 /* Common parameters */
  121.             
  122.                 case 'H':
  123.                 case 'V':
  124.                 case ':':
  125.                 case '\0':
  126.                     error += get_common_args (argv[i], 0, 0, 0, 0);
  127.                     break;
  128.                 
  129.                 /* Supplemental parameters */
  130.                 
  131.                 case 'b':
  132.                     single_precision = !(argv[i][strlen(argv[i])-1] == 'd');
  133.                     if (argv[i][2] == 'i')
  134.                         b_in = TRUE;
  135.                     else if (argv[i][2] == 'o')
  136.                         b_out = TRUE;
  137.                     else
  138.                         b_in = b_out = TRUE;
  139.                     break;
  140.                     
  141.                 case 'F':    /* Filter selection  */
  142.                     switch (argv[i][2]) {
  143.                         case 'B':
  144.                             robust = TRUE;
  145.                         case 'b':
  146.                             filter_type = BOXCAR;
  147.                             filter_width = atof (&argv[i][3]);
  148.                             break;
  149.                         case 'C':
  150.                             robust = TRUE;
  151.                         case 'c':
  152.                             filter_type = COS_ARCH;
  153.                             filter_width = atof (&argv[i][3]);
  154.                             break;
  155.                         case 'G':
  156.                             robust = TRUE;
  157.                         case 'g':
  158.                             filter_type = GAUSSIAN;
  159.                             filter_width = atof (&argv[i][3]);
  160.                             break;
  161.                         case 'M':
  162.                             robust = TRUE;
  163.                         case 'm':
  164.                             filter_type = MEDIAN;
  165.                             filter_width = atof (&argv[i][3]);
  166.                             break;
  167.                         case 'P':
  168.                             robust = TRUE;
  169.                         case 'p':
  170.                             filter_type = MODE;
  171.                             filter_width = atof (&argv[i][3]);
  172.                             break;
  173.                         case 'F':
  174.                             robust = TRUE;
  175.                         case 'f':
  176.                             if ((fp_wt = fopen(&argv[i][3],"r")) == NULL) {
  177.                                 fprintf (stderr, "filter1d: Cannot open file %s\n", &argv[i][3]);
  178.                                 exit(-1);
  179.                             }
  180.                             break;
  181.                         default:
  182.                             error = TRUE;
  183.                             break;
  184.                     }
  185.                     break;
  186.                     
  187.                 case 'D':
  188.                     dt = atof (&argv[i][2]);
  189.                     equidist = FALSE;
  190.                     break;
  191.                 case 'E':
  192.                     use_ends = TRUE;
  193.                     break;
  194.                 case 'I':
  195.                     ignore_val = atof (&argv[i][2]);
  196.                     ignore_test = TRUE;
  197.                     break;
  198.                 case 'L':
  199.                     check_lack = TRUE;
  200.                     lack_width = atof (&argv[i][2]);
  201.                     break;
  202.                 case 'N':
  203.                     sscanf (&argv[i][2], "%d/%d", &n_cols, &t_col);
  204.                     break;
  205.                 case 'Q':
  206.                     q_factor = atof(&argv[i][2]);
  207.                     check_q = TRUE;
  208.                     break;
  209.                 case 'S':
  210.                     check_asym = TRUE;
  211.                     sym_coeff = atof (&argv[i][2]);
  212.                     break;
  213.                 case 'T':
  214.                     sscanf (&argv[i][2], "%lf/%lf/%lf", &t_start, &t_stop, &t_int);
  215.                     out_at_time = TRUE;
  216.                     if (t_int == 0.0) error = TRUE;
  217.                     break;
  218.                 default:
  219.                     error = TRUE;
  220.                                         gmt_default_error (argv[i][1]);
  221.                     break;
  222.             }
  223.         }
  224.         else {
  225.             if ((fp_in = fopen(argv[i],"r")) == NULL) {
  226.                 fprintf (stderr, "filter1d: Cannot open file %s\n", argv[i]);
  227.                 exit(-1);
  228.             }
  229.         }
  230.     }
  231.  
  232.     if (argc == 1 || gmt_quick) {
  233.         fprintf (stderr, "filter1d %s - Time domain filtering of 1-D time series\n\n", GMT_VERSION);
  234.         fprintf (stderr, "usage: filter1d [infile] -F<type><width> [-D<increment>] [-E]\n");
  235.         fprintf (stderr, "\t[-H] [-I<ignore_val>] [-L<lack_width>] [-N<n_cols>/<t_col>]\n");
  236.         fprintf (stderr, "\t[-Q<q_factor>] [-S<symmetry_factor>] [-T<start>/<stop>/<int>] [-V] [-b[i|o]i[d]]\n\n");
  237.         
  238.         if (gmt_quick) exit (-1);
  239.         
  240.         
  241.         fprintf (stderr, "\tfilter1d is a general time domain filter for multiple time series data.\n");
  242.         fprintf (stderr, "\t\tThe user specifies the number of columns of input and which column is the time.\n");
  243.         fprintf (stderr, "\t\t(See N option below).  The fastest operation occurs when the input time series are\n");
  244.         fprintf (stderr, "\t\tequally spaced and have no gaps or outliers and the special options are not needed.\n");
  245.         fprintf (stderr, "\t\tfilter1d has options L,Q,S for unevenly sampled data with gaps and outliers.\n\n");
  246.  
  247.         fprintf (stderr, "\t-F sets Filtertype, type is one of b(oxcar), c(osine arch), g(aussian), m(edian),\n");
  248.         fprintf (stderr, "\t\tor p(maximum liklihood Probability estimator -- a mode estimator),\n");
  249.         fprintf (stderr, "\t\tand specify full filter <width> in same units as time column,\n");
  250.         fprintf (stderr, "\t\tOR, use Ff<name> to give the name of a one-column file of your own coefficients.\n");
  251.         fprintf (stderr, "\t\tUpper case type B, C, G, M, P, F will use robust filter versions:\n");
  252.         fprintf (stderr, "\t\ti.e., replace outliers (2.5 L1 scale off median) with median during filtering.\n");
  253.         fprintf (stderr, "\n\tOPTIONS:\n");
  254.         
  255.         fprintf (stderr, "\t-D<increment> is used when series is NOT equidistantly sampled.\n");
  256.         fprintf (stderr, "\t\tThen <increment> will be the abscissae resolution, i.e. all abscissae\n");
  257.         fprintf (stderr, "\t\twill be rounded off to a multiple of <increment>.\n");
  258.         fprintf (stderr, "\t-E include Ends of time series in output.  Default loses half_width at each end.\n");
  259.         explain_option ('H');
  260.         fprintf (stderr, "\t-I to ignore values; If an input value == <ignore_val> it will be set to NaN.\n");
  261.         fprintf (stderr, "\t-L<width> checks for Lack of data condition.  If input data has a gap exceeding\n");
  262.         fprintf (stderr, "\t\t<width> then no output will be given at that point.  Default does not check Lack.\n");
  263.         fprintf (stderr, "\t-N<n_cols>/<t_col> sets # columns in input and which column contains the independent\n");
  264.         fprintf (stderr, "\t\tvariable (time). The left-most column is # 0, the right-most is # (<n_cols> - 1).\n");
  265.         fprintf (stderr, "\t\tDefault is <n_cols> = 2, <t_col> = 0; i.e. file has t, f(t) pairs.\n");
  266.         fprintf (stderr, "\t-Q<q_factor> assess Quality of output value by checking mean weight in convolution.\n");
  267.         fprintf (stderr, "\t\tEnter <q_factor> between 0 and 1.  If mean weight < q_factor, output is suppressed\n");
  268.         fprintf (stderr, "\t\tat this point.  Default does not check Quality.\n");
  269.         fprintf (stderr, "\t-S<symmetry_factor> checks symmetry of data about window center.  Enter a factor\n");
  270.         fprintf (stderr, "\t\tbetween 0 and 1.  If ( (abs(n_left - n_right)) / (n_left + n_right) ) > factor,\n");
  271.         fprintf (stderr, "\t\tthen no output will be given at this point.  Default does not check Symmetry.\n");
  272.         fprintf (stderr, "\t-T make evenly spaced timesteps from <start> to <stop> by <int>. Default uses input times.\n");
  273.         explain_option ('V');
  274.         fprintf (stderr, "\t-b means binary (double) i/o.  Append i or o if only input OR output is binary\n");
  275.         fprintf (stderr, "\t   Append d for double precision [Default is single precision].\n");        explain_option ('.');
  276.         exit (-1);
  277.     }
  278.  
  279.     /* Check arguments */
  280.  
  281.     if (out_at_time && (t_stop - t_start) < filter_width) {
  282.         fprintf (stderr, "%s: GMT SYNTAX ERROR -T option:  Output interval < filterwidth\n", gmt_program);
  283.         error++;
  284.     }
  285.     if (check_lack && (lack_width < 0.0 || lack_width > filter_width) ) {
  286.         fprintf (stderr, "%s: GMT SYNTAX ERROR -L option:  Unreasonable lack-of-data interval\n", gmt_program);
  287.         error++;
  288.     }
  289.     if (check_asym && (sym_coeff < 0.0 || sym_coeff > 1.0) ) {
  290.         fprintf (stderr, "%s: GMT SYNTAX ERROR -S option:  Enter a factor between 0 and 1\n", gmt_program);
  291.         error++;
  292.     }
  293.     if (filter_type < 0 || filter_type >= N_FILTERS) {
  294.         fprintf (stderr, "%s: GMT SYNTAX ERROR -F option.  Correct syntax:\n", gmt_program);
  295.         fprintf (stderr, "-FX<width>, X one of BbCcGgMmPpFf\n");
  296.         error++;
  297.     }
  298.     if (t_col >= n_cols) {
  299.         fprintf (stderr, "%s: GMT SYNTAX ERROR -N option.  time column exceeds number of columns\n", gmt_program);
  300.         error++;
  301.     }
  302.     if (check_q && q_factor < 0.0) {
  303.         fprintf (stderr, "%s: GMT SYNTAX ERROR -Q option:  Enter a factor between 0 and 1\n", gmt_program);
  304.         error++;
  305.     }
  306.     if (b_in && gmtdefs.io_header) {
  307.         fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data cannot have header -H\n", gmt_program);
  308.         error++;
  309.     }
  310.     
  311.     if (error) exit (-1);
  312.     
  313.     if (filter_type > CONVOLVE) robust = FALSE;
  314.  
  315.     if (b_in) gmt_input = (single_precision) ? bin_float_input : bin_double_input;
  316.     if (b_out) gmt_output = (single_precision) ? bin_float_output : bin_double_output;
  317.  
  318.     if (allocate_space () ) {
  319.         fprintf(stderr,"filter1d: fatal error memory allocation.\n");
  320.         free_space ();
  321.         exit(-1);
  322.     }
  323.  
  324.     if (load_data_and_check_extrema () ) {
  325.         fprintf(stderr,"filter1d: fatal error during data read.\n");
  326.         free_space ();
  327.         exit(-1);
  328.     }
  329.  
  330.     if (set_up_filter () ) {
  331.         fprintf(stderr,"filter1d: fatal error during coefficient setup.\n");
  332.         free_space ();
  333.         exit(-1);
  334.     }
  335.  
  336.     if (do_the_filter () ) {
  337.         fprintf(stderr,"filter1d: fatal error in filtering routine.\n");
  338.         free_space ();
  339.         exit(-1);
  340.     }
  341.  
  342.     free_space ();
  343.     
  344.     gmt_end (argc, argv);
  345. }
  346.     
  347. int load_data_and_check_extrema ()
  348. {    
  349.  
  350.     int    i, more, n_fields;
  351.     double    last_time, new_time, *in;
  352.  
  353.     if (fp_in == NULL) fp_in = stdin;
  354.     
  355.     if (gmtdefs.io_header) for (i = 0; i < gmtdefs.n_header_recs; i++) {
  356.         fgets (buffer, 512, fp_in);
  357.         fputs (buffer, stdout);
  358.     }
  359.  
  360.     n_rows = 0;
  361.     if (robust || (filter_type == MEDIAN) ) {
  362.         for (i = 0; i < n_cols; i++) {
  363.             min_loc[i] = MAXDOUBLE;
  364.             max_loc[i] = -MAXDOUBLE;
  365.         }
  366.     }
  367.     last_time = -MAXDOUBLE;
  368.     
  369.     more = ((n_fields = gmt_input (fp_in,  n_cols, &in)) == n_cols);
  370.     
  371.     while (more) {
  372.     
  373.         for (i = 0; i < n_cols; i++) {
  374.             if (ignore_test && in[i] == ignore_val)
  375.                 data[i][n_rows] = gmt_NaN;
  376.             else {
  377.                 data[i][n_rows] = in[i];
  378.                 if (robust || (filter_type == MEDIAN) ) {
  379.                     if (in[i] > max_loc[i]) max_loc[i] = in[i];
  380.                     if (in[i] < min_loc[i]) min_loc[i] = in[i];
  381.                 }
  382.             }
  383.         }
  384.         
  385.         if (!bad_float( (double)data[t_col][n_rows]) ) {    /* Don't n_rows++ if time is bad this row  */
  386.         
  387.             new_time = data[t_col][n_rows];
  388.             n_rows++;
  389.             if (new_time < last_time) {
  390.                 fprintf(stderr,"filter1d: Error! Time decreases at line # %d\n", n_rows);
  391.                 fprintf(stderr,"\tUse UNIX utility sort and then try again.\n");
  392.                 return(-1);
  393.             }
  394.             last_time = new_time;
  395.         }
  396.         
  397.         if (n_rows == n_row_alloc) {    /* Need more memory */
  398.             n_row_alloc += GMT_CHUNK;
  399.             allocate_more_data_space ();
  400.         }
  401.         
  402.         more = ((n_fields = gmt_input (fp_in,  n_cols, &in)) == n_cols);
  403.     }
  404.     
  405.     if (fp_in != stdin) fclose (fp_in);
  406.     
  407.     if (n_fields == -1) n_fields = 0;    /* Ok to get EOF */
  408.     if (n_fields%n_cols) {    /* Got garbage or multiple segment subheader */
  409.         fprintf (stderr, "%s:  Cannot read line %d, aborts\n", gmt_program, n_rows);
  410.         exit (-1);
  411.     }
  412.     
  413.     for (i = 0; i < n_cols; i++) data[i] = (double *) memory ((char *)data[i], n_rows, sizeof (double), "filter1d");
  414.     if (gmtdefs.verbose) fprintf (stderr, "filter1d: Read %d data lines\n", n_rows);
  415.  
  416.     /* Initialize scale parameters and last_loc based on min and max of data  */
  417.     
  418.     if (robust || (filter_type == MEDIAN) ) {
  419.         for (i = 0; i < n_cols; i++) {
  420.             min_scl[i] = 0.0;
  421.             max_scl[i] = 0.5 * (max_loc[i] - min_loc[i]);
  422.             last_scl[i] = 0.5 * max_scl[i];
  423.             last_loc[i] = 0.5 * (max_loc[i] + min_loc[i]);
  424.         }
  425.     }
  426.  
  427.     return (0);
  428. }
  429.  
  430. int    set_up_filter()
  431. {
  432.     int    i, i1, i2;
  433.     double    t_0, t_1, time, w_sum;
  434.     BOOLEAN    normalize = FALSE;
  435.     
  436.     
  437.     
  438.     t_0 = data[t_col][0];
  439.     t_1 = data[t_col][n_rows-1];
  440.     if (equidist) dt = (t_1 - t_0) / (n_rows - 1);
  441.     
  442.  
  443.     if (fp_wt != NULL) {
  444.         f_wt = (double *) memory (CNULL, n_work_alloc, sizeof(double), "filter1d");
  445.         n_f_wts = 0;
  446.         while (fscanf (fp_wt,"%lf\n", &f_wt[n_f_wts]) != EOF) {
  447.             n_f_wts++;
  448.             if (n_f_wts == n_work_alloc) {    /* Need more memory */
  449.                 n_work_alloc += GMT_CHUNK;
  450.                 allocate_more_work_space ();
  451.             }
  452.         }
  453.         fclose (fp_wt);
  454.         half_n_f_wts = n_f_wts / 2;
  455.         half_width = half_n_f_wts * dt;
  456.         f_wt = (double *) memory ((char *)f_wt, n_f_wts, sizeof (double), "filter1d");
  457.         if (gmtdefs.verbose) fprintf (stderr, "filter1d: Read %d filter weights from file.\n", n_f_wts);
  458.     }
  459.     else if (filter_type <= CONVOLVE) {
  460.  
  461.         get_weight[BOXCAR] = (PFD)boxcar_weight;
  462.         get_weight[COS_ARCH] = (PFD)cosine_weight;
  463.         get_weight[GAUSSIAN] = (PFD)gaussian_weight;
  464.         half_width = 0.5 * filter_width;
  465.         half_n_f_wts = floor (half_width / dt);
  466.         n_f_wts = 2 * half_n_f_wts + 1;
  467.         
  468.         f_wt = (double *) memory (CNULL, n_f_wts, sizeof(double), "filter1d");
  469.         for (i = 0; i <= half_n_f_wts; i++) {
  470.             time = i * dt;
  471.             i1 = half_n_f_wts - i;
  472.             i2 = half_n_f_wts + i;
  473.             f_wt[i1] = f_wt[i2] = ( *get_weight[filter_type]) (time, half_width);
  474.         }
  475.         if (normalize) {
  476.             w_sum = 0.0;
  477.             for (i = 0; i < n_f_wts; i++) w_sum += f_wt[i];
  478.             for (i = 0; i < n_f_wts; i++) f_wt[i] /= w_sum;
  479.         }
  480.     }
  481.     else {
  482.         half_width = 0.5 * filter_width;
  483.     }
  484.  
  485.     /* Initialize start/stop time */
  486.  
  487.     if (out_at_time) {
  488.         
  489.         if (use_ends) {
  490.             while (t_start < t_0) t_start += t_int;
  491.             while (t_stop > t_1) t_stop -= t_int;
  492.         }
  493.         else {
  494.             while ( (t_start - half_width) < t_0) t_start += t_int;
  495.             while ( (t_stop + half_width) > t_1) t_stop -= t_int;
  496.         }
  497.     }
  498.     else {
  499.         if (use_ends) {
  500.             t_start = t_0;
  501.             t_stop = t_1;
  502.         }
  503.         else {
  504.             for (i = 0; (data[t_col][i] - t_0) < half_width; i++);
  505.             t_start = data[t_col][i];
  506.             for (i = n_rows - 1; (t_1 - data[t_col][i]) < half_width; i--);
  507.             t_stop = data[t_col][i];
  508.         }
  509.     }
  510.  
  511.     if (gmtdefs.verbose) fprintf (stderr, "F width: %lg Resolution: %lg Start: %lg Stop: %lg\n",
  512.         filter_width, dt, t_start, t_stop);
  513.     
  514.     return(0);
  515. }
  516.  
  517. int    do_the_filter ()
  518. {
  519.  
  520.     int    i_row, i_col, left, right, n_l, n_r, iq, i_f_wt;
  521.     int    i_t_output, n_in_filter, n_for_call, n_good_ones;
  522.     double    time, delta_time, *outval, wt, val, med, scl;
  523.  
  524.     outval = (double *)memory (CNULL, n_cols, sizeof (double), "filter1d");
  525.     
  526.     if(!out_at_time) {    /* Position i_t_output at first output time  */
  527.         for(i_t_output = 0; data[t_col][i_t_output] < t_start; i_t_output++);
  528.     }
  529.     time = t_start;
  530.     left = right = 0;        /* Left/right end of filter window */
  531.  
  532.     iq = rint(q_factor);
  533.  
  534.     while (time <= t_stop) {
  535.         while ((time - data[t_col][left]) > half_width) left++;
  536.         while (right < n_rows && (data[t_col][right] - time) <= half_width) right++;
  537.         n_in_filter = right - left;
  538.         if ( (!(n_in_filter)) || (check_lack && ( (filter_width / n_in_filter) > lack_width) ) ) {
  539.             if (out_at_time)
  540.                 time += t_int;
  541.             else {
  542.                 i_t_output++;
  543.                 time = (i_t_output < n_rows) ? data[t_col][i_t_output] : t_stop + 1.0;
  544.             }
  545.             continue;
  546.         }
  547.             
  548.         for (i_col = 0; i_col < n_cols; i_col++) {
  549.             n_this_col[i_col] = 0;
  550.             wt_sum[i_col] = 0.0;
  551.             data_sum[i_col] = 0.0;
  552.             if (i_col == t_col) {
  553.                 good_one[i_col] = FALSE;
  554.             }
  555.             else if (check_lack) {
  556.                 good_one[i_col] = !(lack_check(i_col, left, right));
  557.             }
  558.             else {
  559.                 good_one[i_col] = TRUE;
  560.             }
  561.             if (check_asym) {
  562.                 n_left[i_col] = 0;
  563.                 n_right[i_col] = 0;
  564.             }
  565.         }
  566.             
  567.         if (robust || filter_type > CONVOLVE) {
  568.             if (n_in_filter > n_work_alloc) {
  569.                 n_work_alloc = n_in_filter;
  570.                 allocate_more_work_space ();
  571.             }
  572.             for (i_row = left; i_row < right; i_row++) {
  573.                 for (i_col = 0; i_col < n_cols; i_col++) {
  574.                     if (!(good_one[i_col])) continue;
  575.                     if (!bad_float( (double)data[i_col][i_row])) {
  576.                         work[i_col][n_this_col[i_col]] = data[i_col][i_row];
  577.                         n_this_col[i_col]++;
  578.                         if (check_asym) {
  579.                             if (data[t_col][i_row] < time) n_left[i_col]++;
  580.                             if (data[t_col][i_row] > time) n_right[i_col]++;
  581.                         }
  582.                     }
  583.                 }
  584.             }
  585.             if (check_asym) {
  586.                 for (i_col = 0; i_col < n_cols; i_col++) {
  587.                     if (!(good_one[i_col])) continue;
  588.                     n_l = n_left[i_col];
  589.                     n_r = n_right[i_col];
  590.                     if ((((double)abs(n_l - n_r))/(n_l + n_r)) > sym_coeff) good_one[i_col] = FALSE;
  591.                 }
  592.             }
  593.             if ( (filter_type > CONVOLVE) && check_q) {
  594.                 for (i_col = 0; i_col < n_cols; i_col++) {
  595.                     if (n_this_col[i_col] < iq) good_one[i_col] = FALSE;
  596.                 }
  597.             }
  598.                 
  599.             for (i_col = 0; i_col < n_cols; i_col++) {
  600.                 if (good_one[i_col]) {
  601.                     n_for_call = n_this_col[i_col];
  602.                     get_robust_estimates (i_col, n_for_call, robust);
  603.                 }
  604.             }
  605.                 
  606.         }    /* That's it for the robust work  */
  607.  
  608.         if (filter_type > CONVOLVE) {
  609.             
  610.             /* Need to count how many good ones; use data_sum area  */
  611.                 
  612.             n_good_ones = 0;
  613.             for (i_col = 0; i_col < n_cols; i_col++) {
  614.                 if (i_col == t_col) {
  615.                     data_sum[i_col] = time;
  616.                 }
  617.                 else if (good_one[i_col]) {
  618.                     data_sum[i_col] = this_loc[i_col];
  619.                     n_good_ones++;
  620.                 }
  621.                 else {
  622.                     data_sum[i_col] = gmt_NaN;
  623.                 }
  624.             }
  625.             if (n_good_ones) gmt_output (stdout, n_cols, data_sum);
  626.         }
  627.             
  628.         else {
  629.             
  630.             if (robust) for (i_col = 0; i_col < n_cols; i_col++) n_this_col[i_col] = 0;
  631.                         
  632.             for (i_row = left; i_row < right; i_row++) {
  633.                 delta_time = time - data[t_col][i_row];
  634.                 i_f_wt = half_n_f_wts + floor(0.5 + delta_time/dt);
  635.                 if ( (i_f_wt < 0) || (i_f_wt >= n_f_wts) ) continue;
  636.                     
  637.                 for(i_col = 0; i_col < n_cols; i_col++) {
  638.                     if ( !(good_one[i_col]) ) continue;
  639.                     if (!bad_float( (double)data[i_col][i_row])) {
  640.                         wt = f_wt[i_f_wt];
  641.                         val = data[i_col][i_row];
  642.                         if (robust) {
  643.                             med = this_loc[i_col];
  644.                             scl = this_scl[i_col];
  645.                             val = ((fabs(val-med)) > (2.5 * scl)) ? med : val;
  646.                         }
  647.                         else if (check_asym) {    /* This wasn't already done  */
  648.                             if (data[t_col][i_row] < time) n_left[i_col]++;
  649.                             if (data[t_col][i_row] > time) n_right[i_col]++;
  650.                         }
  651.                         wt_sum[i_col] += wt;
  652.                         data_sum[i_col] += (wt * val);
  653.                         n_this_col[i_col]++;
  654.                     }
  655.                 }
  656.             }
  657.             n_good_ones = 0;
  658.             for (i_col = 0; i_col < n_cols; i_col++) {
  659.                 if ( !(good_one[i_col]) ) continue;
  660.                 if ( !(n_this_col[i_col]) ) {
  661.                     good_one[i_col] = FALSE;
  662.                     continue;
  663.                 }
  664.                 if (check_asym && !(robust) ) {
  665.                     n_l = n_left[i_col];
  666.                     n_r = n_right[i_col];
  667.                     if ((((double)abs(n_l - n_r))/(n_l + n_r)) > sym_coeff) {
  668.                         good_one[i_col] = FALSE;
  669.                         continue;
  670.                     }
  671.                 }
  672.                 if (check_q && ((wt_sum[i_col] / n_this_col[i_col]) < q_factor)) {
  673.                     good_one[i_col] = FALSE;
  674.                     continue;
  675.                 }
  676.                 n_good_ones++;
  677.             }
  678.             if (n_good_ones) {
  679.                 for (i_col = 0; i_col < n_cols; i_col++) {
  680.                     if (i_col == t_col)
  681.                         outval[i_col] = time;
  682.                     else if (good_one[i_col])
  683.                         outval[i_col] = data_sum[i_col] / wt_sum[i_col];
  684.                     else
  685.                         outval[i_col] = gmt_NaN;
  686.                 }
  687.                 gmt_output (stdout, n_cols, outval);
  688.             }
  689.         }
  690.  
  691.         /* Go to next output time */
  692.  
  693.         if (out_at_time)
  694.             time += t_int;
  695.         else {
  696.             i_t_output++;
  697.             time = (i_t_output < n_rows) ? data[t_col][i_t_output] : t_stop + 1.0;
  698.         }
  699.     }
  700.     
  701.     free ((char *)outval);
  702.     
  703.     return(0);
  704. }
  705.  
  706. int    allocate_space ()
  707. {
  708.  
  709.     int    i;
  710.  
  711.     data = (double **) memory (CNULL, n_cols, sizeof(double *), "filter1d");
  712.     for (i = 0; i < n_cols; i++) data[i] = (double *) memory (CNULL, n_row_alloc, sizeof(double), "filter1d");
  713.     
  714.     wt_sum = (double *) memory (CNULL, n_cols, sizeof(double), "filter1d");
  715.     data_sum = (double *) memory (CNULL, n_cols, sizeof(double), "filter1d");
  716.     n_this_col = (int *) memory (CNULL, n_cols, sizeof(int), "filter1d");
  717.     good_one = (int *) memory (CNULL, n_cols, sizeof(int), "filter1d");
  718.  
  719.     if (check_asym) n_left = (int *) memory (CNULL, n_cols, sizeof(int), "filter1d");
  720.     if (check_asym) n_right = (int *) memory (CNULL, n_cols, sizeof(int), "filter1d");
  721.     
  722.     if (robust || (filter_type > CONVOLVE) ) {    /* Then we need workspace  */
  723.  
  724.         work = (double **) memory (CNULL, n_cols, sizeof(double *), "filter1d");
  725.         for (i = 0; i < n_cols; i++) work[i] = (double *) memory (CNULL, n_work_alloc, sizeof(double), "filter1d");
  726.         min_loc = (double *) memory (CNULL, n_cols, sizeof(double), "filter1d");
  727.         max_loc = (double *) memory (CNULL, n_cols, sizeof(double), "filter1d");
  728.         last_loc = (double *) memory (CNULL, n_cols, sizeof(double), "filter1d");
  729.         this_loc = (double *) memory (CNULL, n_cols, sizeof(double), "filter1d");
  730.         min_scl = (double *) memory (CNULL, n_cols, sizeof(double), "filter1d");
  731.         max_scl = (double *) memory (CNULL, n_cols, sizeof(double), "filter1d");
  732.         this_scl = (double *) memory (CNULL, n_cols, sizeof(double), "filter1d");
  733.         last_scl = (double *) memory (CNULL, n_cols, sizeof(double), "filter1d");
  734.     }
  735.     return (0);
  736. }
  737.  
  738. int    allocate_more_data_space ()
  739. {
  740.     int    i;
  741.  
  742.     for (i = 0; i < n_cols; i++) data[i] = (double *) memory ((char *)data[i], n_row_alloc, sizeof(double), "filter1d");
  743.     
  744.     return (0);
  745. }
  746.  
  747. int    allocate_more_work_space ()
  748. {
  749.     int    i;
  750.  
  751.     for (i = 0; i < n_cols; i++) work[i] = (double *) memory ((char *)work[i], n_work_alloc, sizeof(double), "filter1d");
  752.  
  753.     return (0);
  754. }
  755.  
  756.  
  757. int    free_space ()
  758. {
  759.     free( (char *)data);
  760.     free( (char *)wt_sum);
  761.     free( (char *)data_sum);
  762.     free( (char *)n_this_col);
  763.     free( (char *)good_one);
  764.     free( (char *)n_left);
  765.     free( (char *)n_right);
  766.     free( (char *)work);
  767.     free( (char *)min_loc);
  768.     free( (char *)max_loc);
  769.     free( (char *)last_loc);
  770.     free( (char *)this_loc);
  771.     free( (char *)min_scl);
  772.     free( (char *)max_scl);
  773.     free( (char *)last_scl);
  774.     free( (char *)this_scl);
  775.     return (0);
  776. }
  777.  
  778. BOOLEAN    lack_check (i_col, left, right)
  779. int    i_col, left, right;
  780. {            
  781.  
  782.     int    last_row, this_row;
  783.     BOOLEAN    lacking = FALSE;
  784.     double    last_t;
  785.     
  786.     last_row = left;
  787.     while (!(bad_float( (double)data[i_col][last_row])) && last_row < (right - 1)) last_row++;
  788.     
  789.     last_t = data[t_col][last_row];
  790.     this_row = last_row + 1;
  791.     while ( !(lacking) && this_row < (right - 1)) {
  792.         while (!(bad_float( (double)data[i_col][this_row])) && this_row < (right - 1)) this_row++;
  793.         
  794.         if ( (data[t_col][this_row] - last_t) > lack_width)
  795.             lacking = TRUE;
  796.         else {
  797.             last_t = data[t_col][this_row];
  798.             last_row = this_row;
  799.             this_row++;
  800.         }
  801.     }
  802.     return (lacking);
  803. }
  804.  
  805. int    get_robust_estimates (j, n, both)
  806. int    j, n, both;
  807. {
  808.     int i, n_smooth, sort_me = TRUE;
  809.     double    low, high, last, temp;
  810.     
  811.     if (filter_type == MODE) {
  812.  
  813.         n_smooth = n / 2;
  814.     
  815.         mode(work[j], n, n_smooth, sort_me, &temp);
  816.     }
  817.     else {
  818.  
  819.         low = min_loc[j];
  820.         high = max_loc[j];
  821.         last = last_loc[j];
  822.     
  823.         median(work[j], n, low, high, last, &temp);
  824.     }
  825.  
  826.     last_loc[j] = this_loc[j] = temp;
  827.     
  828.     if (both) {
  829.         for (i = 0; i < n; i++)
  830.             work[j][i] = fabs(work[j][i] - this_loc[j]);
  831.         low = min_scl[j];
  832.         high = max_scl[j];
  833.         last = last_scl[j];
  834.         median(work[j], n, low, high, last, &temp);
  835.         last_scl[j] = this_scl[j] = temp;
  836.     } 
  837.     return(0);
  838. }
  839.     
  840. double    boxcar_weight (radius, half_width)
  841. double    radius, half_width;
  842. {
  843.     double weight;
  844.     
  845.     if (radius > half_width)
  846.         weight = 0.0;
  847.     else
  848.         weight = 1.0;
  849.         
  850.     return (weight);
  851. }
  852.  
  853. double    cosine_weight (radius, half_width)
  854. double    radius, half_width;
  855. {
  856.     double weight, cosine_constant;
  857.     cosine_constant = M_PI / half_width;
  858.     if (radius > half_width)
  859.         weight = 0.0;
  860.     else
  861.         weight = 1.0 + cos (radius * cosine_constant);
  862.         
  863.     return (weight);
  864. }
  865.  
  866. double    gaussian_weight (radius, half_width)
  867. double    radius, half_width;
  868. {
  869.     double weight, gauss_constant;
  870.     gauss_constant = -4.5 / (half_width * half_width);
  871.     if (radius > half_width)
  872.         weight = 0.0;
  873.     else
  874.         weight = exp (radius * radius * gauss_constant);
  875.         
  876.     return (weight);
  877. }
  878.