Commit d29999b1 authored by Cedrick Ansorge's avatar Cedrick Ansorge
Browse files

Implementation of time series filters in Filter.c

parent d2c662b2
......@@ -17,7 +17,10 @@
/*
This module contains the following operators:
Filter highpass
Filter lowpass
Filter bandpass
*/
......@@ -26,62 +29,158 @@
#include "cdo_int.h"
#include "pstream.h"
#include "stdlib.h"
#define NALLOC_INC 1000
#define PI2 6.2832
#define HALF 0.5
static
void detrend(int nts, double missval1, double *array1, double *array2)
// FAST FOURIER TRANSFORMATION (bare)
void fft2(double *real, double *imag, const int n, const int isign)
{
int j;
int n;
double sumj, sumjj;
double sumx, sumjx;
double work1, work2;
double missval2 = missval1;
sumx = sumjx = 0;
sumj = sumjj = n = 0;
for ( j = 0; j < nts; j++ )
if ( !DBL_IS_EQUAL(array1[j], missval1) )
int nn, mmax, m, j, istep, i;
double wtemp, wr, wpr, wpi, wi, theta, tempr, tempi, tmp;
if ( n < 2 || n&(n-1) ) printf("n must be power of 2\n");
nn = n << 1;
j = 1;
//BIT Reversion of data
for ( i=1; i<nn; i+=2 )
{
if ( j > i )
{
//swap real part
tmp = real[j/2];
real[j/2] = real[i/2];
real[i/2] = tmp;
//swap imaginary part
tmp = imag[j/2];
imag[j/2] = imag[i/2];
imag[i/2] = tmp;
}
m = n;
while ( m >= 2 && j > m )
{
j -= m;
m >>= 1;
}
j += m;
}
//Danielson-Lanzcos algorithm
mmax = 2;
while ( nn > mmax )
{
istep = mmax << 1;
theta = isign*(PI2/mmax);
wtemp = sin(HALF*theta);
wpr = -2.0*wtemp*wtemp;
wpi = sin(theta);
wr = 1.0;
wi = 0.0;
for( m = 1; m<mmax; m+=2 )
{
for ( i = m; i <= nn; i+=istep)
{
j=i+mmax;
tempr = wr*real[j/2]-wi*imag[j/2];
tempi = wr*imag[j/2]+wi*real[j/2];
real[j/2] = real[i/2]-tempr;
imag[j/2] = imag[i/2]-tempi;
real[i/2] += tempr;
imag[i/2] += tempi;
}
wr = (wtemp=wr)*wpr-wi*wpi+wr;
wi = wi*wpr+wtemp*wpi+wi;
}
mmax = istep;
}
if ( isign == -1 )
for( i =0; i<n; i++)
{
sumx += array1[j];
sumjx += j * array1[j];
sumj += j;
sumjj += j * j;
n++;
real[i]/=n;
imag[i]/=n;
}
}
work1 = DIV(SUB(sumjx, DIV(MUL(sumx, sumj), n) ),
SUB(sumjj, DIV(MUL(sumj, sumj), n)) );
work2 = SUB(DIV(sumx, n), MUL(work1, DIV(sumj, n)));
for ( j = 0; j < nts; j++ )
array2[j] = SUB(array1[j], ADD(work2, MUL(j, work1)));
/*include from Tinfo.c*/
void getTimeInc(int lperiod, int deltam, int deltay, int *incperiod, int *incunit);
void create_fmasc(const int nts, const double fdata, const double fmin, const double fmax, int *fmasc)
{
double dimin, dimax;
int i, imin, imax;
dimin = nts*fmin / fdata;
dimax = nts*fmax / fdata;
imin = dimin<0? 0 : floor(dimin);
imax = ceil(dimax)>nts/2? nts/2 : ceil(dimax);
//fprintf(stderr, "%7.5f %7.5f %5i %5i %7.5f\n", fmin, fmax, imin, imax, fdata);
fmasc[imin] = 1;
for ( i = imin+1; i <= imax; i++ )
fmasc[i] = fmasc[nts-i] = 1;
}
void filter(const int nts, const int *fmasc, double *array1, double *array2)
{
int i;
double sum=0;
fft2(&array1[0], &array2[0], nts, 1);
for ( i = 0; i < nts; i++ )
if ( ! fmasc[i] ) array1[i] = array2[i] = 0;
fft2(&array1[0], &array2[0], nts, -1);
return;
}
void *Filter(void *argument)
{
static char func[] = "Filter";
static enum {BAND, HIGH, LOW};
char *tunits[] = {"second", "minute", "hour", "day", "month", "year"};
int iunits[] = {31536000, 525600, 8760, 365, 12, 1};
int operatorID;
int operfunc;
int gridsize;
int nrecs;
int gridID, varID, levelID, recID;
int tsID;
int i;
int nts;
int nts,nts2;
int nalloc = 0;
int streamID1, streamID2;
int vlistID1, vlistID2, taxisID1, taxisID2;
int nmiss;
int nvars, nlevel;
int *vdate = NULL, *vtime = NULL;
int tunit;
int incperiod0, incunit0, incunit, dpy, calendar;
int vdate0=0, vtime0=0, year0, month0, day0, hour0, minute0, vdateold;
double missval;
double *array1, *array2;
double fdata;
FIELD ***vars = NULL;
double fmin, fmax;
int *fmasc;
cdoInitialize(argument);
cdoOperatorAdd("bandpass", BAND, 0, NULL);
cdoOperatorAdd("highpass", HIGH, 0, NULL);
cdoOperatorAdd("lowpass" , LOW, 0, NULL);
operatorID = cdoOperatorID();
operfunc = cdoOperatorFunc(operatorID);
streamID1 = streamOpenRead(cdoStreamName(0));
if ( streamID1 < 0 ) cdiError(streamID1, "Open failed on %s", cdoStreamName(0));
......@@ -92,62 +191,150 @@ void *Filter(void *argument)
taxisID2 = taxisDuplicate(taxisID1);
vlistDefTaxis(vlistID2, taxisID2);
tunit = taxisInqTunit(taxisID1);
calendar = taxisInqCalendar(taxisID1);
dpy = calendar_dpy(calendar);
streamID2 = streamOpenWrite(cdoStreamName(1), cdoFiletype());
if ( streamID2 < 0 ) cdiError(streamID2, "Open failed on %s", cdoStreamName(1));
streamDefVlist(streamID2, vlistID2);
nvars = vlistNvars(vlistID1);
tsID = 0;
tsID = 0;
while ( (nrecs = streamInqTimestep(streamID1, tsID)) )
{
int deltay, deltam;
int julval0, julval, jdelta;
int lperiod, incperiod;
int year, month, day, hour, minute;
if ( tsID >= nalloc )
{
nalloc += NALLOC_INC;
vdate = (int *) realloc(vdate, nalloc*sizeof(int));
vtime = (int *) realloc(vtime, nalloc*sizeof(int));
vars = (FIELD ***) realloc(vars, nalloc*sizeof(FIELD **));
}
{
nalloc += NALLOC_INC;
vdate = (int *) realloc(vdate, nalloc*sizeof(int));
vtime = (int *) realloc(vtime, nalloc*sizeof(int));
vars = (FIELD ***) realloc(vars, nalloc*sizeof(FIELD **));
}
vdate[tsID] = taxisInqVdate(taxisID1);
vtime[tsID] = taxisInqVtime(taxisID1);
decode_date(vdate[tsID], &year, &month, &day);
decode_time(vtime[tsID], &hour, &minute);
if ( tsID )
decode_date(vdate[tsID-1], &year0, &month0, &day0);
if ( month == 2 && day == 29 &&
( day0 != day || month0 != month || year0 != year ) )
{
cdoWarning("filtering of multi-year times series only");
cdoWarning(" works properly with 365-day-calendar.");
cdoWarning(" Please delete the day %i-02-29 (cdo del29feb)",year);
}
vars[tsID] = (FIELD **) malloc(nvars*sizeof(FIELD *));
for ( varID = 0; varID < nvars; varID++ )
{
gridID = vlistInqVarGrid(vlistID1, varID);
missval = vlistInqVarMissval(vlistID1, varID);
nlevel = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
vars[tsID][varID] = (FIELD *) malloc(nlevel*sizeof(FIELD));
for ( levelID = 0; levelID < nlevel; levelID++ )
{
vars[tsID][varID][levelID].grid = gridID;
vars[tsID][varID][levelID].missval = missval;
vars[tsID][varID][levelID].ptr = NULL;
}
}
{
gridID = vlistInqVarGrid(vlistID1, varID);
missval = vlistInqVarMissval(vlistID1, varID);
nlevel = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
vars[tsID][varID] = (FIELD *) malloc(nlevel*sizeof(FIELD));
for ( levelID = 0; levelID < nlevel; levelID++ )
{
vars[tsID][varID][levelID].grid = gridID;
vars[tsID][varID][levelID].missval = missval;
vars[tsID][varID][levelID].ptr = NULL;
}
}
for ( recID = 0; recID < nrecs; recID++ )
{
streamInqRecord(streamID1, &varID, &levelID);
gridID = vlistInqVarGrid(vlistID1, varID);
gridsize = gridInqSize(gridID);
vars[tsID][varID][levelID].ptr = (double *) malloc(gridsize*sizeof(double));
streamReadRecord(streamID1, vars[tsID][varID][levelID].ptr, &nmiss);
vars[tsID][varID][levelID].nmiss = nmiss;
}
{
streamInqRecord(streamID1, &varID, &levelID);
gridID = vlistInqVarGrid(vlistID1, varID);
gridsize = gridInqSize(gridID);
vars[tsID][varID][levelID].ptr = (double *) malloc(gridsize*sizeof(double));
streamReadRecord(streamID1, vars[tsID][varID][levelID].ptr, &nmiss);
vars[tsID][varID][levelID].nmiss = nmiss;
if( nmiss ) cdoAbort("Missing value support for operators in module Filter not added yet");
}
/* get and check time increment */
if ( tsID > 0)
{
julval0 = encode_julval(dpy, vdate[tsID-1], vtime[tsID-1]);
julval = encode_julval(dpy, vdate[tsID], vtime[tsID]);
jdelta = julval - julval0;
lperiod = (long)(jdelta+0.5);
incperiod = (int) lperiod;
deltay = year-year0;
deltam = deltay*12 + (month-month0);
if ( tsID == 1 )
{
getTimeInc(lperiod, deltam, deltay, &incperiod0, &incunit0);
incperiod = incperiod0;
if ( incperiod == 0 ) cdoAbort("Time step must be different from zero\n");
incunit = incunit0;
//fprintf(stderr, "step %i%s\n", incperiod, tunits[incunit]);
fdata = 1.*iunits[incunit]/incperiod;
}
else
getTimeInc(lperiod, deltam, deltay, &incperiod, &incunit);
if ( ! ( incperiod == incperiod0 && incunit == incunit0 ) )
cdoWarning("Time increment in step %i differs from step 0", tsID);
}
tsID++;
}
nts = tsID;
array1 = (double *) malloc(nts*sizeof(double));
array2 = (double *) malloc(nts*sizeof(double));
/* round up nts to next power of two for (better) performance
** of fast fourier transformation */
nts2 = nts-1;
nts2 |= nts2 >> 1; // handle 2 bit numbers
nts2 |= nts2 >> 2; // handle 4 bit numbers
nts2 |= nts2 >> 4; // handle 8 bit numbers
nts2 |= nts2 >> 8; // handle 16 bit numbers
nts2 |= nts2 >> 16; // handle 32 bit numbers
nts2++;
array1 = (double *) malloc(nts2*sizeof(double));
array2 = (double *) malloc(nts2*sizeof(double));
fmasc = (int *) calloc(nts2,sizeof(int));
switch(operfunc)
{
case BAND:
{
operatorInputArg("Lower bound of frequency band:");
fmin = atof(operatorArgv()[0]);
operatorInputArg("Upper bound of frequency band:");
fmax = atof(operatorArgv()[1]);
break;
}
case HIGH:
{
operatorInputArg("Lower bound of frequency pass:");
fmin = atof(operatorArgv()[0]);
fmax = fdata;
break;
}
case LOW:
{
fmin = 0;
operatorInputArg("Upper bound of frequency pass:");
fmax = atof(operatorArgv()[0]);
break;
}
}
create_fmasc(nts2, fdata, fmin, fmax, &fmasc[0]);
for ( varID = 0; varID < nvars; varID++ )
{
......@@ -156,20 +343,35 @@ void *Filter(void *argument)
gridsize = gridInqSize(gridID);
nlevel = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
for ( levelID = 0; levelID < nlevel; levelID++ )
{
for ( i = 0; i < gridsize; i++ )
{
for ( tsID = 0; tsID < nts; tsID++ )
array1[tsID] = vars[tsID][varID][levelID].ptr[i];
detrend(nts, missval, array1, array2);
for ( tsID = 0; tsID < nts; tsID++ )
vars[tsID][varID][levelID].ptr[i] = array2[tsID];
}
}
{
for ( i = 0; i < gridsize; i+=2 )
{
for ( tsID = 0; tsID < nts; tsID++ )
{
array1[tsID] = vars[tsID][varID][levelID].ptr[i];
if ( i < gridsize - 1 )
array2[tsID] = vars[tsID][varID][levelID].ptr[i+1];
}
// zero padding up to next power of to
for ( tsID = nts; tsID < nts2; tsID++ )
{
array1[tsID]=0;
array2[tsID]=0;
}
filter(nts2, fmasc, &array1[0], &array2[0]);
for ( tsID = 0; tsID < nts; tsID++ )
{
vars[tsID][varID][levelID].ptr[i] = array1[tsID];
if ( i < gridsize - 1 )
vars[tsID][varID][levelID].ptr[i+1] = array2[tsID];
}
}
}
}
if ( array1 ) free(array1);
if ( array2 ) free(array2);
......@@ -178,33 +380,33 @@ void *Filter(void *argument)
taxisDefVdate(taxisID2, vdate[tsID]);
taxisDefVtime(taxisID2, vtime[tsID]);
streamDefTimestep(streamID2, tsID);
for ( varID = 0; varID < nvars; varID++ )
{
nlevel = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
for ( levelID = 0; levelID < nlevel; levelID++ )
{
if ( vars[tsID][varID][levelID].ptr )
{
nmiss = vars[tsID][varID][levelID].nmiss;
streamDefRecord(streamID2, varID, levelID);
streamWriteRecord(streamID2, vars[tsID][varID][levelID].ptr, nmiss);
free(vars[tsID][varID][levelID].ptr);
}
}
free(vars[tsID][varID]);
}
{
nlevel = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
for ( levelID = 0; levelID < nlevel; levelID++ )
{
if ( vars[tsID][varID][levelID].ptr )
{
nmiss = vars[tsID][varID][levelID].nmiss;
streamDefRecord(streamID2, varID, levelID);
streamWriteRecord(streamID2, vars[tsID][varID][levelID].ptr, nmiss);
free(vars[tsID][varID][levelID].ptr);
}
}
free(vars[tsID][varID]);
}
free(vars[tsID]);
}
if ( vars ) free(vars);
if ( vdate ) free(vdate);
if ( vtime ) free(vtime);
streamClose(streamID2);
streamClose(streamID1);
cdoFinish();
return (0);
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment