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.

     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

* You can skip a set of instructions conditionally
setenv DOTHIS 0 alt. setenv DOTHIS 1
TSF EDIT target
IF ${






* Targets must not contain blanks.

* You can include command lines that work with the block, for instance examples for use:
TSF EDIT target
;. command line
;. ...
   Run - or examine - the command using script tsex . For usage, issue
   tsex -h

  There is also a method for issuing commands in shell scripts that e.g. may prepare
the environment or shell variables for use in a subsequent command.
* 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)  

    Setting environment variables that are significant for a command blocks:
cf the specific section on assisted setenv

 Some internal variables

      %Statistics variables:
      Computed with the STATISTICS, MEDIAN and SAMPLE commands and retrieved with %SYM symbols.
      Retrieval and formatting:
      Default formats exist but may not be apt to intended use. A format code (Fortran) can be appended with a leading colon.
   NOOP %AVE:f5.1
- if this yields overflow, procedure will fall back to the default.
   NOOP %1         - stored variable 1 is written with default 1p,e20.12
   NOOP %001:i10
- conversion of stored variable 1 (must be coded with three places) to integer with IDNINT.
   NOOP %001:I10   - conversion to integer with IDINT

      Formatting numerical values:

      A default format can be set with
      SET% FMT=fmt

      A couple of numerical values (under  REPDC WMEAN MEDIAN RMS RMSD COMBAVE SAMPLE
      (to be extended)) can be formatted with option
      on virtually any command line (including NOOP or the very command line that issues the printing. The format
      holds until recall/redefinition.
      Reset to default using
      2016-12-03: The current trouble is due to too much ad-hoc programming: each of these commands (REPDC ... ) produces its own output layout.
      And there are already many scripts that process these records. If we bring it into a standard from (which would make tsfedit.f
      a bit shorter), the scripts need revision.  Subroutine sas/p/tsfeditprt.f hasn't been added to libsas.a yet.

      It's not a hash structure, it's because we use the #-character.
      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

      A counter:
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))

A Simple Calculator
   New feature from 2014-08-11, not fully outgrown yet w.r.t. risk of clobber.
   Instead of coding numerical constants, you can use /home/hgs/bin/calc in-line to evaluate arithmetic expressions.
   The expression is evaluated after $environment parameters, #hash and %stat variables have been resolved.

NOOP  F`%RMS *sin(pi/4)`
`%RMS *sin(pi/4) :%10.4f`   
   Expressions evaluate to floating-point numbers if preceded by F` , integer if preceded by I` , (no format option with  `I),
   no other option to calc can be given. The format string must be compatible with both C and Fortran77 (besides the obvious
   rearrangement  %10.4f -> (f10.4) ). And it must be offset with a blankspace after %- and #-variables
   The subroutine that invokes the calc script is ~/sas/p/fcalcs.f .
   calc  is called with options  `-f "%20.12e" -F´ (default format, -F is hard-wired).
   If the expression contains negative values or symbols that evaluate to such, parentheses will prevent problems.

NOOP Median - Average = F`calc ( %MED ) - ( %AVE ):%10.4f`
   The result is written to a file /tmp/ftmp.123... (the extension is a random number; if you need to see it, search by creation date).

An RPN calculator
   Expressions in Reverse Polish Notation
        R`<expr>[ :<fmt>]`
  are evaluated with  subroutine revpolnot(expr) using ~/math/afor/p/rpn.f
  Default format is `*´ . For more, see manual page. Only just begun 2018-11-05.

    Comments are indicated by a ;-sign in the first position, or are preceded 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.
    The NOOP command lets the variable replacements take place and reports the currently defined time-range.

    More versatile comments can be produced with the block
    COMMENT [U=#u] <<target
    text text
    text including $ % and # variables
    The text lines are written to the protocol by default, or to a file opened on unit u. This way you can produce table heads
    with flexible content.

A mechanism to interact with shell scripts, e.g. to set variables
   If the tse-file contains comment lines coded like
   ;:TARGET: command
   then in a csh-script, the following line will execute the command:
   eval `awk '/:TARGET:/{sub(/.:.*:/,"",$0);print}' file.tse`
     Instructions are by default echoed to the stdout.
     This can be suppressed by coding a `/´ into the first column.

Retrieval of what was actually done - environment parameters, internal variables ...
    If you run the job with tsl instead of tslist, a log file is produced in ~/tslist-logs/
    with the symbolic variables resolved. In the tse-file you can add command
    NOOP with text that refers to the symbols you want to inquire
    A   grep <string> <logfile> | tail -1   of such a line will reveal the actual value. Example:


    In the log, the lines with resolved symbols are indicated with <TsfEdi>C> :
    <TsfEdi>C> NOOP setenv PDGLENGTH 36000 ;<-

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.

Under tslist, if you need data (e.g. a column) from the input file, use the open file on unit 11.
(command line option -keep-open ). Use -Llabel],Rwd here or REWIND 11 before the file processing command. 

Also remember that you can open more than one file (option-# filename  and -#L filename -Llabel ) .
Input unit is then incremented +1 (12, 13...). Not guaranteed bug-free though when labels are used.

New syntax, ASCII:

  @iun [+cdir] F=format[]] [T=target[]]] [-k#hms] [TZ=#tz] [MRS=#recmrs] [A=#value] [more options]

New syntax, BIN:

  @iun [+cdir] [L=label[],Rwd]] [A=#value] [more options]

+cdir   directive code must appear in position 5 of the line, else the default (N) is used!
  ASCII: Format can be specified without outer parenthesis. If format or target include blanks, right-delimit the code with a right-bracket.
  BIN-MC: If you need to specify a label with lower-case characters, use  F=l:label[],Rwd]


  target - BIN
  cdir   - N
  format - U
  label  - <none>
  hms    - 0
  tz     - 0
  recmrs - -99999.0
  value  - 1.0

Thus, the most common case of adding the unbiased content of a file to the time-series in memory is simply

  OPEN 21 < the-binary.ts
  @21 DC-

Old syntax:

 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-

The operands:

 iun     - file unit number; though a few symbols are accepted:
??      - will find an unused unit,
 ?i      - up to nine such units can be allocated, see util/afor/p/iunas.f
 O>      - under tslist or any main program that has called
           tsf_edit_fileunit (iun) for an open file, this unit number
           will be used.
           OBS! REWIND
and other file commands have not yet been coded up to also make use of the  O> symbol.

 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    - Directives, char*2
           cdir - meaning                             needs more data [YES] / instructions / options

...........A[!] - multiply with value and append                         NO
                  With !, replace existing samples
           A[!] -
multiply with value and append , adjust the series
                  to the same DC-level in the overlapping section        NO  / / +SP

           X[!] - multiply with value and append with cross-fading       NO  / / XF=length,early
OBS! 'A' with XF= option will switch to cross-fading. 
Specify length of fading interval and a time
offset for starting early w.r.t. Ndata(X)-length
(unit = index). 

           N*   - nothing, multiply with value and add data as is        NO
           S*   - apply a scaling factor and add              1 record:  scale factor sc0
           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) and 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.i/W.i
           C*   - add series cyclically. Length of series determines cycle length.
           *t   - override dt with the one from the data file

           Append calls subroutine append_ts. There are few (currently one) parameter setting
           options, command APPENDDEL ,
           Appending with or without cross-fading: see ~/sas/p/apptss.f

           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 * sc0 ) )
sc0 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.

Gap filling (G):
           OBS! Series must be in synch, x and w must have the same sampling interval.
             Gap parameters:     scale,bias,strategy

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

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

             The options XD or XH request read_fm_exactd('D') or ('H'), respectively.
             In the case of binary files, specify a dot.
and end are indexes of w.

             Option +ADJ requests adjustment of the gap fill's DC-level to the simultaneous part of x

             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
sc0, used with cdir = 'A*' or 'N*' or 'L*' or 'S*' or 'K*' or 'T*' or '+*'

Table:  More options 
(options must be coded with blanks as delimiters)
compatible with cdir
 +O | O+
'A*' - instead of A! , the replacement option can be issued this way.
   +SP | SP+
'A*' - "splice" i.e. adjust the appending series so that there is a zero
  DC-level difference in the overlap
. If SP=#w is specified, a narrower
  range, fraction w of the overlap interval, is used.

 -DC | DC-
- Remove the DC-level of the time-series W
 -DC | DC- 'L*' - the following operation takes place:

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

 -DC | DC- '+*' - 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

 +&& | &&+
'N*' 'K*'
- Use this data series to mask the one in memory (inject missing record
  This masking (still) uses a simplified scheme. The rates of the two
  series are *assumed* equal.

 +SQ | SQ+
'N*' 'L*'
- Add this data series to the one in memory by squares
 +RR | RR+
'N*' 'L*' - Remove this series from the one in memory by regression
- Fill up with zeros beyond file end deprecated
 *X | X*
 *W | W*
'N*' 'L*'
- multiply with
     X.out = sc2
··( [-DC-level] + value) (OBS! +value)

 X/ | /W
'N*' 'L*'
- divide by
     X.out = sc2
· [-DC-level] + value) (OBS! +value)
- divide by scale factor (/value)
  divide by scale factor (/sc0)
  divide by scale factor (/sc1)
