sasm03  - Power Spectrum Estimation with MEM and Bartlett-method

      Program Power_Spectrum

SAS PSP * Bartlett's method; optional confidence intervals from MEM

A number of different issues can be addressed:
(1)  Treat a time series with a filter,
     e.g. a simple difference filter, a Prediction-Error or a Wiener filter

(2)  and estimate the power spectrum using MEM and {Periodogram | Bartlett}
     Display: time series, Filter spectrum,
              filtered time series,
              Spectrum of the filtered time series using
                       Bartlett_ACV and _PSP) + MEM_PSP
                       optionally with the effect of the difference
                       filter undone (useful if the power range is
                       greater than the window sidelobe level).

(3)  Output the PEF filters produced in the MEM step.
              Careful when using it later: the difference filter
              and the PEF must be applied in succession.

The MEM spectrum can be used to draw a confidence region around the
Bartlett spectrum. Parameter LVLLEN controls that.

Prompt at MEM spectrum:

Plotting will go through PEF_SHOW pre-selected filter orders
specified in array PEF_ORDER_SHOW(1..PEF_SHOW) (namelist)
and draw MEM with different colors.
The last one, PEF_ORDER_SHOW(PEF_SHOW) by default, will be used to draw
a backdrop to the power spectrum, the width of it signifying the
confidence interval of the Power Spectrum. During prompting
you can specify another MEM order. 

There is a key_prompt here.
Hit the   Return  key to go to the next curve,
or press  S       to display the rest non-stop
or press  Q       to quit the loop

The last time, at pef_show = +1,
the graphic screen will prompt for an integer:
     &('LPEF ['//buff(ib:4)//'] or < 1 to reset color :',20,470)

implying one of two modes, "O.K." (last filter shown) or "Try" (select a different order)

Positive number:   filter order (TRY step 2: draw in white and output the MEM-spectrum)
  Zero ... -899:   color 1..15  (TRY step 1: change color to something
                                             less conspicuous)
        <= -900:   clear plot area and redraw via text-screen prompt (when
                   diagram has become too clotted)

            ESC:   O.k.
O.K.: colour white, redraw but don't come back
      (continue with non-parametric PSP, but keep the MEM-spectrum
      for drawing the confidence intervals).

TRY:  first with the old colour, redraw, come back without
      prompt, color white, draw according to new pef-order,
      come back for O.K. | TRY.

Output of spectrum+confidence intervals to file:

There are three methods to output spectrum data.
The first will write BIN data that can be retrieved with splist.
The second is from within Sp_Prime_Data_? , the third is during Display_Sp

1. Open a binary file on unit 51 and include in namelist
   The following spectra are written, tagged as
     PSX - The power spectrum
     MEM - The Maximum entropy spectrum
     MEC - Confidence limits for MEM
     PCF - Confidence limits for PSX
   If a file is opened on unit 52, the power spectrum is output again, but this time tagged as
     PSP - ... and without other sections.  

2. Output within Sp_Prime_Data_? is useful in a batch job (qgra=.false.);
   it is activated by opening a file and setting namelist parameter
   iuspout to the unit number. This applies as well to
   the MEM spectrum (iumemspout)
42 B ${ODIR:[o]}/grav-1h.psp
43 B ${ODIR:[o]}/grav-1h.memsp
  iuspout=42, iumemspout=43

3. During Display_Sp
   Two options, interactive of automatic:

3.1 Interactive:
    Use the text prompter
GRDSPS> Power spectrum
Spectral lines, Nyquist: 0 ..  2048  1.2000E+01

        Q Quit or S Skip
        F { bin#range | freq.range }
        D {< Connect connectEmph | Point        n | I- >}
        Y ymin,ymax  [  -2.544E+01   2.936E+01]
        W filename

    to enter W filename, example
    W owf/grav.ra.psp

3.2 Automatic, same example:
    In the namelist, set 
    You can specify environment parameters with defaults, like
First, a file named
owf/grav.ra.psp is written with the Power values only.
    When the confidence is added to the plot, a file named
    will be written, containing both kinds of data.

Files to open: A: ASCII, B: BIN, D: ASCII or BIN depends on target,
                  ASCII input format is specified by FMT in namelist &param

D     21 - signal 'X' input
D     22 ~ signal 'Y' to subtract from 'X'
B     25 ~ PEF filter input, used for power spectrum confidence limits
B     41 - PEF filter output, used for spectral whitening                  1)
A     42 - PSP output
B     51 - Power and MEM spectrum output and confidence, upon
B     52 - Power spectrum, optional
B     $w ~ weights to restore effects of data gaps in covariance          
A iun_wf ~ Wiener filter input
B iun_wh - PEF whitening filter output / re-input                          1)
A     99 ~ diagnostic printing, /dev/pts/#n to a second xterm window
           recommended. Obtain #n from the ps command. 

          ~ open optional
          - open compulsory

