pulsefilter.c

Go to the documentation of this file.

/*
 * This file is part of Titan2Reader
 *
 * Copyright (C) 2006, Sebastien Judenherc <sebastien.judenherc@agecodagis.com>
 *
 * 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
 *
 */
#include <stdio.h>
#include <math.h>
#include <ctype.h>
#include <stdlib.h>
#include <unistd.h>
#include <time.h>
#include <string.h>
#include <assert.h>
#include <sys/types.h>
#include "common.h"

#define MIN_REG_POINTS 10
#define FIRST_PASS_RATIO 10.
#define MIN_WINDOW_LEN (12*3600)

struct xy { double x,y; };

int dtcmp(const void *p1, const void *p2)
{
        struct xy *d1=(struct xy*)p1;
        struct xy *d2=(struct xy*)p2;
        if (d1->x>d2->x)
        {
                return(1);
        }
        else
        {
                return(-1);
        }
}

int LinearRegression(struct xy *data, int n, double *a, double *b)
{
        int i;
        double x0;
        double sx2=0.0,sx=0.0,sy=0.0,sxy=0.0;
        double d;
        double nn=0.;
        for (i=0; i<n; i++)
        {
                if (data[i].x<0.)
                {
                        continue;
                }
                x0=data[i].x-data[0].x;
                sx2+=(x0*x0);
                sx+=x0;
                sy+=data[i].y;
                sxy+=x0*data[i].y;
                nn+=1.0;
        }
        d=sx2*(double)nn-sx*sx;
        if ((d>1.e-6)||(d<-1.e-6))
        {
                *a=(sxy*(double)nn-sx*sy)/d;
                *b=(sx2*sy-sx*sxy)/d;
                return(0);
        }
        else
        {
                return(-1);
        }
}

#define ALLOC_BLOCK 128

