Skip to content

Commit

Permalink
Merge branch 'master' of https://github.com/teuben/nemo
Browse files Browse the repository at this point in the history
  • Loading branch information
@teuben
teuben committed Oct 4, 2023
2 parents dda6569 + f0a0e48 commit 0c6ce44
Show file tree
Hide file tree
Showing 2 changed files with 82 additions and 76 deletions.
50 changes: 38 additions & 12 deletions man/man1/tabsmooth.1
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,10 +7,13 @@ tabsmooth \- smooth columns of a table (Hanning, Savitzky-Golay)
\fBtabsmooth\fP [parameter=value]

.SH "DESCRIPTION"
Quick and dirty program to Hanning (1/4,1/2,1/4) smooth a table.
Quick and dirty program to \fIsmooth\fP the column(s) of a table
without paying attention to the X coordinate (meaning: a constant
step size is assumed).
.PP
Equally quick and dirty is the filter= option, to try some
Savitzky-Golay filters.
The filter= option is used to select the filter, for example
a Hanning and a Savitzky-Golay filter (and some derivates)
can be selected. A pass-through is also allowed.

.SH "PARAMETERS"
The following parameters are recognized in any order if the keyword
Expand All @@ -27,35 +30,58 @@ columns are output to stdout. [1]
Filter option. The way the options are set will likely change, this
is work in progress.
.nf
-1 passthrough
0 hanning (1,2,1)
1 SK 5pt smooth
2 SK 1st derivative
3 SK 2nd derivative
4 SK 7pt smooth
-1 passthrough, the column is not modified
0 hanning (1,2,1)/4
1 SK2 5pt smooth
2 SK2 1st derivative with h=1
3 SK2 2nd derivative with h=1
4 SK4 7pt smooth
5 SK4 1st derivative
.fi
The smoothing versions are normalized to conserve "flux".
Default: 0
.TP
\fBsmooth=\fP
Manually set the smoothing array. By default, the filter= is used to set the smoothing
array. Be sure to enter an odd number
of array elements, for example a Hanning would be
of array elements, for example a normalized Hanning would be
smooth=0.25,0.5,0.25.
Default: not provided.
.TP
\fBnmax=\fP
max size if a pipe [100000] - to be deprecated for new table system

.SH "EXAMPLES"
Here we draw 100000 random values from a normal distribution (mean 0, dispersion 1) and smooth them
with various filters:
.nf

tabgen - 100000 1 2 123| tabsmooth - 1 -1 | tabstat - | grep disp
1.0022
tabgen - 100000 1 2 123| tabsmooth - 1 0 | tabstat - | grep disp
0.6128
tabgen - 100000 1 2 123| tabsmooth - 1 1 | tabstat - | grep disp
0.6973
tabgen - 100000 1 2 123| tabsmooth - 1 4 | tabstat - | grep disp
0.7537

.fi

.SH "CAVEATS"
Near the edge the smoothing array is cut off, thus there is some loss of flux.
Otherwise in general the smoothing filter is normalized and should conserve flux.

.SH "SEE ALSO"
tabmath(1NEMO), tabtrend(1NEMO)
tabmath(1NEMO), tabtrend(1NEMO), tabrows(1NEMO), tabcols(1NEMO), table(5NEMO)

.SH "AUTHOR"
Peter Teuben

.SH "HISTORY"
.nf
.ta +1.5i +5.5i
20-Dec-2010 V0.1 Created PJT
20-dec-2010 V0.1 Created PJT
13-oct-2014 documented smooth=
28-sep-2023 V0.5 added filter= PJT
29-sep-2023 V0.6 converted to table V2 PJT
.fi
108 changes: 44 additions & 64 deletions src/kernel/tab/tabsmooth.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,8 @@
* 24-oct-07 Created quick & dirty PJT
* 12 keyword error
* may-13 smooth= added
* 28-sep-23 filter= added
* 28-sep-23 0.5 filter= added
* 29-sep-23 0.6 converted to table V2 format
*
*/

Expand All @@ -16,23 +17,24 @@
#include <moment.h>
#include <yapp.h>
#include <axis.h>
#include <table.h>
#include <mdarray.h>

/**************** COMMAND LINE PARAMETERS **********************/

string defv[] = {
"in=???\n Input file name",
"xcol=1\n Column(s) to use",
"filter=0\n Select one of the test filters (0,1,2,...)",
"smooth=\n Optional expliciti smoothing array",
"nmax=100000\n max size if a pipe",
"VERSION=0.5\n 28-sep-2023 PJT",
"xcol=1\n Column(s) to use (1=first)",
"filter=0\n Select one of the test filters (-1,0,1,2,...)",
"smooth=\n Optional explicit smoothing array",
"VERSION=0.6\n 29-sep-2023 PJT",
NULL
};

string usage = "(hanning) smooth columns of a table";
string usage = "smooth columns of a table (Hanning, Savitzky-Golay)";

/**************** GLOBAL VARIABLES *****************************/

/**************** LOCAL VARIABLES *****************************/

#ifndef MAXHIST
#define MAXHIST 1024
Expand All @@ -42,21 +44,19 @@ string usage = "(hanning) smooth columns of a table";
#define MAXCOL 256
#endif

#define MAXCOORD 16
#define MAXSM 9

local string input; /* filename */
local stream instr; /* input file */
local table *tptr; /* table pointer */

local int ncol; /* number of columns used */
local int col[MAXCOL]; /* histogram column number(s) */

