/*
htavsd.c Joe Mietus Dec 19 2001
-------------------------------------------------------------------------------
htavsd: compute Hilbert transform averages, standard deviations and threshold
crossings for amplitude and frequency over moving window
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/).
_______________________________________________________________________________
get Hilbert transform averages, standard deviations and threshold crossings
for amplitude and frequency over moving window of length `win'
Usage : htavsd incr win ampthres
Reads stdin of time in secs, amp and freq
outputs time, av, sd and % time within threshold limits
for both amp and freq
incr : output for every 'incr' time steps
win : window length
ampthres : min amp
*/
#include <stdio.h>
#include <math.h>
#define MAXFREQ 0.06
main(argc, argv)
int argc;
char *argv[];
{
int i, j, k, win, start, incr, ydet, zdet;
double ampthres;
double sumy, sumz, sumyy, sumzz, avy, avz, sdy, sdz;
double *x, *y, *z;
double atof();
if (argc < 4) {
usage(argv[0]);
exit(1);
}
if ((incr = strtim(argv[1])) <= 0) {
fprintf(stderr, "%0 : incr must be greater than 0\n", argv[0]);
exit(1);
}
if ((win = strtim(argv[2])) <= 0) {
fprintf(stderr, "%0 : win must be greater than 0\n", argv[0]);
exit(1);
}
if ((ampthres = atof(argv[3])) <= 0) {
fprintf(stderr, "%0 : ampthres must be greater than 0\n", argv[0]);
exit(1);
}
if (incr > win) {
fprintf(stderr, "%s : incr must be <= win\n", argv[0]);
exit(1);
}
if ((x = (double *)calloc(win, sizeof(double))) == NULL ||
(y = (double *)calloc(win, sizeof(double))) == NULL ||
(z = (double *)calloc(win, sizeof(double))) == NULL) {
fprintf(stderr, "%s : insufficient memory\n", argv[0]);
exit(1);
}
ydet = zdet = 0;
sumy = sumz = sumyy = sumzz = 0.;
if (scanf("%lf %lf %lf", &x[0], &y[0], &z[0]) != 3) {
fprintf(stderr, "%s : no data read\n", argv[0]);
exit(2);
}
start = x[0];
sumy = y[0];
sumz = z[0];
sumyy = y[0]*y[0];
sumzz = z[0]*z[0];
if (y[0] >= ampthres)
ydet++;
if (z[0] <= MAXFREQ)
zdet++;
for (i=1; i<win && scanf("%lf %lf %lf", &x[i], &y[i], &z[i]) == 3; i++) {
sumy += y[i];
sumz += z[i];
sumyy += y[i]*y[i];
sumzz += z[i]*z[i];
if (y[i] >= ampthres)
ydet++;
if (z[i] <= MAXFREQ)
zdet++;
}
if (i < win) {
fprintf(stderr, "not enough points in input\n");
exit(2);
}
avy = sumy/win;
avz = sumz/win;
sdy = sqrt((sumyy - sumy*sumy/win)/(win-1));
sdz = sqrt((sumzz - sumz*sumz/win)/(win-1));
printf("%02d:%02d:%02d ", start/3600, (start%3600)/60, start%60);
printf("%f %f %f %f %f %f\n",
avy, sdy, ((double)ydet)/win,
avz, sdz, ((double)zdet)/win);
if (incr == win) {
sumy = sumz =sumyy = sumzz = 0.0;
ydet = zdet = 0;
}
else {
for (j=0; j<incr; j++) {
if (j >= win)
j = 0;
sumy -= y[j];
sumz -= z[j];
sumyy -= y[j]*y[j];
sumzz -= z[j]*z[j];
if (y[j] >= ampthres)
ydet--;
if (z[j] <= MAXFREQ)
zdet--;
}
}
i = 0;
start += incr;
while (scanf("%lf %lf %lf", &x[i], &y[i], &z[i]) == 3) {
while (x[i] < start)
if (scanf("%lf %lf %lf", &x[i], &y[i], &z[i]) != 3)
exit(0);
sumy += y[i];
sumz += z[i];
sumyy += y[i]*y[i];
sumzz += z[i]*z[i];
if (y[i] >= ampthres)
ydet++;
if (z[i] <= MAXFREQ)
zdet++;
if (++i >= win)
i = 0;
for (j=1; j<incr && j<win; i++, j++) {
if (i >= win)
i = 0;
if (scanf("%lf %lf %lf", &x[i], &y[i], &z[i]) != 3)
exit(0);
sumy += y[i];
sumz += z[i];
sumyy += y[i]*y[i];
sumzz += z[i]*z[i];
if (y[i] >= ampthres)
ydet++;
if (z[i] <= MAXFREQ)
zdet++;
}
if (i >= win)
i = 0;
avy = sumy/win;
avz = sumz/win;
sdy = sqrt((sumyy - sumy*sumy/win)/(win-1));
sdz = sqrt((sumzz - sumz*sumz/win)/(win-1));
printf("%02d:%02d:%02d ", start/3600, (start%3600)/60, start%60);
printf("%f %f %f %f %f %f\n",
avy, sdy, ((double)ydet)/win,
avz, sdz, ((double)zdet)/win);
if (incr == win) {
sumy = sumz =sumyy = sumzz = 0.0;
ydet = zdet = 0;
}
else {
for (j=0, k=j+i; j<incr; j++, k++) {
if (k >= win)
k = 0;
sumy -= y[k];
sumz -= z[k];
sumyy -= y[k]*y[k];
sumzz -= z[k]*z[k];
if (y[k] >= ampthres)
ydet--;
if (z[k] <= MAXFREQ)
zdet--;
}
}
start += incr;
}
}
int strtim(buf) /* convert string in [[HH:]MM:]SS format to seconds */
char *buf;
{
int x, y, z;
switch (sscanf(buf, "%d:%d:%d", &x, &y, &z)) {
case 1: return (x);
case 2: return (60*x + y);
case 3: return (3600*x + 60*y + z);
default: return (-1);
}
}
usage(prog)
char *prog;
{
fprintf(stderr, "Usage : %s incr win ampthres\n\n", prog);
fprintf(stderr, " Reads stdin of time in secs, amp and freq\n");
fprintf(stderr, " outputs time, av, sd, and");
fprintf(stderr, " % time within threshold limits for both amp and freq\n\n");
fprintf(stderr, "incr] : output for every 'incr' time steps\n");
fprintf(stderr, "win : window length\n");
fprintf(stderr, "ampthres : minumum amplitude threshold\n");
}