/*
ldetrend.c Joe Mietus Dec 19 2001
-------------------------------------------------------------------------------
ldetrend: Locally detrend signals
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/).
_______________________________________________________________________________
Reads 2 columns from stdin and locally detrends
by subtracting a least squares fitted line
over a sliding window 2*hwin+1 points wide.
Usage: ldetrend hwin
*/
#include <stdio.h>
#define MAXDAT 262144
main(argc, argv)
int argc;
char *argv[];
{
int i, n, hwin, win;
double a, b, sumx, sumy, sumxy, sumx2;
static double x[MAXDAT+1], y[MAXDAT+1];
if (argc != 2) {
usage(argv[0]);
exit(1);
}
if ((hwin = atoi(argv[1])) < 1) {
fprintf(stderr, "%s : hwin must be greater than 0\n", argv[0]);
exit(1);
}
win = 2*hwin+1;
sumx = sumy = sumxy = sumx2 = 0.;
for (n=0; n<MAXDAT && scanf("%lf %lf", &x[n], &y[n]) == 2; n++) {
sumx += x[n];
sumy += y[n];
sumxy += x[n]*y[n];
sumx2 += x[n]*x[n];
}
if (n == MAXDAT && scanf("%lf %lf", &x[i], &y[i]) == 2) {
fprintf(stderr, "%s : input exceeds buffer length, ", argv[0]);
fprintf(stderr, "truncating to %d points\n", MAXDAT);
}
if (win >= n) {
b = (sumxy - sumx*sumy/n) / (sumx2 - sumx*sumx/n);
a = sumy/n - b*sumx/n;
for (i=0; i<n; i++)
printf("%g %g\n", x[i], y[i] - (a + b*x[i]));
}
else {
sumx = sumy = sumxy = sumx2 = 0.;
for (i=0; i<win; i++) {
sumx += x[i];
sumy += y[i];
sumxy += x[i]*y[i];
sumx2 += x[i]*x[i];
}
b = (sumxy - sumx*sumy/win) / (sumx2 - sumx*sumx/win);
a = sumy/win - b*sumx/win;
for (i=0; i<=hwin; i++)
printf("%g %g\n", x[i], y[i] - (a + b*x[i]));
for (i=win ; i<n; i++) {
sumx += x[i]-x[i-win];
sumy += y[i]-y[i-win];
sumxy += x[i]*y[i]-x[i-win]*y[i-win];
sumx2 += x[i]*x[i]-x[i-win]*x[i-win];
b = (sumxy - sumx*sumy/win) / (sumx2 - sumx*sumx/win);
a = sumy/win - b*sumx/win;
printf("%g %g\n", x[i-hwin], y[i-hwin] - (a + b*x[i-hwin]));
}
for (i=i-hwin; i<n; i++)
printf("%g %g\n", x[i], y[i] - (a + b*x[i]));
}
}
usage(prog)
char *prog;
{
fprintf(stderr, "Usage : %s hwin\n\n", prog);
fprintf(stderr, " Reads 2 columns from stdin and locally detrends\n");
fprintf(stderr, " by subtracting a least squares fitted line\n");
fprintf(stderr, " over a sliding window 2*hwin+1 points wide.\n");
}