/*
htmedfilt.c		Joe Mietus		Dec 19 2001

-------------------------------------------------------------------------------
htmedfilt: Separately median filter Hilbert amplitudes and frequencies
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/).
_______________________________________________________________________________

Separately median filter Hilbert amplitudes and frequencies

Usage : htmedfilt win

*/

#include <stdio.h>

#define MAXWIN 2000

static double time[MAXWIN], x[MAXWIN], y[MAXWIN], sx[MAXWIN], sy[MAXWIN];


main(argc, argv)
int argc;
char *argv[];
{
    int i, j, k, win, hwin;
    int dblcmp();

    if (argc != 2) {
        usage(argv[0]);
	exit(1);
    }

    if ((win = atoi(argv[1])) < 3) {
        usage(argv[0]);
        exit(1);
    }
    if (win > MAXWIN) {
	fprintf(stderr, "%s : maximum window size is MAXWIN\n", argv[0]);
	exit(2);
    }

    for (i=0; i<win && scanf("%lf %lf %lf", &time[i], &x[i], &y[i]) == 3; i++)
        ;

    if (i < win) {
        fprintf(stderr, "%s : not enough points in input\n", argv[0]);
	exit(2);
    }
    if (++i >= win)
	i = 0;
    j = hwin = win/2 -1;

    for (k=0; k<win; k++) {
        sx[k] = x[k];
        sy[k] = y[k];
    }

    qsort(sx, win, sizeof(double), dblcmp);
    qsort(sy, win, sizeof(double), dblcmp);
    printf("%g %g %g\n", time[j], sx[hwin], sy[hwin]);

    while (scanf("%lf %lf %lf", &time[i], &x[i], &y[i]) == 3) {

	if (++i >= win)
	    i = 0;
	if (++j >= win)
	    j = 0;

        for (k=0; k<win; k++) {
            sx[k] = x[k];
            sy[k] = y[k];
        }

        qsort(sx, win, sizeof(double), dblcmp);
        qsort(sy, win, sizeof(double), dblcmp);
	printf("%g %g %g\n", time[j], sx[hwin], sy[hwin]);

    }

}


int dblcmp(y1, y2)
double *y1, *y2;
{

    if (*y1 < *y2)
        return(-1);
    else if (*y1 > *y2)
        return(1);
    else
        return(0);

}


usage(prog)
char * prog;
{
	fprintf(stderr, "Usage : %s  win\n\n", prog);
	fprintf(stderr, " Separately median filter Hilbert amplitudes and\n");
	fprintf(stderr, " frequencies with a sliding window hwin points wide.\n");

}