subroutine tsf_edit (trg_opt,iuins,x,ndim,nx,nxr,t0,dt,w,nw,*)

  ... is called from a number of signal processing main programs.

 PURPOSE: Time series * Edit, with commands from instruction-file unit iuins.

 trg_opt   - char**  - trg = target, a search string to find a specific TSF EDIT
                       instruction block on iuins. The instruction block
                       identifier is
.......................'TSF EDIT '//trg
                       The trg search string ought to be right-delimited by ']'

                     - opt = options, can be appended to the search string.
                       Use ',' ' ' or '_' to concatenate.

                       Use '],Skip' to avoid searching the target at all.

                     For finding the instruction block:
                       Add ',Rwd' after ']' to rewind iuins;
                       Add ',NoR' after ']' to never rewind iuins;
                             \___ note: comma, blank or underscore are equivalent.
                       Else: File on iuins will be scanned from the current
                       position, wherever that may be, and be rewound exactly
                       one time if necessary.
                     Add ',P=2' to find a second block with the same target.

                     More advanced trg: Wildcards:
                       Wildcard character in trg is '.', must be given at
                       exact position.
                       Float-length wildcard '*' is only possible as the
                       first character.

                     More advanced trg: List of choices
                       If the first character of trg is any of !:@/|
                       this character is used as a delimiter between
                       trg choices. The choices designate substrings (only here!)
                       The substring may occur anywhere beyond 'TSF EDIT '.
                       The test is case-insensitive.
                       Example: trg '|ABC|DEF|XYZ]' will
                       find blocks marked 'TSF EDIT uvwDEFklm' due to 'DEF'

 iuins     - integer - instruction file unit 
 x(0:nx-1) - real*8  - time series in/out.
 ndim      - integer - dimension of x
 nxr       - integer - length of x at return.
 t0,dt     - real*8  - time of first sample, interval [h], input.
 w(0:nw-1) - real*8  - work array.
 *         - alt exit- taken if instruction block cannot be found.

 Time series editing:
 Works on x by reading+processing+adding w, scaling or filtering x, adding a constant to x
 and other operations on x.

 Unit iuins - Read instruction records from iuins to control the editing.
              Long instruction strings can be continued on the following line;
              indicate that by a lonely / at the end of the line that is

              Instructions comprise two types, file processing and in-core