void PulseFilter(char *dtfile, double dtmax)
{
        FILE *fin=NULL;
        FILE *fout=NULL;
        char fltfile[255];
        struct xy *data=NULL;
        double Tmin=-1.;
        double Tmax=4294967296.0; // 1<<32
        double dt;
        double Treset=-1.;
        int nalloc=0;
        int ndata=0,i;
        char line[255];
        char sdate[64];
        int istart,istop;
        int N,n,cnt;
        double a=0.,b=0.;
        // open fin for read
        fin=fopen(dtfile,"rt");
        if (!fin)
        {
                PrintError("%s: %m\n",dtfile);
                return;
        }
        // open fout for write
        strcpy(fltfile,dtfile);
        i=strlen(dtfile)-2;
        fltfile[i++]='f';
        fltfile[i++]='l';
        fltfile[i++]='t';
        fltfile[i++]=0;
        fout=fopen(fltfile,"wt");
        if (!fout)
        {
                PrintError("%s: %m\n",fltfile);
                goto error;
        }
        // while not at the end of fin
        while (!feof(fin))
        {
                // if ndata>=nalloc
                if (ndata>=nalloc-1)
                // then
                {
                        // nalloc+=ALLOC_BLOCK
                        nalloc+=ALLOC_BLOCK;
                        // realloc x and y
                        assert(NULL!=(data=(struct xy *)realloc(data,nalloc*sizeof(struct xy))));
                // endif
                }
                line[0]=0;
                // read the new line
                fgets(line,254,fin);
                // scan for a time set line and feed x,y
                if (sscanf(line,"%s %lf",sdate,&dt)==2)
                {
                                if (str2time(sdate,&(data[ndata].x))>=5)
                                {
                                        data[ndata].y=dt;
                                        ndata++;
                                }
                }
                else
                {
                        if (strstr(line,"time reset"))
                        {
                                str2time(sdate,&Treset);
                        }
                }
        // end while
        }

        // we have now the complete dt file in memory

        // sort the file
        qsort(data,ndata,sizeof(struct xy),dtcmp);

        // set Tmin as the first time after the time set
        istart=0;
        for (i=0; i<ndata; i++)
        {
                if (data[i].x<Treset)
                {
                        istart=i;
                }
        }
        Tmin=data[istart].x;
        // set Tmax as the last time
        Tmax=data[ndata-1].x;

        // compute N and dt such as Tmax-Tmin=N*dt, 1h<dt<3h
        N=1;
        while (((Tmax-Tmin)/(double)N)>MIN_WINDOW_LEN)
        {
                N++;
        }
        dt=(Tmax-Tmin)/(double)N;
        // dtmax=1e-6*dt/2 for osiris dtmax=1e-4*dt/2 for other stations

        // for each segment
        for (n=0; n<N; n++)
        {
                // count the number of points in this interval
                cnt=0;
                i=0;
                istop=0;
                while (i<ndata)
                {
                        if (data[i].x<Tmin+dt*(double)n)
                        {
                                i++;
                        }
                        else
                        {
                                break;
                        }
                }
                istart=i;
                while (i<ndata)
                {
                        if (data[i].x>Tmin+dt*(double)(n+1))
                        {
                                break;
                        }
                        else
                        {
                                i++;
                                istop=i;
                                cnt++;
                        }
                }
                // if we have enough data for this segment
                if (cnt>=MIN_REG_POINTS)
                // then
                {
                        // first pass: reject points that are obviously false
                        // compute the average
                        double avg1=0.;
                        double avg2=0.;
                        double wg=0.;
                        for (i=istart; i<istop; i++)
                        {
                                avg1+=data[i].y;
                        }
                        avg1/=(double)(i-istart);
                        // recompute a weighted avg
                        for (i=istart; i<istop; i++)
                        {
                                double a;
                                a=1./(data[i].y-avg1);
                                a*=a;
                                avg2+=data[i].y*a;
                                wg+=a;
                        }
                        avg2/=wg;
                        // kill points that are not within avg+-FIRST_PASS_RATIO*0.5*dtmax*dt
                        for (i=istart; i<istop; i++)
                        {
                                if (fabs(data[i].y-avg2)>FIRST_PASS_RATIO*0.5*dtmax*dt)
                                {
                                        PrintWarning("pass1: ignoring %s %f >%f\n",
                                                        time2str(data[i].x),
                                                        data[i].y-avg2,
                                                        FIRST_PASS_RATIO*0.5*dtmax*dt);
                                        data[i].x*=-1.;
                                }
                        }

                        // compute the linear regression for this segment
                        cnt=LinearRegression(&(data[istart]), cnt,&a,&b);
                        // output the data if abs(dt-a*t+b)<dtmax*winlen/2
                        for (i=istart; i<istop; i++)
                        {
                                double x,y;
                                if (data[i].x<0.0)
                                {
                                        data[i].x*=-1.;
                                }
                                x=data[i].x-data[istart].x;
                                y=a*x+b;
                                if (fabs(y-data[i].y)>0.5*dtmax*dt)
                                {
                                        PrintWarning("pass2: rejecting %s %f >%f\n",
                                                        time2str(data[i].x),
                                                        data[i].y-y,
                                                        0.5*dtmax*dt);
                                        fprintf(fout,"#%s %f dt=%f\n",
                                                        time2str(data[i].x),data[i].y,
                                                        fabs(a*data[i].x+b-data[i].y));
                                }
                                else
                                {
                                        fprintf(fout,"%s %f\n",time2str(data[i].x),data[i].y);
                                }
                        }
                }
                // else
                else
                {
                        // dump all the data
                        for (i=istart; i<istop; i++)
                        {
                                fprintf(fout,"%s %f\n",time2str(data[i].x),data[i].y);
                        }
                }
        // end for
        }
error:
        // free allocated memory
        if (data)
        {
                free(data);
        }
        if (fin)
        {
                fclose(fin);
        }
        if (fout)
        {
                fclose(fout);
        }
}

---