1) If unit 25 is not opened, 41 must be opened since a new filter bank
    will be produced.
2) $w is set in namelist: iu_bartlett_cov_reweight = $w , $w > 0

Namelist PARAM:
Parameter   Meaning                                                  [default]
qgra       - .true.: Display graphics

Display    - char*16 - options for displaying time-series and ACV:
             W - weighting
             R - raw data
             E - edited data
             Y - the second time series, 'Y' (from unit 22)
             F - filtered data (Wiener filter, Difference filter)
             L - after linear combination with Y
             I - Residual from Wiener filter
             A - Autocovariance                                    ['REYLIFA']

l_section  - Section length for data segmentation. Length of computed
              autocovariance = 2*l_section;
              power-spectrum length = l_section
l_window   - Window's two-sided length; default= 2*l_section
t_window   - Window type ('KAIS')
p_window   - KAISer window design parameter (1.8D0)
lmems      - (length of spectrum domain for MEM - unclear how to use
              and whether it interferes with l_section)

rec_mrs   &
rec_mrs_y  - Symbolic value for missing records in input file        [2*999.9]
dat_mrs    -    "       "    "     "    internal data

q_tsf_edit - .true.: call tsfedit (after data filtering)
tsf_edit_name        The location target (block to be included
                     the ins-file).

t_filter   - 'WF' - use Wiener filter, ...                              ['  ']
iun_wf     - ... from this file unit                                      [27]
opt_wf     - ... and with this filter option                           ['SU ']

lpef_wh    - > 0 : use whitening filter (PEF) ...                          [0]
iun_wh     - ... from this file unit                                      [41]
iuspout    - output unit for PSP                                          [-1]
iumemspout - output unit for MEM                                          [-1]

qsavesp    - if .true., output BIN spectra to unit 51.               [.false.]
             (Use splist to read and process further.) Output are:           
             1. The PSP, 2. the MEM-sp, 3. the centre of the cfd-interval,
             4. its width. Targets are: PSX, MEM, MEC, PCF

           - file unit of reweighting coefficients, see Experiment 3      [-1]
ntrunc     - truncate data length
tst        - _mrs-test accuracy                                          [1.0]

target    &
target_y   - (char*16) Data transfer starts after this record.
              ] = right target string delimiter                     [2*'>> ]']

fmt       &
fmt_y      - (char*64) Data format

hms_pos   &
hms_pos_y  - (integer) Time record check; 1: HH ... 3: HHMMSS
              FMT must contain the number of HMS_POS integer
              format fields                                              [0,0]

time_zone &
time_zone_y  (integer)                                                   [0,0]

scale     &
scale_y    - (real*8)  Multiply data with scale                        [2*1.0]
start_date(3) - (integer) Data epoch Year,Month,Day
dtf        - if > 0, impose as the sampling interval                   [-1.d0]

lvllen     - Detection of noise level.
              LVLLEN > 0: Length of sliding average,
                       0: Defaults to   "      "     length 6 internally.
              LVLLEN < 0: Read a Prediction-Error filter from unit 25
                          If file not open: Estimate new PE-filter.        [0]
qshlvl     - (logical) Display detected noise level                   [.true.]
                        Needed if lvllen >= 0
cutlvl     - Spectrum display, minimum data

pef_order  - (integer) Order of noise model (PEF-filter length +1)
lpef       - Maximum order of noise model to be estimated (10)

units      - after applying scale
units_y                     scale_y

qrfsp      - restore difference filter in MEM and PSP spectrum       [.false.]
qrpefsp    - restore the PEF filter used to whiten the data          [.false.]
              in MEM and PSP

