00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include <stdio.h>
00023 #include <math.h>
00024 #include <ctype.h>
00025 #include <stdlib.h>
00026 #include <unistd.h>
00027 #include <time.h>
00028 #include <string.h>
00029 #include <assert.h>
00030 #include <sys/types.h>
00031 #include "common.h"
00032
00033 #define MIN_REG_POINTS 10
00034 #define FIRST_PASS_RATIO 10.
00035 #define MIN_WINDOW_LEN (12*3600)
00036
00037 struct xy { double x,y; };
00038
00039 int dtcmp(const void *p1, const void *p2)
00040 {
00041 struct xy *d1=(struct xy*)p1;
00042 struct xy *d2=(struct xy*)p2;
00043 if (d1->x>d2->x)
00044 {
00045 return(1);
00046 }
00047 else
00048 {
00049 return(-1);
00050 }
00051 }
00052
00060 int LinearRegression(struct xy *data, int n, double *a, double *b)
00061 {
00062 int i;
00063 double x0;
00064 double sx2=0.0,sx=0.0,sy=0.0,sxy=0.0;
00065 double d;
00066 double nn=0.;
00067 for (i=0; i<n; i++)
00068 {
00069 if (data[i].x<0.)
00070 {
00071 continue;
00072 }
00073 x0=data[i].x-data[0].x;
00074 sx2+=(x0*x0);
00075 sx+=x0;
00076 sy+=data[i].y;
00077 sxy+=x0*data[i].y;
00078 nn+=1.0;
00079 }
00080 d=sx2*(double)nn-sx*sx;
00081 if ((d>1.e-6)||(d<-1.e-6))
00082 {
00083 *a=(sxy*(double)nn-sx*sy)/d;
00084 *b=(sx2*sy-sx*sxy)/d;
00085 return(0);
00086 }
00087 else
00088 {
00089 return(-1);
00090 }
00091 }
00092
00093 #define ALLOC_BLOCK 128
00094
00099 void PulseFilter(char *dtfile, double dtmax)
00100 {
00101 FILE *fin=NULL;
00102 FILE *fout=NULL;
00103 char fltfile[255];
00104 struct xy *data=NULL;
00105 double Tmin=-1.;
00106 double Tmax=4294967296.0;
00107 double dt;
00108 double Treset=-1.;
00109 int nalloc=0;
00110 int ndata=0,i;
00111 char line[255];
00112 char sdate[64];
00113 int istart,istop;
00114 int N,n,cnt;
00115 double a=0.,b=0.;
00116
00117 fin=fopen(dtfile,"rt");
00118 if (!fin)
00119 {
00120 PrintError("%s: %m\n",dtfile);
00121 return;
00122 }
00123
00124 strcpy(fltfile,dtfile);
00125 i=strlen(dtfile)-2;
00126 fltfile[i++]='f';
00127 fltfile[i++]='l';
00128 fltfile[i++]='t';
00129 fltfile[i++]=0;
00130 fout=fopen(fltfile,"wt");
00131 if (!fout)
00132 {
00133 PrintError("%s: %m\n",fltfile);
00134 goto error;
00135 }
00136
00137 while (!feof(fin))
00138 {
00139
00140 if (ndata>=nalloc-1)
00141
00142 {
00143
00144 nalloc+=ALLOC_BLOCK;
00145
00146 assert(NULL!=(data=(struct xy *)realloc(data,nalloc*sizeof(struct xy))));
00147
00148 }
00149 line[0]=0;
00150
00151 fgets(line,254,fin);
00152
00153 if (sscanf(line,"%s %lf",sdate,&dt)==2)
00154 {
00155 if (str2time(sdate,&(data[ndata].x))>=5)
00156 {
00157 data[ndata].y=dt;
00158 ndata++;
00159 }
00160 }
00161 else
00162 {
00163 if (strstr(line,"time reset"))
00164 {
00165 str2time(sdate,&Treset);
00166 }
00167 }
00168
00169 }
00170
00171
00172
00173
00174 qsort(data,ndata,sizeof(struct xy),dtcmp);
00175
00176
00177 istart=0;
00178 for (i=0; i<ndata; i++)
00179 {
00180 if (data[i].x<Treset)
00181 {
00182 istart=i;
00183 }
00184 }
00185 Tmin=data[istart].x;
00186
00187 Tmax=data[ndata-1].x;
00188
00189
00190 N=1;
00191 while (((Tmax-Tmin)/(double)N)>MIN_WINDOW_LEN)
00192 {
00193 N++;
00194 }
00195 dt=(Tmax-Tmin)/(double)N;
00196
00197
00198
00199 for (n=0; n<N; n++)
00200 {
00201
00202 cnt=0;
00203 i=0;
00204 istop=0;
00205 while (i<ndata)
00206 {
00207 if (data[i].x<Tmin+dt*(double)n)
00208 {
00209 i++;
00210 }
00211 else
00212 {
00213 break;
00214 }
00215 }
00216 istart=i;
00217 while (i<ndata)
00218 {
00219 if (data[i].x>Tmin+dt*(double)(n+1))
00220 {
00221 break;
00222 }
00223 else
00224 {
00225 i++;
00226 istop=i;
00227 cnt++;
00228 }
00229 }
00230
00231 if (cnt>=MIN_REG_POINTS)
00232
00233 {
00234
00235
00236 double avg1=0.;
00237 double avg2=0.;
00238 double wg=0.;
00239 for (i=istart; i<istop; i++)
00240 {
00241 avg1+=data[i].y;
00242 }
00243 avg1/=(double)(i-istart);
00244
00245 for (i=istart; i<istop; i++)
00246 {
00247 double a;
00248 a=1./(data[i].y-avg1);
00249 a*=a;
00250 avg2+=data[i].y*a;
00251 wg+=a;
00252 }
00253 avg2/=wg;
00254
00255 for (i=istart; i<istop; i++)
00256 {
00257 if (fabs(data[i].y-avg2)>FIRST_PASS_RATIO*0.5*dtmax*dt)
00258 {
00259 PrintWarning("pass1: ignoring %s %f >%f\n",
00260 time2str(data[i].x),
00261 data[i].y-avg2,
00262 FIRST_PASS_RATIO*0.5*dtmax*dt);
00263 data[i].x*=-1.;
00264 }
00265 }
00266
00267
00268 cnt=LinearRegression(&(data[istart]), cnt,&a,&b);
00269
00270 for (i=istart; i<istop; i++)
00271 {
00272 double x,y;
00273 if (data[i].x<0.0)
00274 {
00275 data[i].x*=-1.;
00276 }
00277 x=data[i].x-data[istart].x;
00278 y=a*x+b;
00279 if (fabs(y-data[i].y)>0.5*dtmax*dt)
00280 {
00281 PrintWarning("pass2: rejecting %s %f >%f\n",
00282 time2str(data[i].x),
00283 data[i].y-y,
00284 0.5*dtmax*dt);
00285 fprintf(fout,"#%s %f dt=%f\n",
00286 time2str(data[i].x),data[i].y,
00287 fabs(a*data[i].x+b-data[i].y));
00288 }
00289 else
00290 {
00291 fprintf(fout,"%s %f\n",time2str(data[i].x),data[i].y);
00292 }
00293 }
00294 }
00295
00296 else
00297 {
00298
00299 for (i=istart; i<istop; i++)
00300 {
00301 fprintf(fout,"%s %f\n",time2str(data[i].x),data[i].y);
00302 }
00303 }
00304
00305 }
00306 error:
00307
00308 if (data)
00309 {
00310 free(data);
00311 }
00312 if (fin)
00313 {
00314 fclose(fin);
00315 }
00316 if (fout)
00317 {
00318 fclose(fout);
00319 }
00320 }