- divide by sc0, multiply with sc1
- multiply with sc0, divide by sc1
- divide by sc0, divide by sc1
 +++ 'N*'

- Stack this data series onto X cut to the length of
  (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='+*'
 +CY | CY+ 
- cyclically: W is wrapped over X

 +++CY 'N*'
- stack cyclically
- trim-option with BIN input: aligns the input with the currently valid
  time origin t0

 !T=[X][u:]#v (any)
- move time origin by v, units of hours is default
      u=S - #v in seconds
      u=D - #v in days
  X - move time origin to the one of series X and add #v

 !T={+|-}ST (any)
- Summertime adjustment, e.g. stacking days synchronous with civil time.
- adjust time origin to the next integer v seconds (r=s)
  or minutes (
r=m) or  hours (r=h and default).
  With T, adjust with respect to the master time series, else adjust
  without offset.


Protocoll: Read carefully, eventually use the DEBUG command!
W) -
Works with other cdir's ?

Application example for !RT:
You have a bunch of ts-files that internally keep the same time interval, e.g. 10s, but the origins are shifted by a random number of seconds.
You would use tslist-app to concatenate them into one file, tolerating the offset biases. In order to align the bunch with the first file,
   setenv TSAPP_ROPT '\!RTs=10.0'
   tslist-app -I + AGSG_201502??b*.mc 

Application example for  tslist-app with cross-fading
Two files with a considerable overlap but suspected transients at the end.
Cross-fading will take place 24 minutes before the end of the first file, cross-over interval is 6 minutes long.

   setenv TSAPP_ROPT XF=360,1440
   tslist-app -o tmp/tmpxfade.ts -I + G2015.055/GCF.2015:055:02:00:00.3U93Z2-ZACC-1Hz.ts \

                     In-core processing:
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. 
  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
  COMMAND  Around:-30,30  2014 08 10 23 50 50 000
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.
 DEL To YMD+14 12 0 0 0 0

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

##k   is equivalent to   For ##(k-j+1)
  PERIODOGRAM /N From #740 For #4230
Sample numbers are in the range
Introducing a search margin (or a shift):
#m1 #m2     - m1 and m2 is in number of samples. The margins can be incremented:
 TRMARGSTEP [#ms]   -
adds ms to both m1 and m2 , default ms = 1
COMMAND Time-range

Example: apply an operation to the sample immediately before 2011 07 04 12   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]   
                                  B L O C K S

 DO Time-range                         In a DO - DONE block time-range expressions in lines below the
 DONE [+EP]                            DO will not be evaluated. There is no stack of ranges (yet)
                                       The purpose of this command is to avoid repetitive time-range
                                 +EP - Set a new epoch. Useful when an invasive command is carried out,
                                       like FFTINT which moves the result if the operation to sample #0
                                       (in general, commands that change the sampling interval).

 BEG[IN] RESCALE #f                    Within this block amplitudes in "additional arithmetic commands"
 commands                              are scaled with Is reset to 1.0 at entry.

 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] [F:fmt] <label>     Loop while i#n (start at #a and step by #s, all integers).
                                       If i > #n the loop content will be skipped.
                                       Defaults: #s=1, #a=1.
 ENDLOOP <label> [(condition)]         A Label is compulsory, must be keyed with brackets < >
                                       The loop variable is ###LOOP# =
condition - N> ###INDEX#  - continue loop while the index variable is within
                                                       the scope of the time series

                                          LOOP ${DOXXX:[1]} <XXX>
                                          ENDLOOP <XXX>
                                       is a way to skip a section using environment parameter DOXXX = 0.
             OBS! Changes 2014-03-29