The tsf_edit instruction file (`tse´)


TSF EDIT target

* You can jump into a block:
TSF EDIT target
TSF EDIT anothertarget

* Targets must not contain blanks.

* Long lines can be continued on the next line if indicated by a trailing ' /'

* Instructions may indicate a temporal scope (Keyword Time Range below.)
  Coding is either by putting the Time Range information at the end of the line or 
  using the DO instruction, e.g.

      DO From 2009 11 04 12 30 00 To
2009 11 04 12 40 00
      SCALE 0.987

(scope persists until DONE)

      SCALE 0.987
From 2009 11 04 12 30 00 To 2009 11 04 12 40 00

(scope comprises only the actual command line).
However, the DO instruction MUST be followed by an instruction that
   is neither a comment nor has a continuation line.

* A script may be constructed with the tse block in the file:
   set f=$1
   tslist $f -E$0,E $*

Control from environment
      Instructions may contain shell environment parameters, e.g.
       Don't use too keen tricks though. Prefer the ${ } form of quoting.
       There is more information on the REPLACE_ENVARS_L subroutine (~/util/afor/p/repenvl.f)

Some internal variables

      %Statistics variables:
      Computed with the STATISTICS, MEDIAN and SAMPLE commands and retrieved with %SYM symbols.

      At present (Nov. 2011) the procedure can handle some symbolic integer variables that can be written into the command line string. There is also a blunt internal counter.
      See the special section below.
      Just shortly, two examples.

Named variable SHIFTS:
OPEN 51 < file###SHIFTS#.ts
      would generate a file named file+015.ts

OPEN 51 < file###.ts
OPEN 52 < file###.ts
   would generate two files, file+001.ts and file+002.ts

  (It would be nice if we could) generate dates, e.g. Julian date as a counter
  ###JULD# and the translation into dates using ###YMD# / #YYMD#
  Then, e.g. OUTTS U=31 < data###YMD#.ts From ... To ...
  would first analyse From-To, then calculate JULD and then translate to YMD  
  sas_kalendar (teff,jul,idate,ihms)
  replace ### with either jul (ihashv(4)) or idate (ihashv(5))

     Comments are indicated by a ;-sign in the first position, or are preceeded by a !-sign
     later on the instruction line.
     If the !-sign is used, keep one character each to the left and to the right blank.
     The comments are printed to stdout.
     Quiet comments are indicated by a blank first position.

     Instructions are by default echoed to the stdout.
     This can be suppressed by coding a `/´ into the first column.

X and W
      X is the array passed through the call line and processed along the way.
      W is an array that may be input, see File Processing immediately below. In some instances it is used as a work area.

               FILE PROCESSING:

 Upon instruction, tsfedit calls read_fuf for time-series file reading.
 A file instruction is assumed if the first word on the line is not resolved as a command.
 Let designate the time series currently in core, the time series that is to be
 read, and x.out the result.
 Instruction consists of the following parameters:

     iun, target, recmrs, fmt, khms, itz, cdir, value  [more options]

 Example (to subtract a filed time-series from the one currently in memory):

 22, 'BIN]',-99999.0,'U', 0,   0,  'N', -1.d0  DC-

 iun...itz are standard parameters the meaning of which is documented 
 also at other places. Special to tsfedit is cdir and value

 iun     - file unit number.

 target  - max 8 characters: Position file at "target";
           to rewind before positioning, specify:        "target] Rwd"
           "TOP]" - read ascii file from the top.           "TOP] Rwd"
           "BIN"  - binary data (getts family)               "BIN Rwd"
           "SAC"  - seismological SAC files  (read_sac) 

           "BRX"  - Bruxelles format                         "BRX Rwd"
                    also: ETERNA format
                    the target specifies an END-OF-FILE mark.
                    c.f. readbrx.f.

           C.f. proc_f_on_rec  or read_fuf for more options that can be
           passed through target.

 recmrs  - Missing record symbolic value (in formatted data sets only).

 fmt     - Case ASCII data: The time and data format in FORTRAN, for
retrieving one data column.
           Case Binary data, MultiComp, the label of the requested segment

                  "l:text" "L:text" "*:text"

           Case Bruxelles: Sub-format for samples, e.g. 'F9.0'
                  or 'F10.0:ETERNA' for ETERNA files

           other (ASCII): The data format in FORTRAN.
           This argument is passed to read_fuf.Thus,
           can be appended to describe format conversion for the date part.

           C.f. proc_f_on_rec or read_fuf for more options that can
           be passed through fmt.

 khms    - depth of time information (only ASCII files non BRX).
           The three items of the date are mandatory. The additional time items
           are hour minute second and millisecond. This is designated by khms,
           a value between 0 and 4.

 itz     - Time zone w.r.t UTC.
           The time records in the file are possibly referring a zone's time while
           the program wants t refer to UTC.

 cdir    - char*2 - meaning                          needs more data / instructions
...........A* - multiply with value and append                     NO
           N* - nothing, multiply with value and add data as is    NO
           S* - apply a scaling factor and add        1 record:  scale factor
           F* - apply a filter and add                1 record:  filter params.
                                                    + n values:  filter coeffs.
           L* - bias and scale, and add               1 record:  bias1,sc1,bias2,sc2
           +* -
bias and scale, and stack             1 record:  bias1,sc1,bias2,sc2

           G* - use as gap filling data               1 record:  gap parameters
           J* - use as jump data                                 NO
...........K* - scale, (mask) and Keep in work array W
...........KT - scale, (mask), filter (with the SAVE'd one),
                and Keep in work array W
           T* - Filter (using the SAVE'd filter) and keep in work array W
           P* - Compute phase from arctan X/W
           C* - add series cyclically.
           *t - override dt with the one from the data file

           The options N S F L T and C request editing of the internal time-series
           by means of linear combination with the time-series input here.
           One effect of the editing may be an undesired truncation to the
           length of the most recently input time-series; specify lower-case
           characters n s f and l to keep the length unchanged.
           Scaling factor (S):
              x <- ( x + ( w * value * sc ) )
                              where sc is the read-in scale factor
           Cyclical addition (C):
  x.i <- ( x.i + ( w.j * value ) ) for all i
                              with j = mod(i,lw), where lw is the length of w.
(itz may be used for an origin shift (really?))

           Bias and scale:   [bias1 [, sc1 [, bias2 [, sc2]]]] [-DC]
                               x  <- ( x + ( w + bias1 )*sc1+ bias2 )*sc2 )

                              Default: 0.0, 1.0, 0.0, 1.0

           Filter parameters:  scale,nfminus,nfplus,symmetry

..............................Example: 1.0, 0, 16, 'S'
                              for a symmetric filter with 16 unique coefficients
                              i.e. 31 coefficients when expanded.

           OBS! No interpolation if series aren't in synch.
           Gap parameters:     scale,bias,strategy

           strategy:           '[X{D|H}]A]'                  to replace missing and append all.

                               '[X{D|H}]A [begin end]'       to replace missing and append,
...................................................         .begin end are sample numbers
                               '[X{D|H}]R [begin end]'       replace missing, no appending, default = all
                               '[X{D|H}]F [begin end]'       force over, no appending,        "        "

           The options XD or XH request read_fm_exactd('D') or ('H'), respectively.

           Defaults for scale and bias are not available; record is read with
           The following operation takes place
...........x.out = ( if valid and not Force )
                   ( else if valid ) * value * scale + bias

 value     - ampl.factor, used with cdir = 'A*' or 'N*' or 'L*' or 'S*' or 'K*' or 'T*' or '+*'

more options
  -DC       - Remove DC-level of the time-series W

more options pertaining to cdir='L*' or '+*'
           With cdir='L*' and the -DC option (? need to check that one),

           the following operation takes place

..........   .   X.out = sc2· + value·sc1·( - DCLvl + bias1) + bias2

           With cdir='+*' and the -DC option
           the following operation takes place

...........      X.out = X.out + sc2· + value·sc1·( - DCLvl + bias1) + bias2

                 X.out(i) = X.out(i) / Nstacks(i)     at TSF EDIT END

 +&&       - Use this data series to mask the one in memory (inject missing record symbols).
             This masking (still) uses a simplified scheme. The rates of the two series are
             *assumed* equal.
 +SQ       - Add this data series to the one in memory by squares.
 +RR       - Remove this series from the one in memory by regression
 M>0       - Fill up with zeros beyond file end, as much as available.

 *X [-DC]  - with  cdir='L*'   multiply with
                X.out = sc2
··( [-DC-level] + value) (OBS! +value)

 X/ [-DC]  - with  cdir='L*'   divide by
                X.out = sc2
· [-DC-level] + value) (OBS! +value)

         - divide by scale factor (cdir = 'N*' or 'G*')
 /.        - divide with sc0, multiply with sc1
 ./        - multiply with sc0, divide by sc1
 //        - divide by sc0, divide by sc1
 +++       - Stack this data series (at return the series will be normalised by the
             number of stackings at each point. Stack with equal weight, no bias.
             (If you need to rescale and bias, use cdir='+*'
      - cyclically (affects only stacking so far): W is wrapped over X
     - stack cyclically.

!T=[X][u:]#v  - move time origin by v, units of hours except when
              u=S - #v in seconds
              u=D - #v in days
                X - move time origin to the one of series X and add #v)
{+|-}ST - Summer time adjustment for e.g. stacking days synchronous with civil time.

           Note that `/´ as the last character on the line is interpreted as the
           a continuation mark. To specify the `/´ `/.´ or `//´ options last on
           the line, add ` -´ or any string that is not an option.

                     In-core processing:
Parameter symbols are shown in  italic, compulsory text in  bold, optional text in [brackets].
A #-sign before a symbol indicates that the system expects a numerical value.
Time range definition                      Subroutine time_range                     File  ~/sas/p/timerange.f

By dates

Let time-range :=

 { [ From #YYYY #MM #DD #HH #MM #SS #FFF ] [ To #YYYY #MM #DD #HH #MM #SS #FFF ]  |
   At #YYYY #MM #DD #HH #SS #FF  |
   ALL }

  An empty time range means (usually) all. 
  If the command has no other parameters than a time range, the At codeword can be dropped. 
  Example:  DEL 2001 1 1 0 0 1 0

#p,#q #YYYY...  means that the range extends from  p  samples to  q  samples centred on the given date.
 Example: For a symmetric interval specify  Around:-30,30
Obs! No blanks in this expression.

The year month and day can be replaced by  YMD[{+|-}#n]   to mean the day of the epoch ±, optionally,  n  days.
Example: DEL To YMS+14 12 0 0 0 0

By sample numbers
{ [ From ##j ] [ { To | For } ##k ] | [ At ##j ] | Around:#p,#q ##j }

##k   is equivalent to   For ##(k-j+1)
 Example: PERIODOGRAM /N From #740 For #4230
 Sample numbers are in the range 0 .. NXR-1

Introducing a search margin (or a shift):
#m1 #m2     - m1 and m2 is in number of samples. The margins can be incremented:
TRMARGSTEP [#ms]   - add ms to both m1 and m2

COMMAND Time Range

Example: apply operation to the sample immediately before 2011 07 04 12 00 00 00
TRMARG -1 -1
SCALE  9. At 2011 07 04 12 00 00 00  
  N.B. range by dates:  After  FromTo  and  At the seven time fields must all be specified, else a reading failure will result.
  The read command used is 
        READ (INSTR(I:), *,ERR=99,END=99)  IDATE, IHMS
This remark is probably obsolete

 Keyword                         Meaning                                        [default]   

 DO Time-range                         In a DO - DONE block time-range expressions beyond the
 DONE                                  DO will not be evaluated. There is no stack of ranges (yet)

 TRMARG #m1 #m2                        Margin for time range, m1 and m2 is in number of samples.
                                       The margins can be incremented:
 TRMARGSTEP [#ms]                      add ms to both m1 and m2

 TRMARGSTEP [#m1s #m2s]                add m1s to m1 and m2s to m2

 LOOP #n [#s] [#a] <label>             Loop #n times (step #s and add #a, all integers).
<label>                       A Label is compulsory, must be keyed with brackets < >

 DEBUG [END] [ L={#n|D} ]              Diagnostic printing. Toggled;
.......................................END can be stated for clarity.
                                       The message level in some ts-routines can be set with n
                                       (0=debug, 1=much, 2=normal, 3=important-only, 5=quiet)
                                       L=D for the level in effect at call time.

 ECHO                                  Echoes the commands and reports time-range
                                       analysis to STDOUT. ECHO is the default.
 NOECHO                                In the case of samples that become deleted
                                       a report of the delete command and the number
                                       of deleted samples will be issued.


...................................F I L E   C O M M A N D S

 OPEN                                  The next lines contain an Open-a-file block.
                                       Last line = `b b b

 OPEN s                                String s = `ii o filename´

                                  ii - unit number or `*n´ (an asterisk and a number between 1 and 9).
                                       If the latter form is used, an unused file unit number will
                                       be searched for. In a unit number specification of the kind U=
the symbolic value can be given in the `*n´ form.
                                       After this eventual replacement, string s is passed to

                                 ### - File name may contain `###´, which will be replaced
                                       by a running signed three-digit number. One occurance
                                       per file name only!

                                       Some special variables can be accessed:
                          ###TRMARG# - will be replaced by the value of TRMARG,
                                       the "time-range margin"
                          ###SHIFTS# - the accumulated value of origin shift a,
                                       from SHIFT S=s : a = sum s
                                         or SHIFT T=s : a = sum s/dt

 CONT[QUIET] #n                        next line will be an Open-a-file block.
                                       Instruction input will continue
                                       from unit n (only in the scope
                                       of the current TSFEDIT invocation).
                                       Cannot be nested, just chained.
                                       The QUIET version will suppress ECHO until
                                       the end on unit #n.

 CONT[QUIET] s                         String s is passed to openfs, specifying
                                       file unit, option, and file name. Otherwise
                                       like 'CONT #n'  above.

 CLOSE #n                              Close file unit #n. Any file may be closed as needed.
                                       If the TSFEDIT instruction file unit is
                                       closed, input will continue from the
                                       file unit defined on the call list.
                                       That one cannot be closed.

 DUMP  #u [L=label]] [O:opt] [ time-range ]    Write (binary) the time series to the file on unit u.
[L=label]] [O:opt]                   Write (binary) the work array to the file on unit u.
                                 opt - Options:
                                  +E - Strip off leading Missing-records and adjust epoch data written to
                                       the file header.          

 OUTTS [W] U=#u < filename] [O:opt] [ time-range
                                       Open the file, write the time series
(binary), and close.

 OUTTS #u [O:opt] [ time-range ]      
OUTTSL #u L:label [O:opt] [ time-range ]      
                                       Short forms to write series X to an already opened file.

 OUTTSL [W] U=#u < filename] L:label [O:opt] [ time-range ]
Open the file, write labeled multi-column file, keep open.

 OUTTSL [W] U=#u L:label [O:opt] [ time-range ]
                                       Append another column, keep open.

                                       The part `U=#ub<bfilename´ has a strict format,
                                       one blank inbetween. Acceptable options besides < are > % +

                                       see openf.f
. Note that the file name is right-delimited by `]´
                                       File name may contain `###´, which will be replaced
                                       by a running three-digit number.

                             L:label - can be coded    PART1,PART2   or   PART1_PART2  or  PART1______PART2.

                                   W - Output the work array.
                                       N.B.: The time-range as well as the origin of W are assumed
                                       equal to the master time-series.

                                opt - Options:
                                  +E - Strip off leading Missing-records and adjust epoch data written to
                                       the file header.

                          time-range - not valid with W.          

..............................I N F O R M A T I O N   C O M M A N D S

 MISS [ time-range ]                   Report missing samples

 PRINT #u [F=fmt] [time-range]         Print time series X in free format or format fmt on unit u
[F=fmt]                     Print time series W ...
 PRINTF #u [F=fmt]                     Print filter ...
                                fmt -  One integer field followed by one real*8 field.

 REPDC [opt] [ time-range ]            Report DC-level to protocol (unit 6)
                            opt  +TY - include full date information and decimal year
                                  +Y - ... decimal year
                                  +T - ... full date
[/]->S[-] - set scale parameter (the the inverse value if /-S ),
                                       invert the sign if appended by `-´,
for use in
                                       a range of arithmetic commands

                                       Example for +TY:
                                       REPDC +TY  From 2009   0 187 13 35 18 00 To 2009   0 191 13 25 13 59
                                       ` <TsfEdi>>> REPDC 2009 11 06 10 28 18 000  2009.860195 -2.844177E+01´

 RMS  [ time-range ]                   (ignored, obsolete:  Report RMS )

 RMSD   [opt] [ time-range ]           Report RMS about mean.

 RMSDW  [opt] [ time-range ]           ... do that on the work array W. However, here the time-range
RMSD W [opt] [ time-range ]           may be a pitfall (it relates to the origin and rate of time-series X);
                                       prefer to use array indices instead.

                            opt  +TY - include full date information and decimal year
                                  +Y - ... decimal year
                                  +T - ... full date
                           [/]->S[-] - set the scale parameter to the result (or 1/result if /->S ),
                                       invert the sign if appended by `-´,
for use in
                                       a range of
arithmetic commands.
                                       The reported dates are derived from the origin and rate of the
                                       respective time-series.

 MEDIAN [s-opt] [ time-range ]         Report median
                               s-opt - [/]->S[-] store the result as scale value for use in
                                       a range of 
arithmetic commands.

 WDR #t [C=text]] [U=#u] [M[c]=#ml,#mu[,c]] [ -DC[=#dc]] [ time-range ]
                                       Write records for a given signal threshold t.
By default the
                                       records are written as tsfedit DEL commands..
                             C=text] - text will replace DEL.

                                       Write a record for each group of samples x with | x - dc | > t
                                       If t<0, the threshold is |t| times the result of the most
                                       recent RMSD command.
                                 -DC - Determine dc, the series' DC-level, else use dc=0
                             -DC=#dc - Calculate deviation around dc.
                                U=#u - Write to file unit u. Default is 6
