/****************************************************************************/ /* */ /* Spike.h */ /* */ /****************************************************************************/ /****************************************************************************/ /* This is a C language header file that defines routines for searching and */ /* removing spikes and similar structures from magnetometer data files. */ /****************************************************************************/ /****************************************************************************/ /* Lasse Hakkinen */ /* Finnish Meteorological Institute */ /* Department of Geophysics */ /* P.O.Box 503 */ /* FIN-00101, Helsinki, Finland */ /* e-mail: Lasse.Hakkinen@fmi.fi */ /* phone : (+358)-9-19294634 */ /* fax : (+358)-9-19294603 */ /* */ /* version 1.02 14.08.2006 */ /****************************************************************************/ /****************************************************************************/ /* Version history: */ /* */ /* 1.02 14.08.2006 Fixed the handling of end time in the situation where */ /* missing data was encountered. */ /* 1.01 06.11.1998 Fixed a bug where a spike went unnoticed if the average */ /* of field variations before the spike rounded to zero */ /* 1.0 04.11.1998 First official release */ /****************************************************************************/ #include #include /*--------------------------------------------------------------------------*/ /* Function prototypes */ /*--------------------------------------------------------------------------*/ static long ComputeDifferences(StationPtr s,char Comp,Time_sec Time,long Count); static long Max2(long N1,long N2); long LocateSpike(StationPtr s, char Comp, Time_sec CurrTime, Time_sec *SpikeStart, Time_sec *SpikeEnd, long TriggerValue,long AmplBefore, long AmplAfter,float SlopeParam); void RemoveData(StationPtr s, char Comp,Time_sec StartTime, Time_sec EndTime); void WriteSpikeInfo(StationPtr s,char Comp,Time_sec StartTime,Time_sec EndTime, long FieldResolution,long GapLength); /*--------------------------------------------------------------------------*/ /* Compute the average of successive field differences for given period */ /*--------------------------------------------------------------------------*/ static long ComputeDifferences(StationPtr s,char Comp,Time_sec Time,long Count) { long Sum = 0; long Good = 0; long i,Value; long X0,X1; X0 = GetDataValue(s,Time,Comp); Time += s->TimeStep; if (Count == 0) return (0); for (i=0;iTimeStep; } if (Good == 0) return (0); if (Sum/Good == 0) return (1); /* smallest possible field variation */ else return (Sum/Good); } /*----------------------------------------------------------------------------------*/ /* Find the largest difference of two successive field values for the given period */ /*----------------------------------------------------------------------------------*/ static long FindLargestVariation(StationPtr s,char Comp,Time_sec Time,long Count) { long Sum = 0; long Good = 0; long i,Value; long X0,X1; X0 = GetDataValue(s,Time,Comp); Time += s->TimeStep; if (Count == 0) return (0); for (i=0;iTimeStep; } return (Sum); } /*--------------------------------------------------------------------------*/ /* Find the larger one of two longs. */ /*--------------------------------------------------------------------------*/ static long Max2(long N1,long N2) { if (N1 > N2) return (N1); else return (N2); } /*--------------------------------------------------------------------------*/ /* Locate a spike structure and return the start and end times. The spike */ /* is located by first finding a jump in the field value whose magnitude */ /* is larger than TriggerValue and larger (by a factor of AmplBefore) than */ /* the average difference of N (=8) previous field values. */ /* Then the spike end is located by finding two successive field values */ /* that are reasonably close to the last good field value before the spike */ /* and whose difference is less than AmplAfter times the field variation */ /* before the spike start. */ /* If no spike is found then SpikeStart and SpikeEnd are both set to */ /* s->EndTime and 0 value is returned. */ /*--------------------------------------------------------------------------*/ long LocateSpike(StationPtr s, char Comp, Time_sec CurrTime, Time_sec *SpikeStart, Time_sec *SpikeEnd, long TriggerValue, long AmplBefore, long AmplAfter,float SlopeParam) { Time_sec t,t0; long X0,X1,Xbefore,Xup,Xdown; long MaxStep; long FieldDiff,Var; long N = 8; /* --- Find the start of the spike structure --- */ X0 = GetDataValue(s,CurrTime,Comp); t = CurrTime + s->TimeStep; while (t < s->EndTime) { X1 = GetDataValue(s,t,Comp); if ((X0 != MissingValue) && (X1 != MissingValue)) { /* Var = ComputeDifferences(s,Comp,t-N*s->TimeStep,N); */ Var = FindLargestVariation(s,Comp,t-N*s->TimeStep,N); if ((Var != 0) && (abs(X0-X1) > Max2(TriggerValue,AmplBefore*Var))) break; } X0 = X1; t += s->TimeStep; } if (t == s->EndTime) { /* No spikes found */ *SpikeEnd = s->EndTime; *SpikeStart = s->EndTime; return 0; } /* Now X0 is the last good data point before the spike */ /* and X1 is the first bad data point of the spike */ /* and t is the time corresponding to X1. */ *SpikeStart = t; /* --- Spike start has been located, now find the end --- */ /* The end of spike structure is found when two conditions are met: */ /* 1. The difference between two successive field points must be */ /* about same (by a factor of AmplAfter) as field variations */ /* before the spike. */ /* 2. The first point after spike should be reasonably close to the */ /* last point before spike. The interval where the first point */ /* should be increases linearly with time. */ Xbefore = X0; /* Last good field value before spike */ FieldDiff = AmplAfter*Var; /* Max diff of two successive points */ X0 = X1; t0 = *SpikeStart; t += s->TimeStep; while (t < s->EndTime) { X1 = GetDataValue(s,t,Comp); if (X1 != MissingValue) { Xup = Xbefore + SlopeParam*(t - t0); /* Max field value */ Xdown = Xbefore - SlopeParam*(t - t0); /* Min field value */ if ((abs(X1-X0) < FieldDiff) && (X1 < Xup) && (X1 > Xdown)) break; } else break; X0 = X1; t += s->TimeStep; } if (t == s->EndTime) *SpikeEnd = s->EndTime; else { if (X1 == MissingValue) *SpikeEnd = t - s->TimeStep; else *SpikeEnd = t - 2*(s->TimeStep); } return 1; } /*--------------------------------------------------------------------------*/ /* Fill data values from StartTime to EndTime with MissingValue. */ /*--------------------------------------------------------------------------*/ void RemoveData(StationPtr s, char Comp,Time_sec StartTime, Time_sec EndTime) { long *DataPtr; Time_sec T; DataPtr = GetDataPtr(s,StartTime,Comp); T = StartTime; while (T <= EndTime) { *DataPtr++ = MissingValue; T += s->TimeStep; } } /*--------------------------------------------------------------------------*/ /* Write information about the spike structure */ /*--------------------------------------------------------------------------*/ void WriteSpikeInfo(StationPtr s,char Comp,Time_sec StartTime,Time_sec EndTime, long FieldResolution,long GapLength) { char TimeStr[20]; long SpikeLength; long X0,X1,FieldStep; SecsToStr(StartTime,TimeStr); fprintf(stderr," %s %c : %.6s %s - ",s->StationID,Comp,TimeStr,TimeStr+6); if (EndTime == s->EndTime) fprintf(stderr,"Spike end not found. Nothing done.\n"); else { X1 = GetDataValue(s,StartTime,Comp); X0 = GetDataValue(s,StartTime-s->TimeStep,Comp); FieldStep = X1-X0; SecsToStr(EndTime,TimeStr); SpikeLength = EndTime-StartTime+s->TimeStep; fprintf(stderr,"%.6s %s%6d s %7.1f nT ",TimeStr,TimeStr+6, SpikeLength,(1.0*FieldStep)/FieldResolution); if (EndTime-StartTime < GapLength) fprintf(stderr,"yes\n"); else fprintf(stderr," no\n"); } }