(an idea: a FIND command, a value, an MRS, a gap of a minimum length, returning ###FOUND# - should need arise)

 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.

 NOOP[c] text and expressions [ time-range ]
                                       A test and debug feature. Expressions are evaluated and
                                       start and end indices of the time-range are reported;
                                       however, nothing is executed.
                                   c - If blank, the response is prepended with `<TsfEdi>N> ´
                                       and `NOOP´ is not echoed.
                                       Else, the command is echoed `<TsfEdi>C> ->NOOP [...] ;<-´


...................................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  [T=target]            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  [T=target]             String s is passed to openfs, specifying
                                       file unit, option, and file name. Otherwise
                                       like 'CONT #n'  above.

 REWIND {#n|O>}                        Rewind file unit #n. Any file may be rewound as needed.
 CLOSE  {#n|O>}                        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  
_PART2   or 

                                   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.       
APPENDDEL #b #e                        Sets a delete range at the start and the end of a series
                                       that is appended with a file-reading command, directive 'A*' or 'A!'
                               #b #e - Delete this much samples at the beginning and at the end, respectively.
                                       Specify a negative number to preserve the data. Is in effect
                                       until reset with APPENDDEL -1 -1

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

 COMMENT U=#u <<target                 Write comments to file unit u. Suppress if u < 0
                                       All text between this line and the next line that contains only the
                                       word `target´ is first subjected to variable replacement
                                       (# -"hash", % -statistics, $ -environment) before it is printed. 

 MISS [ time-range ]                   Report missing samples

 FAIL [ time-range ]                   Abort TsfEdit if a sample is missing within the time range.

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

 REPMISS [ time-range ]              - Reports number of missing samples in the time range
                                       (calls subroutine remdc_ct) and sets ###NMISS# e

 MODE #n [ time-range ]              - Reports the mode of the series (the most often recurring value
                                       in a 5-rms interval to both sides of the mean divided into n bins).
                                  #n - n ≤ 1000. n < 0 uses the previous number; default n = 1000.
                                       (No option available to copy the result to the scale parameter yet.
                                       If you call STATISTICS, you can use the %MOD variable.)

 REPDC [opt] [ time-range ]            Report DC-level to protocol (unit 6)
                                 opt - options may include a blank-free string of print-tsf options

                                N=#n - n = 2 requests also printing of the valid samples' "centre of gravity"
[/]->S[-] - set scale parameter, for use in a range of arithmetic commands.
the inverse value if /->S ),
invert the sign if appended by `-´

                                       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    [opt] [ time-range ]           Report RMS    
  [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 - options may include a blank-free string of print-tsf options

                               /SQRN - (simplistic) standard deviation of an arithmetic mean 
                           [/]->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.

 MEDIAN [s-opt] [opt] [ time-range ]   Report median to file unit 6, and save value in the statistics
                                       (retrieve with %MED). Optionally reports the time of the median in
                                       kalendar form.

                                 opt - options may include a blank-free string of print-tsf options

                               s-opt - ->S[-] store the result as scale value (invert sign with S-) for use
                                       in  a range of 
arithmetic commands.

                                  +T - in a special format:
                                       [<TsfEdi>>> ]MEDIAN date time MJD #imed :: median STDEV DMEAN stdev mean
^^^^^- if u = 6 (defunct)

 WMEAN [opt] [O=wmopt]] [
s-opt] [ time-range ]
                                       Report weighted mean
to file unit 6, and save values in the statistics
                                       (retrieve with %XWM ).

                                       The weights are expected in array X starting immediately after
                                       X(ndim). This is the case with  tslist -L'G|VAL' -L'G|SIG'

                    wmopt - passed on to subroutine weighted_mean (sas/p/wmeans.f):
                                +DBG - debug option for the scope of the subroutine. TsfEdit-wide DEBUG
                                       is obeyed.   
                              +SQ[D] - compute weighted mean sum-square.
                                       With D , on the deviations from the weighted mean.
                             +RMS[D] -
compute weighted RMS.
                                       With D , on the deviations from the weighted mean. 
                                  +N - By default, the results are stored in the Statistics.
                                       +N inhibits this.

                                 opt - options may include a blank-free string of print-tsf option    

                               s-opt - Further option:
                              ->S[-] - store the result as scale value (invert sign with S-) 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.

                                   I N F O R M A T I O N   T O   F I L E  

 This is a new project (Dec. 2016) with the ambition to harmonize the options to the few information commands
 that write dated values to file. This concerns the following commands:


 (STATISTICS)  PRTSF  (see the example) 

 The subroutine called for the purpose is PRINT_TSF in sas/p/tsfeditprt.f

 COMMAND [SET|RESET] [U=iun] +options [time-range]

 The +options part may have to be changed to e.g. PO=options
                         U=#iun  - Output to file (must have been opened on unit iun)
                                   The file unit is private to each COMMAND and remembered.
                                   Default is STDOUT.  

 Output consists of (default order):

 Date-information : Value(s) : COMMAND [:: commentary ]

                        +options - Common options, a string free from blanks. We shall refer to
                                   them as "+options".
                                   An option starts with `+´ or `-´ .
                                   The leading `+´ is mandatory 
                                   (e.g. +-D to mute the Value(s) delimeter `:´).

                                   +options are remembered until redefined or
                                   expressly reset.

                         Default - Taken from environment $TSFEDIT_PROPT or ' '

                              +  - as the only option: use the bare default.

                           +ENV  - as the only option: reset to environment $TSFEDIT_PROPT
                                   i.e. the default at startup.

                             +V  - Switches debug option on. Evaluates the same control
variable that's set by TsfEdit's DEBUG command.

                            +XC  - Date information is the centre of the time range:
                                   At yyyy mm dd  hh mm ss fff         
                            +XT  - Date information is the time range
yyyy mm dd  hh mm ss fff To yyyy mm dd  hh mm ss fff 
                             +I  - Date information is index numbers of time range:
                                From #iiiiiii To #jjjjjjj 
                            +IC  - Date information is index number central to the time range:

                             +Y  - simple format: Decimal year and Value, nothing else.

                          +Oabc  - where abc is a permutation of `D´ `V´ and `C´ , defines in what
                                   order Date, Value and Command are printed. DVC is the default.

                        +C:text  - replaces COMMAND with text . This option must be given last!
                                   The text string must be free of blanks;
Use `_´ as a placeholders.
                             -C  - Remove
COMMAND from the output.        
                             -A  - Remove the `At´ word normally ahead of centred date- (index-) information.
                            +Dd  - Use d as the delimiter surrounding the Value(s) (default: `:´); one character only!
                             -D  - Use a blank instead of a delimiter. 
                             -N  - If output is to STDOUT, the printed records are prepended with
                                   ` <TsfEdi>P> ´ by default. Use -N to mute this.    

                             +F  - Use the numerical format (optionally specified with FMT=format-code] )
                                   verbatim. By default, if a command generates more than one value, the
                                   format code is multiplied to make the result string fit on one line.
                                   A missing outer pair of parenthesis is always compensated.

                         +P[:d- If the command line contains d text , this text is added as a commentary
                                   (prepended with `:: ´ as the unmutable commentary delimiter).
(Thus, result values from one command can be preserved in the comment
                                   part of a new command, see the example.)
                             -P  - Don't print the comment part at all.

                            SET  - Parameters and options are defined, but the task is not executed.
                                   Therafter, a sequence of COMMAND Time-range lines may follow that don't
                                   need to repeat the parameters/options. The defaults set is put into
                                   effect before the specifications are processed.

                          RESET  - Like SET, but also carrying out the task of the command.

 SETPROPT +options [FMT=fmt]     - Sets the +options for all COMMANDs in this chapter (however not PRTSF)
                                   and provides an occasion to specify the numeric format fmt .

 PRTSF [SET|RESET] [U=iun] N=n [+options]  %statvar_1 ...  %statvar_n
                                 - Prints selectable results from the previous STATISTICS command, including
                                   its time range. See the example.

 An easy way to inspect what's going on is
 in the tse-command file.

(1)   Averages and Medians section-wise and simultaneous in two files.

First tse-file is produced by  ~/TD/a/bin/tp4sets : setcmds.tse
OPEN 31 B tmp/201405a-set-ra-means.dat
$COMMAND From #0 To #51535  ::  201405a
$COMMAND From #51536 To #51585  ::  201405a
$COMMAND From #51586 To #51635  ::  201405a

The output shall be: $COMMAND time-range the-results-from-averaging  (REPDC)
so that $COMMAND can subsequently be MEDIAN )
$ tslist o/scg-cal-merged-O-sdr-expf.ra.ts -Esetcmds.tse,CAMP01 -I > ! tmp/tslist-log
$ m tmp/201405a-set-ra-means.dat
$COMMAND From #0000000 To #0051535  :  -1.673780E-03 :   ::  201405a
$COMMAND From #0051536 To #0051585  :  -4.842953E-01 :   ::  201405a

We could make the next step a little nicer, like reporting to a file again, but we'll
be content with STDOUT. Thus, let's make a tse-file from the output simplest way:
$ wraptse -o tmp/201405a-set-ra-means.tse R tmp/201405a-set-ra-means.dat
$ tslist o/scg-cal-merged-O-sdr-expf.wr.ts -Etmp/201405a-set-ra-means.tse,R -I | fgrep '>P>'
 <TsfEdi>P> MEDIAN :   1.474900E-03 :  At #0019542  :: -1.673780E-03 : 201405a
 <TsfEdi>P> MEDIAN :  -4.393862E-01 :  At #0051575  :: -4.842953E-01 : 201405a
 <TsfEdi>P> MEDIAN :  -4.016674E-01 :  At #0051612  :: -4.092216E-01 : 201405a


(2)   Look at TD/a/bin/tp4sets ,
$ tp4sets -p <campaign> [<campaign...>]

which exercises a.o.
      awk -v c=$c -v g=$gap 'BEGIN{o=0} {if(sqrt($2*$2)>g){print "$COMMAND From #"o,"To #"$1-1," :: ",c;o=$1}}' > ! tmp/tmp-$n.tsi
        '/^ <DumMer>>> /{if ($5=="::" && substr($4,1,2)=="SL" && $6 ≥ b && $7 ≤ e ){print "$COMMAND From #"$6,"To #"$7," :: ",$11}}' \
  setenv COMMAND "REPDC U=31 N=1 +IC-P"
  setenv COMMAND "WMEAN U=31 +IC-P"
  setenv COMMAND "RMS U=31 /SQRN +IC"

to produce a plottable series (altough not without further awks),
here for campaign c = 201405a  what = prj  type = ra

$ cat tmp/$c-$what.$type.dat
1792.512905 0054805 -3.982267E-01 5.230419E+01 6.454223E-01

1793.744387 0060511  2.856266E-02 5.205854E+01 7.613090E-01
1794.791840 0065253 -4.307186E-01 5.146370E+01 7.558082E-01
1796.334433 0072397 -1.136964E-02 5.049295E+01 5.163716E-01

central time of prj
day f.start  index     Average    MeanDropStd   UncertAve

$ cat tmp/$c-$what.$type.txt
1792.512905 0054805 -3.982267E-01 5.230419E+01 6.454223E-01 Onsala_AC_20140527a.drop.txt FG5-233 north
1793.744387 0060511  2.856266E-02 5.205854E+01 7.613090E-01 Onsala_AC_20140529a.drop.txt FG5-233 south
1794.791840 0065253 -4.307186E-01 5.146370E+01 7.558082E-01 Onsala_AA_20140530a.drop.txt FG5-233 south
1796.334433 0072397 -1.136964E-02 5.049295E+01 5.163716E-01 Onsala_AA_20140531a.drop.txt FG5-233 north

$ cat statcmds.tse

OPEN 31 B tmp/201405a-stats.dat
STATISTICS From #51486 To #58124  ::  Onsala_AC_20140527a.drop.txt
PRTSF +IC+C:STATISTICS_MIN_MAX_AVE U=31 N=3 : %MIN %MAX %AVE :: Onsala_AC_20140527a.drop.txt

$ tslist o/scg-cal-merged-O-sdr-expf.ra.ts -E statcmds.tse,C -I
$ cat tmp/201405a-stats.dat
At #0054805  :  -5.120114E+00  4.006673E+00 -3.982267E-01 :  STATISTICS MIN MAX AVE  :: Onsala_AC_20140527a.drop.txt

Explanations:  +C:string  replaces the command name normally been printed. Underscore is replaced with blankspace.
N=3 :  value value value  in PRTSF reads the stat-vars after their conversion to numbers

(4)   very simple task: Print the DC-level without any decorations

An arithmetic option:
In commands that set the scale variable, the value can be stored using
and called back using  %V#i
Not proven yet whether a sign change or a reciprocal quote of the value (-%V#i  or /%V#i) will work.
 REPDC ->V01
 ADD %V01
This does not collide with the save option ->S

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

 X->W                               - Copy series to work space
 W->X                                  ... vice versa, carrying over the dimension and the timing.

 V#i = v                             - Store this value in a variable that can be quoted back with %V#i
0 <= i < 100 

 CUT #n                              - Cut away n samples from the end
 RETAIN time-range                     ... with adjustment of starting time and length.
                                       For a suite of retains (delete-inbetween) see below

 CLEARW [ time-range ]               - Set up a work area the size and timing of X, filled with
                                       MRS according
to the specified time-range.

 X->W  time-range                    - Copy a segment of the series to work space
  time-range                      ... vice versa, without carrying over the dimension and the timing.
                                       CLEARW, a suite of X->W time-range, and W->X in sequence provide a means
                                       to preserve segments of a file, i.e. the opposite of DEL or RJC

 WCDEL +P                            - Similar to CLEARW to be followed by one or more X->W time-range
time-range                      but here, the data in X is deleted.

 SUPPLW [+A] [ time-range          - Supplements data from W into X where X has gaps, and optionally
                                  +A - appends X with W if W ends later
(subroutine sas/p/supplts.f).
                                       The mean of W is adjusted to the mean of X .
                                       (The difference of the means is computed from the simultaneous subsets.)                                           
                                       OPEN 31 ^ ${SUPPLFILE:[/home/hgs/SMHI/o/tgg-rin-100301.ts]}
                                       @31 +K
                                       SUPPLW +A
                                       supplements Ringhals tide gauge values into gaps of the OSO Bubbler

 X*W O=s [+P] [ time-range ]         - Uses mult_ts_ts (sasu061.f) to perform X op W, result in X
 W*X O=s [+P]
[ time-range ]         - Uses mult_ts_ts (sasu061.f) to perform W op X, result in W
                                   s - Operator and options, starting with / for divide or * for multiply,
for squaring the second series before the operation,
                                       Q for square-root,
for Let (don't synchronize),
for replacing missing records with DC-level.
                                       Default operation is multiplication (`*´)
                                  +P - Talk option = .true. in mult_ts_ts.

 XWMIN                               - Retain in X the minimum or ...
... maximum, respectively, of Xj, Wj , 0 < j < N-1
                                       An example is in ~/Ttide(SCG/stackextremes.tse
                                       where the min/max of a stack of files is computed.
                                       (A new feature, not born out yet. For instance there is no
                                       adjustment nor test on synchronicity.)

 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][R=#m] [+S]  [ time-range ]  
                                       Fit Legendre polynomials up to degree n and ...
                                   I - return the polynomial prediction (default: the residual)
                                R=#m - evaluate the polynomial up to an order m < n
                                   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
                                  +S - Print slope results per hour and per year

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                         (source file sas/p/sasu06.f)                                      
                                   n - polynomial order
                       p0 p1 ... p#n - n+1 coefficients in the following lines.
                                         X.out <- SUM ( p(j)^j , j=0..#n )

 CBS [+W] At { YYYY MM DD MM HH SS FFF | #j }
                                       Connect at a jump, "Connect Boxes and Slopes"
                                  +W - After X->W the jump height in X is reduced by the jump in W
                                       (useful to connect added time-limited slopes and other
                                       synthetic signals). The compensating offset is added to all
                                       ordinates from the indicated time to the end.

 REPLACE #v [|] [X<#xmin] [X>#xmax] [ time-range ]
                                     - replace X.i with v if xminX.i xmax
                                       With `|´ take the clause with |X.i|
                                       Default xmin=-1.e20, xmax=1.e20

 CLIP #a #b   [ time-range ].........-.Take min(max(, a),b)
ABS          [ time-range ].........-.Take absolute value of
 X**#n        [ time-range ].........-.Take power n of
 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 font-weight: bold;">                                        0.5 0.5
                                       SHIFT T=-0.5r

  SYMM [+A]  [ time-range ]           - Make series symmetric in time (in the sense of Fourier domains,
                                       i.e. the first sample is kept unaltered, the sample at or directly
                                       after the centre is kept unaltered too.
                                  +A - make anti-symmetric (first sample will be set to zero).

............................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.

[+UG] 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.
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
                                     - Apply a saved ¨(n=0) or stored (n>0) filter                        
                                     - Fetch a saved filter                                               
- Fetch, alter, and save a filter                                     [n=1]
                                       N.B.: the operation is done on the actual filter; a change
                                       of code to operate on the stored copies only might be desirable.
FILTERCOMMAND [+C] [+UG] options      Compute or read, and apply instantly.
                                  +C - Apply on a circularly connected domain.
                                 +UG - Rescale to unity-gain (changes coefficients).

 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

 FILTER ... [APPLY...] +K            - (executing; only with this command!) Useful with 2-coefficient
filters: Keep time origin. Example:
                                       FILTER 0 1 ' ' +K
                                       1.0 -1.0

 FILTER ... APPLY... ON:{W|#c}       - (executing; only with this command!) Apply the filter on ...
                                   W - on the work array
                                   c - under tslist: on a data column;
                                       c is the number+1 of a column read in e.g. using labels.
                                       The filter should be applied on X last (no `ON:´ option)
                                       if time origin and length is to be preserved.

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) [STORE->#n]
                                       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)
                                       Default: don't rewind.
                                       +N  - unit-normalize the filter
                            M=string - input the filter that has been marked with string
                                       If string contains blanks, right-delimit it with `)´
                                       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 programs
                                       fslist options filename

 TRUNCF #m,#n                        - Truncate a filter
 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..), i.e. only 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  [ SAVE | STORE->#n ]
                                     - convert time series X or, respectively, W to a filter.
                        SAVE | STORE - like in the FILTER commands
                                 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)

  FILTS U=#u [ F=f ] P=#m #n 's'] [-0ENDS]   
                                       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

#m #n 's' [+N]                 Input coefficients ASCII

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

.............................#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 NLIN=,                
  f(0),f(1),...,f(n)                  - autoregressive non-linear system. See Example 19
- mixing coefficient for all past y ("output")
                                  #ξ - mixing coefficient for x ("input")
                                       where yi = xi f0
(1 + ξ yi-1) + sum(k=1..n) yi-k fk (1 + ζ yi-k)

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

                                       The filter block must begin with one record supplying #m #n 's'
                                       and the search string, if used, must match it exactly (!)
                                  +N - normalize filter gain during convolution
                                 +UG - adjust the filter to yield unity gain at f=0

                              -0ENDS - (applies to all versions of FILTER with file input):
                                       Cut the ends of the filter if these coefficients are 0.

 WIENER U=#u F=f L=#lw  [-0ENDS]       Process a time-domain transformed gain function
                                       (a.k.a. impulse response), for instance, use output during
                                       sasm03 impulse response display; you can taper it here
                                       using a HANN window of duration #lw
                                   f - a file name, ASCII 

                                       File record with filter coefficients: READ (IU,*) I,F(I) I=1,LW
                                       i f(i)


 SIGFILTER #n #m [->#t]                Convert a filter (SAVE'd in storage n) into one for uncertainties
                                       using an autocorrelation series (SAVE'd in storage m) and
                                       save the product in storage t. Default t=1
                                       The autocorrelation series may be read in with GETFS ... SAVE->m
                                       The filter is most often a decimation filter (low-pass). The idea:

                                       FILTER D:WD ...          STORE->2
                                       GETFS 41 < o/mem-acv.fs] STORE->3
                                       SIGFILTER 2 3            STORE->1
                                       FILTER                   APPLY<-1
                                       SQRT +F
                                       DECIMATE ...
                                       (new command, introduced 2017-05-10 and yet untested)

 PEF [U=#u F=f] L=#l O=opt [+F]        Use a Prediction 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

.......................................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] [C@0] [U] [+TR] [{ P | SAVE }] [ time range ]    

                                       Set up a window and, if ` P ´ ("prepare") is missing
                                       (SAVE is synonymous with P), apply it on the data.
                                       The length of the window is determined from the time-range.    
                                 +TR - The series outside the time range is zeroed by default.
                                       +TR prohibits that.
                                   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].
                                 C@0 - Center window on start of series (useful for AUTOCOV)  
                                                                                    [centre of time-range]
[:v] type RESTORE [C@0] [U] {[{ P | SAVE }] | [ time range ]}   - (SAVE not with time range!)
                                       Uses the parameters (incl. length) of the most recent

 WIN[:v] [U] [+TR] [ 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.

 MOVXC[ORR] #m #l O=opt              - Moving (sliding) cross-correlation, X * W
                                       If a Window has been declared (with SAVE), it will be used
                                       in the sums. 
                                  #m - index increment for shift of interval, must be > 0.
                                  #l - length of interval
                                 opt - `XG´ to output  sum X*W/sum W2,  `YG´ to output sum X*W/sum X2
                                       else sum X*W/sqrt(sum X2 sum W
                                       Returns the result in X.

 XCORRU  [R=#kb,#ke] [P=#iun] [KU=#iui] [+LET] [S=#apshift] [ time range ]
                                       Cross-correlation between two data sets in X and W 
                                       with different time origins and sample rates as long as
                                        - they have an overlap
                                        - their sampling rates are commensurate
                                        - W has the longer time step
                                +LET - is liberal on sampling window, but does not replace
                                       series X with the XCORR-values. Thus, many XCORRU commands
                                       can be executed in succession.

                                       Without +LET , returns in X the time series in synch with the
                                       one in W and delayed with the highest cross-correlation
                                       (incl. the effect of the time-range). See Example 22. .
                             #kb,#ke - the desired lag range centred around apshift
apshift - an a-priori shift that can be set using using S=value     [0]

                                #iun - file unit for results. The correlation series is
                                       only printed; it is not saved in memory as of current.
                                #iui - the file unit from which W was filled; is used to
                                       generate tslist / tsfedit code on iun for output of
                                       time series in synch with each other.

 CCV [L=#l] [M=#m] [D=#d] [+CORR] [+SLIDE]
                                       Computes cross-covariance between X and Y
                               +CORR - Computes the correlation.
                                       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.
                                       If not +SLIDE , the covariance or correlation
                                       -l ... l series is returned

                              +SLIDE - The centre covariance or correlation is returned
from a moving interval with stride m

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

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

                                  +W - Keep result in W else return the series in X with dt=1

                                       (Hints: Write out X from tslist with options -Ni -nr
                                               Write out W from tsfedit with PRINTW )

 AUTOCOR[RELATION] #k [+W] [+BW] [FRA>#f] [time-range]
                                       Compute auto-correlation up to lag k
AUTOCOV[ARIANCE#k [+W] [+BW] [FRA>#f] [time-range]
                                       Compute auto-covariance up to lag k
                                       Both commands use W for intermediate storage of result!
Specify an absurdly big lag to compute all terms.
                                  +W - Keep result in array W 
                                 +BW - Apply the Bartlett window
FRA>#f - Require at least a fraction f of possible value pairs      [0.0]

 ALLANV[ARIANCE] [L=#m] [O=opt [time-range]
                                       Compute Allan variance up to ...
                                M=#m - interval 2m                                                  [20]
                                 opt - Options (code without intervening blanks):
                                       +P[:#u] - print to file unit u                                [6]
                                                 (the results are retained in array W(0..m-1) )
                                       +ADEV   - compute the Allen deviation                  [variance]
                                       /DT     - divide by time interval in seconds              [don't]

                                       Results, format:

    i         2i    count   value
    0         1     57268  4.0873E-01
    1         2     56944  1.0750E+00
    2         4     56466  2.2024E+00
    m         2m        

 MEMSP #n    [R:s] [L,][U[=#h,]][Q][{/N|/SQRN}][dB] [OMEMSP=#u[,name]]
                                       After producing or reading a PEF, compute the
                                       Maximum Entropy spectrum, domain length n.
                                       Options see PERIODOGRAM below.
                                       See Example 24 how to compute a series of MEM-spectra with
                                       a selection of filter orders.

 PERIODOGRAM [R:s] { [L,][{R,||R||Q,|I,}][U[=#h,]][{/N|/SQRN}][,dB][,+DC][+R][,+CW0][,M<#m] | , }
                   [/Hz] [/MISS] [OPSP
=#u[,name]] [ time range ]

                                 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
                                       The spectrum series will obtain dt=1
                                       The frequency interval and Nyquist is reported through
                                       COMMON /CPDGRAM/ in units of cyc/s cyc/h and cyc/d.
                                       See the -nFrq option of tslist for suitable frequency
                                       tagging (options for automatic scaling are largely untested yet).

                                       "Report" = Print the Nyquist-frequency and the spectral bin
                                       interval on the protocol.
                                       Read on about the other options...

 PERIODOGRAM +W[SQR][,+DC][,M<#m] OCSP=#u[,name] [ time range ]

                                  +W - With +W, array X is unchanged (except editing and
                                       DC-level removal); the complex spectrum will be scaled with
                                       /N or ...
                               +WSQR - Like +W, /sqrt(N) however.
                               OCSP= - The complex spectrum is written to BIN sp-file using OUTSP (ufiosps.f).
                                  #u - unit number of a file that must have been opened.
                               ,name - max 16 characters, string is written to the label of the
                                       sp-file. Read back using GETSP (main: splist).
                          Without +W - Return the periodogram in X. If called from tslist, the
                                       output series should be printed using -N -nFrqu .
                                       Before the periodogram is computed, the time series
                                       is repaired (linear interpolation of missing data)
                                       and the DC-level is subtracted.

                               OPSP= - The spectrum (in the X-array) is written to BIN sp-file using OUTSP.
                                       #u and ,name like under OCSP 

                                       The following options must be coded without
                                       intervening blankspace; prefer comma:

                                   , - use all defaults
                                   Q - return the power                                   [Amplitude]

|R| I - return real-part, abs(real-part) or imag-part,
                                       respectively                                       [Amplitude]

                                   L - return the logarithm (base 10) of the spectrum.       [Linear]

                                  /N - scale power by the number of samples                   [unity]
                                       Automatically respected in Amplitude mode:
                                       scale with the square-root. BUT!...

                               /SQRN - scale power by the square root of the number of
                                       samples (needed with dB
because of clumsy programming, sorry!)

                                  dB - return decibels         
                             U[=#h,] - Compensate the effect of a filter (1,-#h)   [don't, else h=-1]

                                +CW0 - compensate window's mainlobe leakage (due to DC-removal)
                                       to the lowest-frequency bin (a WINDOW must be declared and
                                       applied; the unity-response option is recommended).
                                       This option is not recommended unless a REMDC is used.

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

                                       The following options must be coded after 
                                       intervening blankspaces:
                                 /Hz - Scale output to spectral power density (per Hz)

                                 +DC - keep the DC-level                                     [remove]

                                  +R - Reset DFFTRI-work array
                                       (applicable in loops / repeated calls)                  [keep]

                               /MISS - adjust power for missing data

{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[C] [{-XT|+XT}]] [-SPP] /
[+D] [+H[-B] [-F|+J] [+STEP|+DIFF] [+P] [-P] [[dB]] [SU=#u] [DU=#u] /
[FSPO=#u[:label] /
                [context-options] [ time range ]
                                       Complex filter applied in the spectrum domain.
                                       Several models can be applied (polynomials in z, pole-zero,
                                       Bessel, Butterworth, 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.

                                       B.t.w. band limitation can be used to spectral z-filter
                                       a time series with a sampled impulse response of an analog
                                       filter with TRIVIAL
                                       (the Impulse invariant method of Oppenheim & Schafer Chap. 5.1.1).
                                       The crucial point here is the time scale of sampling.
                                       See Mathematica/Bessel-filter-analog.nb on ORE).

                                       Context-dependent options: The context is given by
                                       the type of the filter operation.

                                  /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
                                       (creates the analytical complement)                    [don't]
                                       Simplest application: ZSPF +H

                                 +BL - band-limit with a recently stored FILTER
                                +BLC - like +BL , apply complex-conjugate, however
                                       (needed with asymmetric filters)                    [straight]
                             +XT -XT - With +BL : Crop the time domain to avoid end-effects.
                                       +XT sets the new origin via spectral_zf_time
                                           which implies extension before processing,
                                           which only works with CPAD and a specified offset;
                                       -XT sets the time origin and length to fit the
                                           shrunken series.
                                       By default, time domain is like filter-in-place.
                                  -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]

                     CSPO=#u[,label] - Output complex spectrum to file (unit u must be open)
                                       using subr. outsp (ufiosps.f)
                                       Use PERIODOGRAM +W,+DC
OCSP=#u[,label] instead!
                     FSPO=#u[,label] - Output filter spectrum (incl. band-limiting filter) to
                                       pre-connected file unit u using
subr. outsp (ufiosps.f) 
                               SU=#u - output spectrum to ascii file on unit number u
                                       (instead of the print to stdout)                       [don't]
                                [dB] - output and print (!) spectrum with amplitudes in dB 
                              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]
                                -SPP - suspends print of spectrum                                [do]
                              FF=#fs - a scaling factor
in the output on dimensionless
frequency, in effect if [F]= is not used               [1.0d0]
                            [F]=unit - In the output list (+P), the frequency that
                                       is displayed is i/N * tscale * fs, i=0..
                                       where tscale is derived from the [F]=expression
                                       (N is at the Nyquist). 
   [dimensionless so that frq=k/N≤0.5]

                               +DISC - discrete version; poles and zeroes in variable exp(iω)
                                       The analog (continuous) version employs the bilinear
                                       transformation and is the default.
                                       iw = 2 wscale/dt (1-exp(2pi i k/N)/
(1+exp(2pi i k/N) 

                        [w]=rad/unit - unit = { s | h | d } scaling of angular velocity in
                                       analog filter designs.
It is ignored with option +DISC
                                       Cases POWER P&Z ZPOLY RPOLY :
                                       the unit is significant in the bilinear transform,
                                       parameter wscale in spectralzspf.f
                                       Cases BESSEL and BUTTERW : The bilinear transform is
                                       carried out in the subroutines besselfs.f and
                                       butterws.f, respectively (unless +DISC is given).           [h]

                                type - TRIVIAL
                                       The rest of the options is listed under P&Z
                                       TRIVIAL is the special case of P&Z:0,0
                                       Example (12.2) shows how this type can be used to fast-filter
                                       a long time series with a long filter.

                          type:order - BUTTERW:#n   Butterworth design, where n is the
                                       filter order.
                              FC=#fc - design parameter, scaling time, give a unit in terms
                                       of time, e.g.  0.2 s
                                       Most of the other parameters of BESSEL can be used.

                          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

                                 +SR - OBS! more than just show: calculate the coefficients
                                       the roots way.
                                       (For increased functionality this might have to be changed:
                                       parameter-setting call sequence, default handling and
                                       tsfedit command-options.)

                                       With +SR there's a possibility to vary the roots individually:
                                       Command VARBFR preceeding ZSPF BESSEL:...

                                  +T - Trace filter action (frequency and response) from
                                       within Bessel-filter subroutine
                               +GRPD - also print the group delay
                               RL=#p - In ZSPF BESSEL**-1   crop the filter gain top at 10p
LL=#p - In ZSPF BESSEL (any) crop
the filter gain bottom at 10p
                                       Default is no cropping.

                              Example: ZSPF BESSEL:8 SU=44 FC=0.0174 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
                                       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

                    type:poles,zeros - P&Z:code
with the same options as above.
                                       Internally defined
sets of poles and zeroes can be
                                       retrieved by code:
                                code - STS2G - The Streckeisen STS2 Generic 120s Seismograph
                                       CMGT3 - The Guralp 120s Seismograph
                                       L1HZ  - The Lennartz 1Hz Seismograph
                                       (DATA statement in sas/p/spectralzspf.f
                                       Subroutine get_zeroes_and_poles) 

 VARBFR #n #v1 ... #v              - Preceding ZSP BESSEL: ... +SR ...
ary root i by 1+vi , all are read in at once.
 VARBFR [+R] [N=#n] C=#i,#vi         - alternatively: n is 2 times the number of roots.
                                       This version is cumulative.
                                  +R - will reset the array v
                                       n must be given at least once, else the default 32 is used.
blanks in the parameter pair after C= cannot be avoided,
                                       use `]´ as the right delimiter.

 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  [time-range]         Computes a series of length n, where each sample
(i) = SUM {j /. MOD(j,n)==i : }
n - cycle length

 MODSET #m #n #v [+EQ] [+T] [time-range]   
                                       Set x(j) = v  if (( MOD(j,m)==n ) XOR ( not +EQ ))

                                  +T - Normally j is counted from the start of x (i.e. x(0)),
                                       not from the start of the time range. If j is to be
from the start of the time range, use this option.

                                       To introduce zeroes between RESAMPLED values in a 10s-series:
                                       RESAMPLE 0 0 1.0 [s]
                                       MODSET 10 0 0.0
                                       To set every 10'th value to 1.0:

                                       RESAMPLE 0 0 1.0 [s]
                                       MODSET 10 0 1.0 +EQ

 AUTOSEGMENT [U=#u] [C=cmd>] [M=#m] [D=d] [+L]
                                       Writes Time-range records for detected segments prepended by cmd .
                                       A segment is defined as starting and ending with a valid sample
                                       and containing at maximum #m missing values.
                                 cmd - A command string, delimited by `>´ . A different delimiter
                                       can be specified.              Default cmd is ${COMMAND:[COMMAND]}
                                   d - delimiter for the command string                               [>]
                                  #u - output file unit. File must have been opened before.           [6]
                                  #m - tolerance limit for missing samples, which also is the
                                       criterion for the minimum gap length that defines two segments.
                                  +L - import lines written before (B:) or after (A:) each
                                       Time-range record, or once and first (F:) or once and last (L:):
                                       This is useful if segments are to be exported to different files.
                                       #HASH and $ENVIRONMENT variables will be replaced. Example:

                                       OPEN 32 B makesegs.tse
                                       AUTOSEGMENT U=32 C=OUTTS 31> +L
                                       F:TSF EDIT MAKESEGMENTS

                                       B:OPEN 31 B ${DESTDIR:[.]}/segment-###.dat
                                       A:CLOSE 31


 AUTOSLMEAN [O=subopt] [U=#u] [{ +W | ->W [-W] }] [G=#n] [{+DOY|+AG|+TSF}] [time-range] [:: comment text]
  Auto-detected sliding mean by gap width exceeding n samples.
                                       A method to reduce Absolute gravity drop-data to set-means,
                                       optionally weighted means. For weighted means under tslist,
                                       the next column must contain standard deviations. 

                                       This function exploits a peculiarity of the input array X:
                                       Normally, a signal channel ("column") is handed over one at a time.
                                       However, the later channels are available at multiples of the storage
                                       offset NDIM. So here, weights can be imported in the channel next
                                       to the signal. The +W option uses these weights. The work space
                                       (W(i)) is a separate space, so the result can still be put there.
                                       The notations  +W  and ->W might be a little confusing.

                              subopt - option passed on to sliding_mean_auto (sas/p/autoslidemeans.f)
                                 +SQ - compute the weighted sum-squares, return it together with the
                                       sum of weights and the number of accepted data values.
                                   n - decides when to output the mean over the recent bunch of samples.
                                  +W - Use weights in the next signal channel (under tslist).
                                 ->W - like +W, return results in array W (default is X).
                              ->W -W - no weights,
return results in array W.
                                   u - output unit (ASCII data, tsf-style). Default = no ASCII output.
                                       Without +SQ, the data columns are:
                                       (weighted) mean, the standard deviation (weighted),
                                       the number of samples accepted;
                                       and optionally  `:: comment text´.
                                +DOY - ASCII output: date format: with month number,
                                 +AG - ASCII-output similar to g-soft, uGal
                                +TSF - Standard TSF format (also the default), nm/s2
                               +JTSF - Standard TSF format + float MJD, nm/s2

                                   n - decimation factor
                                   s - initial shift

 DECISYNC #m [#n [#u]]                 In anticipation of DECIMATE, synchronize sample origin to integer
                                   m - m dt. E.g. m=100 to align
a 100 Sps series on integer seconds.
                                       Thus, you don't need to anticipate the MOVE parameter.
                                       Operation is a SHIFT with k; however simple as that would be, the
DECIMATION section of the code is used, employing W as a workspace. 
                               n [u] - If specified, DECIMATE n k u  will take place.   

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

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

 DECI[MATE] #n #s [#u] G=#f,#m         ...with Gaussian-bell filter
                                   m - half-length.
Total length = 2m+1.
                                   f - Filter's -3dB point is at f = fraction of Nyquist.
                                       Small m will lead to cutoff effects.

 DECI2W #n #s [#u] [time-range]        Decimation (undersampling) to work array

  #n #s [#u] [O=options]          
                                       Three kinds of decimations,
taking the max or min or median
value of the data segments.
- decimation length
                                   s - initial shift
                                   u - a "hidden" initial shift (no effect on time origin)
                            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 - options string, may contain
T:#lt,#s[,+P]  - with r>1, to taper the spectrum at the
                                                        suture using the positive half of a
                                                        Gaussian bell
                                       S:#sh          - to shift the time series

                                       lt - length of the taper, number of bins.
                                        s - width of Gaussian bell, argument (i/(lt*s))2, 0≤ilt
                                       +P - print tapering action on protocol.
                                       sh - A real-valued shift, positive values shift ahead,
                                            in sampling intervals

 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 #dtr [T!] [s] [FILL=#v]     
                                       Resample with polynomial interpolation, using subsample_ts
                                       Uses W as a work array.
                                   o - polynomial order (sorry, we can only do linear if dtr < dtin),
                                   t - shift of origin time
                                 dtr - sampling interval

                                       origin time is with respect to the previously first sample,
with respect to epoch (T!)
.......................................t and dtr are in hours or in seconds ([s]).
v - If o=-1 the fill-value is substituted at interstitial times.
                                       If o=0 a sample-and-hold value is substituted.

 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] [N=#n] [O=#t] [+A] [D=delim] [ [[+Q[Q] ]->#j] At #Y #M #D [#h [#m [#s [#f]]]] ]
[U=#u] [N=#n] [O=#t] [+A] [D=delim] [ [[+Q[Q] ]->#j] At {##s |#N |#N-#s }   ]
                                       Sample by polynomial interpolation. The At-part can be left out; then
                                       the output options are set for the subsequent SAMPLE commands.
                                       (+Q , if desired, must be given in each line).
                                   t - (integer) offset in seconds w.r.t. the date given after `At´.
                                       Must be reset expressly.
                              At ... - 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.
                                   n - Polynomial order +1. Default: Linear.
                              -A  +A - Print output even if the sampled value is invalid.
+A can be negated by -A .
                                 ->j - save value in a storage (retrieve with %j; for details, c.f. STATISTICS)
->j - do this quietly (overrides +A).
                                 +QQ - don't write at all
                              -Q -QQ - negate a previous +Q or +QQ 
                                  At - this keyword must be given if an option is specified.
                            +C[=del] - copy the comment, the remainder of the line beyond the delimiter del.
                                       Default delimiter `;´)
                             D=delim - The delimeter in front of the sampled value, up to 16 chars long. Default `:´
#s - The second method uses the sample index. An index relative to
                                       the end of the series is specified by  #N  (meaning N-1)
                                       or  #N-#s (meaning N-s-1). See Example (20).

 STATISTICS [+S] [options] [time-range]       Save some key characterizing values in a storage:

                                       %MIN, %MAX, %AVErage, %RANge, %RMS-dev, %MODe,
                                       %NVS = N.Valid Samples, %RVS = sqrt(%NVS) .

                                       %MEDian is saved when then
                                       MEDIAN command is issued.

                             BINS=#n - an option for the MODe calculation: The 5*rms ranges centred on the mean
                                       is divided into n bins. Default (= maximum) is 1000. 

                                       %XWM (weighted mean), %
WSM (sum of weights), %XSQ = %WST (WRMS), 
                                       and %WSQ (sum of weights-squared) are saved when the
                                       WMEAN command is issued (without the RMS option!).

                                       Up to 10 samples of the series can be stored with the
                                       SAMPLE command with option ->#j  
                                       (e.g. SAMPLE ->1 At 2014 01 28 22 37 00 000)

                                       Retrieve in any command line using, respectively,
                                       %MIN %MAX %AVE %MOD %RAN %RMS %NWM %NVS %RVS
                                                                      - computed  upon STATISTICS      command
                                       %XWM %WST %WSM %NWM            - computed  upon WMEAN           command
                                       %XSQ %WSQ %NSM                 - WSM = WSQ upon WMEAN O=+SQ     command
                                       %MED                           - computed  upon MEDIAN          command
                                       %j (1 ≤ j ≤ 10)                - set by commands that support the ->#j option
                                       To invert the sign and/or taking the reciprocal value, prepend
                                       the `%´-symbol with ` -´ , ` 1/´ or ` 1/-´ 
                                       e.g. ` 1/-%1 ´  or ` -%MED´
  or ` -%XWM´ )

                                       The values can be modified with a trailing arithmetic expression
                                       (one basic operator and one constant for each %VAR) like in
                                       %AVE+122  %RMS*5 or %RMS*%RMS . This option has no full-fledged arithmetic parser!
                                       There is no analysis of parenthesis and no operator precedence list.

                                       The values, if floating-point, can be formatted with a trailing expression
                                :fmt - Fortran format code, e.g. %MIN:f10.4
                                       This code is not remembered. The default 1p,e20.12 can be modified with ...
 SET[ ]%FMT=fmt                        and, with fmt specified as `default´ , reset to the "factory" default. 

                                       The consistency of prepended signs and reciprocals with the appended
                                       arithmetics has not been tested yet (2014-06-24)
                                       You can try and use the F`arithm.expr.` mechanism, e.g. F`28*%MED`

                                       See Example (18)
                                  +S - Show

 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] [/yr]
                                  /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. 1) 

 PRB  #v [u] [ Iref=#n | Tref yyyy mm dd hh mm ss fff ] [ time-range ]     
                                  add a parabola; same syntax as SLO

                                  However, all units for time scaling must be coded
                                  with an appended `^2´, e.g. `[h^2]´ or

 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.
                                  with Iref= :   a exp((j-n) p dt 
                                  with Tref= :   a exp(j p dt - tref) 1)
                              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] [Iref={B|C|E|#n} | Tref yyyy mm dd hh mm ss fff ] [ 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] [Iref={B|C|E|#n} | Tref yyyy mm dd hh mm ss fff ] [ 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)

                        Iref=#n - Set the time reference of the feature to time series index n
Iref={B|C|E} - Set the time reference of the feature to the time range's begin, centre or end.
                           Tref - Set a reference point in time.

                           Tref   Caveat: The date and time is tentatively read down to the millisecond level.
                                  Thanks to an ERR and an END clause, the operation may retrieve what you intended.
                                  A right-delimiter that is not a valid numerical symbol can be used
                                  to make do without hours ... milliseconds (all default to 0)

#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 [Iref={B|C|E|#n}] [ time-range ]
add a rectangular wave
offset+scale*{(MOD(i+ip,nf) ≤ iduty?) +1 else -1}

 TRI #a #nf #ip #iduty #offset [Iref={B|C|E|#n}] [ 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 ] [ C=#c ] [GMP=#g] [ time-range ]
                              a - amplitude
                              s - random seed
                              g - Gauss-Markov parameter: y <- g y + ran(s), x += a y

                 type = PEF     - Uses a Prediction Error filter transformed to the spectral
                                  domain. The filter must be set up using command
                                  PEF [options] SAVE
(new on 2018-05-06, so far untested)

                 type = FRACTAL - with power-law parameter p and a constant c
                                  P = a
/(c + f-p)
                                  Default c = 0, no default for 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 is indefinite)
..............REPAIR L≤#n        repairs breaks of length n by linear interpolation
..............REPAIR              finds all breaks (REP ALL)
..............(REP L ALL is perhaps too difficult for the program.)
                             +P - to print out the places to repair.

[ time-range ]  - Fill gaps in X with W, using interpolation if necessary
                             #m - polynomial order, 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-squares 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; beware 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 and 2: linear-varying sine amplitude, 3 and 4: linear ramp)
                               5: also the frequency parameter (not recommended,
                                  however the default is 5 (!)

 DESPIKE #t [+D] [ time-range ] - Remove spikes with the criterion
                                  |x(j) - Xmean| > t * Xrmsdev
+D - delete the samples, else fill in Xmean
+R - threshold t specifies the factor on probability 1-1/N
DESPIKE #t +S W=#w1,#w2
[ time-range ]  - A simple version with the criterion
x(j) - x(j-1)| > t
                                  replaces x(j) with the mean of x from x(j+w1) to
                                  Note, if w2=0 and w1<0, the mean is taken from the already
                                  despiked part. 

                                  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?

 DEL  [ {<|>} #v ]  [ time-range ]   
                                - delete (criterion abs(x(i)) < or >
v , repectively)

 DEL  [ {<<|>>} #v ]  [ time-range
                                - delete (criterion x(i) < v
or > v , repectively)
                                  If no operator and value is given, the scope for DEL is
                                  the 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.)

 DELBYPOS U=#u[ < filename]]    - using a *.po.ts-file, makes elements if X missing.
                                  Such files are produced by urtapt.
                             #u - File unit number. If opening option and file name
                                  aren't specified, the file must have been opened before.
                                  Doesn't work with W yet.
                                  N.B.: Uses work space directly behind X(n). Should be moved
                                  far away.

 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. 

 REMOVEJUMP [ +L ][W=#w] [O=#l,#r] { From | At } YYYY MM DD HH MM SS FFF
                                - removes a jump by forcing DC-level of left side of a
                                  breakpoint (From-date) to equal right side.
                                  At and From are synonymous.
                             +L - The bias is added to the left of the breakpoint
                             #w - width of the section on which the DC-level is calculated
                                  Default is one sample.
                          #l,#r - offset from the breakpoint (number of samples)
                                  to the left and the right, respectively. Default = 0,0                                                                    

 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

 TRIG[GER] #v #d [R=#w] [ time-range ]             
                                - Create a series in W consisting of const. w in triggered intervals
                             #v - tigger level
                             #d - +1 or -1 for sign of slope
                             #w - default = 1.0
                                  see Example 23 (trigger at sub-zero temperatures)       

 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-#i                     - truncate #i samples from the end.
 TRUNC #n                       - should also work like N=#n
EXTEND[W] #n #v
XTEND[W] #n [V={DC|LAST|MRS|#v}]
 EXTEND[W] by #m
[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 MRS
                                  or repeat the last value or
mark as missing, the latter by default.  
                                  Case DC: 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.
                       +I       - requests unpadding in ZSPF and AMDEMOD (subroutines
                                  spectralzspf.f and amdemod.f)

 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.
                                  Except O= the options are the same as for PAD

                           O=#l - integer l: Add l additional samples to avoid contamination from an eventual
                                  jump at the the start of the pad.
O=FH - Add the length of the positive half of the recently saved FILTER
O=F  - Add the total off-zero length of the recently saved FILTER.  

                                  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.

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

 RSYNCUNPAD +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).
RSYNCUNPAD [#n]                  However, if CPAD is called repeatedly in a loop, RSYNCUNPAD +A
                                  becomes uncertain. Specify the sample number n manually.
                                  Perhaps there's demand for obtaining n from an expression
RSYNCUNPAD 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 (Sep. 2015) 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=9 of 20 possible.

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
LOOP is the index to the data array before incrementation. It is initialized with the LOOP command
NMISS is set by commands REPDC and REPMISS
is reset to the last index of the primary series after every command

RESET# var = #val
can be used to set an internal variable.
An earlier REPMISS command for instance might be irrelevant; in that case, use
or issue a new REPMISS command 

Many extensions can be imagined. The procedure is simple:
In tsfedit.f, define a name in NHASHV(k), a default value in IHASHV(k), a format in FHASHV(k)
and increment KHASHV. Store the variable of your choice in
However, remember that only integers can be handled. Floating-point values can be handled like
those under STATISTICS.  

For quoting these numbers in the instruction block,
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)

Formats can be changed:

Thus, the loop index is formatted to yield `  1´ ... `999´ until it overflows.
Sometimes you'll need i3.3, e.g. when using the index in file names.

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 rgb(0, 153, 0);"> 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
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 000
41 < o/rt.pef


Example (7)
; samples to a sequence of files, each time shifting the
; sampling instance one time step ahead, 10 times
TRMARG  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#)

Examples (12)

Example (12.1)

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.
The ZSPF command is called with the band-limit option (+BL), which applies the spectrum of the saved filter
(FILTER ... SAVE) on the spectrum of the time-series segment.
The domain is the usual discrete domain (+DISC), the z-polynomial is trivial, so we simply use z0 (POWER:0).
We'll preserve the DC-level (+D)

; 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 (12.2)
Like before, but with a different way to preserve the DC-level
It's assumed the input is augmented with
tslist-app <opt> ++ -Efastlpf.tse,LPF + <files ...>
xxso, e.g. doing decimation, a final set of files can be concatenated.
FILTER D:${FWD:[WD KAIS 401 2.4 0.0 0.002 0.5 0.0]} SAVE
LOOP 78 <122>
DO From # ###INDEX# For #60000
ENDLOOP <122> (N> ###INDEX#)
${DECIMATE:[; ]}
; e.g. setenv DECIMATE "DECISYNC 60 100 0"
$ setenv DECIMATE "DECISYNC 3600 100 0"
$ tslist-app -I -o tmp/glpf.ts ++ -E fastlpf.tse,L + ~/TD/d/G1_garb_18070?

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=B 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 `%´

Example (19)

IIR filter as a nonlinear dynamic system: We create nonlinear tides

tslist o/ -C3 -E modulate.tse,N -o tmp/modulated.ts -I

with modulate.tse
FILTER 0 1 ' '
1.d0 -1.d0
IIR 1 NLIN=0.0,0.001
1. 0.9
We send a theoretical tide into the system, differentiate and integrate. The nonlinear parameter Xi causes higher influence of the input when the output is high and vice-versa.
The system is seen to create a sequence of harmonic overtones. When both non-linearity parameters are zero, the output is identical to the input except for rounding errors.

Example (20)
Cure an end-effect, e.g. the last 15 samples are replaced with the one before.
SAMPLE +Q ->1 At #N-16
POKE %1 From  #N-15 To #N
Example (21)
A smart thing, listing the ts-file names in an instruction file, and replacing environment variables in passing.
Here:    Ttide/SCG/urtap-openend.ins
which can be sourced in tcsh since it has a header with setenv commands and an exit command.
source urtap-openend.ins
grep -e '^.. [<B^R] .*\.ts$' -e '^.. [<B^R] .*\.mc' urtap-openend.ins | awk '{print "NOOP",$3}' > ! tmp.tmp

tslist _ -B -r1 -E `wraptse L tmp.tmp` | awk '/NOOP/{print $3}'

Example (22)
Using the  XCORRU  command:
tslist $seisf -I -Eanycorr.tse,CE -D -O:`label SEIS,ZACC` $sfn
set cmd = ( `fgrep tslist $CORRF` -o  )
setenv CEOUT o/ag-in-synch.ts
where  $seisf  is a seismogram file (Z, acceleration, 0.2 sampling time),
$sfn (= o/ ) is the output in synch (with the lag detected)
with  $AGF (= d/scg-cal-201502p10.dtr.ra.ts ), a residual drop series from absolute gravity
The cross-correlation results are printed in ASCII to  $CORRF (= o/seis-Z-ag-p18.corr ).
A subset will be correlated if  $FROM  and  $FOR  are set. In this example they were
$FROM = 10800
$FOR = 10800
OPEN 31 < ${AGF}
DO From # ${FROM:[0]} For # ${FOR:[###NDATA#]}
31, 'BIN',-99999.0,'${AGLAB:[U]}',0, 0, 'K', 1.0d0
where  $CORRL = 200.  Note especially  KU=31 , which results in the following lines in $CORRF:
tail $CORRF
Lag  +195 <Corr2U>>> CORR: -0.02484, CXY CXX CYY= -1.8525E+02  4.3185E+01  1.7271E+02, ndata=  597

Lag  +196 <Corr2U>>> CORR: -0.03625, CXY CXX CYY= -2.6894E+02  4.2952E+01  1.7271E+02, ndata=  597
Lag  +197 <Corr2U>>> CORR: -0.04520, CXY CXX CYY= -3.3315E+02  4.2676E+01  1.7271E+02, ndata=  597
Lag  +198 <Corr2U>>> CORR: -0.04984, CXY CXX CYY= -3.6478E+02  4.2378E+01  1.7271E+02, ndata=  597
Lag  +199 <Corr2U>>> CORR: -0.04890, CXY CXX CYY= -3.5547E+02  4.2087E+01  1.7271E+02, ndata=  597
Lag  +200 <Corr2U>>> CORR: -0.04211, CXY CXX CYY= -3.0444E+02  4.1863E+01  1.7271E+02, ndata=  597
tslist d/scg-cal-201502p10.dtr.ra.ts -I -Eo/seis-Z-ag-p18.corr,OUTCE
OUTTS U=31 < ${CEOUT:[tmp.ts]} From #2160 To #4189
KU=31 assists the subroutine (~/sas/p/corr2uncp.f) to find the file name.
The data sets in synch are then:
Seismograph:    -O:`label SEIS,ZACC` $sfn
Gravimeter:         CEOUT = o/ag-in-synch.ts
Finish up with
tslist $CEOUT -I -O:`label AG,VAL` $sfn
and list with
tslist $CEOUT -LS -LA -MQ -C3 -qqq

Example (23)
TRIG command:
cd ~/wx/ARC
rm -f
tslist ORTH2016.ts -I -E trigger.tse,TRIGG -O:`label TEMP,OSO`
tslist -LTE -LTR],Rwd -C3                  # 1)
tslist -LTR -LTE -C3                       # 2)
tslist -LTR -LTE -C3 -Ft36,f10.0,t36,f10.2 # 3)
1) Temperatures first, MRS will show as `-1.0000E+05´
2) prints only at subzero temperatures, two data columns
3) overprints the trigger-on signal with temperature

; create tigger-signal in workspace
; save in labeled MC-file
TRIG 0.0 -1 R=0.0

; delete the original series outside trigger-on
; never use a time-range here!
TRIG 0.0 -1 R=1.0
Example (24)
A hash-variable, here ###LOOP#, together with environment definitions,
gaining from the order in which replacement is carried out: first hash-variables, then environment parameters.
setenv PARAM002 83
setenv PARAM003 73
setenv PARAM001 90

tslist _16348 -B2017,1,1 -r1. -E manymemsp.tse,M -I
   with manymemsp.tse
OPEN 31 < o/manymemsp.sp
LOOP 3 <101>
PEF U=41 F=o/tmp.wr-1h.pef L=${PARAM###LOOP#}
  A spectrum can be retrieved with e.g.
splist o/manymemsp.sp X X
splist -n/92 o/manymemsp.sp RSP

Example (25)
Compute the complex spectrum of a Wiener filter, save it to binary file:
setenv WFLEN 48
setenv WFFILE o/
tslist _`calc "2*$WFLEN+2**15"` -r1. -B2017,1,1 -E wf.tse,B -Nf16.9 -I
  with  wf.tse
ADD 1.0 At # I`2*${WFLEN:[048]}`
GETFS U=31 < ${WFFILE:[o/]}] M=WF_00${WFLEN:[048]})
OPEN 32 < tmp/wf.csp
ZSPF POWER:0 [F]=cpd +BL +D -P [dB] -SPP CSPO=32:WF${WFLEN:[048]}
 Note that sasm06 saved the filter in a binary file (see ~/Ttide/SCG/sasm06-ra-ochy.ins ). Use fscat to show file content.
  Print the phase spectrum:
splist -cpd -Pd -n/048 tmp/wf.csp CSP

Setting environment variables for a specific block

For obtaining assistance from a tcsh source-utility /home/hgs/bin/setenv4tse (for comfort, define an alias), the command block may contain lines like

  setenv PDGFROM   "$<"    # starting index of time-series
  setenv PDGDUR    "$<"    # duration
  setenv PDGDB     "$<"    # dB if you want decibel
  setenv PDGCAL    "$<"    # ;  if you want to skip
  setenv PDGCALAMP "$<"    # calibration sinusoid calibration amplitude
  setenv PDGCALFRQ "$<"    #... frequency [mHz]

each having a blank first column. setenv4tse converts this to a temporary tcsh-script and executes it.
The user is prompted with the text after the comment mark ` to set values.
(In responses, sets of values or strings containing blankspace must not be enclosed with quotes. Instead, the setenv-lines in the tse-file must enclose `$<´ in double-quotes if such strings are allowed.)
Since the GETENV subroutine in Absoft F77 does not distinguish between empty or undefined environment parameters, empty answers will engage default values.
Default values in command lines are coded as follows:


The assistent is called as follows
source ~/bin/setenv4tse tse-file,target

Main programs that call tsf_edit :

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

Some subroutines at least: