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.
%Statistics variables:
Computed with the STATISTICS, MEDIAN and SAMPLE commands
and retrieved with %SYM symbols.
#Hash-variables:
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 51
OPEN 51 < file###.tswould generate two files, file+001.ts and file+002.ts
DUMP 51
FILTER ...
OPEN 52 < file###.ts
DUMP 52
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.
Instruction consists of the following parameters:
iun, target, recmrs, fmt, khms, itz, cdir, value [more options]
Example (to subtract a filed time-series from the one currently in memory):
22, 'BIN]',-99999.0,'U', 0, 0, 'N', -1.d0 DC-
iun...itz are standard parameters the meaning
of which is documented
also at other places. Special to tsfedit is cdir
and value
iun - file unit number.
target - max 8 characters: Position file
at "target";
to rewind before positioning,
specify: "target]
Rwd"
"TOP]"
- read ascii file from the
top.
"TOP] Rwd"
"BIN"
- binary data (getts family)
"BIN Rwd"
"SAC"
- seismological SAC files (read_sac)
"BRX"
- Bruxelles
format
"BRX Rwd"
also:
ETERNA format
the
target specifies an END-OF-FILE mark.
c.f. readbrx.f.
C.f. proc_f_on_rec
or read_fuf for
more options that can be
passed through target.
recmrs - Missing record symbolic value (in formatted data sets only).
fmt - Case ASCII data: The
time and data format in FORTRAN, for
retrieving
one data column.
Case Binary data, MultiComp, the label of the requested segment
"l:text" "L:text" "*:text"
Case Bruxelles: Sub-format for samples, e.g. 'F9.0'
or 'F10.0:ETERNA' for ETERNA files
other (ASCII): The data format in FORTRAN.
This argument is passed to read_fuf.Thus,
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.
cdir - char*2 -
meaning
needs more data / instructions
------------------------------------------------------------------------
...........A*
- multiply with value and append
NO
N*
- nothing, multiply with value and add data as is
NO
S*
- apply a scaling factor and
add 1 record:
scale factor
F*
- apply a filter and
add
1 record: filter params.
+ n values: filter coeffs.
L*
- bias and scale, and add
1 record: bias1,sc1,bias2,sc2
+*
- bias and scale, and stack
1 record:
bias1,sc1,bias2,sc2
G* - use as gap filling
data
1 record: gap parameters
J*
- use as jump
data
NO
...........K* - scale, (mask) and Keep in work array W
...........KT - scale, (mask), filter
(with the SAVE'd
one),
and Keep in work array W
T* - Filter (using the SAVE'd
filter) and keep in work array W
P* - Compute phase from arctan X/W
C* - add series
cyclically.
*t
- override dt with the one from the data file
The options N S F L
T and C request editing of the
internal time-series
by means of linear combination with the time-series input here.
One effect of the editing may be an undesired truncation to the
length of the most recently input time-series; specify
lower-case
characters n s f and l to keep the length
unchanged.
Scaling factor (S):
x <- ( x + ( w * value * sc ) )
where sc is the
read-in scale factor
Cyclical addition (C):
x.i <- ( x.i + ( w.j *
value ) ) for all i
with j
= mod(i,lw),
where lw is the length
of w.
(itz may be used for an
origin shift (really?))
Bias and scale: [bias1 [, sc1 [,
bias2 [, sc2]]]] [-DC]
x
<- ( x + ( w + bias1 )*sc1+
bias2 )*sc2 )
Default: 0.0, 1.0, 0.0, 1.0
Filter parameters: scale,nfminus,nfplus,symmetry
..............................Example:
1.0, 0, 16, 'S'
for a symmetric filter with 16 unique coefficients
i.e. 31 coefficients when expanded.
cdir='G*'
OBS! No interpolation if series aren't in
synch.
Gap parameters: scale,bias,strategy
strategy:
'[X{D|H}]A]'
to replace missing and append all.
'[X{D|H}]A [begin
end]' to replace
missing and append,
...................................................
.begin end are sample numbers
'[X{D|H}]R [begin
end]' replace missing,
no appending, default = all
'[X{D|H}]F [begin
end]' force over,
no appending,
" "
The options XD or XH request
read_fm_exactd('D') or ('H'), respectively.
Defaults for scale and bias are not available;
record is read with
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, used
with cdir = 'A*' or 'N*'
or 'L*' or 'S*' or 'K*' or 'T*' or '+*'
more
options pertaining to cdir='L*'
or '+*'
With cdir='L*' and the -DC option ,
the following operation takes place
.......... .
X.out = sc2·X.in + value·sc1·(W.in
- DCLvl + bias1) + bias2
...........
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
+&&
-
Use
this
data
series to mask the one in memory (inject missing record
symbols).
This masking (still) uses a simplified scheme. The rates of the
two series are
*assumed* equal.
+SQ
- Add this data series to the one in memory by squares.
+RR
-
Remove
this series from the one in memory by regression
M>0
- Fill up with zeros beyond file end, as much as available.
*X [-DC] - with
cdir='L*' multiply X.in with
W.in:
X.out = sc2·X.in·(W.in
[-DC-level] + value) (OBS! +value)
X/ [-DC] - with
cdir='L*' divide X.in by W.in:
X.out = sc2·X.in/(W.in
[-DC-level] + value) (OBS! +value)
/
- divide by scale factor (cdir
= 'N*' or 'G*')
/.
- divide with sc0,
multiply with sc1
./
- multiply with sc0,
divide by sc1
//
- divide by sc0,
divide by sc1
+++ - Stack this data
series (at return the series will be normalised by the
number of
stackings at each point. Stack with equal weight, no bias.
(If you need to rescale and bias, use cdir='+*'
+CY -
cyclically (affects only stacking so far): W is
wrapped over X
+++CY - stack cyclically.
!T=[X][u:]#v - move time origin by v, units of hours except
when
u=S - #v in seconds
u=D - #v in days
X - move time origin to
the one of series X and add #v)
!T={+|-}ST - Summer time adjustment for e.g. stacking
days synchronous with civil time.
Note that `/´ as the last character
on the line is interpreted as the
a
continuation mark. To specify the `/´ `/.´ or `//´ options last
on
the
line, add ` -´ or any string that is not an option.
In-core
processing:
Parameter symbols are shown in italic,
compulsory text in bold, optional text in
[brackets].
A #-sign before a symbol indicates that the system
expects a numerical value.
| Time
range definition
Subroutine
time_range
File ~/sas/p/timerange.f |
By dates Let time-range := { [ From #YYYY #MM #DD #HH #MM #SS #FFF ] [ To #YYYY #MM #DD #HH #MM #SS #FFF ] | At #YYYY #MM #DD #HH #SS #FF | Around:#p,#q #YYYY #MM #DD #HH #SS #FF | ALL } An empty time range means (usually) all. If the command has no other parameters than a time range, the At codeword can be dropped. Example: DEL 2001 1 1 0 0 1 0 Around:#p,#q #YYYY... means that the range extends from p samples to q samples centred on the given date. Example: For a symmetric interval specify Around:-30,30 Obs! No blanks in this expression. The year month and day can be replaced by YMD[{+|-}#n] to mean the day of the epoch ±, optionally, n days. Example: DEL To YMS+14 12 0 0 0 0 By sample numbers { [ From ##j ] [ { To | For } ##k ] | [ At ##j ] | Around:#p,#q ##j } To ##k is equivalent to For ##(k-j+1) Example: PERIODOGRAM /N From #740 For #4230 Sample numbers are in the range 0 .. NXR-1 Introducing a search margin (or a shift): TRMARG #m1 #m2 - m1 and m2 is in number of samples. The margins can be incremented: TRMARGSTEP [#ms] - add ms to both m1 and m2 COMMAND Time Range Example: apply operation to the sample immediately before 2011 07 04 12 00 00 00 TRMARG -1 -1 SCALE 9. At 2011 07 04 12 00 00 00 |
| N.B. range by dates: After 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]
-----------------------------------------------------------------------------------------------------------
DO Time-range
In a DO - DONE block time-range expressions beyond the
DONE
DO will not be evaluated. There is no stack of ranges (yet)
TRMARG #m1 #m2
Margin for time range, m1 and
m2 is
in number of samples.
The margins can be incremented:
TRMARGSTEP [#ms]
add ms to both m1 and m2
TRMARGSTEP
[#m1s #m2s]
add m1s to m1 and m2s to m2
LOOP #n [#s] [#a] <label>
Loop #n times
(step #s and add #a, all integers).
ENDLOOP <label>
A Label is compulsory, must be keyed with brackets
< >
DEBUG [END] [ L={#n|D}
]
Diagnostic printing. Toggled;
.......................................END
can be stated for clarity.
The message level in some ts-routines can be set with n
(0=debug, 1=much, 2=normal, 3=important-only, 5=quiet)
L=D for the level in
effect at call time.
ECHO
Echoes the commands and reports time-range
analysis to STDOUT. ECHO
is the default.
NOECHO
In the case of samples that become deleted
a report of the delete command and the number
of deleted samples will be issued.
...................................F I L E C O M M A N D S
OPEN
The next lines contain an Open-a-file
block.
Last line = `b
b 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!
Some special variables can be accessed:
###TRMARG# - will be
replaced by the value of TRMARG,
the "time-range margin"
###SHIFTS# - the
accumulated value of origin shift a,
from SHIFT S=s : a = sum s
or SHIFT T=s : a = sum s/dt
CONT[QUIET]
#n
next line will be an
Open-a-file block.
Instruction input will continue
from unit n (only in the scope
of the current TSFEDIT invocation).
Cannot be nested, just chained.
The QUIET version will suppress ECHO until
the end on unit #n.
CONT_QUIET or CONTNOECHO or CONT_NOECHO is also possible
CONT[QUIET]
s
String s is passed to openfs, specifying
file unit, option, and file name. Otherwise
like 'CONT #n' above.
CLOSE #n
Close file unit #n.
Any file may be closed as needed.
If the TSFEDIT instruction file unit is
closed, input will continue from the
file unit defined on the call list.
That one cannot be closed.
DUMP #u
[L=label]] [O:opt] [ time-range
] Write (binary) the
time series to the file on unit u.
DUMPW #u [L=label]]
[O:opt]
Write (binary) the
work array to the file on unit u.
opt - Options:
+E - Strip off leading Missing-records and adjust epoch
data written to
the file header.
OUTTS
[W] U=#u < filename] [O:opt]
[ time-range ]
Open the file, write the time series (binary),
and close.
OUTTS #u [O:opt]
[ time-range ]
OUTTSL #u L:label [O:opt]
[ time-range ]
Short forms to write series X to an already opened file.
OUTTSL
[W] U=#u < filename] L:label
[O:opt] [ time-range
]
Open the file, write labeled multi-column file,
keep open.
OUTTSL
[W] U=#u L:label [O:opt]
[ time-range ]
Append another column, keep open.
The part `U=#ub<bfilename´
has a strict format,
one blank inbetween. Acceptable options besides < are > % +
see openf.f
. Note that the file name is
right-delimited by `]´
File
name
may
contain
`###´,
which
will
be
replaced
by a running three-digit number.
L:label - can be coded PART1,PART2
or PART1_PART2
or PART1______PART2.
W - Output the work array.
N.B.:
The time-range as well as the origin of W are assumed
equal to the master time-series.
opt - Options:
+E - Strip off leading Missing-records and adjust epoch
data written to
the file header.
time-range - not
valid with W.
MISS [ time-range ]
Report missing samples
PRINT #u [F=fmt] [time-range]
Print time series X in free format or format fmt on unit u
PRINTW #u [F=fmt]
Print time series W ...
PRINTF #u [F=fmt]
Print filter
...
fmt
- One integer field followed by one real*8 field.
REPDC [opt] [ time-range
]
Report DC-level to protocol (unit 6)
opt +TY - include full date
information and decimal year
+Y - ... decimal year
+T - ... full date
[/]->S[-] - set scale
parameter (the the inverse value if /-S ),
invert the sign if appended by `-´, for
use in
a range of arithmetic commands.
Example for +TY:
REPDC +TY From
2009 0 187 13 35 18 00 To 2009 0 191
13 25 13 59
` <TsfEdi>>>
REPDC 2009 11 06 10 28 18 000 2009.860195 -2.844177E+01´
RMS [ time-range ] (ignored, obsolete: Report RMS )
RMSD [opt] [ time-range ] Report RMS about mean.
RMSDW [opt] [ time-range
] ... do that on the work array W.
However, here the time-range
RMSD W [opt] [ time-range
]
may be a pitfall (it relates to the origin and rate of
time-series X);
prefer to use array indices instead.
opt +TY - include full date
information and decimal year
+Y - ... decimal year
+T - ... full date
[/]->S[-] - set the scale
parameter to the result (or 1/result if /->S ),
invert the sign if appended by `-´, for use
in
a range of arithmetic
commands.
The reported dates are derived from the origin and
rate of the
respective time-series.
MEDIAN [s-opt]
[ time-range ]
Report median
s-opt - [/]->S[-] store the
result as scale value for use in
a range of arithmetic
commands.
WDR #t [C=text]]
[U=#u] [M[c]=#ml,#mu[,c]] [ -DC[=#dc]]
[ time-range ]
Write records for a given signal threshold t. By default the
records are written as tsfedit DEL commands..
C=text]
- text will replace DEL.
Write a record for each
group of samples x
with | x - dc | > t
If t<0, the
threshold is |t| times
the result of the most
recent RMSD command.
-DC - Determine dc, the series' DC-level,
else use dc=0
-DC=#dc - Calculate deviation around dc.
U=#u - Write to file unit u. Default is 6
M=#ml,#mu - Output time ranges with a lower and
upper margin of +ml +mu
of the sampling interval. E.g. -1,+1 to extend the deletion
window by one sample at both ends. Default
-0.5,+0.5
MH=#ml,#mu - Margin of +ml +mu
in units of hours ( M=#ml,#mu,h should also work)
MS=#ml,#mu - Margin of +ml +mu
in units of seconds ( M=#ml,#mu,s should also work)
Consider to use WDR on
the result of MOVRMS
(don't opt for -DC
then)
WDR MISS [C=text]]
[U=#u] [M[c]=#ml,#mu[,c]] [ time-range ]
Like WDR, writes DEL records for currently missing samples.
CUT #n
Cut away n samples from the end
RETAIN time-range
... with adjustment of starting time and length
X->W
Copy series to work space
W->X
... vice versa
X*W O=s [+P] [
time-range ]
- Uses mult_ts_ts (sasu16.f) to multiply the series, result in
X
W*X O=s [+P] [
time-range ]
- Uses mult_ts_ts (sasu16.f) to multiply the
series, result in W
s - Operator and options, / for divide, Q
for squaring the second
series before the operation, L for Let, 0 for
replacing missing
records with DC-level.
SUMTS[X|W] [[/]->S[-]]
[ time-range ]
Sum the series X or W,
respectively. Default: X
Set scale parameter (the reciprocal value if
`/S´ , append `-´
to reverse the sign) for use in a range of arithmetic
commands
REMDC [ time-range ]
Remove DC level. The DC-level of the time range is
subtracted within the time range.
ADJDC [ time-range ]
Remove DC level. The DC-level of the time range is
subtracted within the full scope of the time series.
GETDC [ time-range
]
Saves DC level internally. See ADD,
SUB, etc. further below.
DETREND [P|I] [ time-range ]
Detrend (slope and bias)
P - Poor man's detrend
(defined by start and end point)
I - Inverse (unweighted)
using DETREND_UNW
POLYRED #n [O=[P][I][A] [ time-range ]
Fit Legendre polynomials up to degree n and ...
I - return the
polynomial prediction (default: the residual)
P - diagnostic printing
A - Analyse-only mode.
Option I is not applicable then.
Notably, with AP, the analysis will show the Akaike
criterion
for increasing model order 0...nr
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
(sasu06.f)
n - polynomial order
p0
p1 ... p#n -
coefficients
X.out <- SUM ( p(j)X.in^j
, j=0..#n )
ABS [ time-range
].................Take absolute
value of X.in.
X**#n [ time-range ].................Take power n of X.in
(inject MRS if negative).
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
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
0 1 ' '
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 SAVE
- Compute or get a filter, save but
don't apply it yet.
Allows a filter of max length 200,000, uses storage location 1.
FILTERCOMMAND
options STORE[->#n]
- Like SAVE, however specifically in storage location n
[n=1]
Allows 5 filters of max length 40,000 or 1
filter of length 200,000
FILTERCOMMAND [+C]
APPLY[<-#n]
- Apply a saved filter
[n=1]
FILTERCOMMAND
[+C] FETCH[<-#n]
- Fetch a saved filter
[n=1]
FILTERCOMMAND [+C]
options
- Compute or read, and apply instantly.
+C - Apply on a circularly connected domain
FILTER FORGET[->#n]
- Maybe useful after FILTER APPLY
(works only with those filters that use COMMON /CFILTER/
F(0:nfdim)...
ZSPF in the first case.
FORGET->#n will wipe only the filter stored at n.
FORGET will wipe the filter stored at 1 and the
active filter
NORMFILTER [+N][+H][+U]
+N - Give filter unity-response at frequency 0. Default
is +N
HIGHPASS [+N][+H][+U]
+H - Convert to High-pass.
Default is +H
+N+U - Give filter unity-response at Nyquist
frequency
can be applied between SAVE and APPLY
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 program
fscat filename
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 from an ASCII file
f(m),f(m+1),...,f(n)................- filter with
explicitly supplied parameters
.......................................Filters
the time series. Read filter coefficients
......................................from instruction file.
.............................#m
#n s - filter length (begin end) and symmetry.
+N
- normalise coefficients to yield unity response
Symmetry codes:
'S' - symmetric,
#m=0, #n=n,
=> f(1)..f(n) is mirrored to
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
IIR #n NONLIN=ζξ
-
uses
the
coefficients
that
have been input under FILTER or GETFS
FILTS U=#u [ F=f ] P=#m #n 's']
- input a time series from unit #u and use as a filter.
P= - Same rules apply to #m #n 's' as with >loc>
in FILTER
except for the coding of the command line:
either commas are used and blanks are
avoided
or string is
right-delimited with ]
#m - Sets the
anticausal start point. The most common
filter type would have #m=0
and 's'=' '
#n - Any big value will be
accepted. The actual parameter
is calculated from
#m+L-1
where L is the length of the time series.
f - if a file name is given, the file will be
opened
else an open file
is expected on unit #u
FILTER U=#u F=f [+N]
[>loc> [O=opt]]
.....................................-
file instruction
.......................................File
unit #u and name f are passed to openfs,
location string loc and option string opt passed
to
.......................................qloc_in_file.
E.g.
FILTER U=23 F=bandpass_1_2.fpf >-16
16 'A'> O=AB
The file must begin with one record supplying #m #n 's'
and the search string, if used, must match exactly (!)
+N - normalize filter
gain
WIENER U=#u
F=f L=#lw
Process a time-domain transformed gain function
(a.k.a. impulse response), for instance output during
sasm03 step response display and time-limit it
using a HANN window of full duration #lw
File record: read (iu,*)
i f(i)
PEF [U=#u
F=f] L=#l
O=opt [+F]
Use a Predicition Error Filter, length l
PREDF
[U=#u
F=f] L=#l
O=opt [+F] SAVE
Prepare the corresponding prediction filter for use
under IIR APPLY
Padsasdsd.............Use a Predicition Error Filter, length l
.......................................If
options U= and|or
F= are used, read
from filter bank
from file unit u,
file name f. If U= is not used, a unit
is assigned automatically.
Else the filter bank prepared by a preceeding BURG command
is used.
Option opt is passed to read_pef
or take_pef, respectively
(see ~/sas/p/pefs.f)
+F -
Swap the filter type (PEF <-> PRF) by brute force.
If not used, the type is assumed as detected from the file.
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 [O=opt] [ time
range ]
Creates a prediction-error filter bank, max order l. You need
[ #u < filename
]
the PEF command to apply a filter from the bank.
Options opt
O - enable output; one
more instruction line is needed to open
the file (binary!) using openfs (start in column 1):
#u < filename
File
name
may
contain
`###´,
which
will
be
replaced
by a running three-digit number.
S - scale filter for
unity power output n
UPDFILTER{A|M|R} #i
#v
Update filter coefficient i:
A - add v
M - multiply with v
R - replace with v
MEDIANFI[LTER] #l
[M=#m] [STORE]
Filters by running median, segment length #l
Can be moved foreward by m (default 1) for
decimation
STORE - Will not apply the
filter. The parameters are set aside
for subsequent MEDIANFILTER
APPLY or
file import to array W under option 'T.' ]
MEDIANFI[LTER] APPLY
WINDOW[:v] type
[P=#p] [T=#t] [ P
] [U] [ time
range ]
Set up a window and, if ` P
´ ("prepare") is missing,
apply it on the data.
The length of the window is determined from the time-range.
U - unit-normalise the
response (window shape and missing
data will modify the response).
:v = :W -
Apply it on the work space, else on the primary data.
:v = :F - Apply it on the
filter currently stored away with SAVE
type - 4-character code, see ~/sas/p/window.f
P=#p - Design parameter if needed. Default is good
for KAIS [1.5]
T=#t - Delayed taper, 0 < t <
1 [0.0].
WIN[:v] [U] [ time range ]
Use
the
currently
defined
window
and
apply
it
on
the data.
:v = :W - Apply it on the work
space, else on the primary data.
:v = :F - Apply it on the
filter currently stored away with SAVE
U - unit-normalise the
response (window shape and missing
data
will
modify
the
response).
CORR [ X
][ W ][ C ]
- Report the central cross-covariance coefficient in the
protocol
along with the RMS of X and W.
Y is synonymous with W.
Options X W (Y) and C request that the DC-level
in
the computation of the parameters is to be kept in.
CCV [L=#l] [M=#m] [D=#d]
Compute a running cross-covariance between X and Y
Use a file processing command with
option cdir = K*
to read and store Y.
Output is normalized with the number of valid data
pairs.
The two-sided domain length is #l and the center
is shifted by #m. Y is shifted with the
amount #d.
The t0 and dt
is adjusted to the new conditions.
Defaults: l=101,
m=1, d=0.
STRUCTURE[FUNCTION] [L=#l] [N=#q] [A]
[ time range ]
Return the structure function up to lag 2#l-1
[10]
#q - The norm to be used (floating
point)
[2.0]
A - Analyse, print the
result in the protocol:
[log(s(k))-log(s(1)]/(k log 2), k=2..l
where
s(k) = Average<|x(i+2k-1)-x(i)|#q>
Returns the series s along with dt=1
AUTOCOR[RELATION] #k [+W] [time-range]
Compute
auto-correlation up to lag k
Uses W for intermediate
storage of result!
AUTOCOV[ARIANCE] #k [+W]
[time-range]
Compute auto-covariance up to lag k
Specify an absurdly big lag to compute all terms.
+W - Keep result in array W
MEMSP #n [R:s] [L,][U[=#h,]][Q][/N][dB]
After producing or reading a PEF, compute the
Maximum Entropy spectrum, domain length n.
Options see PERIODOGRAM below.
PERIODOGRAM [R:s] [L,][U[=#h,]][Q][/N][dB][,M<#m] [ time range ]
Return the periodogram. If called from tslist, the
output series should be printed using -NI8 -n1-1
Before the periodogram is computed, the time series
is repaired (linear interpolation of missing data)
and the DC-level is subtracted.
R:s - Report Nyquist
parameters in units of:
R:s - Hz and seconds
R:h - cyc/h and hours
R:d - cyc/day and days
Returns the spectrum along with dt=1
Prints the Nyquist-frequency and the spectral bin
interval on the protocol.
The
following options must be coded without
intervening blankspace:
Q - return the
power
[Amplitude]
L
- return the logarithm (base 10) of the
spectrum. [Linear]
/N - scale power by the
number of
samples
[unity]
(amplitude: the square root)
dB - return decibels
U[=#h,] - Compensate the effect
of a filter (1,-#h)
[don't, else h=-1]
M<#m - Repair option: Maximum
#m repairs, linear
interpol-
ation is
used.
[all but 2]
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
Computes a series of length n,
where each sample
x.out(i) = SUM
{j /. MOD(j,n)==i : x.in(j) }
n - cycle length
DECI2W #n
#s [time-range]
Decimation (undersampling) to work array
n - decimation factor
s - initial
shift
DECI[MATE] #n #s
Decimation (undersampling). The offset is applied
before. Thus, x_out(i)=x_in(i*n+s), i=0,1,...
n - decimation factor
s - initial shift
DECI[MATE] #n #k S=#m ...with moving-average filter of length m (odd)
DECI[MATE] #n #k G=#f,#m
...with Gaussian-bell filter of length 2m+1.
Filter's -3dB point is at
f = f(-3dB)/f(Nyquist)
DECIMIN | DECIMAX | DECIMED
Three kinds of decimations:
DECI... #n #s [O=options]
Decimate by taking the max or min or median value
of the data segments.
n - decimation length
s - initial shift
option P - diagnostic
print
FOURIERINT[ERPOLATION] R=#r [O=opt]
(using ~/sas/p/spectralf.f ENTRY Fourier_Interp)
Fourier interpolation with optional taper
r - rate factor > 1
opt - string #lt,#s[,+P]
lt - length of the taper, number of bins.
s - width of Gaussian bell, argument (i/(lt*s))2,
0<=i<=lt
+P - print tapering action on protocol.
FFTINT[ERPOLATE] I=#r
[ Time-range ] (using ~/sas/p/fftintps.f)
Oversample with Fast-Fourier. If a time range is specified,
the resulting series will be shifted and t0 (the origin)
be adjusted. The main purpose is to interpolate through gaps.
r - integer rate r > 1
RESAMP[LE] #o #t #dt
[T!] [sec]
Resample with polynomial interpolation.
o - polynomial order,
t - origin time
dt - sampling interval
origin time is with respect to the previously first sample,
or with respect to epoch (T!)
.......................................t
and dt are in hours or in seconds (sec).
RESAMPLE TT
The following records
constitute a time-tie, terminated with
END S[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.
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.
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]
/N: the ramp goes from 0 to +v over the
length of the series
/TR: the ramp goes from 0 to +v over the
current time range
The zero-point can be changed:
Iref=#n - by time-series
index
Iref=C - at the center of the time-range
(zero-mean)
Tref - at the specified point of time
EXP #a #p u [
Iref=n | Tref yyyy mm dd hh mm
ss fff ] [ time-range ]
add an exponential a exp(j p dt), p in [1/h] real-valued.
Iref= Tref - see SLO
u - the units in which p is given
(default: [1/h]):
[1/s] [1/d] [1/y] [s]
[h] [d] [y]
ADD {#v | =DC | -DC } [RSQ]
[ time-range ] add value v to valid
samples.
RSQ - x <- sqrt(x2 + v2)
if v>0
x <- sqrt(max(x2 - v2,0))
else.
ADDRSQ is possible as a command.
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]
[ time-range ]
add a sine wave with amplitude a, frequency f [cyc/h],
and
phase p [deg]
COS #a #f #p [u]
[F'=d]
[FM=s] [AM=s] [
time-range ]
add a cosine wave...
u - measurement unit.
Specify [cpd] [cph]
[Hz] [mHz]
for frequency
or [d] [h] [s] for period,
respectively. Default: cph.
For cycles per time step, specify [DT]
For
time steps per cycle, specify [/DT]
PULSE #a -1 time-range
- add a zero-mean pulse of amplitude a and off-centre
half-width w.
-1 - Full width is given by time-range.
a < 0 - Amplitude is
computed by RMS-normalisation multiplied by |a|.
w = 0 - The pulse is a single spike of amplitude |a|.
MOD #a #nf #ip #iduty #offset [
time-range ]
add
a rectangular wave
offset+scale*{(MOD(i+ip,nf)
<= iduty?) +1 else -1}
TRI #a #nf #ip #iduty #offset [
time-range ]
add a triangular
wave, rising-flank length = iduty.
Based on the same
MOD-equation as above, but signal
varies between 0 and 1 when offset=0
and a=1
RDC [ time-range
]
remove a DC-level inside the window
ADDNOISE #a,#s type [ P=#p
] [GMP=#g] [ time-range
]
a - amplitude
s - random seed
g - Gauss-Markov parameter: y <- g y + ran(s),
x += a y
type FRACTAL - with power-law parameter p
type GAUSS -
type UNIFORM -
This function uses workspace W . Consider
dumping it
with DUMPW #u for a check.
RSM
{I|0} [ time-range
] call
repair_single_misses (isolated dropouts)
Mandatory option 'I' or '0' (zero).
REP [L] [ time-range ]
REPAIR [ opt
] [+P] [ time-range
] Options: V=#v
| L |
L<=#n | H#n
..............REPAIR
repairs breaks by substituting DC level.
..............REPAIR V=#v
repairs breaks
by substituting numerical value v.
..............REPAIR L
repairs breaks by linear interpolation.
..............REPAIR H[#n]
"sample-and-hold",
repeat last valid value #n
times at maximum
into a gap (default #n:
indefinite)
..............REPAIR
L<=2 repairs
breaks
of length 1 or 2 by linear interpolation
..............REPAIR
finds all breaks (REP ALL)
..............(REP L ALL
is perhaps too difficult for the program.)
+P - additional option to print out the places to repair.
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.
GAPFINT [P=#m] [ time-range ] - Fill
gaps in X by interpolation of W,
polynomial order m
Default m=2
WGAPF #a,#c [SQR]
[ time-range ] - Fill gaps with a half-sine-square of
amplitude a2,
add constant c2 everywhere. This function is
used to
generate
weights for interpolated tides in gaps.
The
constant a can be computed with
urtipgt @tintpolerr.ins; fgrep 3600
o/tintpolerr.dat # ( -> 22.122 )
SQR - Return the square root
HARMONICRESTORE
#f1 #f2 #fny #thrl #thru
[ time-range ]
Interpolation of missing samples or samples out of
range thrl < thru by Least-square fit of
a
band-limited spectrum, f1
< f < f2. The frequency
scale is
given by fny, the Nyquist frequency.
Mediocre results, close to useless so far.
No
defaults for frequency range, -1040 to 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
excessive printing.
#tol - tolerance for solution misfit,
default 1.d-6
#itr - maximum number of iterations, default
1000
#p
{4|5} - Number of parameters to
solve.
4: all the linear
dependencies
(1+2: linear-varying sinus amplitude, 3+4: linear ramp)
5: also the frequency
parameter (not recommended,
however the default is 5 (!)
.
DEL [ > #v ] [ time-range ] - delete (criterion abs(x(i)) > v )
DEL [ < #v ] [ time-range ] - delete (criterion x(i) < v )
DEL [ >> #v ] [ time-range
] - delete (criterion x(i) > v )
DEL [ => #v ] [ time-range
] - synonymous
If no value is given, scope for DEL is due to time-range.
DEL == #v +-#d [
time-range ] - delete if |x(i)-v|
< d . OBS! no blank between `+-´ and number.
(A complementary command /= is still missing.)
RJC directive [ time-range ]
Copy the missing-value symbols in X and W
according to the following directives
X=W - if x(i) or w(i)
is missing, mark both as missing
X<W -
Missing values in W are copied to X
X>W - ...
or vice-versa
RJC [ > #v
] [ -DC=#v ][ time-range ]
- delete (criterion abs(x(i) - DClevel) > v
)
The other options known from DEL exist too,
but might be less useful.
REJECTSHS #n [
time-range ] - Reject short segments,
amount of valid samples in each
segment will be > n.
INJ[ECT] { #v | MRS | ITEM } [ N=#n ] [ At time ]
Shifts the series towards the end by n samples, extending its
length, starting at the specified time. The shift may be
negative. If positive, the value v is filled into the gap.
Symbolic values are MRS (missing) or ITEM (the value at the
emerging gap). Default n=1
TRUNC [ time-range
] -
truncate beyond that time. Pefer to specify
At #YYYY #MM #DD #HH #MM #SS #FFF
TRUNC N=#n
- truncate to retain #n samples.
TRUNC #n
- should also
work like N=#n
EXTEND[W] #n #v
EXTEND[W] #n [V={DC|#v}]
EXTEND[W] by #m [V={DC|#v}]
EXTEND[W] [V={DC|#v}]
To #YYYY #MM #DD #HH [#MM [#SS [#FFF]]]
Extend the series X resp. W to
obtain n samples or add m or to
the specified date. Fill with value v or DC level or
mark as missing, the latter by default. Issue REPDC first, in
order to apply the actually valid DC-level.
TRIM time-range
- retain indicated section. Time origin will be adjusted.
MIXADJUST [O:[+R[+T]][+S]]
- Use W (low rate) to adjust X
by adding a linear curve at and
between the support points. Useful to combine Atmacs data (W)
with local pressure loading (X).
+R
- Replace W with the difference W - X
+R+T - ... and
trim to fit time scope of X
+S - If a supporting sample of X
is missing, delete the samples
of the whole interval; else try and find a sample in the
neighbourhood. See
TD/Atmacs/HOW.TO
Zero-padding for spectral filtering,
periodogram etc.
PAD -> spectral filter or periodogram or ... -> UNPAD
is
recommended in order to avoid effects of circular correlation
PAD [{N|#n}*{DC|#v}]
[+I] Extend the
series with amount n injecting value v. Default
N*0.0d0
PAD ={DC|#v}
If default for n is o.k.
N - to double the length of the series
DC - to substitute the series' DC-level
CPAD [{N|#n}*{DC|#v}]
[O=#l] [+I] time-range
- For use with ZSPF: Copy a padded
segment into a work area.
Add l
additional samples to avoid contamination from an eventual
jump
at the the start of the pad.
The
work area is allocated at the end of the signal array X,
so
that the calling program must take this extra
space into account.
tslist, for instance, should not be called with more than one
data
column in order to prevent invasion.
[C]PAD
+I
- For
time-integrating ZSPF applications, refresh the pad space
(equivalent to a zero constant of integration).
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.
RESYNCUNPAD
After CPAD with nonzero overlap, resets the time of the first
sample due to the overlap.
RESYNCUNPAD +A
In a block
DO time-range // CPAD // ... // DONE
with a time-range not starting at zero we
store this offset and apply it now to adjust the time of the
first
sample. See example (12).
RESYNCUNPAD [#n]
However, if CPAD is called repeatedly in a loop, RESYNCUNPAD +A
becomes uncertain. Specify the sample number n manually.
Perhaps there's demand for obtaining n from an
expression
RESYNCUNPAD At time-range
T S F - E D I T B L O C
K S T A R T A N D E N D:
TSF EDIT codeword
target word = TSF EDIT (capitals!)
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 00 00
DEL From 1993 10 12 20 48 00 00 To 1993 10 12 20 48 00 00
DEL From 1993 10 12 18 48 00 00 To 1993 10 12 19 12 00 00
DEL From 1993 10 13 00 48 00 00 To 1993 10 13 00 48 00 00
DEL From 1993 10 12 16 48 00 00 To 1993 10 12 16 48 00 00
DEL From 1993 10 07 20 48 00 00 To 1993 10 07 20 48 00 00
Example (6)
Take a segment of a signal, construct a prediction-error
filter, and apply it to the whole signal
Save the filter for later use.
TSF EDIT BURG
BURG L=20 O=O From 2009 10 02 12 00 00 00 To 2009 10 02 13 00 00 00
41 < o/rt.pef
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 0 0
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 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) ~/Seismo/gcf/resmp.tse
The input has 100 S/s and is loooong. We have to start filtering
at offset 20 samples
to obtain output on precise seconds. That's done in the DO
command.
Because of the overlap option in CPAD we must call RSYNCUNPAD
after the loop.
It's also the place to reset the start time given to the first DO
command.
TSF 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
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 0.0 SAVE
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=AB 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 0.0 SAVE
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
./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