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,*)
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
continued.
Instructions comprise two types, file processing and in-core
processing.
setenv DOTHIS 0 alt. setenv DOTHIS 1
TSF EDIT target* Targets must not contain blanks.
instruction
IF ${DOTHIS:[0]}
instruction
...
ENDIF
...
END
TSF EDIT targetRun - or examine - the command using script tsex . For usage, issue
;. command line
;. ...
instruction
instruction
...
END
* Long lines can be continued on the next line if indicated by a trailing ' /'
Setting environment variables that are significant for a command block:
¤String variables:
Assigned with E`expression`
->¤name in
NOOP (and planned: a few other) instructions,
and the result when an F` ,
I` or R`-expression`
->¤name
is evaluated (one per line; OBS!
not D` yet).
¤name is evaluated
at instances of ¤name
anywhere.
%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.
Examples:
NOOP %AVE:f5.1 -
if this yields overflow, procedure will fall back to the default.
NOOP %AVE:(f10.1)
- given in with parenthesis blankspace will be removed.
OPEN 21 < file-%AVE:(f10.1).dat
- will tightly embed the
number in the file name.
NOOP %AVE:(R)
- resets to the default format.
Undefined variables' codes are
left unchanged and so are %% %- %_ a.m.o.
Also, some arithmetic commands
store a value in a variable that can be retrieved with a simple %
WMEAN ...
->S
Retrieval and storing is done in
subroutine Replace_Statvars (sas/p/tsstats.f).
POLYRED ...
+S->O&S
stores the estimates of offset and
slope so they can be recalled with %OFF
and %SLO, respectively.
A time-range expression
(also in a DO instruction) stores the indexes of
begin and end in
%TRB and %TRE,
respectively. To quote them in a conversion to kalendar date
, e.g.
in wide,
blank-separated form, use D`%TRB`
in tight,
comma-separated form, use D`%TRB_:3,`
the latter in order to
defer the format option away from the parsing as a statistics
variable.
The date at the begin, end or
centre of the current interval can also be inquired using
D`B` D`E` or D`C`,
respectively
The NOOP (or
| ) command can be used to set the numbered
parameters %1 .. %100
(the so-called sample array).
For instance, if you want to include numeric values into file
names,
consider
| %->%1
OPEN 21 B abcd-%1:(f10.1).dat
OPEN 21 B
abcd-%1:(I10.3).dat - conversion to integer
and padded with zeros
The sample symbols can be coded
with one, two or three decimal places, %1 %01 %001
Formatting numerical values of
kind % :
A default format for floating-point
results 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
FMT=code]
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 with
FMT=
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 form (which would make tsfedit.f
a bit shorter), the scripts need revision.
The subroutine print_tsf (tsfeditprt.o) has been added to the code.
#Hash-variables:
It's not a hash structure, it's
only 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:would generate a file named file+015.ts
SHIFT S=15
OPEN 51 < file###SHIFTS#.ts
DUMP U=51
OPEN 51 < file###.tswould generate two files, file+001.ts and file+002.ts
DUMP U=51
FILTER ...
OPEN 52 < file###.ts
DUMP U=52
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.
Example:
NOOP F`%RMS *sin(pi/4)`Expressions evaluate to floating-point numbers if preceded by F` , integer if preceded by I` , (no format option with `I),
NOOP F`%RMS *sin(pi/4) :%10.4f`
MEDIANThe 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).
STATISTICS +S
NOOP Median - Average = F`calc ( %MED ) - ( %AVE ):%10.4f`
COMMENT [U=#u] <<targetThe text lines are written to the protocol (STDOUT) by default, or to a file opened on unit u. This way you can produce table heads
text text
text including ¤ $ % and # variables
target
;:TARGET: commandthen in a csh-script, the following line will execute the command:
eval `awk '/:TARGET:/{sub(/.:.*:/,"",$0);print}' file.tse`Executable shell script code inside a TSF EDIT block
e.g. when a Unix command line's option requires a value consistent with a
specific section in a tse-file, like command -shift $shft
eval `awk '/:2M20HZ:/{sub(/.:.*:/,"",$0);print}' /Seismo/gcf/pdgram.tse`
where
;:1HW: set shft=3600
;:3000HW: set shft=15000
;:2M10HZ: set shft=120
;:10M10HZ: set shft=600
;:5M10HZ: set shft=300
;:ANY: set shft=$PDGLENGTH
appear in the tse-file.
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 x.in designate the time series
currently in core, w.in 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 (options
-# 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] |
|
Binary 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- |
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,
DF:date_format
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.
value - a numerical value, purpose
depends on context, mostly multiplication
of
the input series
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). With !, replace existing
samples
N* - nothing, multiply with value and add data
as is NO
NK -
like N, but adds data of non-aligned series
in synchronicity; no MRS symbols are introduced,
and series X is not cropped.
Directive +KL is an alternative.
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
) )
where 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.
begin and end are indexes of w.
Option +ADJ requests adjustment of the gap filling's
DC-level to the simultaneous part of x
Defaults for scale and bias are not available;
record is read with
READ (INSTR, *) SCALE,
BIAS, STRATEGY
The following operation takes place
...........
x.out = ( if valid x.in and not Force ) x.in
( else if valid w.in )
w.in * value * scale + bias
value - ampl.factor sc0,
used with cdir = 'A*' or 'N*'
or 'L*' or 'S*' or 'K*' or 'T*' or '+*'
| Option |
compatible
with cdir directives |
Function |
Check! Remarks
|
| +O
| O+ |
A* | - instead of A! , the replacement option can be issued this way. | |
| +SP
| SP+ SP=#w
|
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- |
(any) |
- Remove the DC-level of the time-series W | |
| -DC | DC- | L* | - the following operation
takes place: . X.out = sc2·X.in + value·sc1·(W.in - DCLvl + bias1) + bias2 |
|
| -DC | DC- | +* | - the following operation takes
place X.out = X.out +
sc2·X.in + value·sc1·(W.in
- 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 symbols). 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 | |
| M>0 | - Fill up with zeros beyond file end | deprecated |
|
| *X
| X* *W | W* |
N* L* S* |
- multiply X.in with W.in: X.out = sc2·X.in·(W.in [-DC-level] + value) (OBS! +value) |
Protocol DEBUG |
| X/
| /W |
N* L* S* |
- divide X.in by W.in: X.out = sc2·X.in/(W.in [-DC-level] + value) (OBS! +value) |
Protocol DEBUG |
| /
|
N* G* S* L* |
- divide by scale factor (/value) divide by scale factor (/sc0) divide by scale factor (/sc1) |
Protocol DEBUG |
| /. |
L* |
- divide by sc0, multiply with sc1 | |
| ./ |
L* |
- multiply with sc0, divide by sc1 | |
| // |
L* |
- divide by sc0, divide by sc1 | |
| +++ | N* |
- Stack this
data series onto X cut to the length of X
(at return the series will be normalised by the stack height at each point). Stack with equal weight, no bias. (If you need to rescale and bias, use cdir=`+*´ |
W) |
| +CY
| CY+ |
+* |
- cyclically: W is wrapped
over X |
W) |
| +++CY | N* |
- stack cyclically |
W) |
| +KL |
N* S* L* |
- keep length of series X,
and it's origin time; cdir=`{N|S|L}K´ would do the same. |
|
| !T=T |
(any) |
- 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. | |
| !R[T][r]=#v
|
(any) |
- 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. |
|
| xxxxxxxxxxxxxx |
________________________________________________________________________
Protocol: Read carefully, eventually use the DEBUG command! W) - Works with other cdir's ? |
In-core
processing:
Time-range
| 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 | Around:#p,#q #YYYY #MM #DD #HH #SS #FF | ALL } OBS! The keywords From, At, etc. must be coded with initial upper-case followed by lower-case (case conversion was dismissed 2020-02-24). An empty time range means (usually) all. If the command has no other parameters than a time range, the At codeword can be dropped. Example: DEL 2001 1 1 0 0 1 0 Around:#p,#q #YYYY... means that the range extends from p samples to q samples relative to 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. Example: DEL To YMD+14 12 0 0 0 0 By sample index { [ From ##j ] [ { To | For } ##k ] | [ At ##j ] | Around:#p,#q ##j } To ##k is equivalent to For ##(k-j+1) Example: PERIODOGRAM /N From #740 For #4230 Introducing a search margin (or a shift): TRMARG #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 A sample time or index can be stored in an internal variable if a single index is specified (At). Two variants exist: YMD->¤variable The full date, year ... milliseconds, is stored as a character string MJD->%n The modified Julian date (floating point) is stored in numeric variable %n The time range can be printed out using NOOP +S[ M=comment]] |
| OBS! In order to avoid parsing errors, right-end delimit the time-range expression with ` ; ´ Example DEL At 2020 02 26 16 33 00 000 ; YMD->¤DELDATE (YMD is not intended here as an instruction for epoch reference). N.B. range by dates: After From, To 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 [+R] [+EP]
DO will not be evaluated. There is no stack of ranges (yet)
The purpose of this command is to avoid repetitive time-range
specifications.
+R - To keep the end of the time-range of DO beyond DONE;
else, the end is set to the one
at subroutine entry.
+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 f Is reset to 1.0 at END.
END RESCALE
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
These variables can be accessed in other commands:
###TRMARG#
###TRMARG2# - will be
replaced by the value of m1 TRMARG, and m2
TRMARG2,
respectively, the "time-range margins"
###SHIFTS# - the
accumulated value of origin shift a,
from SHIFT S=s : a = sum s
or SHIFT T=s : a = sum s/dt
Applies to time-range expressions and SAMPLE RESAMP
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# = i
condition - N> ###INDEX#
- continue loop while the index variable is within
the scope of the time series
Thus,
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 ]
|
with an almost synonymous `| ´ that
suppresses the printing
of the NOOP command line <TsfEsi>C> with expressions
resolved;
the <TsfEdi>N> stripped of the command is
printed anyway.
A test and debug feature. Expressions are evaluated and
start and end indices of the time-range are reported;
however, nothing is executed. Can be used to assign
assign values to internal variables.
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 bQ´
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 openfs
### - File name may contain `###´, which will be replaced
by a running signed three-digit number. One occurance
per file name only!TR
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 or CONTNOECHO or CONT_NOECHO is also possible
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=#u [L=label]] [O:opt] [ time-range
] Write (binary) the
time series X or
DUMPW U=#u [L=label]]
[O:opt]
the work array W to the file
on unit u.
opt - Options:
+E - Strip off leading Missing-records and adjust epoch
data written to
the file header.
OUTTS
[W] U=#u < filename] [O:opt]
[ time-range ]
Open the file, write the time series (binary),
and close.
OUTTS #u [O:opt]
[ time-range ]
OUTTSL #u L:label [O:opt]
[ time-range ]
Short forms to write series X to an already opened file.
OUTTSL
[W] U=#u < filename] L:label
[O:opt] [ time-range
]
Open the file, write labeled multi-column file,
keep open.
OUTTSL
[W] U=#u L:label [O:opt]
[ time-range ]
Append another column, keep open.
The part `U=#ub<bfilename´
has a strict format,
one blank inbetween. Acceptable options besides < are > % +
see openf.f
. Note that the file name is
right-delimited by `]´
File
name
may
contain
`###´,
which
will
be
replaced
by a running three-digit number.
L:label - can be coded PART1,PART2
or
PART1_PART2 or
PART1______PART2
W - Output the work array.
N.B.:
The time-range as well as the origin of W are assumed
equal to the master time-series.
opt - Options:
+E - Strip off leading Missing-records and adjust epoch
data written to
the file header.
time-range - not
valid with W.
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
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
PRINTW U=#u [F=fmt] [time-range]
Print time series W ...
PRINTXW U=#u [F=fmt] [time-range]
Print time series X and W
(supposed synchronous; MRS suppressed) ...
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.
set 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´
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
file.mc -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
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 Output consists of (default order): Date-information :
Value(s) : COMMAND [::
commentary ]
+options
- Common options, a string free from blanks. We
shall refer to
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
+V
- Switches debug option on. Evaluates the same control
+XC
- Date information is the centre of the time range:
+Y - simple format: Decimal
year and Value, nothing else.
+C:text
- replaces COMMAND with text .
This
option must be given last!
+F - Use the numerical format (optionally
specified with FMT=format-code] )
+P[:d]
- If the command line contains d
text , this text is added as a commentary
SET
- Parameters and options are defined, but the task is
not executed.
RESET - Like SET, but also carrying out
the task of the command. PRTSF [SET|RESET] [U=iun]
N=n [+options]
: %statvar_1
... %statvar_n
An easy way to inspect what's going on is EXAMPLES: First tse-file is produced by ~/TD/a/bin/tp4sets
: setcmds.tse The output shall be: $COMMAND time-range
the-results-from-averaging (REPDC) We could make the next step a little nicer, like
reporting to a file again, but we'll (2) Look at TD/a/bin/tp4sets , which exercises a.o. $ cat tmp/$c-$what.$type.dat (3) STATISTICS and
PRTSF $ tslist o/scg-cal-merged-O-sdr-expf.ra.ts -E
statcmds.tse,C -I Explanations: +C:string
replaces the command name normally been printed.
Underscore is replaced with blankspace.
|
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
W->X
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
WCDEL 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.)
Example:
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, first character / for
divide or * for multiply.
Optionally add
2 for squaring the second series before the operation,
Q for square-root,
L for Let (don't synchronize),
M for Mending, replacing missing records in the second
series
with its DC-level;
0 for removing a DC-level from the first series at
return.
Here, the second series is treated as weights on the first,
so the operation is DC-level = Sum X(i)*W(i)/Sum W(i)
Default operation is multiplication (`*´)
+P - Talk option = .true. in mult_ts_ts.
XWMIN
- Retain in X the minimum or ...
XWMAX
... 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[->O&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 0...nr
+S - Print slope results per hour and per year
+S->O&S
and optionally store the result in parameters %OFF
and %SLO
+S->#s or in sample array element s.
POLYREDJ 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.in - x)/s)·s
I - Keep integer result.
MODT #s #x [K] [ time-range
]
Transform data by modulus
X.out <- MOD(X.in - 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)X.in^j
, j=0..#n )
CLIP #a
#b [ time-range
].........-.Take
min(max(X.in, a),b)
ABS
[ time-range ].........-.Take
absolute value of X.in
X**#n [ time-range
].........-.Take
power n of X.in
EXPTR
[ time-range ].........-.Take
exp of X.in
LOGTR
[ time-range ].........-.Take
log10 of X.in (inject MRS if
negative).
DB [ R=#v ] [
time-range ].........-.Convert
X.in/v to dB (inject MRS if
negative). Default v=1.
MAX [ > #v ] [ time-range
].........-.Take
max(v, X.in).
Default v=0.
SQRT
[ time-range ].........-.Take
square root of X.in (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
.......................................Cf.downws.f
DCLOCK #n #t0' #dt' [opt]
[ time-range ]
- Interpolation over time simulating a clock drift using
n - a polynomial of order n. The time series is
assumed in W,
the resulting one is returned in X.
t0' dt' - The start time is t0', slightly
different from the nominal
value, and dt' is the slightly lengthened sampling
interval
(a positive value for a slow clock).
opt = [s/h] - alternative parametrisation:
the start time is t0+t0'/3600 h,
the sampling interval is (1+dt'/3600)dt
In practice, read the signal into X, then
X->W
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# (IHASHV(3))
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
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.
FILTERCOMMAND options [+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.
FILTERCOMMAND
options [+UG] STORE[->#n]
- Like SAVE, however specifically in storage
location n
[n=1]
Allows 5 filters of max length 40,000 or 1
filter of length 200,000
FILTERCOMMAND [+C]
[+UG] APPLY[<-#n]
[HIPASS]
- Apply a saved ¨(n=0) or stored (n>0)
filter
[n=0]
FILTERCOMMAND
[+C] [+UG]
FETCH[<-#n]
- Fetch a saved filter
[n=1]
FILTERCOMMAND [+C]
[+UG]
UPD[<-#n]
[HIPASS]
- 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.
FILTERCOMMAND may be
FILTER FILTS PEF AR IIR
WIENER BURG
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
INPLACE anyother
preserves the signal (switches the option off).
See example 11 near the bottom of this page.
To convert a filter into a time-series, use
GETFS
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
or
fscat 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.
*-*-*-*-*-*-*-*-*
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
WDZF[,P][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
#p - filter order (e.g. 8)
#fc - corner frequency
{mHz|Hz|cpd}
- units for #fc, Default is 1/dt
P=#p2,#p3... - additional design
parameters (not used); #p1 is identical to #fc.
(not tested yet)
FILTER #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
f(-1)..f(-n)
'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
f(0),f(1),...,f(n)
- autoregressive (infinite impulse response) filter
with #n+1
coefficients.
IIR -1
-
uses
the
coefficients
that
have been input under FILTER or GETFS
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.
E.g.
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
X**2.
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.
O=P
- under PEF, convert to a prediction filter already at the
reading instance.
FILTHR [S] #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=]#l [O=opt] [U=#u] [ time range
]
Creates a prediction-error filter bank, max order l. You need
the PEF command to apply a filter from the bank. Writes the
filter bank to a binary file and closes it thereafter.
U=#u - if given, the output file on unit u
is expected open.
opt - Options:
[ #u
< filename ]
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. Also statvars %... and
envars ${...} will be resolved.
S - scale filter for
unity power output.
Q - quiet, though minimal printing is to be expected.
GAPFPEF [S=#seed]
[P=#pwr] [D=#dff] [SPC=#crit,#wdth[,#mxit]]
[{+F|>X}] [ time range
]
After urtap analysis:
Reconstruct a gap-free residual (.ra.ts)
from the whitened residual
(.wr.ts) using the PEF and difference
pre-filter (coefficient -dff
(minus!) of the job (urtap*.ins).
The PEF must have been retrieved with PREDF ... SAVE
before.
See a commented example
and supplemental
material.
Specification
of a time range is discouraged until further notice.
S=
seed - (integer) a seed for the Gaussian variates
[-1]
P= pwr - the
prediction error power
[taken from the PEF]
The
RMS-dev of the gap filling (seen in the protocol)
should
match the RMS-dev of the gappy .ra.ts series.
D= dff - the
negative of the difference filter's
coefficient
at lag 1
[0.0]
Despiking:
If crit >
0
[-1.0 = don't]
Is recommended as there may be
spikes introduced on the
sides of a gap. Despiking will
delete the neighbouring samples.
A statistics of spike
indicators is printed to the protocol.
SPC= crit - after the
prediction stage and before the diff,
abs double difference > crit indentifies a spike
wdth - half-width of
interval around the spike where
samples
are invalidated
(integer)
[1]
mxit - maximum number
of attempts to despike and reprocess [10]
+F - the series of gap fillings is
returned in X
[kept in Wn+i]
the restored series in W.
>X - the restored series is returned in X,
the difference-filtered in W.
Neither +F nor >X - X
will contain the Gauss variates in the gaps;
for all else, X
is unchanged.
+F and >X - >X takes
precedence.
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.' ]
MEDIANFI[LTER] APPLY
WINDOW[:v] type
[P=#p] [T=#t] [C@0] [U] [+TR] [{
P | SAVE }] [
time range ]
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 ] [ time range ]
- 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
- comma-separated, recognized in movxcorr.f
XG to output sum X*W/sum W2,
YG to output sum X*W/sum X2
else sum X*W/sqrt(sum X2 sum W2)
Returns the result in X.
-DC to remove DC-levels in each segment
-AI to ignore an incomplete last segment
- recognized in tsfedit:
-IGNT0 to ignore a time offset between the series.
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
where
s(k) = Average<|x(i+2k-1)-x(i)|#q>
+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
FFTRC [R->X] [OSCP=#iun[,comment]]
[time-range] -
Fast-Fourier forward (IMSL DF2TRF) of X -> W(0,...)
OCSP= - Write complex spectrum using subroutine outsp to
pre-connected
unit iun labelled with comment.
R->X - The real part of the
spectrum is copied back to X .
MEMSP #n [R:u] [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.
R:N
- If the sampling interval is to be kept
(e.g. when the MEMSP command is in a loop);
else the sampling interval will be converted
into a
frequency interval, see R:u below. An example is in
~/Ttide/SCG/burg.tse,BURGSEGM
PERIODOGRAM [R:u] { [L,][{R,||R||Q,|I,}][U[=#h,]][{/N|/SQRN}][,dB][,+DC][+R][,+CW0][,M<#m] | ,
}
[/Hz] [/MISS] [OPSP=#u[,name]]
[ time range ]
R:u - Report Nyquist
parameters in units of:
R:s - Hz and seconds
R:h - cyc/h and hours
R:d - cyc/day and days
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)
Default is R:S (!).
"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 miss.data 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
, - use all defaults
Q - return the
power
[Amplitude]
R
|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
AMDEMOD[ULATE] {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
[don't]
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)kn )/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
RMS-computation.
MOVAVE #m #w [M=Q] [time-range]
- Compute a time series of running averages of X.in
where #m is a decimation factor and #w
the full(!)
width of the averaging window. If values are
missing in x.in, the output samples will be
located
in the "centre of gravity" of the valid samples in
each averaging interval. This feature differs from
MOVRMS.
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
x.out(i) = SUM
{j /. MOD(j,n)==i : x.in(j) }
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
counted from the start of the time range, use this
option.
Example:
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:):
F:text
F:text
B:text
B:text
A:text
...
EOL
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
L:END
EOL
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, nm/s2
+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
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
DECIMIN | DECIMAX | DECIMED
#n #s
[#u] [O=options]
Three kinds of decimations, taking the max or min or
median
value of the data segments.
n - 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≤i≤lt
+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] [O=opt] [ time-range
]
Resample with polynomial interpolation, using subsample_ts
Uses W as a work array. Normally (no or empty O=
specified),
the resampled series replaces X. A time-range
inside the
range of X will upset series X.
If re-use of X is intended, e.g. in a LOOP,
opt=K or W>iun must be given.
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,
or 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.
opt - K - keep
result in W
W>#iun - write to BIN-ts file on pre-connected
unit iun.
File will be closed.
RESAMPLE TT
The following records
constitute a time-tie, terminated with
END S[AMPLE]
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.
MONTHLYMEANS [ Time-range ] - Produce monthly means per calendar months
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.
%MIN %MAX %AVE %MOD %RAN %RMS %NWM %NVS %RVSTo invert the sign and/or taking the reciprocal value, prepend
- computed upon STATISTICS
%XWM %WST %WSM %NWM - computed upon WMEAN
%XSQ %WSQ %NSM - WSM = WSQ upon WMEAN O=+SQ
%MED - computed upon MEDIAN
%j (1 ≤ j ≤ 10) - set by commands that
support the ->#j option
DT #t (Re-)defines sampling interval t
; text Comment; text is completely ignored
;; text Comment; text will be printed to protocol
| 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 `/N^2´
EXP #a #p u [
Iref=#n
| Tref yyyy mm dd hh mm ss fff ] [ time-range
]
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.
SUB {#v |
=DC }[
time-range ] substitute v for missing samples
POKE {#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]
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=#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.
GAPFINT [P=#m] [
time-range ] -
Fill gaps in X with W, using
interpolation if necessary
#m - polynomial order, default m=2
INTERPEF [S=#i]
[ time-range
] - After importing a PEF with SAVE, a gappy, noisy
series
is repaired using the PEF as a model.
#i - (integer) random seed, default -1. See Example 26
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 1040
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 x(j+w2)
Note, if w2=0 and w1<0,
the mean is taken from the already
despiked part.
DELETE's
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)) < v 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=#v
EXTEND[W] #n [V={DC|LAST|MRS|#v}]
EXTEND[W] by #m [V=...]
EXTEND[W] [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.
TRIM -MRS
- trim off MRS-records of the end.
MIXADJUST [O:[+R[+T]][+S]]
- (Highly special)
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 inside ZSPF and AMDEMOD after the action
is done
(subroutines spectralzspf.f and amdemod.f); UNPAD is redundant
then.
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).
UNPAD [+K]
Undo padding. +K will keep the length of the series
so that
the result of circular correlation / overlap management can be
inspected.
CUNPAD
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!)
instructions
END
instruction word = END (capitals!)
Try this:
TSF EDIT XX
LOOP 10 <102> #-#
ADD 1.0 From YMD+### 0 0 0 0 To YMD+#.# 0 0 30 0
ENDLOOP <102>
END
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.
TSF EDIT KIRU
CONT 31
31 B $TAPDIR/KIRU.tse
Q
The file $TAPDIR/KIRU.tse contains TSF EDIT commands; a
block END mark
is not needed. Notice that we can enter environment
parameters.
TSF EDIT CTH grav
22,'BRX 9999]',-99999.d0,'F6.0',0,0,'G',1.0d0
-1.4462,1973.0,'F]'
22,'BRX 9999]',-99999.d0,'F6.0',0,0,'G',1.0d0
-1.4462,533.0,'F]'
22,'BRX 9999]',-99999.d0,'F6.0',0,0,'G',1.0d0
-1.4462,533.0,'F]'
22,'BRX 9999]',-99999.d0,'F6.0',0,0,'G',1.0d0
-1.4462,535.0,'F 17 999]'
23,'TOP]',-99999.d0,'(i2,3i3,f8.1)',1,0,'J',1.0d0
SCALE=-0.094878d0
24,'TOP]',-99999.d0,'(i2,3i3,f9.2)',1,0,'J',1.0d0
END
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
Q
TSF EDIT AIRP
FILTER U=42 F=../sas/r/${section}_wtr_ap.inf O=AB >-16 16 'A'>
END
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 P
for trace printing and N 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)
TSF EDIT BEGIN
REPAIR
SCALE .178856
OPEN 23 < $TAPDIR/METS-resid.mc
23,'BIN]',-99999.0,'L:RAD_VAL]_M>0',0,0,'N', -.065634 -DC
CLOSE 23
OPEN 23 < $TAPDIR/UMEA-resid.mc
23,'BIN]',-99999.0,'L:RAD_VAL]_M>0',0,0,'N', .346502 -DC
CLOSE 23
OPEN 23 < $TAPDIR/SUND-resid.mc
23,'BIN]',-99999.0,'L:RAD_VAL]_M>0',0,0,'N', .198229 -DC
CLOSE 23
OPEN 23 < $TAPDIR/MART-resid.mc
23,'BIN]',-99999.0,'L:RAD_VAL]_M>0',0,0,'N', .135893 -DC
CLOSE 23
OPEN 23 < $TAPDIR/LEKS-resid.mc
23,'BIN]',-99999.0,'L:RAD_VAL]_M>0',0,0,'N', -.016194 -DC
CLOSE 23
OPEN 23 < $TAPDIR/SVEG-resid.mc
23,'BIN]',-99999.0,'L:RAD_VAL]_M>0',0,0,'N', .113417 -DC
CLOSE 23
OPEN 23 < $TAPDIR/OSTE-resid.mc
23,'BIN]',-99999.0,'L:RAD_VAL]_M>0',0,0,'N', .068454 -DC
CLOSE 23
OPEN 23 < $TAPDIR/VILH-resid.mc
23,'BIN]',-99999.0,'L:RAD_VAL]_M>0',0,0,'N', .047797 -DC
CLOSE 23
OPEN 23 < $TAPDIR/ARJE-resid.mc
23,'BIN]',-99999.0,'L:RAD_VAL]_M>0',0,0,'N', .081275 -DC
CLOSE 23
OPEN 23 < $TAPDIR/KIRU-resid.mc
23,'BIN]',-99999.0,'L:RAD_VAL]_M>0',0,0,'N', .008231 -DC
CLOSE 23
REMDC
END
Example (5)
Signal conditioning (differencing), cleaning and information
TSF EDIT COND
FILTER 0 1 ' '
1. -1.
REMDC
MISS
RMSD
RMSD To 2000 02 24 0 0 0 0
RMSD From 2000 02 24 0 0 0 0
DEL > $cond_limit ALL
MISS
RMSD
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
$cond_cont
END
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 DEL From 1993 10 13 00 48 00 00 To 1993 10 13 00 48 DEL From 1993 10 12 16 48 00 00 To 1993 10 12 16 48 DEL From 1993 10 07 20 48 00 00 To 1993 10 07 20 48
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.
TSF EDIT BURG
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
PEF L=4
END
TSF EDIT SAMPLE
; 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
CLOSE 51
TRMARGSTEP 1
ENDLOOP <101>
END
TSF EDIT SV
; samples to a sequence of files, each time shifting the
; sampling instance one time step ahead, 201 times.
SHIFT S=-101
LOOP 201 <101>
SHIFT S=+1
OPEN 51 > svdata/sv-hpf-s###SHIFTS#.mc
DECI2W 500 0 From 2011 06 15 00 20 00 000 To 2011 06 15 12 09 05 000
DUMPW U=51 L=SV_____________${SVC}]
CLOSE 51
ENDLOOP <101>
END
Example (8)
The spectrum of a finite impuls response
filter, here
a prediction filter
firspect.tse:
TSF EDIT S
ADD 1.0 At #${PRFL}
PEF U=41 F=${PRF} L=${PRFL} O=R
PERIODOGRAM /N From #0 To #4095
END
with
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
prfspect.spa
(4300: we need a margin for the filter cutoff
effects)
Example (9)
Undo the phase effects of the Guralp
seismometer and compute acceleration
TSF EDIT 10
PAD
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 )
UNPAD
DECIMATE 10 0
FILTER 0 1 ' '
1.0 -1.0
SCALE 0.2
END
Example (10)
Apply the GWR Bessel-8 filter
TSF EDIT B
PAD
ZSPF BESSEL:8 FC=0.125 Hz -P
UNPAD
END
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
TSF EDIT J
SETINDEX 60000
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
ADDINDEX 60000
ENDLOOP <101> (N> ###INDEX#)
END
Example (12.1)
The input has 100 S/s and is loooong. We have to start filtering at offset 20 samples
Example (12.2) ~/Seismo/gcf/fastlpf.tseTSF EDIT SDECI100
; 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
SETINDEX 20
LOOP 78 <103>
DO From # ###INDEX# For #60000
CPAD O=2880
ZSPF POWER:0 +DISC +BL +D -P
CUNPAD
DONE
ADDINDEX 60000
ENDLOOP <103> (N> ###INDEX#)
RSYNCUNPAD 20
DECIMATE 100 0
END
TSF EDIT LPF
FILTER D:${FWD:[WD KAIS 401 2.4 0.0 0.002 0.5 0.0]} SAVE
REPDC ->S
REMDC
SETINDEX 0
LOOP 78 <122>
DO From # ###INDEX# For #60000
CPAD O=FH
ZSPF TRIVIAL +BL -P -SPP
CUNPAD
DONE
ADDINDEX 60000
ENDLOOP <122> (N> ###INDEX#)
RSYNCUNPAD
ADD %
${DECIMATE:[; ]}
; e.g. setenv DECIMATE "DECISYNC 60 100 0"
END
$ setenv DECIMATE "DECISYNC 3600 100 0"
$ tslist-app -I -o tmp/glpf.ts ++ -E fastlpf.tse,L + ~/TD/d/G1_garb_18070?-1s.mc
Shows the use of the OUTFS command with conversion from
time series to filter.
TSF EDIT FDECI100
; 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 ZSPF POWER:0 +DISC +BL +D -P
OPEN 31 < fdeci100.fs
DECIMATE 100 0
OUTFS U=31 X> +S From #0 To #72
END
TSF EDIT A2D
; acceleration to displacement on 3-s or 10-s data
; Scale must be set accordingly.
;
; 1. restore SCG's anti-aliasing filter
REMDC
OPEN 31 < tmp/dump.mc
CPAD O=1200 +I
ZSPF BESSEL:8 /F FC=0.0174 Hz +S -P DU=31
RSYNCUNPAD -1200
CUNPAD
;
; 2. Integrate by Spectral Z-filter
CPAD +I
FILTER D:WD KAIS 1200 2.0 0.002 0.5 0.5 0.5 SAVE
ZSPF POWER:-1 -P +BL DU=31
ZSPF POWER:-1 -P +BL DU=31
UNPAD
REMDC
; This is for mm from nm/s^2:
SCALE -1.e-6
END
Example (15)
~/TD/rms.tse,EQ & ~/TD/weq.tse,W
Delete measurements affected by earthquakes,
interpolate and provide weights for the interpolated
ordinates.
Input = 1s.mc-data, 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
rms.tse:
TSF EDIT EQ
; 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
DECIMATE 20 ${MOVE} 0
FILTER 0 1 ' '
1.0 -1.0
MOVRMS 1 10
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
END
weq.tse:
TSF EDIT W
; 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
DECIMATE 30 ${MOVEKWD30} 27
END
TSF EDIT WFTGGRetrieves and applies a Wiener filter, and lists it in ASCII on file filter.coeff
REPAIR L ALL
GETFS U=11 < o/tgwr-grwr-1h.wwf] M=WIENER INF: 32) O=B SAVE
FILTER APPLY
PRINTF F=I4,1p,e12.4 U=41 B filter.coeff]
END
tslist owf/g100201-130515-1h-lapww-M6.ra.ts -E amdemod.tse,AMS -Un23240 -I -o tst.tswith amdemod.tse:
TSF EDIT AMSNote the big loss of data because of the length of the low-pass filter.
FILTER D:WD KAIS 2880 2.0 0.0 0.008 1.0 CPAD 2880*DC O=2880
AMDEMOD +D +BL FD=1.9322736 FS=12.0
CUNPAD
RSYNCUNPAD
END
TSF EDIT SNote: MEDIAN must be called separately while RMS is computed in the STATISTICS command.
MEDIAN
STATISTICS +S
ADD -%MED
SCALE /%RMS
END
IIR filter as a nonlinear dynamic system: We create
nonlinear tides
with modulate.tsetslist o/g090701-140131-1h-boa-nk.pt.ts -C3 -E modulate.tse,N -o tmp/modulated.ts -I
TSF EDIT NONLINWe 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.
REPAIR L
FILTER 0 1 ' '
1.d0 -1.d0
IIR 1 NLIN=0.0,0.001
1. 0.9
END
TSF EDIT CURE_END_EFFECTExample (21)
SAMPLE +Q ->1 At #N-16
POKE %1 From #N-15 To #N
END
source urtap-openend.insExample (22)
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}'
d/g100202-OPNEND-1h.mc
d/PMG100202-OPNEND-1h.ts
d/barlap100202-OPNEND-1h.ts
d/bagl100202-OPNEND-1h.ts
d/bagn100202-OPNEND-1h.ts
d/bp-rin-oso.ra-1h.ts
d/gnt100202-OPNEND-1h.ts
tslist $seisf -I -Eanycorr.tse,CE -D -O:`label SEIS,ZACC` $sfnwhere $seisf is a seismogram file (Z, acceleration, 0.2 sampling time),
set cmd = ( `fgrep tslist $CORRF` -o )
setenv CEOUT o/ag-in-synch.ts
cmd
TSF EDIT CEwhere $CORRL = 200. Note especially KU=31 , which results in the following lines in $CORRF:
OPEN 31 < ${AGF}
OPEN 41 B ${CORRF}
DO From # ${FROM:[0]} For # ${FOR:[###NDATA#]}
31, 'BIN',-99999.0,'${AGLAB:[U]}',0, 0, 'K', 1.0d0
XCORRU R=-${CORRL},${CORRL} P=41 KU=31
DONE
CLOSE 31
END
tail $CORRFKU=31 assists the subroutine (~/sas/p/corr2uncp.f) to find the file name.
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
TSF EDIT OUTCE
OUTTS U=31 < ${CEOUT:[tmp.ts]} From #2160 To #4189
END
tslist $CEOUT -I -O:`label AG,VAL` $sfnand list with
tslist $CEOUT -LS -LA -MQ -C3 -qqqExample (23)
cd ~/wx/ARC1) Temperatures first, MRS will show as `-1.0000E+05´
rm -f trig.mc
tslist ORTH2016.ts -I -E trigger.tse,TRIGG -O:`label TEMP,OSO` trig.mc
tslist trig.mc -LTE -LTR],Rwd -C3 # 1)
tslist trig.mc -LTR -LTE -C3 # 2)
tslist trig.mc -LTR -LTE -C3 -Ft36,f10.0,t36,f10.2 # 3)
TSF EDIT TRIGGExample (24)
; create tigger-signal in workspace
; save in labeled MC-file
TRIG 0.0 -1 R=0.0
OUTTS W U=31 < trig.mc] L:TRIG
END
TSF EDIT TRIGX
; delete the original series outside trigger-on
; never use a time-range here!
TRIG 0.0 -1 R=1.0
X*W
END
setenv PARAM002 83with manymemsp.tse
setenv PARAM003 73
setenv PARAM001 90
tslist _16348 -B2017,1,1 -r1. -E manymemsp.tse,M -I
TSF EDIT MEMPSPA spectrum can be retrieved with e.g.
OPEN 31 < o/manymemsp.sp
LOOP 3 <101>
PEF U=41 F=o/tmp.wr-1h.pef L=${PARAM###LOOP#}
REWIND 41
MEMSP 16384 OMEMSP=31,ORDER${PARAM###LOOP#}
ENDLOOP <101>
END
splist o/manymemsp.sp X XExample (25)
splist -n/92 o/manymemsp.sp RSP
setenv WFLEN 48with wf.tse
setenv WFFILE o/ra-ecco1-lpt-1h.wf
tslist _`calc "2*$WFLEN+2**15"` -r1. -B2017,1,1 -E wf.tse,B -Nf16.9 -I
TSF EDIT BFSPNote that sasm06 saved the filter in a binary file (see ~/Ttide/SCG/sasm06-ra-ochy.ins ). Use fscat to show file content.
ADD 1.0 At # I`2*${WFLEN:[048]}`
GETFS U=31 < ${WFFILE:[o/bpra-gra-1h.wf]}] M=WF_00${WFLEN:[048]})
FILTER +UG UPD
OPEN 32 < tmp/wf.csp
ZSPF POWER:0 [F]=cpd +BL +D -P [dB] -SPP CSPO=32:WF${WFLEN:[048]}
END
splist -cpd -Pd -n/048 tmp/wf.csp CSP
TSF EDIT PEFREPAIR
REMDC
OPEN 41 ^ d/grav-1h.pef
PEF U=41 L=9 SAVE
INTERPEF S=-6
END
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:
${PARAMETER:[default]}
The assistent is called as follows
source ~/bin/setenv4tse tse-file,target
./Apload/m/apltsm.f
./Oload/smhi/m/apltsm.f
./sas/p/mg/elder/sasm03.f
./sas/p/mg/sasm01.f
./sas/p/mg/sasm03.f
./sas/p/mg/sasm06.f
./sas/p/mg/xyd.f
./sas/p/mt/pastescans.f
./sas/p/mt/tsl-x.f
./sas/p/mt/tslist.f
./sas/p/mt/twostepm.f
./sas/tslist.sub/main.f
./tap/m/old/urtap-old.f
./tap/m/urtap.f
./tap/m/urtaphw.f
./tap/t/m/urtapt.f
./sas/p/downws.f
./sas/p/tsapps.f
./tap/p/urtapu9.f
./tap/t/m/urtapt.f
.bye