Do_Psp     - (char*1) 'Y' - compute PSP of time series Y
                            (after eventual filter and data reduction
                      'X' or any other char: compute PSP of time series X

Sp_File   - (char*64) - file name for spectrum output                   [' ']
Acv_File   - (char*64) - file name for autocov. output                   [' ']

logxax     - .true.: Spectrum abscissa logarithmic
time_scale - controls the frequency units in display
              > 0: cycles/(time_scale*hours)  e.g. 24.0 for cyc/d
              < 0: cycles*|time_scale|/hours  e.g. -3.6 for mHz

1. Experiment with PEF's
The data series can be described with PEF's of varying order.
After the (one) call to BURG with max filter order lpef we use the
filter bank to explore by graphics. The PEF with order per_order will
always be implemented at the end.


pef_show   - integer - the number of different filters (except pef_order)
pef_order_show()     - the PEF filter orders to show. It is advisable to
                        start with the highest number, since the plot range
                        is determined from the first call, which therefore
                        should contain most of the spectral detail.
                        pef_order_show(pef_show) will be used for confidence
                        estimates: a MEM-spectrum with PSP confidence limits.

2. Experiment with spectral power index
See the call to psp_powindex. Associated subroutines in sas/p/psppowix.f

mpix       - integer - If greater than zero, estimate the parameters
                        of a spectral model. Currently implemented
                        is p = a f**b + c, where f is frequency.
                        Estimated parameters are pix(1..3) = { a b c }
                        Thus, if mpix=2, log p = a + b * log f
pix(*)     - real*8 -  A starting solution (e.g. 1.,-1.,1.)
iclr_powindex ?

3. Experiment with reweighting for data gaps
First run:
Make a series with zeroes where the data to be examined has data gaps, one's elsewhere.

Namelist parameters with identical spectrum design parameters, except:
Output the ACF from the graphic window using W-command

Second run with the original data series:
For reweighting, include the file in the open-file block
59 B onoff.acv

In the namelist, add

and keep your fingers crossed.


Example input file: Estimate MEM- and Power Spectrum

Note that the file open block is scanned for environment parameters, e.g.
setenv sasm03_in on/
to override the default.

21 ^ ${sasm03_in:[owf/g100201-130515-1h-lapww-M6.ra]}.ts
25 * ^ ${sasm03_pef}                                             (for by-passing MEM, open this file)           
40 * ^
${sasm03_in:[owf/g100201-130515-1h-lapww-M6.ra]}.pef      (to open if iun_wh=40: remove "* ")
41 < ${sasm03_in:[owf/g100201-130515-1h-lapww-M6.ra]}.pef
42 B o/grav-1h.psp
43 B o/grav-1h.memsp
99 B /dev/pts/1
  lvllen=-1, qshlvl=.false.
  rec_mrs=-99999.9, dat_mrs=-99999.9
  l_section=1024, logxax=.true.
  target='BIN', fmt='$noiseexp'
  mx_miss=0, dff1=0.d0, ibeg_burg=24, qrfsp=.false.
  lpef=202, pef_order=200, lpef_wh=-1, iun_wh=-1, qrpefsp=.false.
  pef_show=11, pef_order_show=200,180,40,20,12,9,7,6,5,4,3
  t_window='KAIS', p_window=2.2
  mpix=-1, pix(1)=1.d0, pix(2)=-1.d0, pix(3)=0.001d0
  q_tsf_edit=.false., tsf_edit_name='LOPASS'
  qonly_mem=.false., nmems=16384, iumemspout=43

FILTER  D:WD KAIS 100 2.1 0.0 0.02 0.5 0.0

Example: generated by   pspd +T -Emodulate.tse,R + o/ -w -r
The time-series display options and ACV display are switched off.
Auto-write of the power spectrum and confidence data is switched on. Output to
tsf-edit is switched off here; it was invoked in the pre-tslist step.

21 ^ /data1/hgs/scratch/ts/tsftmp.ts
41 < ${ODIR:[o]}/
42 B ${ODIR:[o]}/
43 B ${ODIR:[o]}/
99 * nothing
   C setenv OTHERTERM 'B /dev/pts#'
  Display=' '
  lvllen=-1, qshlvl=.false.
  rec_mrs=-99999.9, dat_mrs=-99999.9
  l_section=2048, logxax=.true.
  target='BIN', fmt='U'
  dff1=0.0d0, qrfsp=.false.
  lpef_wh=-1, iun_wh=-1, qrpefsp=.false.
  lpef=202, pef_order=200
  pef_show=9, pef_order_show=200,150,110,80,50,30,20,10,5
  t_window='KAIS', p_window=2.2
  mpix=-1, pix(1)=1.d0, pix(2)=-1.d0, pix(3)=0.001d0
  qgra=.true., q_shut_graphics=.false.
  qsp_autowf=.true., wfname='o/'
  q_tsf_edit=.false., tsf_edit_name='LOPASS'
  iuspout=42, iumemspout=43
  qonly_mem=.false., nmems=16384