SASM06 - Cross-spectrum analysis
System input = X from
file unit iux
" output = Y from file unit iuy
Y = hx X + hz Z, partial coherence X with Z, requires four jobs:
1 - X/Z 2 - Y/X 3 - (Y/Z)|X 4 -
(Y/X)|Z
T F T
F T
T F
T = qsavesp, qpartial
Computes power spectra, cross-covariance, coherence, gain, phase
Wiener filters
A disadvantage by construction:
Spectrum data output to file together with confidence limits can
only be issued from within graphic display.
splist
salvages the flaw to some extent.
Yet, you can plot the results with GMT
Missing documentation:
Partial cross-spectrum; examples.
There are a few environment parameters that supersede namelist
definitions:
SASM06_NDATA ndata
SASM06_INTERACTIVE [Y[ES]]
SASM06_FDATE date
SASM06_FHOUR start
SASM06_MARKINF mark_inf
Namelist parameters:
NAME
- TYPE - DIM
MEANING
[Default]
---------------------------------------------------------------------------------------------------------------
qbatch -
logical - Batch mode, no graphics
[.false.]
Run-time deselection: setenv SASM06_INTERACTIVE YES
About data input:
No time series input
- See examples
qinput_cxy - logical
- Input CXY from
unit 55, PXX from 56, and PYY from 57.
Please specify fny_1h, time_scale, ls, and dtx in
namelist [.false.]
qinput_acv
- logical
- Input CXY from
unit 55, ACVX from 56 and ACVY from 57.
Please specify fny_1h, time_scale, ls, and dtx in
namelist [.false.]
dtx
- real*8 -
Sampling interval in hours. If dtx<0, 1/dtx ->
dtx
[0.0d0]
Time series input
ndata
- integer
- Input data
length, default =
maximum
[Parameter NTDIM]
Run-time setting supersedes:
setenv SASM06_NDATA ndata
data_mrs, tst - real*8 -
Internal Missing Record Symbol and test
accuracy [-99999.9d0, 1.0d0]
date
- integer - 3 Year
Month Day of epoch; default: let input
define [-1,
-1, -1]
start -
real*8 - Starting hour from
epoch to cut input, default = don't cut
[-9999.9d0]
Run-time setting supersedes:
setenv SASM06_FTIME date
setenv SASM06_FHOUR start
time_scale - real*8
- Time scaling for
spectrum frequency
units
[24.0d0]
time_scale < 0: use 1/|time_scale|
e.g. -3.6 will give units in mHz, -3600. in Hz,
1.0 in cyc/h, 24.0 in cyc/day
fny_1h - real*8
- Scaling of the frequency
axis.
[12.0d0]
The displayed Nyquist-frequency is
fny = (fny_1h/24) * (time_scale/dt) where dt is the
sampling
interval of the input. Thus the defaults will give
cyc/days
and a Nyquist frequency of 12. if the data is hourly
sampled.
Specify 180.0 to obtain deg/h.
fny
- real*8 -
Overriding the Nyquist frequency that would be
calculated
[0.0d0]
ufny
- char*10 - Units
for the Nyquist, e.g.
[Hz]$
['[?]$']
Some automatics will try to determine it.
iux, iuy -
integer - Input
file units for series X and
Y
[21,22]
fx_mrs, fy_mrs - real*8
- ASCII input
files' Missing Record
Symbols
[2*-99999.9d0]
trgx, trgy - char*16 -
ASCII file location targets or 'BIN'
for binary files [2*'TOP]']
fmtx, fmty - char*64
- ASCII file
formats
or
[2*'(3i2,1x,i2,5x,f15.0)']
'U' for unlabeled binary ts-files or
'L:label' for multi-component binary
files.
'L@' will take the string `L:label´
from the comment
in the open-file block.
khmsx, khmsy - integer
- ASCII file time
record depth, 1 - 4
[1]
4 means all, i.e. Hours Minutes Seconds
Milliseconds
itzx, itzy - integer
- Files' time
zones, default = UTC
[0,0]
dtx, dty -
real*8 -
Sampling interval to impose. Default=take from
input
[2*0.0d0]
scale_x, scale_y - real*8
- Input scale
factors
[2*1.0d0]
q_edit_x
q_edit_y -
logical - Call tsfedit
with instructions at location target ...
[2*.false.]
edit_block_x - char*16
- in the
instruction file
['X]', 'Y]']
edit_block_y
q_repair_x
- logical
- Allow to repair
missing samples. Subroutine sas/p/repair.f
q_repair_y
is called with option 'I' still hard-wired
[2*.false.]
fw_mrs, trgw, fmtw, khmsw,
itzw Like the X and Y
equivalents. This data may be added to Y
wadd
- real*8 -
Scale factor on
W
[1.0d0]
Prefer to use the TSF EDIT feature!
qdisplay -
logical - Display
time series
[.true.]
qdisplay_raw - logical
- Display also the
raw time series
[.true.]
mx_miss
-
integer
-
How
many
missing
samples
are
permitted (they may be repaired)
[0]
If
more
than
mx_miss
samples
are
missing,
the record segment
is rejected. Contemplate the Bartlett method!
qdownw
-
logical
-
Downweight
the
data,
istructions
from
a `Downweight´ block in
the instruction file.
xstat
-
real*8
-
Remove
xstat*X
from
Y
(correct
a static, known contamination) [0.0d0]
Alignment of series X and Y
shift_x
- real*8
- Lead of X with
respect to Y, initial guess, in units of dt
[0.d0]
shiftc
-
real*8
-
While
shift_x
can be changed at prompter, shiftc is always added [0.d0]
A positive number means that the start of X is earlier than
Y.
Decimation: We might require
decimation of one or the other series. The following parameters
control the process
len_df
- integer - Length
of the
filter
[48]
wtype_df - char*4
- Window
code
['KAIS']
wdp_df
- real*8 -
Window design
parameter
[2.5d0]
qff_df
-
logical
-
Use
the
float-frequency
method
(Subroutine
sas/p/fwwds.f)
[.false.]
safety
-
real*8
-
Distance
between
Nyquist
and
upper
pass-band corner
as a fraction of the Nyquist, i.e.
s=(f-fny)/fny
[0.0d0]
qfix_t00, qfix_dt - logical
- In alignment,
the t0 and dt of both series can be forced
[2*.false.]
to
be
the
same
by brute force, namely...
t00,
dt00 -
real*8 -
Spectral estimation: A
prediction-error filter is read from unit 25 (if open). The
purpose of the filter is to get a white series X.
A simple 1-parameter differencing
filter may be applied first
method_cov -
char*1 - 'B' for Bartlett, any other
character: explicit
(useful
if
series
have
many and long gaps)
['B']
Specify 'X' if spectra or covariance series are input
(qinput_cxy=.true. or qinput_acv=.true.
pef_order
-
integer
-
Filter
order
(number
of
off-zero
coefficients)
[4]
dff
-
real*8
-
2
Difference
filter,
e.g. 1,-1 to whiten
Brownian noise [1.0d0, 0.0d0]
ls
- integer -
Segmentation: Length of data
segment
[128]
maxsect
- integer -
Maximum number of segments
[9999]
trunc
-
real*8
-
By
default,
the
covariance
window
will have ls
as
the
truncation
point.
This
can
be
narrowed in with trunc < 1 [1.0d0]
wtype_sp -
char*4 - The
window code for the spectrum
estimation
['KAIS']
wdp_sp
- real*8 -
Window design
parameter
[2.5d0]
lsmooth_psp_lvl - integer
- For the PSP
level around which the confidence interval is
drawn, how many bins to each side are included in the
average [ls/10]
Partial cross-spectra input uses file units 61 (power) and 62
(cross), output 51 (both power and cross)
qpartial
-
logical
-
Partial
cross-spectra:
Input
the
third
agent's spectrum
Spectrum output:
qsavesp
- logical - Save
real-valued spectra in binary, unit 51,
for later use as a third agent
or printing/manipulating with splist
PSX PSY CSP
[.false.]
iun_zadmit - integer -
If > 0, save complex valued
spectra in binary to this file unit:
CAD ADC - Complex admittance and confidence limits
CPH APH - Complex phase and confidence limits.
File must have been opened. qignore_cfd is in
effect.
[-1]
iun_cohc -
integer - If > 0,
save coherence spectra in binary to this unit:
KOH KOC - Coherence and confidence limits.
File must have been opened. qignore_cfd is
in effect.
[-1]
If iun_cohc > 0 and not 51, KOH and KOC are written to
51 too.
qignore_cfd - logical -
If .false., only values where coherence
> conf.cohc
are
written
[.false.]
adm_file
pha_file cohc_file (function is currently disabled) file
names for spectrum output
infa_file infp_file
Wiener
filter assembly parameters
The system prompts (in graphic mode) for explicit filter length. If
no value is entered, the system automatically increases or decreases
the filter length in a loop, halting either when length reaches 2 or
the domain length:
The logic behind this is perhaps too spaghettish. See code at
label 402.
len_wf
-
integer
-
Starting
length
[32]
Recursive variation of filter length:
len_wfi_prc
- integer
- if .ne.0,
percentage increase (or decrease if negative) of length
[0] ?? working ??
len_wfi
- integer
- if .ne. 0, add this
value to the length or ...
[-4]
len_wfd
- real -
if .ne. 1.0, multiply length with this
value
[1.0]
qinf_binary - logical
- Output Wiener-filters
to binary file (unit
26)
[.false.]
mark_inf
- char*40 - Mark the binary filter
labels with 'WF_nnnnn'//gully//mark_inf
[' ']
Run-time specification supersedes:
setenv SASM0_MARKINF mark_inf
Wiener filter response tuning:
q_tune_by_ratio - logical
- If itune_wf=0, tuning
is by adding the loss. If this is
not
desired,
set
q_tune_by_ratio=.true.
See itune_wf.
[.false.]
itune_wf
- integer
- Spectral bin number
where resulting Wiener filter is to
be amplitude-adjusted to
the nonparametric gain spectrum.
If
itune_wf=0,
the
missing response is added unless
q_tune_by_ratio=.true. Then, the whole filter is
up-scaled.
Default
-1
means,
don't
tune
[-1]
frq_tune_wf - real*8 -
Alternative
to
specify
the
tuning
frequency
(fraction of Nyquist)
itune=frq_tune_wf*domain_length
[-1.0d0]
Other
arlx nrlx - need to
check the meaning of these two guys. Traces lead to SKB KAS03
project.
arlx is used in zprelax (sas/p/prelaxs.f),
nrlx is the number of elements in nrlx. Probably a spectrum tuner,
e.g. how borehole pressure relaxes over time
after a surface pressure impulse.
Graphics: Titles etc.
series_x series_y series_yl series_w
-
titles;
_yl:
after
eventual
linear
combination
(xstat)
E.g. 'GWPress'
series_x
series_y Special value '@' -
take the string from the respective file comment
xunits
yunits - measurement units
(short forms) before scaling, e.g. 'hPa'
units_x units_y - ...
after scaling, bracketed forms
units_psp_x units_psp_y - e.g. '[hPa]$'
npal
- ordinal number in the palette file
q_shut_graphics
- .false. keeps graphic screens, .true. closes and opens.
q_grey
- grey scale plots
n_xw,
m_xw
- graphics window size [640,480]
Spectrum graphics:
q_unfilter
-
Show
power
spectra
as
they
are
after filtering
(dff
and
npef)
if
.false.,
else
restore
filter action
cutpha
- 180 to show phases in
range -180 to +180
0 to show phases in
range 0 to 360
cutlvl
-
Discard
spectral
power
below
this
level
when setting
up plotting ranges.
q_psp_logax - Power spectrum, logarithmic abscissa
q_cohc_logax
- Coherence spectrum,
logarithmic abscissa
q_adm_logax - Gain spectrum,
logarithmic abscissa
q_pha_logax - Phase spectrum, logarithmic abscissa
q_ifs_logax - Information
Filter spectrum, logarithmic abscissa
pha_min pha_max -
ordinate limits phase graphics
adm_min adm_max - ordinate limits gain graphics
adm_add
- In the cross-spectrum, adds adm_add * X to Y/X
opt_display_sp_pwr opt_display_sp_adm
opt_display_sp_chc opt_display_sp_pha
opt_display_fs_cor opt_display_fs_cov
opt_display_fs_rsp
-
char*64
-
Options
to
display
and/or
write spectrum:
G+ G- -
show screen graphics + yes - no
W+ W-
- write series to file + yes - no
W+:filename - write series to a specific file
answer_wf_prompt - Batch mode: answer the
Wiener filter prompter
answer_shift_prompt - Batch mode: answer the
covariance-shift prompter
nrjd
- number of harmonic overtones of a diurnal to be
cleaned out (subroutine sas/p/rjdius.f)
Unused for a long time.
Batch
mode
The following namelist parameters must be set in
order to completely skip graphics:
qbatch=.false., qdisplay=.true.,
qdisplay_raw=.false.
answer_wf_prompt='answer',
len_wf=length
answer - what otherwise is entered
at the prompter for transfer model
control, e.g. 'IMP L=-1'
For skipping the transfer model and Wiener filter assembly,
answer='E'
length - The starting order
("length") for the Wiener filters. The series is assembled
in a loop
in which the
lengths are computed as described under Wiener
filters above.
Examples for batch mode:
/home/hgs/Ttide/SCG/sasm06-rios-atmacs.ins
/home/hgs/Ttide/SCG/sasm06-rios-atmacs.ins
qbatch=.true.,
qdisplay=.false., qdisplay_raw=.false.
answer_wf_prompt='IMP L=-1', len_wf=128
q_stop_after_cxy=.false., iu_cov=31
Files
File unit1) type2) C3) I/O Content
---------------------------------------------------------------------------------------------------
N: iux[21] A B *
I Time series `X´
N: iuy[22] A
B * I Time series
`Y´
N: iuw[23] A
B I Time series `W´
25
B
I PEF-filter. Use sasm03 with a
diff-filter and use diff and PEF
here in the same way
26
A
O Wiener filter bank. Must be opened in
the open-file block
if the coefficients are to be saved.
51
B
O Power spectra, suggestedly real-valued
optional; e.g. according to namelist parameters
qsavesp=.true., iun_cohc=51
52
B
O Cross-spectra, preferably complex-valued
optional; e.g. according to
namelist parameters
iun_zadmit=52, iun_zphase=52
55
A
I j CXY(j) // RMS(X) RMS(Y)
56
A
I j f PSP(X)
57
A
I j f PSP(Y)
Files 56 and 57 must have a header:
[ dB ][ C=P ]
"dB" signifies that the subsequent power records are in dB
"C=#value" specifies a constant power level. Subsequent
records will be ignored.
Should be extended to a range of power spectral models:
fractal: P*f^alpha,
Brownian: P*cos(pi f),
Gauss-Markov: 1 + alpha^2 - alpha cos(pi f) ...
MEM (AR): C/|Poly(p_0 ... p_n ,
z)|^2
51
B O
Input, partial cross-spectrum estimation,
PSP(X) PSP(Y) and CSP(X,Y)
(changed - but to what? - controlled how?)
61 B
I Partial spectra, [PSP(X) PSP(Z) skipped]
CSP(X,Z)
62
B
I Partial spectra, [PSP(Y)
skipped] PSP(Z) CSP(Y,Z)
8
A i I
Palette. File must exist (link!) in the working directory.
9
A i O Used
for spectrum output within display_sp (sas/p/grdsps1.f &
grdsps3.f)
99
A i O
some printing useful for debug esp. graphic routines,
e.g. /dev/pts/3 (a second xterm window; use ps to find out
pts-number)
(Some ominous files are created but nothing is written to
them: sp-hh:mm:ss.dat . Some sort of spectrum output is not
performed.)
1) N -
Namelist parameter and [default]
2) A - ASCII,
B - binary. Select using N: trgx trgy trgw, see Namelist
description.
3) `*´ - Compulsory, unless
qinput_cxy=.true.
`i´ - Used internally ´.
We write spectrum data from
within display_sp and display_sp_cfl
File names are to be entered at the Spectrum plot leader page
(text mode)
If a file name e.g. gravity_pressure.adm
is entered, a name for the corresponding confidence file is
suggested in the graphic window.
gravity_pressure.adm.cfd
for the example. The name can be edited if desired.
Prompts
The Shift prompter after the
cross-covariance plot:
Enter a shift, a positive number to move the largest peak to the
right of the mid point closer to zero lag.
Before time-series graphics appear ... (grdtss.f, needs to
be written)
Before filter, auto-, cross-covariance ... (grdfss.f,
needs to be written)
Before spectrum graphics appear, e.g. (emphasis added):
GRDSPS> PSPX zacc
Spectral lines, Nyquist: 0 .. 512
5.0000E-01
Q Quit
or S Skip
F {
bin-range(int) | freq.range(float)}
D {Connect
connectEmph | Point | pIn [ D- I- P-]}
Y
ymin,ymax [ 0.000E+00 0.000E+00]
W
filename
of which W filename prepares file
output from within the graphic subroutine, spectrum estimate
and confidence into separate files - after confirming prompts.
The STEP or
IMPULSE RESPONSE menu with Transfer model / Information Filter
assembly
Q
- Quit and proceed with Wiener filter assembly
E
- End the program
?O
- Help functions
?R
?U
In batch processing, the
following options can be coded in the answer string
answer_wf_prompt in the namelist
E - End the program
I options
- Enter all options on the same line. Use ?-help
functions. C=T
causality options yields strange results. Bug or bad idea??
Information filter (Wiener filter) must be based in the Impulse
response. Enter
I option
where option L=-1 is
suggested (weighting the cospectrum by coherence - the preferred
method).
I L=-1 FCUT=f
where 0<f<1 specifies the lower fraction of
the spectrum domain where coherence is left unchanged; above,
coherence will be sharply decreased using a Gaussian half-bell.
Other examples:
I L=10 dB O=Min
sets the floor at 10 dB above the minimum of the cospectrum.
I L=-20 dB O=Max
sets the floor at 20 dB below the maximum of the cospectrum.
I L=0 dB
O=RMS
sets the floor at the RMS of the cospectrum
S L=-1
Step response, coherence-weighted.
Q
(automatically in batch mode) Next prompter (after impulse response
graphics): choose a starting length for the window.
Plotting the results with GMT
Try splist-plot-cxy
-h
(works with splist - if you have binary
psp/csp-files)
Also: plot-xsp , uses ASCII spectra. May be best to
show an example:
- Devise a project name indicative of the data types, here
`ap´ for air pressure and `tg´ for tide gauge.
- Create/rename the instruction file to sasm06 accordingly:
`sasm06-$prj.ins´ . sasm06 will pick up the file name
and replace extension `.ins´ with one indicative of
the type of spectrum - PSP for power, ADM for
admittance....).
- Devise a subdirectory with a name that's specific for
the parameter setup. Here, `270´ reflects the section
length and `-0.8´ the coefficient of the difference
filter (anything that you consider could be varied in your
project).
set prj = ap-tg-1d
sasm06 @
sasm06-$prj.ins :U
Save spectra in the graphics display using the W command.
When sasm06 has completed, save the files in a subdirectory:
mkdir sasm06-270-0.8
cp sasm06-${prj}.* sasm06-270-0.8/
Make the plot:
plot-xsp [
-N ] -p +++-%%%.ps -I
sasm06-270-0.8 $prj [ go ]
The ps- and png-file names will reflect the $prj and the subdir
name as follows:
sasm06-270-0.8-ap-tg-1d.png
i.e. `+++´ is replaced with the subdir name and `%%%´
with $prj
Other means: Using binary output of the suite of spectra,
splist
and splist-plot
will make life easy.
Examples for no
time-series input, i.e. ...
1. CXY cross-covariance and Power spectra
cd
/data1/hgs/Seismo/gcf/AG/
cat sasm06-*s*-ag.ins
2. CXY cross- and
auto-covariances
cd ~/TD
cat sasm06-ap-bp.ins
Data format for cross-covariance (sas/p/inputcxy.f):
Lag
CXY
CXX
CYY NDATA
-144
1.5646e-02 7.0912e-01
2.3097e+00 8423
-143 1.9396e-02
7.0943e-01 2.3015e+00 8424
-142 6.1095e-03
7.0956e-01 2.3049e+00 8425
(output from tslist/tsfedit XCORRU)
Data format for auto-covariance (sas/p/inputacv.f):
0 5.0415E-01
1
4.5989E-01
2
4.9618E-01
3
4.7452E-01
(output from tslist/tsfedit AUTOCOV, tslist
printing with -Ni -n-1 -p1 -qqq -w filename )
Data format for power spectra (sas/p/inputcxy.f):
2.1 - A PEF-spectrum, specfies a PEF-filter, here /data1/hgs/Seismo/gcf/AG/ag.psp
(sorry, needs more explanation)
C=2000.
PEF
0 3 '
' 1.388889E-03
2.115008E+02
1.00000000000000 7.082633060317301E-002
-2.006348203850501E-002 3.626794396894208E-002
C=225.
2.2 - A Power spectrum line by line specifying decibels, here
/data1/hgs/Seismo/gcf/AG/2HzsvE-250.psp
i f(i) Pwr(i)[dB]
0 0.0000E+00
7.0360E+01
1 4.0000E-03
6.7335E+01
2 8.0000E-03
6.4000E+01
3 1.2000E-02
6.0428E+01
4 1.6000E-02
5.8624E+01
5 2.0000E-02
5.7392E+01
(Output from tslist (?))
3. Flexible input from binary files (
~/Ttide/SCG/sasm06-ra-rain.ins )
We have an mc-file thanks to a range of
IIR-parameters, integrating 1-h rain height in order to mock runoff.
Open-file block, file 21:
21 ^
${SASM06_INPUT:[~/TD/tmp/ORH2009-OPNEND-iir-1h.ts]}
L:P${IIR1:[0.9]}
Namelist:
fmtx='L@'
Run:
tslql ~/TD/tmp/ORH2009-OPNEND-iir-1h.mc
P0.9 P0.95 P0.975
P0.09875 P0.99375
setenv SASM06_INPUT ~/TD/tmp/ORH2009-OPNEND-iir-1h.mc
setenv IIR1 0.99375
sasm06 @ sasm06-ra-rain.ins
:U
.bye