Category : C Source Code
Archive   : SMDATA.ZIP
Filename : SMOOTH.C
| smooth - smooth X,Y data
|----------------------------------------------------------------
| Reads X,Y pairs and performs a rolling average of any one of
| several types, uniform distribution assumed.
|----------------------------------------------------------------
| Author Bill Davidsen ([email protected]) June 1989
| v1.9 last change 8/15/89
****************************************************************/
#include
#define R 1
#if 111
#define version "v1.9"
#else
#define version "beta test version"
#endif
#ifndef MAXINT
#define MAXINT (((unsigned)~0) >> 1)
#endif
/* floating absolute difference */
float T_float;
#define FAD(x,y) ((T_float=(x)-(y))<0 ? -T_float : T_float)
/* getopt stuff */
extern int optind, opterr;
extern char *optarg;
/* global flags */
int verbose = 0, debug = 0;
int minflag = 0, maxflag = 0, weightflag = 0, fixrange = 0;
int endX = 0;
int ncycles = 2;
float minval, maxval;
/* pointers for reading data */
int pointsread = 0, arraysize = 0;
float *X, *Y, *Xstart, *Xend;
/* pointers for weighted average */
int *weightpts, weightlen;
float weightrange;
float *weightsum, weightlow = 1e6, weighthigh = -1e6;
float setweight();
FILE *datafile = stdin;
/* algorithm for selecting points to average */
enum { OFFSET, RANGE } R_select = OFFSET;
/* smoothing routines */
void smoothu(), smoothw(), smootha(), smoothi(), smooth2(), smoothp();
void (*smoothptr)() = smoothu;
/* data input routines */
void readXY(), readXXY();
void (*readptr)() = readXY;
/* useful routines */
char *malloc(), *realloc();
FILE *fopen();
double atof();
main(argc,argv)
int argc;
char *argv[];
{
int n;
float sfact = 3;
char *Options = "r:s:A:a:vx:m:M:i:I:2wWhE?";
while ((n = getopt(argc, argv, Options)) != EOF) {
switch (n) {
case 's': /* smoothing factor */
sfact = atof(optarg);
break;
case 'r': /* smooth by range */
sfact = atof(optarg);
R_select = RANGE;
break;
case 'A': /* average in N width strips */
fixrange = 1;
weightrange = atof(optarg);
readptr = readXXY;
smoothptr = smootha;
break;
case 'a': /* average over an area */
weightlen = atoi(optarg);
weightsum = (float *)malloc(sizeof(float) * weightlen);
weightpts = (int *)malloc(sizeof(int) * weightlen);
readptr = readXXY;
smoothptr = smootha;
break;
case 'E':
endX = 1;
break;
case 'v': /* increase verbose */
++verbose;
break;
case 'x': /* debug */
debug = atoi(optarg);
break;
case 'm': /* set minimum */
++minflag;
minval = atof(optarg);
break;
case 'M': /* set maximum */
++maxflag;
maxval = atof(optarg);
break;
case 'i': /* interval data */
sfact = atof(optarg);
smoothptr = smoothi;
break;
case 'I': /* number of intervals */
ncycles = atoi(optarg);
break;
case '2': /* average every two points */
smoothptr = smooth2;
break;
case 'w': /* linear weight the average */
weightflag = 1;
smoothptr = smoothw;
break;
case 'W': /* weight on square of X distance */
weightflag = 2;
smoothptr = smoothw;
break;
case 'h': /* hi-lo average */
smoothptr = smoothp;
break;
case '?': /* invalid option */
explain();
exit(1);
}
}
/* if there is a filename, use it */
if (optind < argc) {
datafile = fopen(argv[optind], "r");
if (datafile == NULL) {
fprintf(stderr, "Can't open data file %s\n", argv[optind]);
exit(1);
}
}
/* read the data into the global X and Y arrays */
(*readptr)();
if (verbose) {
fprintf(stderr, "\nRead %d points, smooth %f\n",
pointsread, sfact);
}
/* smooth uniform data */
(*smoothptr)(sfact);
exit(0);
}
/*****************************************************************
| explain - give usage hints
|----------------------------------------------------------------
| This does not replace the man page, just lists the options
****************************************************************/
explain() {
printf("smooth / %s\n\nUsage:\n smooth [options] file\n\n",version);
printf( "Where options are:\n");
printf( " /s N smooth over +/-N points\n");
printf( " /r N smooth over points within +/- N\n");
printf( " /a N accept wide data of the form Xstart, Xstop, Y\n");
printf( " and distribute over N equal points\n");
printf( " /A W accept wide data and break into W wide strips\n");
printf( " /E duration date is Xstart, Xend, Y format\n");
printf( " /h average the high and low values within +/-N\n");
printf( " /2 average every 2 points\n");
printf( " /v verbose\n");
printf( " /w linear weighting of data\n");
printf( " /W pseudo-normal weighting of data\n");
printf( " /x N set debug level to N\n");
printf( " /i N sample at X intervals of N\n");
printf( " /I N sample N intervals (def 2)\n");
printf( " /m N set minimum Y value to N, ignore lower\n");
printf( " (used with /r and /i only)\n");
printf( " /M N set maximum Y value to N, ignore higher\n");
printf( " (used with /r and /i only)\n");
printf( "\n");
}
/*****************************************************************
| inrange - test if a Y value is inrage for averaging
****************************************************************/
int
inrange(val)
float val;
{
return (!minflag || val >= minval) && (!maxflag || val <= maxval);
}
/*****************************************************************
| setweight - return weight to apply to a point at this X offset
****************************************************************/
float
setweight(offset, max)
float offset, max;
{
float w;
switch (weightflag) {
case 0: /* linear */
w = 1;
break;
case 1: /* simple linear */
w = (max - offset)/max;
break;
case 2: /* square linear */
w = (max - offset)/max;
w *= w;
break;
}
if (debug > 2) {
fprintf(stderr,
"Type %d weight %f for %f and %f\n",
weightflag, w, offset, max
);
}
return w;
}
/*****************************************************************
| setlimits - set the start and end subscripts
|----------------------------------------------------------------
| Pick points based on either range or points offset from the
| current subscript or X value.
*****************************************************************/
static void
setlimits(crntss, range, start, end)
int crntss, *start, *end;
float range;
{
int n, width = range, maxpt = pointsread-1;
float limit;
switch (R_select) {
case OFFSET: /* all within N points */
/* set low, not less than zero */
if ((n = crntss-width) < 0) n = 0;
*start = n;
/* set high, not more than data available */
if ((n = crntss+width) > maxpt) n = maxpt;
*end = n;
break;
case RANGE: /* all within a given X range */
limit = X[crntss]-range;
/* scan down until limit found */
for (n = crntss; n >= 0 && X[n] >= limit; --n);
*start = n + 1;
/* scan up until limit found */
limit = X[crntss]+range;
for (n=crntss; n <= maxpt && X[n] <= limit; ++n);
*end = n - 1;
break;
}
}
/*****************************************************************
| readXY - read the X and Y data from the data file
|----------------------------------------------------------------
| The arrays are expanded as needed
****************************************************************/
void
readXY() {
float tx, ty;
int n;
while (fscanf(datafile, "%f %f", &tx, &ty) == 2) {
/* got more data, be sure there's room */
if (pointsread == arraysize) {
/* allocate or grow the buffers */
if (arraysize == 0) {
/* allocate the initial buffers */
arraysize = 1024;
X = (float *)malloc(1024 * sizeof(float));
Y = (float *)malloc(1024 * sizeof(float));
if (verbose) fprintf(stderr, "Reading");
}
else {
/* grow the buffers */
arraysize += 1024;
X = (float *)realloc(X, arraysize * sizeof(float));
Y = (float *)realloc(Y, arraysize * sizeof(float));
}
/* check the status */
if (X == NULL || Y == NULL) {
fprintf(stderr, "Can't alloc buffer!\n");
exit(1);
}
if (verbose) putc('.', stderr);
}
/* add the data to the buffer */
X[pointsread] = tx;
Y[pointsread] = ty;
++pointsread;
}
}
/*****************************************************************
| readXXY - read Xstart, Xend, and Y data from the data file
|----------------------------------------------------------------
| The arrays are expanded as needed
****************************************************************/
void
readXXY() {
float txstart, txend, tx, ty;
int n;
while (fscanf(datafile, "%f %f %f", &txstart, &txend, &ty) == 3) {
/* got more data, be sure there's room */
if (pointsread == arraysize) {
/* allocate or grow the buffers */
if (arraysize == 0) {
/* allocate the initial buffers */
arraysize = 1024;
Xstart = (float *)malloc(1024 * sizeof(float));
Xend = (float *)malloc(1024 * sizeof(float));
Y = (float *)malloc(1024 * sizeof(float));
if (verbose) fprintf(stderr, "Reading");
}
else {
/* grow the buffers */
arraysize += 1024;
Xstart = (float *)realloc(Xstart, arraysize * sizeof(float));
Xend = (float *)realloc(Xend, arraysize * sizeof(float));
Y = (float *)realloc(Y, arraysize * sizeof(float));
}
/* check the status */
if (Xstart == NULL || Xend == NULL || Y == NULL) {
fprintf(stderr, "Can't alloc buffer!\n");
exit(1);
}
if (verbose) putc('.', stderr);
}
/* see if the 2nd dat point is an ending value */
if (endX) txend -= txstart;
/* add the data to the buffer */
Xstart[pointsread] = txstart;
Xend[pointsread] = txend;
Y[pointsread] = ty;
/* check the limits */
if (txstart < weightlow) weightlow = txstart;
if ((txstart += txend) > weighthigh) weighthigh = txstart;
/* final incr of the points read */
++pointsread;
}
}
/*****************************************************************
| smoothu - smooth uniform data
|----------------------------------------------------------------
| Average N points and output average, just the points at ends
****************************************************************/
void
smoothu(fwidth)
float fwidth;
{
int n, m, npoints, width = fwidth;
int hi, lo, oldhi, oldlo;
float sum = 0, avg;
/* output the fully averaged points */
for (n = 0; n <= pointsread-1; ++n) {
/* save the current limits and get the new ones */
oldhi = hi;
oldlo = lo;
setlimits(n, fwidth, &lo, &hi);
if (debug) {
fprintf(stderr, "lo %4d, hi %4d\n", lo, hi);
}
/*
* This is a timesaver... drop the points no longer in the average
* and add the new ones, rather than recalc the average.
*/
if (n > 0) {
/* don't do this the first time */
for (m=oldlo; m < lo; ++m) sum -= Y[m];
for (m=oldhi+1; m <= hi; ++m) sum += Y[m];
}
else {
/* do this the first time instead, build the initial sum */
for (sum = 0, m = lo; m <= hi; ++m) sum += Y[m];
}
npoints = hi-lo+1;
/* calculate the average value */
avg = sum/npoints;
/* output the point */
if (debug) {
fprintf(stderr, "pt %-4d np %-2d sum %-5.0f\n",
n, npoints, sum);
}
printf("%f %f\n", X[n], avg);
}
}
/*****************************************************************
| smoothw - smooth data or a given range
|----------------------------------------------------------------
| Average all points within +/- width of the current point
****************************************************************/
void
smoothw(width)
float width;
{
int n, m;
int hi, lo;
float sum = 0, avg, limit, pweight, npoints;
/* output the fully averaged points */
for (n = 0; n <= pointsread-1; ++n) {
sum = Y[n];
npoints = 1;
/* set the limits and create the weighted sum */
setlimits(n, width, &lo, &hi);
sum = npoints = 0;
for (m = lo; m <= hi; ++m) {
/* get the weight factor and apply it */
pweight = setweight(FAD(X[n],X[m]), width);
sum += Y[m] * pweight;
npoints += pweight;
}
avg = sum / npoints;
/* output the point, then slide the average up */
if (debug) {
fprintf(stderr, "pt %-4d np %-6.2f sum %-5.0f\n",
n, npoints, sum);
}
printf("%f %f\n", X[n], avg);
}
}
/*****************************************************************
| smootha - average data by time into points
|----------------------------------------------------------------
| This is used on data with a duration, in the form Xstart,
| Xens, Yval. The range of all X values is broken into a number of
| slices, and the value for any slice is the average of all Y
| values who's X value range includes this slice.
****************************************************************/
void
smootha() {
int n, m;
int first, last;
float tx, ty;
float weightincr;
float wt_tmp1;
/* check for fixed range specified */
if (fixrange) {
/* determine the number of strips */
weightlen = 1 + (weighthigh - weightlow) / weightrange;
weightsum = (float *)malloc(sizeof(float) * weightlen);
weightpts = (int *)malloc(sizeof(int) * weightlen);
weightincr = weightrange;
}
else weightincr = (weighthigh - weightlow) / (weightlen - 1);
/* init the arrays */
for (n=0; n
}
/* process a point at a time loop */
for (n=0; n < pointsread; ++n) {
/* map into starting and ending indices */
wt_tmp1 = (Xstart[n] - weightlow);
first = wt_tmp1 / weightincr;
last = (wt_tmp1 + Xend[n]) / weightincr;
/* add the values to the vectors */
for (m = first; m <= last; ++m) {
weightsum[m] += Y[n];
weightpts[m]++;
}
}
/* form the averages here */
for (n = 0; n < weightlen; ++n) {
if (m = weightpts[n])
weightsum[n] /= m;
}
/* output the first point */
printf("%f %f\n", weightlow, weightsum[0]);
/* loop to the next to last point */
for (n = 1; n < (weightlen-1); ++n) {
printf("%f %f\n", weightlow+n*weightincr, weightsum[n]);
}
/* output the last point */
if (!fixrange)
printf("%f %f\n", weighthigh, weightsum[weightlen-1]);
}
/*****************************************************************
| smoothi - smooth by averaging points at intervals
|----------------------------------------------------------------
| This applies to data at uniform (and assumed cyclic) intervals.
| The data at the +/- N interval of X is averaged with the current
| value.
****************************************************************/
void
smoothi(interval)
float interval;
{
int n, m, npoints, tp;
float sum, testval, limit, tick;
/* set the limit of closeness of match */
limit = X[1] - X[0];
tick = interval/limit;
/* loop thru the values */
for (n=0; n < pointsread; ++n) {
sum = npoints = 0;
if (debug) {
fprintf(stderr, "n=%-4d base X %15f\n", n, X[n]);
}
/* scan +/-2 intervals (for now) */
for (m = -ncycles; m <= ncycles; ++m) {
tp = n + (m * tick);
if (0 <= tp && tp < pointsread && inrange(Y[tp])) {
sum += Y[tp];
++npoints;
if (debug) {
fprintf(stderr, " used X=%15f\n", X[tp]);
}
}
}
/* output what we have */
printf("%f %f\n", X[n], sum/npoints);
}
}
/*****************************************************************
| smooth2 - just average every two points
|----------------------------------------------------------------
| This reduces the number of points by one
****************************************************************/
void
smooth2(sfact)
float sfact;
{
int n;
float sum, Xavg, Yavg;
/* just a simple loop thru the points */
for (n=1; n < pointsread; ++n) {
printf("%f %f\n",
(X[n-1]+X[n])/2,
(Y[n-1]+Y[n])/2
);
}
}
/*****************************************************************
| smoothp - smooth average of peaks within +/- N points
|----------------------------------------------------------------
| Finds the highest and lowest point within the specified range
| and takes the average of just those points.
****************************************************************/
void
smoothp(range)
float range;
{
int width = range;
int start, end;
int n, m;
float high, low, crntY;
/* loop through all points */
for (n=0; n < pointsread; ++n) {
/* set the starting and ending subscripts */
setlimits(n, range, &start, &end);
/* now scan for the peak values */
low = MAXINT;
high = -low;
for (m=start; m <= end; ++m) {
crntY = Y[m];
if (crntY < low) low = crntY;
if (crntY > high) high = crntY;
}
/* output the X and average */
printf("%f %f\n", X[m], (high+low)/2);
}
}
Very nice! Thank you for this wonderful archive. I wonder why I found it only now. Long live the BBS file archives!
This is so awesome! 😀 I’d be cool if you could download an entire archive of this at once, though.
But one thing that puzzles me is the “mtswslnkmcjklsdlsbdmMICROSOFT” string. There is an article about it here. It is definitely worth a read: http://www.os2museum.com/wp/mtswslnk/