M=#ml,#mu - Output time ranges with a lower and upper margin of +ml +mu
                                       of the sampling interval. E.g. -1,+1 to extend the deletion
                                       window by one sample at both ends. Default -0.5,+0.5
                          MH=#ml,#mu - Margin of +ml +mu in units of hours    ( M=#ml,#mu,h  should also work)
                          MS=#ml,#mu - Margin of +ml +mu in units of seconds  ( M=#ml,#mu,s  should also work)

                                       Consider to use WDR on the result of MOVRMS (don't opt for -DC then)

 WDR MISS [C=text]] [U=#u] [M[c]=#ml,#mu[,c]] [ time-range ]
                                       Like WDR, writes DEL records for currently missing samples.

...................................A R I T H M E T I C   C O M M A N D S

 CUT #n                                Cut away n samples from the end
 RETAIN time-range                     ... with adjustment of starting time and length

 X->W                                  Copy series to work space
 W->X                                  ... vice versa

 X*W O=s [+P] [ time-range ]         - Uses mult_ts_ts (sasu16.f) to multiply the series, result in X
 W*X O=s [+P]
[ time-range ]         - Uses mult_ts_ts (sasu16.f) to multiply the series, result in W
                                   s - Operator and options, / for divide, Q for squaring the second
                                       series before the operation, L for Let, 0 for replacing missing
                                       records with DC-level.

 SUMTS[X|W] [[/]->S[-]] [ time-range ]
                                       Sum the series X or W, respectively. Default: X
                                       Set scale parameter
(the reciprocal value if `/S´ , append `-´
                                       to reverse the sign) for use in a range of
arithmetic commands

 REMDC [ time-range ]                  Remove DC level. The DC-level of the time range is
                                       subtracted within the time range.
 ADJDC [ time-range ]
                  Remove DC level. The DC-level of the time range is
                                       subtracted within the full scope of the time series.
 GETDC [ time-range ]                  Saves DC level internally. See ADD, SUB, etc. further below.

 DETREND [P|I] [ time-range ]          Detrend (slope and bias)
                                   P - Poor man's detrend (defined by start and end point)
                                   I - Inverse (unweighted) using DETREND_UNW

 POLYRED #n [O=[P][I][A] [ time-range ]   Fit Legendre polynomials up to degree n and ...
                                   I - return the polynomial prediction (default: the residual)
                                   P - diagnostic printing
                                   A - Analyse-only mode. Option
I is not applicable then.
                                       Notably, with AP, the analysis will show the Akaike criterion 
                                       for increasing model order
From time                    Prepares POLYRED to insert (i.e. solve for) a jump at the
From´ time. Time-range syntax applies.
                                       The POYRED time-range should be narrow around the jump time. (?)
                                       The jump will be applied to all data at and after the jump.
                                       POLYREDJUMP should be followed by a POLYRED command;
                                       since analysis mode is necessary, this is set automatically.
 POLYREDJUMP #n M=#j,#k [O=[P][I][A]]
From time  
                                       All-in-one command, if the time range for the polynomial
                                       can be specified as j samples before and k samples after the
                                       `From´ time.

 SCALE [/] #x  [ time-range ]          Multiplies series with ...
                                   x - scaling value
                                   / - Divides series with x
 SCALEM [/] #x  [ time-range ]         Scale and map to ordinate interval [-x,+x]
                                       ( with / to [1/x,1/x] )
 SCALER [/] #x  [ time-range ]         Scale to x·RMS-deviation
                                       ( with / to 1/x·RMS-deviation )

 SCALE [D] N={#p|00} [ time-range ]   
                                       Scale to unity norm
                                   p - type of norm, 2 yields the usual unity-variance.
                                       N=00 means infinity-norm (maximum=1)
                                   D - "Deviate", take norm with respect to
mean value.

 TDGAIN [ D | D=#v ] [O:opts]        - xi <- xi f(i,wj) where j and i are near each other in time
                                       and f() interpolates linearly between j and j+1.
                                   D - delete the value if j is out of range
                                D=#v - wj = v if j is out of range.
                              O:opts - Options to subroutine apply_gain_ts (~/sas/p/shrtadms.f). 

 NINT #s #x [I] [ time-range ]         Transform data by rounding
                                         X.out <- NINT(( - x)/s)·s
I - Keep integer result.

 MODT #s #x [K] [ time-range ]         Transform data by modulus
                                         X.out <- MOD( - x, s) + x
K - Don't add the x back. 

 POLYTR #n [ time-range ]              Polynomial transformation using subroutine poly_nonlin_ts.
 p0 p1 ... p#n                         (sasu06.f)
                                   n - polynomial order
                       p0 p1 ... p#n - coefficients 
                                         X.out <- SUM ( p(j)^j , j=0..#n )

 ABS    [ time-range ].................Take absolute value of
X**#n  [ time-range ].................Take power n of (inject MRS if negative).
 EXPTR  [ time-range ].................Take exp of
 LOGTR  [ time-range ].................Take log10 of (inject MRS if negative).
 DB  [ R=#v ]
[ time-range ]...........Convert to dB (inject MRS if negative). Default v=1.
 MAX [ > #v ] [ time-range ]...........Take max(v, Default v=0.
 SQRT   [ time-range ].................Take square root of (inject MRS if negative).

 PDWB [U=#iu] T=target                           Process a downweight block from the specified file.
                                  iu - file unit
target - target string
.......................................Command discarded, would lead to recursive calls  

 SHIFT #n                              Shifts time series origin using sasu0111.f shift_ts  
                                   n - number of samples

 SHIFT S=#n [+R]                       Shifts the series by n samples without affecting the
                                       origin (using sasu0111.f dodge_ts). Pads with MRS.
                                       Shifts are accumulated in ###SHIFTS# (
                                   n - number of samples
                                  +R - resets
###SHIFTS# first.

 SHIFT T=#t[s|r] [+R]                  Adjusts the time series origin t0 -> t0 + #t
###SHIFTS# action as above, incl. +R
                                       The string
T=#t[s] must be coded without blanks.
                                  #t - additive adjustment in hours or, if s is given, in seconds
                                       or, if r is given, in units of sampling rate

                                       Example: Interpolate a series half ways between given samples
                                       FILTER 0 1 ' '
                                       0.5 0.5
                                       SHIFT T=-0.5r

............................F I L T E R   A N D   D S P   C O M M A N D S

 This concerns the following commands:  FILTER FILTS PEF WIENER IIR REGRES:
 Note that filters and regression don't work with time-range. Read in to exactly match first sample and TRUNCate before use.

 REGRES                              - After reading a time series with cdir = K or T (to keep it
                                       in W, unfiltered or filtered), W will be subtracted from X
                                       by regression (CD-levels are removed in the process).
                                       W may be longer and start earlier than X; W will be tailored.

SAVE          - Compute or get a filter, save but don't apply it yet.
                                       Allows a filter of max length 200,000, uses storage location 1.
FILTERCOMMAND options STORE[->#n]   - Like SAVE, however specifically in storage location n               [n=1]
                                       Allows 5 filters of max length 40,000 or 1 filter of length 200,000
FILTERCOMMAND [+C] APPLY[<-#n]      - Apply a saved filter                                                [n=1]
FILTERCOMMAND [+C] FETCH[<-#n]      - Fetch a saved filter                                                [n=1]
FILTERCOMMAND [+C] options          - Compute or read, and apply instantly.
                                  +C - Apply on a circularly connected domain

 FILTER FORGET[->#n]                 - Maybe useful after FILTER APPLY
                                       (works only with those filters that use COMMON /CFILTER/ F(0:nfdim)...
                                       ZSPF in the first case.
                                       FORGET->#n will wipe only the filter stored at n.
                                       FORGET will wipe the filter stored at 1 and the active filter

 NORMFILTER [+N][+H][+U]          +N - Give filter unity-response at frequency 0. Default is +N
 HIGHPASS   [+N][+H][+U]          +H - Convert to High-pass
. Default is +H
                                +N+U - Give filter unity-response at Nyquist frequency
                                       can be applied between SAVE and APPLY

FILTERINPLACE FILTSINPLACE ...        In the following, the forms without INPLACE are described. The
                                       INPLACE variants have exactly the same arguments and options.

 FILTERINPLACE [time-range]            is possible. However, the work array W will have to be used.

                                       The INPLACE variants apply the filter such that the ends of the
                                       input series are preserved (length out = length in).
                                       This may lead to stark jumps in high-pass filter, but may
                                       return useful results in the case of low-pass and averaging filters.
                                       In high-pass filters the ends should be set to zero. Command ...
 INPLACE 0                             ... provides the corresponding option.
 INPLACE H[IPASS]                     
alternate form
anyother                      preserves the signal (switches the option off).
                                       See example 11 near the bottom of this page.

  GETFS U=#u[ < filename]] O=opt M=string]
                                       Input from a BIN file (written with OUTFS, sas/p/filterios.f)
                                       Implies the SAVE function. To apply the filter on the data,
                                       add a line FILTER APPLY. To disable, use FILTER FORGET
                                 opt - B   - rewind first
                                       R   - allow one rewind (works only with search string)

string - input the filter that has been marked with string

                                       See Example 16

                                       To convert a filter into a time-series, use
                                       ADD 1.0 At #0
                                       ZSPF POWER:0 +DISC +D +BL -P [F]=Hz

                                       To get a list of filters stored in the file, use program

 OUTFS [X> opt] [+D] U=#u[ < filename]] M=string [ time-range ]
                                       Save a filter to a file. Open by file name. To open for append,
                                                                 use `>´ instead of `<´ .
                                       With only `U=#u´, use an already opened file.
                            M=string - Make filter searchable with the string.

                              X> opt - Convert the time series into a filter before output. See Example (13)
                                 opt - additional options for conversion:
                                  +S - symmetric, copy the indicated time range as f(0..), the positive side.
                                  +Q - antisymmetric
                              +A:m,p - asymmetric, indicated time range goes into f(m..p)
                              +B:m,p - "back-end", x(0..p) goes into f(0..p) and
                                       x(N-1..N+m) goes into f(-1..m) (Obs: m<0).
                                  +D - Show filter coefficients.  

 OUTFS [X> opt] U=#u M=string          Add a filter to the file already open on unit #u.
                                       "A filter" = the currrently defined filter.

 TS2FS {X> | W>} opt                 - convert time series X or, respectively, W to a filter.
                                 opt - as under OUTFS. Set DEBUG to print the coefficients.


 FILTER D:mmp string                 - Create a filter
                                  mm - a short string designating the method. Currently implemented:

                                  SW - a simple extreme-low-pass filter using a window
                        WD[,P][C][D] - a bandpass filter using the window-design method
[C] - WD with a fix to get rejection at zero-frequency for
                                       a high- or band-pass filter.

              SW[H|L]   SIMPLE WINDOW (= weighted) AVERAGING
                                   L - applies the coefficents straight away.
H - for the  complementary filter
                                           f(i) <- -
f(i) for |i| .ne. 0, and
f(0) <- 1 - f(0)
                              string - WTYP #tpoint [d]
Window type, truncation point, design parameter(s)
                                WTYP - char*4 Window type:
                                       RECT BART HANN HAMM KAIS DOCH HANQ KBTA KBTB
cf. ~/sas/p/window.f

              GB   GAUSSIAN BELL   cf. sas/p/fgauss.f                                       
                              string - #L #f [#fny]
                                  #f - -3dB-frequency as a fraction of the ...
                                #fny - Nyquist
                                  #L - filter's positive half-length (filter is symmetric)

              WD[,P][C|D] or WDZF[,P]    WINDOW-DESIGN METHOD   cf. sas/p/wffs.f
                                   P - specify periods instead of frequencies. To specify
                                       infinite period, give a 0.
                                   C - return the hi-pass complement of the low-pass
                                       declared in string
D - modify to exactly suppress DC signal:
                                       subtract h = (sum f(i))/Nf
                                       (degrades rejection at high frequencies!)
                                       This is different from the ZF-method which uses
                                       window times h

                              string - WTYP [I] #tp #d #Flo #Fhi #Fny #Fcal

                                WTYP - char*4 Window type:
                                       RECT BART HANN HAMM KAIS DOCH HANQ KBTA KBTB
cf. ~/sas/p/window.f
                                   I - use band limits rounded to integers                [float]
                                 #tp - truncation point = half-length of filter
                                  #d - design parameter
                                #Flo - lower band limit
                                #Fhi -
upper band limit
                                #Fny - Nyquist frequency for scaling the following parameters
(or -#f for scaling relative to 0.5/dt/f  N.B.: [dt]=hours)
                               #Fcal - calibration point
                                       E.g. if you want to provide the band parameters in Hz,
                                       specify #Fny=-3600.

                                       Example: A free-oscillation filter independent of the actual
                                       sampling rate. The first line specifies frequencies in Hz,
                                       the second line periods in seconds.

                                       FILTER D:WD   KAIS 121 2.0 0.0002 0.0025 -3600. 0.001
                                       FILTER D:WD,P KAIS 121 2.0 3800. 400. -3600. 1000.

                These filters can be explored with sasm02:
                Examples:              sasm02 -f WD HANN 1200 2.0 0.0 0.01469 5.0 0.0
                                       sasm02 -f WDC HANN 1200 2.0 0.0 0.01469 5.0 0.0
(unfortunately not with the P option!)

                BF       BESSEL FILTER (transformed to the time domain)
OBS! This is a one-sided IIR filter; truncation problems may arise.
                         In ZSPF you can try with its spectrum-domain equivalent.
                              string - #L #p #fc [{mHz|Hz|cpd}] [P=#p2,#p3...]

                                  #L - filter length
- filter order (e.g. 8)
- corner frequency
                        {mHz|Hz|cpd} - units for #fc, Default is 1/dt
- additional design parameters (not used); #p1 is identical to #fc.

                                       (not tested yet)

#m #n 's' [+N]                 Input from an ASCII file

  f(m),f(m+1),...,f(n)................- filter with explicitly supplied parameters
.......................................Filters the time series. Read filter coefficients
 ......................................from instruction file.

.............................#m #n s - filter length (begin end) and symmetry.
                                +N   - normalise coefficients to yield unity response

                                       Symmetry codes:
                                       'S' - symmetric,
                                             #m=0, #n=n, =>  f(1)..f(n) is mirrored to

                                       'Q' - antisymmetric,
                                             #m=0, #n=n, =>  f(1)..f(n) is mirrored to
                                                           -f(-1)..-f(-n) and f(0)=0

                                       any other code: general filter f(m)..f(n)

 IIR #n                             
- autoregressive (infinite impulse response) filter
                                       with #n+1 coefficients.
 IIR -1                              - uses the coefficients that have been input under FILTER or GETFS

 IIR #n NONLIN=ζξ                              - uses the coefficients that have been input under FILTER or GETFS

 FILTS U=#u [ F=f ] P=#m #n 's']     - input a time series from unit #u and use as a filter.
                                 P=  - Same rules apply to #m #n 's' as with
>loc> in FILTER
                                       except for the coding of the command line:
                                          either commas are used and blanks are avoided
                                          or string is right-delimited with ]
                                  #m - Sets the anticausal start point. The most common
                                       filter type would have #m=0 and 's'=' '
                                  #n - Any big value will be accepted. The actual parameter
                                       is calculated from #m+L-1
                                       where L is the length of the time series.
                                   f - if a file name is given, the file will be opened
                                       else an open file is expected on unit #u

 FILTER U=#u F=f [+N] [>loc> [O=opt]]
.....................................- file instruction
.......................................File unit #u  and name f are passed to openfs,
                                       location string loc and option string opt passed to
                                       FILTER U=23 F=bandpass_1_2.fpf >-16     16 'A'> O=AB

                                       The file must begin with one record supplying #m #n 's'
                                       and the search string, if used, must match exactly (!)
                                  +N - normalize filter gain

 WIENER U=#u F=f L=#lw                 Process a time-domain transformed gain function
                                       (a.k.a. impulse response), for instance output during
                                       sasm03 step response display and time-limit it
                                       using a HANN window of full duration #lw
                                       File record: read (iu,*)
                                       i f(i)


 PEF [U=#u F=f] L=#l O=opt [+F]        Use a Predicition Error Filter, length l
 PREDF [U=#u F=f] L=#l O=opt [+F] SAVE Prepare the corresponding prediction filter for use
                                       under IIR APPLY

Padsasdsd.............Use a Predicition Error Filter, length l
.......................................If options U= and|or F= are used, read
from filter bank
                                       from file unit u, file name f. If U= is not used, a unit
                                       is assigned automatically.

                                       Else the filter bank prepared by a preceeding BURG command
                                       is used.
                                       Option opt is passed to
read_pef or take_pef, respectively
                                       (see ~/sas/p/pefs.f)
                                  +F - Swap the filter type (PEF <-> PRF) by brute force.
                                       If not used, the type is assumed as detected from the file.
#p [ #p ... ]              Test the currently defined filter's harmonic response
                                       at period #p given in hours (max. 10 values).
                                   S - Periods are given in seconds.

 FILTHR TAPRL=#u[ o filename]]         Process an urtap result sheet (urtap.prl), read from unit u.
(Use is for information only).
                                       File must be declared in an OPEN command unless option o and
                                       filename are given here.

 FLIPF                                 Flip orientation of filter from Forward to Backward or vice versa.

 BURG #l [O=opt] [ time range ]        Creates a prediction-error filter bank, max order l. You need
 [ #u < filename ]                     the PEF command to apply a filter from the bank.
                                       Options opt
                                       O - enable output; one more instruction line is needed to open
the file (binary!) using openfs (start in column 1):
                                           #u < filename
                                           File name may contain `###´, which will be replaced
                                           by a running three-digit number.

                                       S - scale filter for unity power output n 

 UPDFILTER{A|M|R} #i #v                Update filter coefficient i:
                                       A - add v
                                       M - multiply with v
                                       R - replace with v

 MEDIANFI[LTER] #l [M=#m] [STORE]      Filters by running median, segment length #l
                                       Can be moved foreward by m (default 1) for decimation
                               STORE - Will not apply the filter. The parameters are set aside
                                       for subsequent MEDIANFILTER APPLY or
                                       file import to array W under option 'T.' ]


 WINDOW[:v] type [P=#p] [T=#t] [ P ] [U] [ time range ]     
                                       Set up a window and, if ` P ´ ("prepare") is missing,
                                       apply it on the data.
                                       The length of the window is determined from the time-range.
U - unit-normalise the response (window shape and missing
                                       data will modify the response).
                             :v = :W - Apply it on the work space, else on the primary data.
                             :v = :F - Apply it on the filter currently stored away with SAVE
                                type - 4-character code, see ~/sas/p/window.f
P=#p - Design parameter if needed. Default is good for KAIS         [1.5]
T=#t - Delayed taper, 0 < t < 1   [0.0].

 WIN[:v] [U] [ time range ]            Use the currently defined window and apply it on the data. 
                             :v = :W - Apply it on the work space, else on the primary data.
                             :v = :F - Apply it on the filter currently stored away with SAVE

                                   U - unit-normalise the response (window shape and missing
                                       data will modify the response).

 CORR [ X ][ W ][ C ]                - Report the central cross-covariance coefficient in the protocol
                                       along with the RMS of X and W. Y is synonymous with W.
Options X W (Y) and C request that the DC-level in
                                       the computation of the parameters is to be kept in.

 CCV [L=#l] [M=#m] [D=#d]              Compute a running cross-covariance between X and Y
                                       Use a file processing command with option cdir = K*
                                       to read and store Y.
                                       Output is normalized
with the number of valid data pairs.
                                       The two-sided domain length is #l and the center
                                       is shifted by #m. Y is shifted with the amount #d.
                                       The t0 and dt is adjusted to the new conditions.
                                       Defaults: l=101, m=1, d=0.

 STRUCTURE[FUNCTION] [L=#l] [N=#q] [A] [ time range ]        
                                       Return the structure function up to lag
2#l-1               [10]
                                  #q - The norm to be used (floating point)                      [2.0]
                                   A - Analyse, print the result in the protocol:

                                       [log(s(k))-log(s(1)]/(k log 2), k=2..l
                                       s(k) =

                                       Returns the series s along with dt=1

 AUTOCOR[RELATION] #k [+W] [time-range]
                                       Compute auto-correlation up to lag k
Uses W for intermediate storage of result!
 AUTOCOV[ARIANCE#k [+W] [time-range]
                                       Compute auto-covariance up to lag k
Specify an absurdly big lag to compute all terms.
                                  +W - Keep result in array W 

 MEMSP #n    [R:s] [L,][U[=#h,]][Q][/N][dB]
                                       After producing or reading a PEF, compute the
                                       Maximum Entropy spectrum, domain length n.
                                       Options see PERIODOGRAM below.

 PERIODOGRAM [R:s] [L,][U[=#h,]][Q][/N][dB][,M<#m] [ time range ]
                                       Return the periodogram. If called from tslist, the
                                       output series should be printed using -NI8 -n1-1
                                       Before the periodogram is computed, the time series
                                       is repaired (linear interpolation of missing data)
                                       and the DC-level is subtracted.
                                 R:s - Report Nyquist parameters in units of:
                                 R:s - Hz and seconds
                                 R:h - cyc/h and hours
                                 R:d - cyc/day and days
                                       Returns the spectrum along with dt=1
                                       Prints the Nyquist-frequency and the spectral bin
                                       interval on the protocol.

                                       The following options must be coded without
                                       intervening blankspace:
                                   Q - return the power                                   [Amplitude]
                                   L - return the logarithm (base 10) of the spectrum.       [Linear]
                                  /N - scale power by the number of samples                   [unity]
                                       (amplitude: the square root)                                         
                                  dB - return decibels         
                              U[=#h,] - Compensate the effect of a filter (1,-#h)  [don't, else h=-1]

                                M<#m - Repair option: Maximum #m repairs, linear interpol-
                                       ation is used.                                     [all but 2]

{FD=#fd | PD=#pd} [{+USB|+LSB|+CC}] [M<#m]
               {FS=#fs | [F]={Hz|mHz|cph|cpd} | [P]={s|h|d]} }
[+BL] [+D] [ time range ]
                                       Amplitude demodulation with carrier frequency fd
                                       (alt.: period pd) See example 17.

                   +USB | +LSB | +CC - Upper or lower side-band or complex conjugate [both & straight]

                     FD=#fd | PD=#pd - Demodulation frequency or period, respectively
FS=#fs - Scaling frequency for fd, Nyquist in units if cyc/h
                              [F]=u  - A standard unit for fd instead of a specific fs.

                              [P]=u  - A standard unit for pd instead of a specific fs.

                                       If no unit nor fs is specified, fd is interpreted as
                                       a fraction of the Nyquist frequency of the actual
                                       sampling interval

                                M<#m - repair max. m missing samples                          [don't]
                                  +D - preserve DC-level                                     [remove]
                                 +BL - band-limit with a recently stored FILTER

  ZSPF type:order [/F][,M<#m] [FC=#fc] [[F]={Hz|mHz|cph|cpd}] [+BL] /
[+D] [+H[-B] [-F|+J] [+STEP|+DIFF] [+P] [-P] [SU=#u] [DU=#u] /
                [context-options] [ time range ]
                                       Complex filter applied in the spectrum domain.
                                       Several models can be applied (polynomials in z, pole-zero,
                                       Bessel, powers of frequency)

                                       Consider to pad the time domain using the CPAD command!

                                       If a time-domain FILTER has been SAVEd, it can be
                                       applied here on the resulting spectrum for band-limiting,
                                       primarily to avoid near-Nyquist jitter.

                                       Order of operations:
                                       First scaling (P&Z only), then /F inversion,
                                       then +STEP or +DIFF, then band limitation

                                  /F - divide by the filter                                [multiply]
                                M<#m - repair max. m missing samples                          [don't]
                                  +D - preserve DC-level                                     [remove]
                                  +H - do a Hilbert transform on the input spectrum           [don't]
                                       Simplest application: ZSPF +H
                                 +BL - band-limit with a recently stored FILTER
                                  -B - do not do backtransform                                   [do]
                            -F or +J - do not do forward transform, use unity spectrum instead   [do]
                        +STEP  +DIFF - divide, respectively multiply, spectrum with iw        [don't]
                               SU=#u - output spectrum to file on unit number u
                               DU=#u - output time series before and after ZSPF action to
                                       BIN files and
ZSPF-x.out on unit number u    [don't]
                                  +P - force print of spectrum        [no print unless ECHO or DEBUG]
-P - suspends tsf_edit's debug mode (but not +P)            [DEBUG]
                             SPS=#fs - "samples-per-second", a scaling  frequency factor      [1.0d0]
                            [F]=unit - In the output list (+P and U=#u) the frequency that      [cph]
                                       is displayed is i/N * tscale * fs, i=0..
                                       where tscale is derived from the [F]=expression and
                                       fs eventually modified with the SPS option
                                       (N is at the Nyquist).

                          type:order - BESSEL:#n   where n is the filter order (e.g. 8)
FC=#fc - corner frequency
                                       If no unit is given, fc is assumed in units of the         [?]
                                       sampling frequency.
                                  +S - Show coefficients from within Bessel-filter subroutine
                                  +T - Trace filter action (frequency and response) from
                                       within Bessel-filter subroutine
                              Example: ZSPF BESSEL:8 U=44 FC=0.125 Hz +D +S
                                       for the GWR anti-aliasing filter.
Note that "Hz" is recognized by tsf_edit and used to
                                       set the corner frequency's absolute value.
                                       In the call to spectral_zfilter (and cbessel_f therein)
                                       the value fc in Hz becomes
                                       3600*fc*dt where dt is in hours. [F]=Hz would do
                                       something else.

                          type:order - POWER:#p   where p is the (floating-point) exponent
FC=#fc - scaling frequency (under +DISC only)
                                       With +DISC the spectral function is (2π i f/fc)p
where 0 <= f <= 0.5
                                       Thus, if e.g. 3 s is your sampling interval and the
                                       result of a time-integration (p=-1) is to be computed in
                                       units of seconds, the corner frequency is three times
                                       the sampling frequency, so FC=3.0
                                       Without +DISC, FC is not used.
                                       Example 14, not discrete.

                          type:order - RPOLY:#n   where n is the polynomial order, like ZPOLY,
                                       with real-valued coefficients, however.

                          type:order - ZPOLY:#n   where n is the polynomial order
                                       Filter F(z) =
Π(1+p(k)z)k, k=1..n,
where z = exp(2πi j/N), j=0...N/2

                                       Filters can be input in a number of ways:
                   FC={p(1),..,p(n)} - from the command line
                               FC=P  - from the calling parameters (not applicable here)
                               UC=#u - from file unit #u. From the next lines of the tse-file
                                       if #u=-1
  [defaults p(k) = -1]
                                       This is an integrating filter: ZSPF ZPOLY:1 /F

type:poles,zeros - P&Z:#p,#z   with p poles and z zeroes
                               +DISC - discrete version; poles and zeroes in variable exp(iω)   [
                                       Input poles and zeros, format

                                       ZP(#i) ( #zr, #zi )
                                       ZZ(#i) ( #zr, #zi )

                               UC=#u - from file unit u. From the next lines of the tse-file
                                       if u=-1
                                +Z=0 - set all coefficients to zero before reading in new ones.
                                 +UG - Unity gain but phase shift.

                               SG=#g - scale spectrum with factor g
                                       Example for two seismometers in

 FRACTAL SPF [M<#m] P=#p1[,..,#p10]
[ time range ]

                                       Spectral filter. Return the data after a Fourier
                                       transform forward-backward pair with a spectral
                                       operation inbetween. This is under construction.
                                       Functional at the moment: multiply with f^#p1
(using function POWERLAW, included in tsfedit.f)

 FOG2U SPF [M<#m] P=#p2[,..,#p10] [ time range ]

                                       Spectral filter. Return the data after a Fourier
                                       transform forward-backward pair with the spectral
                                       operation inbetween that converts acceleration into
                                       displacement; suitable for slow processes like
                                       Free Oscillations. p2 = Love number combination
                                       p2 = ( 2hn - (n+1)k
n )/hn
                                       p3 = g/a  ~ ( 9.81 / 6.371e6 )
                                       A tricky issue and a dirty method!
                                       We would actually need the frequency spectrum
                                       of the Love numbers and the associations of degree
                                       and frequency. For high frequencies the asymptotic
                                       values would suffice.

                                       (p1 is used for the sampling interval.)

                                       ...SPF: Developer, extend functionality by either branching
                                       in POWERLAW or FO_G2U (sas/p/spectralf.f)

                                       or create new if-blocks and keywords in TSF_EDIT!

 MOVRMS #m #w [M=meth]                 Compute and return the running RMS of the series
                                       where #w is the off-centre half-width of the window
                                       #m the shift (subsample ratio).
                                       In every window the DC-level is subtracted before
                                       the RMS-calculation.
                                meth - Method:
                      { DC | ME }[W] - derive the DC-level by arithmetic mean (DC) or
                                       by median (ME)                                         [DC]
                                   W - to apply a raised-cosine window in the

 MOVAVE #m #w  [M=Q] [time-range]    - Compute a time series of running averages of
                                       where #m is a decimation factor and #w the full(!)
                                       width of the averaging window. If values are
                                       missing in, the output samples will be located
                                       in the "centre of gravity" of the valid samples in
                                       each averaging interval. This feature differs from
                                       The work array W will be used.
                                 M=Q - Instead of the mean, compute the RMS (not the RMS-dev!).
                                       Sliding averages of a time series together with its
                                       standard deviations can thus be computed.

 MODAVE[RAGE] #n                       Computes a series of length n, where each sample
(i) = SUM {j /. MOD(j,n)==i : }
n - cycle length

 DECI2W #n #s [time-range]             Decimation (undersampling) to work array
                                   n - decimation factor
                                   s - initial shift    

DECI[MATE] #n #s                      Decimation (undersampling). The offset is applied
                                       before. Thus, x_out(i)=x_in(i*n+s), i=0,1,...
                                   n - decimation factor
                                   s - initial shift

 DECI[MATE] #n #k S=#m                 ...with moving-average filter of length m (odd)

 DECI[MATE] #n #k G=#f,#m              ...with Gaussian-bell filter of length 2m+1.
                                       Filter's -3dB point is at
                                       f = f(-3dB)/f(Nyquist)

 DECIMIN | DECIMAX | DECIMED           Three kinds of decimations:
 DECI... #n #s [O=options]             Decimate by taking the max or min or median value
of the data segments.
- decimation length
                                   s - initial shift
                            option P - diagnostic print

 FOURIERINT[ERPOLATION] R=#r [O=opt]   (using ~/sas/p/spectralf.f ENTRY Fourier_Interp)
                                       Fourier interpolation with optional taper
                                   r - rate factor > 1
                                 opt - string #lt,#s[,+P]
                                  lt - length of the taper, number of bins.
                                   s - width of Gaussian bell, argument (i/(lt*s))2, 0<=i<=lt
                                  +P - print tapering action on protocol.

 FFTINT[ERPOLATE] I=#r  [ Time-range ] (using ~/sas/p/fftintps.f)
                                       Oversample with Fast-Fourier. If a time range is specified,
                                       the resulting series will be shifted and t0 (the origin)
                                       be adjusted. The main purpose is to interpolate through gaps.
                                   r - integer rate r > 1

 RESAMP[LE] #o #t #dt [T!] [sec]       Resample with polynomial interpolation.
                                   o - polynomial order,
                                   t - origin time
                                  dt - sampling interval

                                       origin time is with respect to the previously first sample,
with respect to epoch (T!)
.......................................t and dt are in hours or in seconds (sec).

 RESAMPLE TT                           The following records constitute a time-tie, terminated with
                                       END S
  RESAMPLE U=#u                         Reads a time-tie from file unit u (`END S[AMPLE can be omitted)
                                       Records are time-range data
                                       [ From YYYY MM DD hh mm ss ff ] [ To
YYYY MM DD hh mm ss ff ]
                                       [ At
YYYY MM DD hh mm ss ff ]
                                       `At´ can be omitted.

 REGULA[RIZE] DT=#dt[S] [T0=#t0[S]] [TA=#ta[S]] [F=#k] [ +D ] [ +R ] [ time-range ]
                                       Interpolation of a sparse and irregular series to a lower rate
                                       regular series. Parabolic interpolation.

                                       Irregular sampling in this context means that a series is resolved
                                       on a regular basis but with much higher sampling rate, whereby a
                                       varying number of missing samples are located between valid samples.

                                  +D - move starting point t0 such that the first sample is a valid sample.
                                  +R - round output time origin with respect to dt.
                                   S - with time parameters: the value is in seconds
                                  dt - resampling interval.
                                  t0 - start time from epoch. Default is as early as possible.
                                  ta - add ta to start time of input series.
                                   k - Maximum number of missing samples allowed between the suspension points
                                       of the polynomial, else a missing-sample will be injected.

                            Example:   REGULAR T0=10S DT=10S F=600 +D From #0 To #86420
                                       for Hannover ZLS Burris gravimeter data, oversampled at 1s.

  SAMPLE [[U=#u] [->#j] At] #Y #M #D [#h [#m [#s [#f]]]]
                                       Sample by linear interpolation.
                                       OBS! The first TRMARG parameter applies. Interpolation is thus
                                       carried out as if the series was translated in time.
                                  #u - output to a file that has been opened with unit number u.
                                  #j - save value in a storage (retrieve with %j; for details, c.f. STATISTICS)
                                  At - this keyword must be given if an option is specified.

 STATISTICS [+S] [time-range]          Save some key characterizing values in a storage:
                                       MIN, MAX, AVErage, RANge, RMS-dev, N.Valid Samples
                                       Median is saved when MEDIAN command is issued.
                                       Up to 10 samples of the series can be stored with the
                                       SAMPLE command with option ->#j   (e.g. SAMPLE ->1 2014 01 28 22 37 00 000)
                                       Retrieve in any command line using, respectively,
                                       %MIN %MAX %AVE %RAN %RMS %NVS %MED %j (1 <= j <= 10)
                                       To invert the sign and/or taking the reciprocal value, prepend
                                       the `%´-symbol with `-´ , `/´ or `/-´  e.g. /-%1  or  -%MED )
                                       See Example (18)

 DT #t                                 (Re-)defines sampling interval t

 ; text                                Comment; text is completely ignored

 ;; text                               Comment; text will be printed to protocol

..............A D D I N G   o r   S U B S T I T U T E - M I S S I N G
...        - 

     value                       A numeric value or the symbolic value ` % ´ (with surrounding blanks!),
                                 applies to ADD SUB MUL POKE DEL RJC commands below
                                 Maybe not implemented consistently yet.

                                 Some commands set the %-value upon the [/]->S[-] option:
                                 MEDIAN  SUMTS REPDC   RMS..

Achilles heel: We might have to change strategy with the % symbol: String replacement of % with scale value. Don't expect too much, it's still not fully tested.

 MUL  #v [ time-range ]           multiply. Try 1/value to divide, maybe it works.

 SLO  #v [u] [ Iref=n | Tref yyyy mm dd hh mm ss fff ] [ time-range ]      
                                  add a slope.
v - if no measurement unit is given, v is rate/sampling interval
                              u - DT or [/h] : v is in units of 1/h.
                                  Other units: [/s] [/d]
                                  /N: the ramp goes from 0 to +v over the length of the series
                                  /TR: the ramp goes from 0 to +v over the current time range
                                  The zero-point can be changed:
                        Iref=#n - by time-series index
                        Iref=C  - at the center of the time-range (zero-mean)
                           Tref - at the specified point of time                      

 EXP #a #p u [ Iref=n | Tref yyyy mm dd hh mm ss fff ] [ time-range ]
                                  add an exponential  a exp(j p dt), p in [1/h] real-valued.
                     Iref= Tref - see SLO
u - the units in which p is given (default: [1/h]):
                                  [1/s] [1/d] [1/y] [s] [h] [d] [y]

 ADD  {#v | =DC | -DC } [RSQ] [ time-range ]   add value v to valid samples.
                            RSQ - x <- sqrt(x2 + v2) if v>0
                                  x <- sqrt(max(
x2 - v2,0)) else.
                                  ADDRSQ is possible as a command.
{#v | =DC }[ time-range ]   substitute v for missing samples
{#v | =DC }[ time-range ]   substitute v anywhere.

                        =DC -DC - Use the result from the most recent REPDC command
                                  (it determined the DC-level, optionally within a time-range for v.
-DC - the negative DC value (with ADD only).

 SIN #a #f #p [u] [F'=d] [FM=s] [AM=s] [ time-range
                                  add a sine wave with amplitude a, frequency f [cyc/h],
                                  and phase p [deg]
 COS #a #f #p [u] [F'=d] [FM=s] [AM=s[ time-range
                                  add a cosine wave...
                              u - measurement unit.
                                  Specify [cpd]
[cph] [Hz] [mHz] for frequency
                                  or [d] [h] [s] for period,
respectively. Default: cph.
                                  For cycles per time step, specify [DT]
For time steps per cycle, specify [/DT]

                           F'=d - Linear frequency ramp, f(t) = f (1 + d t)
                           FM=s - Phase modulation with the signal in W. Additional p = s pi/180 wj
                           AM=s - Amp. modulation with the signal in W. amp = (a + s wj)

#a #w At time

 PULSE #a -1 time-range         - add a zero-mean pulse of amplitude a and off-centre half-width w.
                             -1 - Full width is given by time-range.
                          a < 0 - Amplitude is computed by RMS-normalisation multiplied by |a|.
                          w = 0 - The pulse is a single spike of amplitude |a|.

 MOD #a #nf #ip #iduty #offset [ time-range ]
add a rectangular wave
offset+scale*{(MOD(i+ip,nf) <= iduty?) +1 else -1}

 TRI #a #nf #ip #iduty #offset [ time-range ]
                                  add a triangular wave, rising-flank length = iduty.
                                  Based on the same MOD-equation as above, but signal
                                  varies between 0 and 1 when offset=0 and a=1

 RDC [ time-range ]               remove a DC-level inside the window

 ADDNOISE #a,#s type [ P=#p ] [GMP=#g] [ time-range ]
                              a - amplitude
                              s - random seed
                              g - Gauss-Markov parameter: y <- g y + ran(s), x += a y
                 type   FRACTAL - with power-law parameter p
type   GAUSS   -
type   UNIFORM -
                                  This function uses workspace W . Consider dumping it
                                  with  DUMPW #u   for a check.

 RSM {I|0} [ time-range ]         call repair_single_misses (isolated dropouts)
                                  Mandatory option 'I' or '0' (zero).

 REP [L] [ time-range ]

 REPAIR [ opt ] [+P] [ time-range ]    Options: V=#v | L | L<=#n | H#n    

..............REPAIR              repairs breaks by substituting DC level.
..............REPAIR V=#v         repairs breaks by substituting numerical value v.
..............REPAIR L            repairs breaks by linear interpolation.
..REPAIR H[#n]        "sample-and-hold", repeat last valid value #n times at maximum
                                  into a gap (default #n: indefinite)
..............REPAIR L<=2         repairs breaks of length 1 or 2 by linear interpolation
..............REPAIR              finds all breaks (REP ALL)
..............(REP L ALL is perhaps too difficult for the program.)
                             +P - additional option to print out the places to repair.

                                  If the time-range for delete includes the first
                                  sample, the origin will be moved physically.
                                  This fixed behavior might not be the best choice - maybe on option.

[ time-range ]  - Fill gaps in X by interpolation of W, polynomial order m
                                  Default m=2

 WGAPF #a,#c [SQR] [ time-range ] - Fill gaps with a half-sine-square of amplitude a2,
                                  add constant c2 everywhere. This function is used to
                                  generate weights for interpolated tides in gaps.
                                  The constant a can be computed with
                                    urtipgt @tintpolerr.ins; fgrep 3600 o/tintpolerr.dat # ( -> 22.122 )
                            SQR - Return the square root

 HARMONICRESTORE #f1 #f2 #fny #thrl #thru [ time-range ]
                                  Interpolation of missing samples or samples out of
                                  range thrl < thru by Least-square fit of a
                                  band-limited spectrum, f1 < f < f2. The frequency
                                  scale is given by fny, the Nyquist frequency.
                                  Mediocre results, close to useless so far.
                                  No defaults for frequency range, -1040 to 10
40 for
                                  signal. (Subroutine sas/p/sigrestore.f)

 UNCLIP [O=opt] [T=#tol] [I=#itr] [S=#p] [ time-range ]
                                  Restoration of a clipped signal. (Subroutine sas/p/unclip.f)
                                  `Unclip´ bridges a gap by inserting
                                  a linear ramp 
                                  plus a half-period sine
                                  that starts and ends in the points of fastest signal change
                                  immediately before and after the gap. The locations of these
                                  points determine the frequency. A linear ramp is added.
                                  Only the missing values are replaced; however OUTTS W can be
                                  used to save the predicted series (which may have gaps where
                                  the input signal has long stretches of continuity.
                                  Works ~o.k.

             opt [P+|P-][D+|D-] - Print or debug excessive printing.
                           #tol - tolerance for solution misfit, default  1.d-6
                           #itr - maximum number of iterations, default 1000
                       #p {4|5} - Number of parameters to solve.
                               4: all the linear dependencies
                                  (1+2: linear-varying sinus amplitude, 3+4: linear ramp)
                               5: also the frequency parameter (not recommended,
                                  however the default is 5 (!)                 .

 DEL  [ > #v ]  [ time-range ]  - delete (criterion abs(x(i)) > v )

 DEL  [ < #v ]  [ time-range ]  - delete (criterion x(i) < v )

 DEL  [ >> #v ] [ time-range ]  - delete (criterion x(i) > v )
 DEL  [ => #v ] [ time-range ]  - synonymous

                                  If no value is given, scope for DEL is due to time-range.

 DEL == #v +-#d [ time-range ]  - delete if |x(i)-v| < d . OBS! no blank between `+-´ and number.
                                  (A complementary command /= is still missing.)

 RJC  directive [ time-range ]    Copy the missing-value symbols in X and W
                                  according to the following directives
                          X=W   - if x(i) or w(i) is missing, mark both as missing
                          X<W   - Missing values in W are copied to X
                          X>W   - ... or vice-versa

 RJC  [ > #v ] [ -DC=#v ][ time-range ]
                                - delete (criterion abs(x(i) - DClevel) > v )

                                  The other options known from DEL exist too,
                                  but might be less useful.

 REJECTSHS #n  [ time-range ]   - Reject short segments, amount of valid samples in each
                                  segment will be > n.

 INJ[ECT] { #v | MRS | ITEM } [ N=#n ] [ At time ]
                                  Shifts the series towards the end by n samples, extending its
                                  length, starting at the specified time. The shift may be
                                  negative. If positive, the value v is filled into the gap.
                                  Symbolic values are MRS (missing) or ITEM (the value at the
                                  emerging gap). Default n=1

 TRUNC [ time-range ]           - truncate beyond that time. Pefer to specify
                                  At #YYYY #MM #DD #HH #MM #SS #FFF
 TRUNC N=#n                     - truncate to retain #n samples.
 TRUNC #n                       - should also work like N=#n
EXTEND[W] #n #v
XTEND[W] #n [V={DC|#v}]
 EXTEND[W] by #m
[V={DC|#v}] To #YYYY #MM #DD #HH [#MM [#SS [#FFF]]]
                                  Extend the series X resp. W to obtain n samples or add m or to
                                  the specified date. Fill with value v or DC level or
                                  mark as missing, the latter by default. Issue REPDC first, in
order to apply the actually valid DC-level.

 TRIM time-range                - retain indicated section. Time origin will be adjusted.

 MIXADJUST [O:[+R[+T]][+S]]     - Use W (low rate) to adjust X by adding a linear curve at and
                                  between the support points. Useful to combine Atmacs data (W)
                                  with local pressure loading (X).
                             +R - Replace W with the difference W - X
                           +R+T - ... and trim to fit time scope of X
                             +S - If a supporting sample of X
is missing, delete the samples
                                  of the whole interval; else try and find a sample in the
                                  neighbourhood. See TD/Atmacs/HOW.TO                           

                     Zero-padding for spectral filtering, periodogram etc.
                                  PAD -> spectral filter or periodogram or ... -> UNPAD
                                  is recommended in order to avoid effects of circular correlation

 PAD [{N|#n}*{DC|#v}] [+I]        Extend the series with amount n injecting value v. Default N*0.0d0
 PAD =
{DC|#v}                     If default for n is o.k.
                              N - to double the length of the series
                             DC - to substitute the series' DC-level

 CPAD [{N|#n}*{DC|#v}] [O=#l] [+I] time-range  
                                - For use with ZSPF: Copy a padded segment into a work area.
                                  Add l additional samples to avoid contamination from an eventual
                                  jump at the the start of the pad.
                                  The work area is allocated at the end of the signal array X,
so that the calling program must take this extra space into account.
                                  tslist, for instance, should not be called with more than one
                                  data column in order to prevent invasion.

[C]PAD +I                      - For time-integrating ZSPF applications, refresh the pad space
                                  (equivalent to a zero constant of integration).
NPAD [+K                      Undo padding. +K will keep the length of the series so that
                                  the result of circular correlation / overlap management can be
Undo CPAD padding.

 RESYNCUNPAD                      After CPAD with nonzero overlap, resets the time of the first
                                  sample due to the overlap.

 RESYNCUNPAD +A                   In a block
                                  DO time-range // CPAD // ... // DONE
                                  with a time-range not starting at zero we
                                  store this offset and apply it now to adjust the time of the
                                  first sample. See example (12).
RESYNCUNPAD [#n]                 However, if CPAD is called repeatedly in a loop, RESYNCUNPAD +A
                                  becomes uncertain. Specify the sample number n manually.
                                  Perhaps there's demand for obtaining n from an expression
RESYNCUNPAD At time-range

             T S F - E D I T  B L O C K  S T A R T  A N D  E N D:
 TSF EDIT codeword                target word = TSF EDIT (capitals!)
 END                              instruction word = END (capitals!)

Some internal variables
At present (Jan. 2013) the procedure can handle some integer variables that can be written into the command line string.

The variables must be introduced in the Fortran source code, and set appropriately as parts of the different functions.

Names:  nhashv(k)
Values: ihashv(k)
Number of parameters enabled: khashv; 1 <= k <= khashv
At present, khashv=6,
TRMARG and TRMARG2 are set with the TRMARG command
SHIFTS is set with the SHIFT command.
JUL and YMD are set with the Time-Range From or At or Around:#p
INDEX is set with the SETINDEX and ADDINDEX commands
Many extensions can be imagined. 

For access to these numbers, code ###name#, e.g. ###SHIFTS#
There is also an internal counter for every ### that appears in a command string.
The counter is stepped up at encounter before the value replaces ###.
To preserve the value of the ### counter, code #.# instead.
To reset the counter, code #-# #0# or #1# to start with -1 0 or 1, respectively
To reset it to a different number, use SET ### n  (and don't forget that the value of ###
will be incremented by 1 BEFORE the next encounter!). Use SET #-# n  to set the counter to n-1.

The strings ###[name#] (and variants) will be replaced in the command string
with a string that is produced with    WRITE (RPLCMSTR, FMT) NUM
The numbers are written with at least three decimals (leading zeros as required) and a plus sign in front of non-negative numbers if there is space. (Needs verification, needs comments as to FMT, needs method to change FMT)

Thus, you can e.g. generate file names or labels like in

PEF U=41 F=${PFFE} L=${PFLE} O=R
LOOP ${SHIFTL} <102>
OPEN 51 > prdata/v${VER}/p-f${PFLE}.mc
DECI2W 500 0 From 2011 06 15 00 20 00 000 To 2011 06 15 12 09 05 000
DUMPW 51 L=SPAG_______E###SHIFTS#]
(This example is from Seismo/gcf/AG/predict.tse)

Try this:

LOOP 10 <102> #-#
ADD 1.0 From YMD+### 0 0 0 0 To YMD+#.# 0 0 30 0

Purpose: to add a step lasting 30 seconds at every midnight the first 10 days.
The #-# leaves the counter with a value -1
The first time that ### is encountered, the value is stepped up by 1 before it is evaluated, thus giving 0.


Example (1)

 The file $TAPDIR/KIRU.tse contains TSF EDIT commands; a block END mark
 is not needed. Notice that we can enter environment parameters.

Example (2)
Gravimeter data editing
22,'BRX 9999]',-99999.d0,'F6.0',0,0,'G',1.0d0
22,'BRX 9999]',-99999.d0,'F6.0',0,0,'G',1.0d0
22,'BRX 9999]',-99999.d0,'F6.0',0,0,'G',1.0d0
22,'BRX 9999]',-99999.d0,'F6.0',0,0,'G',1.0d0
-1.4462,535.0,'F 17 999]'

 with the following files (/geo/hgs/TGrav:)

   C instructions for sasm01: Open files
21 B goadc84.dat    observed ADC data
22 B gog84m.dat     merge manual data, 4 sections
23 B gog84j.dat     jumps1
24 B gog84jg.dat    jumps2 after scaling
31 < gog84.ts       output

Example (3)
First part of this example uses a Wiener filter. The filter has been written by  sasm06.
It is used in  urtap:
FILTER U=42 F=../sas/r/${section}_wtr_ap.inf O=AB >-16    16 'A'>

where  $section  is defined in the environment
setenv section ka1131b_1_A
and the filter file contains a section

   -16    16 'A'
   .00000000D+00 -4.22510640D-05  4.25726030D-05 -6.99387654D-04  1.26891712D-03
  7.17114170D-04  2.42920550D-03  3.40108467D-03  8.16553083D-03  2.70715305D-03
  1.47260492D-02  1.74710095D-02  2.97779297D-02  3.79714458D-02  5.97417858D-02
  5.49289255D-02  3.90932678D-02  4.27229613D-02  2.31071835D-02  5.51845096D-03
  4.89410900D-03  9.15132112D-03  3.30929002D-03  3.77412568D-03  2.28769855D-04
 -1.19480978D-03 -6.20019744D-03 -1.40937889D-03 -2.23349367D-03 -5.52621497D-04
 -2.75012143D-04 -4.77684169D-05   .00000000D+00 

This section is located by qloc_in_file. using  -16    16 'A'  as the target. The option   O=AB implies:  maximum one rewind (by default), scan anywhere on the lines, backspace after successful locate. You might also consider  for trace printing and  for no rewind.
The zeroes at the end points of the filter are due to the HANN window. Some windows don't taper to zero.

Other filters:
Window design, can be produced with   sasm02
At the prompter below the spectrum plot, enter U (capital U). The filter record will be written to the terminal.
From there you can cut and paste into a filter bank file like above.

If your filter is of a running average type (smoothing, low-pass filter), you can specify option +N in the
FILTER command line. This will avoid the long breaks at the expense of more noise near the breaks.
The output will also be normalized (sum of filter coefficients normalized to unity)

Example (4)
Computing common mode after EOF analysis.
This file is produced automatically by ~tap/Postp/wr-plot
SCALE .178856
23,'BIN]',-99999.0,'L:RAD_VAL]_M>0',0,0,'N', -.065634  -DC
23,'BIN]',-99999.0,'L:RAD_VAL]_M>0',0,0,'N', .346502  -DC
23,'BIN]',-99999.0,'L:RAD_VAL]_M>0',0,0,'N', .198229  -DC
23,'BIN]',-99999.0,'L:RAD_VAL]_M>0',0,0,'N', .135893  -DC
23,'BIN]',-99999.0,'L:RAD_VAL]_M>0',0,0,'N', -.016194  -DC
23,'BIN]',-99999.0,'L:RAD_VAL]_M>0',0,0,'N', .113417  -DC
23,'BIN]',-99999.0,'L:RAD_VAL]_M>0',0,0,'N', .068454  -DC
23,'BIN]',-99999.0,'L:RAD_VAL]_M>0',0,0,'N', .047797  -DC
23,'BIN]',-99999.0,'L:RAD_VAL]_M>0',0,0,'N', .081275  -DC
23,'BIN]',-99999.0,'L:RAD_VAL]_M>0',0,0,'N', .008231  -DC

Example (5)
Signal conditioning (differencing), cleaning and information

FILTER  0 1 ' '
1. -1.
RMSD To 2000 02 24 0 0 0 0
RMSD From 2000 02 24 0 0 0 0
DEL > $cond_limit ALL
RMSD To 2000 02 24 0 0 0 0
RMSD From 2000 02 24 0 0 0 0
;      A command like
;        setenv cond_cont "CONT 24 R tap/r/${section}_wapsl.tse"
;      will include the specified file for additional editing

which employs the environment parameter $cond_limit. Example:
setenv cond_limit 0.3
Reports of RMS and number of missing samples are requested before and after the outlier cleaning.

The following example considers a file tap/r/${section}_wapsl.tse

DEL From 1993 09 17 14 48 00 00 To 1993 09 17 16 48 00 00
DEL From 1993 10 12 20 48 00 00 To 1993 10 12 20 48 00 00
DEL From 1993 10 12 18 48 00 00 To 1993 10 12 19 12 00 00
DEL From 1993 10 13 00 48 00 00 To 1993 10 13 00 48 00 00
DEL From 1993 10 12 16 48 00 00 To 1993 10 12 16 48 00 00
DEL From 1993 10 07 20 48 00 00 To 1993 10 07 20 48 00 00

Example (6)
Take a segment of a signal, construct a prediction-error filter, and apply it to the whole signal
Save the filter for later use.

BURG L=20 O=O  From 2009 10 02 12 00 00 00 To 2009 10 02 13 00 00 00
41 < o/rt.pef

Example (7)
; samples to a sequence of files, each time shifting the
; sampling instance one time step ahead, 10 times
LOOP 10 <101>
OPEN 51 B sampled###.asc
SAMPLE U=51 2011 06 15 01 40 01 000

A similar case, now incrementing SHIFT in a loop and output to an MC-file with label
; samples to a sequence of files, each time shifting the
; sampling instance one time step ahead, 201 times.
SHIFT S=-101
LOOP 201 <101>
OPEN 51 > svdata/
DECI2W 500 0 From 2011 06 15 00 20 00 000 To 2011 06 15 12 09 05 000
DUMPW 51 L=SV_____________${SVC}]

Example (8)
The spectrum of a finite impuls response filter, here a prediction filter

ADD 1.0 At #
PEF U=41 F=${PRF} L=${PRFL} O=R
PERIODOGRAM /N From #0 To #4095

  setenv PRF solnew/svE-ag-CYCZEN.prf
  setenv PRFL 77
  tslist _4300,0.0 -Efirspect.tse,S -N -n1-1/4096 -B2011,1,1 -w

(4300: we need a margin for the filter cutoff effects)

Example (9)
Undo the phase effects of the Guralp seismometer and compute acceleration

FILTER D:WD KAIS 240 2.0 0.0 0.005 1.0 0.0 SAVE
ZSPF P&Z:5,2 /F [F]=rad/s UC=-1 -P +UG +BL
ZP(1) (-3.7e-2 ,  3.7e-2)
ZP(2) (-3.7e-2 , -3.7e-2)
ZP(3) (-5.03e2 ,  0.0   )
ZP(4) (-1.01e3 ,  0.0   )
ZP(5) (-1.13e3 ,  0.0   )
ZZ(1) ( 0.0    ,  0.0   )
ZZ(2) ( 0.0    ,  0.0   )
FILTER 0 1 ' '
1.0 -1.0

Example (10)
Apply the GWR Bessel-8 filter

ZSPF BESSEL:8 FC=0.125 Hz -P

Example (11)
A series spliced together at interval 60000 samples contains small end effects at each joint.
Each interval is also biased with a DC-level. We adjust the levels by repeatedly calling
POLYREDJUMP. We attenuate the end effects by smoothing within a short interval around the joints.
Also note that sample numbers in time-range expressions don't need to be abutted to the # symbol.
Also note the termination criterion of the loop: The index variable becoming greater than the
length of the time series.

tslist d/tmp.ts -Ejump.tse,J -I -o d/tmpj.ts

LOOP 73 <101>
POLYREDJUMP ${POLYORD:[3]} M=50,50 O=P From # ###INDEX#
FILTERINPLACE 0 1 'S'   Around:-3,3 # ###INDEX#
0.5 0.25
ENDLOOP <101> (N> ###INDEX#)

Example (12)

The input has 100 S/s and is loooong. We have to start filtering at offset 20 samples
to obtain output on precise seconds. That's done in the DO command.
Because of the overlap option in CPAD we must call RSYNCUNPAD after the loop.
It's also the place to reset the start time given to the first DO command.

; resample a long time-series with low-pass
; filtering it first, using a spectral filter
FILTER D:WD KAIS 2880 2.0 0.0 0.008 1.0 0.0 SAVE
LOOP 78 <103>
DO From # ###INDEX# For #60000
CPAD O=2880
ENDLOOP <103> (N> ###INDEX#)

Example (13)

Shows the use of the OUTFS command with conversion from time series to filter.

; calculate the low-pass mixer filter down-decimated to 1 S/s
; the filter that is to be used with the inverse Bessel in
;  ~/TD/besself.tse,BSFM
;  tslist _60000 -qqq -rs0.01 -B2011,1,0 -Eresmp.tse,FDECI100 -w tmp/tmp.fsf
; To examine the spectrum, use
;  sasm02 -F fdeci100.fs
ADD ${PULSE:[1.0]} At #${PULSET:[0]}
FILTER D:WD KAIS 7200 2.0 0.0 0.0004 1.0 0.0 SAVE
OPEN 31 < fdeci100.fs
OUTFS U=31 X> +S From #0 To #72

Example (14) ~/TD/kamchatka.tse
Double time integration for displacements from accelerations

; acceleration to displacement on 3-s or 10-s data
; Scale must be set accordingly.
; 1. restore SCG's anti-aliasing filter
OPEN 31 < tmp/
CPAD O=1200 +I
ZSPF BESSEL:8 /F FC=0.0174 Hz +S -P DU=31
; 2. Integrate by Spectral Z-filter
FILTER D:WD KAIS 1200 2.0 0.002 0.5 0.5 0.5 SAVE
; This is for mm from nm/s^2:
SCALE -1.e-6

Example (15) ~/TD/rms.tse,EQ & ~/TD/weq.tse,W
Delete measurements affected by earthquakes, interpolate and provide weights for the interpolated ordinates.
Input =, working on 20s-data, output = 600s-data.

Find disturbances:
  eq2wdr 090701 1415 # uses rms.tse,EQ
Interpolate gravity and create weights: use script  ~/TD/weights4gaps , see SCG-HOW.TO


; Note: Reasonable defaults so you can run tslist directly
;       although this section is dedicated to ~/TD/eq2wdr
FILTER D:WD KAIS 40 2.1 0 0.05 1.0 0.0
FILTER 0 1 ' '
1.0 -1.0
OPEN 31 ${EQTSEO:[A]} ${EQTSE:[eq.tse]}
; You might like to use the margin option, e.g. MS=-600,600
WDR ${EQTHR:[20]} U=31 MS=${EQWDRM:[0,0]} From #0 To #4520
; 4520: we need overlap, too much doesn't hurt
; 70 minutes for 600s-data plus the 120s due to MOVRMS
; calc '(86400+3660)/20+10' -> 4513


; Purpose: generate weights inside gaps. Used by ~/TD/weights4gaps
; 1. Delete ordinates for the very day and the day after

CONT 31 B d/eq_${YMD}.tse
CONT 31 B d/eq_${YMDP}.tse
; 2. Insert function

WGAPF 22.122,${WGAPFAC:[0.0]} SQR
; 3. convolve with filter and decimate

FILTER D:WD KAIS 58 2.1 0 0.0333 1.0 0.0

Example (16)
Filter from a binary file .fs:  (Ttide/SCG/aprestap.tse)
GETFS U=11 < o/tgwr-grwr-1h.wwf] M=WIENER INF:    32] O=AB SAVE

PRINTF F=I4,1p,e12.4 U=41 B filter.coeff]
Retrieves and applies a Wiener filter, and lists it in ASCII on file filter.coeff

Example (17)
Demodulating an M2 tide residual:
Imagine you fit tides to a time series by least-squares, yet a power spectrum of the residual shows
remaining peaks at and around tide frequencies.
tslist owf/g100201-130515-1h-lapww-M6.ra.ts -E amdemod.tse,AMS -Un23240 -I -o tst.ts
with amdemod.tse:
FILTER D:WD KAIS 2880 2.0 0.0 0.008 1.0 0.0 SAVE
CPAD 2880*DC O=2880
AMDEMOD +D +BL FD=1.9322736 FS=12.0
Note the big loss of data because of the length of the low-pass filter.

Example (18)
Remove a bias using the median and scale with the RMS-deviation
Note: MEDIAN must be called separately while RMS is computed in the STATISTICS command.
The inverse and reciprocal values are replaced already in the call to subroutine replace_statvars (sas/p/tsstats.f);
reason: no blanks between operator and `%´

Calling mains:

(beware calls from subroutines not mentioned in the mains!)

Some subroutines at least: