/*
filt.c Joe Mietus Dec 19 2001
-------------------------------------------------------------------------------
filt: moving average signal filter
Copyright (C) 2002 Joe Mietus
This program is free software; you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free Software
Foundation; either version 2 of the License, or (at your option) any later
version.
This program is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with
this program; if not, write to the Free Software Foundation, Inc., 59 Temple
Place - Suite 330, Boston, MA 02111-1307, USA.
You may contact the author by e-mail (joe at physionet dot org) or postal mail
(MIT Room E25-505A, Cambridge, MA 02139 USA). For updates to this software,
please visit PhysioNet (http://www.physionet.org/).
_______________________________________________________________________________
filt filt hwin [options]
Reads 2 columns from stdin and filters outliers by
deleting those intervals outside of `filt' range of
the average within a window hwin' distance on either
side of the current interval.
options are :
[-x min max] : exclude date outside min - max
*/
#include <stdio.h>
#include <stdlib.h>
#define MAXWIN 1024
main(argc, argv)
int argc;
char * argv[];
{
int hwin, win, i, j, k;
int xflag;
double filt, filtmin, filtmax, min, max, mtmp, sum, av;
static double x[MAXWIN+1], y[MAXWIN+1], lastx, lasty;
FILE *fp;
xflag = 0;
sum = 0.0;
if (argc < 3) {
usage(argv[0]);
exit(1);
}
i = 1;
if ((filt = atof(argv[i])) <= 0) {
fprintf(stderr, "%s: filt must be greater than 0\n",
argv[0]);
exit(2);
}
if ((hwin = atoi(argv[++i])) <= 1) {
fprintf(stderr, "%s: hwin must be integer greater than 1\n",
argv[0]);
exit(2);
}
if ((win = 2*hwin) > MAXWIN) {
fprintf(stderr, "%s: max hwin = %d\n", argv[0], MAXWIN/2);
exit(2);
}
for ( ; ++i < argc && *argv[i] == '-'; ) {
switch(argv[i][1]) {
case 'x': if (i+3 > argc) {
usage();
exit(1);
}
min = atof(argv[++i]);
max = atof(argv[++i]);
if (min > max) {
mtmp = min;
min = max;
max = mtmp;
}
xflag = 1;
break;
default: usage(argv[0]);
exit(1);
}
}
for (i=0; i<=win; i++) {
if (scanf("%lf %lf", &x[i], &y[i]) != 2) {
fprintf(stderr, "%s : insufficient data\n", argv[0]);
exit(2);
}
if (xflag) {
while (xflag && (y[i]<min || y[i]>max)) {
if (scanf("%lf %lf", &x[i], &y[i]) != 2) {
fprintf(stderr, "%s : insufficient data\n", argv[0]);
exit(2);
}
}
}
}
i = 0;
j = hwin;
for (i=0; i<=win; i++) {
sum += y[i];
}
sum -= y[j];
av = sum/win;
sum += y[j] - y[0];
filtmax = filtmin = filt * av;
if (y[j] <= av+filtmax && y[j] >= av-filtmin) {
printf("%.9g %g\n", x[j], y[j]);
lastx = x[j];
lasty = y[j];
}
i = 0;
while (scanf("%lf %lf", &x[i], &y[i]) == 2) {
if (xflag) {
while (y[i]<min || y[i]>max) {
if (scanf("%lf %lf", &x[i], &y[i]) != 2) {
exit(0);
}
}
}
if (++j > win)
j = 0;
sum += y[i] - y[j];
av = sum/win;
if (++i > win)
i = 0;
sum += y[j] - y[i];
filtmax = filtmin = filt * av;
if (y[j] <= av+filtmax && y[j] >= av-filtmin) {
printf("%.9g %g\n", x[j], y[j]);
lastx = x[j];
lasty = y[j];
}
}
exit(0);
}
usage(prog)
char *prog;
{
fprintf(stderr, "Usage : %s filt hwin [options]\n\n", prog);
fprintf(stderr, " Reads 2 columns from stdin and filters outliers by\n");
fprintf(stderr, " deleting those intervals outside of `filt' range of\n");
fprintf(stderr, " the average within a window hwin' distance on either\n");
fprintf(stderr, " side of the current interval.\n\n");
fprintf(stderr, " options are :\n");
fprintf(stderr, " [-x min max] : exclude date outside min - max\n");
}