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 withing graphic display.
Yet, you can plot the results with GMT

Missing documentation:
Partial cross-spectrum; examples.

There is one environment parameter that optionally (if set) overwrites a namelist definition:
( => ndata)

Namelist parameters:

NAME              - TYPE    - DIM     MEANING                                                         [Default]
qbatch            - logical -         Batch mode, no graphics                                         [.false.]

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

ndata             - integer -         Input data length, default = maximum                            [P ntdim]

start             - real*8  -         Starting hour from epoch to cut input, default = don't cut    [-9999.9d0]
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)']
                                      'L:label' for multi-component binary files.

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_y          - logical -         Call tsfedit with instructions at location target ...          [2*.false.]
edit_block_x      - char*16 -         in the instruction file                                       ['X]', '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

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
qsavesp           - logical -         Save spectrum for later use as a third agent.

Spectrum output:
adm_file pha_file cohc_file  (function is currently disabled) file names for spectrum outout
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 ??
- integer -        if .ne. 0, add this value to the length or ...                       [-4]
- 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      [' ']

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

- real*8 -         Alternative to specify the tuning frequency (fraction of Nyquist)
                                      itune=frq_tune_wf*domain_length                                  [-1.0d0]

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:

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

 qbatch=.true., qdisplay=.false., qdisplay_raw=.false.
 answer_wf_prompt='IMP L=-1', len_wf=128
 q_stop_after_cxy=.false., iu_cov=31


File unit1) type2) C
3) 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.

   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)
   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. for the example. The name can be edited if desired.

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.

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

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 the coherence). 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.

(automatically in batch mode) Next prompter (after impulse response graphics): choose a starting length for the window.

Plotting the results with GMT

May be best to show an example:
  1. Devise a project name indicative of the data types, here `ap´ for air pressure and `tg´ for tide gauge.
  2. 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....).
  3. 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 -I sasm06-270-0.8 $prj [ go ]

The ps- and png-file names will reflect the $prj and the subdir name as follows:


i.e. `+++´ is replaced with the subdir name and `%%%´ with $prj

Examples for no time-series input
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)

  0    3 ' '       1.388889E-03  2.115008E+02
   1.00000000000000  7.082633060317301E-002 -2.006348203850501E-002  3.626794396894208E-002


2.2 - A Power spectrum line by line specifying decibels, here /data1/hgs/Seismo/gcf/AG/2HzsvE-250.psp
       i    f(i)        Pwr(i)
       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 (?))