real *coldat[MAXCOL];
local int nmax; /* lines to use at all */
local int npt; /* actual number of points */

local int nsm;
local int xcol[MAXCOL]; /* histogram column number(s) */
local int nxcol; /* actual number of columns used */
local mdarray2 x; /* x[col][row] */
local int npt; /* actual number of points (rows) in table */

local real sm[MAXSM]; /* smoothing array */
local int nsm; /* actual length of smoothing array */


/* Savitzky-Golay tables for fixed stepsize */
Expand All @@ -68,21 +68,21 @@ typedef struct cst {
} cst, *cstptr;


cst cst0 = {3, 4, {1, 2, 1}}; // Hanning
cst csts[] = {
{3, 4, {1, 2, 1}}, // Hanning
{5, 35, {-3, 12, 17, 12, -3}}, // SK4 - smooth
{5, 12, { 1, -8, 0, 8, -1}}, // SK4 - 1st der
{5, 7, { 2, -1, -2, -1, 2}}, // SK4 - 2nd der
{7, 21, {-2, 3, 6, 7, 6, 3, -2}},
NULL,
local cst csts[] = {
{3, 4, { 1, 2, 1}}, // 0: Hanning (3pt)
{5, 35, {-3, 12, 17, 12, -3}}, // 1: SK2 - smooth (5pt)
{5, 12, { 1, -8, 0, 8, -1}}, // 2: SK2 - 1st der
{5, 7, { 2, -1, -2, -1, 2}}, // 3: SK2 - 2nd der
{7,231, { 5,-30, 75,131, 75,-30, 5}}, // 4: SK4 - smooth (7pt)
{7,252, {22,-67,-58, 0, 58, 67,-22}}, // 5: SK4
};

/* forward references */

local void setparams(void);
local void read_data(void);
local void build_filter(int);
local void smooth_data(void);
local void smooth_data_old(void);



Expand All @@ -100,16 +100,18 @@ local void setparams()
real sum;
int i;
int filter = getiparam("filter");

int nrows, ncols;

input = getparam("in");
ncol = nemoinpi(getparam("xcol"),col,MAXCOL);
if (ncol < 0) error("parsing error xcol=%s",getparam("xcol"));
nxcol = nemoinpi(getparam("xcol"),xcol,MAXCOL);
if (nxcol < 0) error("parsing error xcol=%s",getparam("xcol"));

nmax = nemo_file_lines(input,getiparam("nmax"));
if (nmax<1) error("Problem reading from %s",input);

instr = stropen (input,"r");

tptr = table_open(instr, 0);
nrows = table_nrows(tptr);
ncols = table_ncols(tptr);
dprintf(1,"Table: %d x %d\n", nrows, ncols);

nsm = nemoinpr(getparam("smooth"),sm,MAXSM);
if (nsm < 0) error("smooth=%s parsing error",getparam("smooth"));
if (nsm == 0)
Expand All @@ -118,24 +120,16 @@ local void setparams()
if (nsm % 2 != 1) error("smooth=%s needs an odd number of values",getparam("smooth"));

for (i=0, sum=0.0; i<nsm; i++) sum += sm[i];
dprintf(0,"Smooth sum= %g\n",sum);

dprintf(1,"Smooth sum= %g\n",sum);
}



local void read_data()
{
int i,j,k;

dprintf(0,"Reading %d column(s)\n",ncol);
for (i=0; i<ncol; i++)
coldat[i] = (real *) allocate(sizeof(real)*nmax);
npt = get_atable(instr,ncol,col,coldat,nmax); /* read it */
if (npt == -nmax) {
warning("Could only read %d data",nmax);
npt = nmax;
}
x = table_md2cr(tptr, nxcol, xcol, 0, 0);
npt = tptr->nr;
if (nxcol == 0) nxcol = tptr->nc;
}

local void build_filter(int filter)
Expand All @@ -147,30 +141,16 @@ local void build_filter(int filter)
sm[0] = 1;
return;
}
int maxfilter = sizeof(csts)/sizeof(cst);
if (filter >= maxfilter) error("Only %d filters available",maxfilter);

cst *c = &csts[filter];
nsm = c->n;
for (i=0; i<nsm; i++)
sm[i] = c->coeff[i]/c->norm;
}


local void smooth_data_old(void)
{
int i,j;

for (j=0; j<ncol; j++)
printf("%g ",0.667*coldat[j][1] + 0.333*coldat[j][0]);
printf("\n");

for (i=1; i<npt-1; i++) {
for (j=0; j<ncol; j++)
printf("%g ",0.25*coldat[j][i-1] + 0.5*coldat[j][i] + 0.25*coldat[j][i+1]);
printf("\n");
}
for (j=0; j<ncol; j++)
printf("%g ",0.667*coldat[j][npt-2] + 0.333*coldat[j][npt-1]);
printf("\n");
}

local void smooth_data(void)
{
Expand All @@ -179,12 +159,12 @@ local void smooth_data(void)
int nsmh = (nsm-1)/2; /* nsm better be odd */

for (i=0; i<npt; i++) {
for (j=0; j<ncol; j++) {
for (j=0; j<nxcol; j++) {
sum[j] = 0.0;
for (k=-nsmh; k<=nsmh; k++) {
ipk = i+k;
if (ipk<0 || ipk>=npt) continue;
sum[j] += coldat[j][ipk] * sm[k+nsmh];
sum[j] += x[j][ipk] * sm[k+nsmh];
}
printf("%g ",sum[j]);
}
Expand Down

0 comments on commit 0c6ce44

Please sign in to comment.