#include "stdio.h" #include "string.h" #include "astring.h" #include "amag.h" #include "flatfile.h" #include "stdlib.h" #include "math.h" int* remove_spike(char *filein, char *fileout, float maxv, float maxdbdt, int zerof) { char *vn, *vu, *vs, *vt, no[1000]; char stime[24],etime[24], outstime[24], outetime[24]; double suttime, euttime, uttime[1], *oldtime; int inunit, outunit, ncol[1], csize[1], error, i, found, *nptsm, nc; long srow, erow, nrow[1], orow, irow; float *data, *olddata, *startdata, *dbdt, missing, temp; float *avedata, tot, dt; missing = -1.0e32; avedata = ave_value(filein); inunit = ffopen(filein); vn = calloc(sizeof(char),ff_ncol[inunit]*16+2); vu = calloc(sizeof(char),ff_ncol[inunit]*16+2); vs = calloc(sizeof(char),ff_ncol[inunit]*31+2); vt = calloc(sizeof(char),ff_ncol[inunit]*2+2); data = calloc(sizeof(float),ff_ncol[inunit]); dbdt = calloc(sizeof(float),ff_ncol[inunit]); olddata = calloc(sizeof(float),ff_ncol[inunit]); startdata = calloc(sizeof(float),ff_ncol[inunit]); oldtime = calloc(sizeof(double),ff_ncol[inunit]); nptsm = calloc(sizeof(int),ff_ncol[inunit]+1); error = rdhead(inunit,vn,vu,vs,vt,no, stime,etime,nrow,ncol,csize); ncol[0] = n_columns(inunit); for (i=1;i<=ncol[0]-1;i++) nptsm[0] = 0; nptsm[0] = ncol[0]-1; strcat(no,",Spikes removed from file : "); strcat(no,filein); strcat(no," (max deviation : "); strcat(no,c_int_str((int) maxv,2)); strcat(no," ; max dbdt : "); strcat(no,c_int_str((int) (maxdbdt*100),2)); strcat(no,"e-2)"); error = mkhead(fileout,ncol[0],vn,vu,vs,vt,no); outunit = ffcreate(fileout); if (outunit != -1) { srow = 0; erow = nrow[0]-1; orow = 0; nc = n_columns(inunit); for (i=0;i<=nc-2;i++) { startdata[i] = avedata[i*6+1]; } for (i=0; i<=nc-2; i++) olddata[i] = missing; for (irow=srow; irow<=erow; irow++) { error = fdget(inunit,irow,uttime,data); for (i=0; i<=nc-2; i++) { if ((olddata[i] == missing) && (data[i] != missing)) { olddata[i] = data[i]; oldtime[i] = uttime[0]; dbdt[i] = 0.0; } else { if ((olddata[i] != missing) && (data[i] != missing) && (uttime[0] != oldtime[i])) { dt = (float) (uttime[0] - oldtime[i]); dbdt[i] = fabs(data[i] - olddata[i])/dt; olddata[i] = data[i]; oldtime[i] = uttime[0]; } else { dbdt[i] = 0.0; } } } found = 0; for (i=0; i<=nc-2; i++) { if (dbdt[i] > maxdbdt ) { found = 1; } } if (!found) { for (i=0; i<=nc-2; i++) { if (data[i] != missing) { if (fabs(data[i]-startdata[i]) > maxv ) { found = 1; } } } } if (!found) { tot = 0.0; for (i=0;i= 2) { for (i=2; i<=argc; i++) { if (strpos(argv[i-1],"-")==0) { flag = 0; if (strpos(argv[i-1],"z")>0) { zerof = 1; flag=1; } if (strpos(argv[i-1],"r")>0) { removef = 1; flag=1; } if (strpos(argv[i-1],"m")>0) { maxf = 1; maxv = atof(strmid(argv[i-1],2,strlen(argv[i-1])-2)); flag=1; } if (strpos(argv[i-1],"d")>0) { dbdtf = 1; maxdbdt = atof(strmid(argv[i-1],2,strlen(argv[i-1])-2)); flag=1; } if ((strpos(argv[i-1],"h")>0) || (!flag)) { helpf = 1; printf("usage : \n"); printf("aveflat [-rz] [-m#] [-d#] [flatfilein] [flatfileout]\n"); printf(" -r rm new datafile if npts removed is 0\n"); printf(" -z rm datapoint if all components are 0 nT\n"); printf(" -m# # is the max allowable deviation from the mean\n"); printf(" -d# # is the max allowable db/dt\n"); } } else { if (!ifilef) { ifilef = 1; strcpy(filein,argv[i-1]); } else { ofilef = 1; strcpy(fileout,argv[i-1]); } } } } if (!helpf) { if (!ifilef) { printf("Enter flat file to copy : "); gets(filein); } if (!ofilef) { printf("Enter output file name : "); gets(fileout); } if (!maxf) { printf("Enter maximum data perturbation (in nT) : "); scanf("%f", &maxv); } if (!dbdtf) { printf("Enter maximum dB/dt (in nT/s) : "); scanf("%f", &maxdbdt); } npts = remove_spike(filein, fileout, maxv, maxdbdt, zerof); totpts = 0; for (i=1; i<=npts[0]; i++) { totpts = totpts + npts[i]; } printf("points removed : %d\n",npts[1]); if ((totpts == 0) && (removef)) { printf("No data points removed - removing new data file.\n"); strcpy(tempf,"/bin/rm "); strcat(tempf,fileout); strcat(tempf,".dat"); system(tempf); strcpy(tempf,"/bin/rm "); strcat(tempf,fileout); strcat(tempf,".hed"); system(tempf); } } }