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 PEF or a Wiener filter and
     estimate the power spectrum using MEM and {Periodogram | Bartlett}
     Display: time series, Filter spectrum,
              filtered time series,
              Spectrum of the filtered t.s. using
                       Bartlett_ACV and _PSP) + MEM_PSP

(2) Treat the time series with a simple difference filter and
     estimate the power spectrum using  MEM and {Periodogram | Bartlett}
     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.
     Display: time series and
              filtered time series
              Spectrum of the filtered t.s. using
                       Bartlett_ACV and _PSP) + MEM_PSP
                       optionally with the effect of the difference
                       filter undone.

The MEM spectrum is always used to draw a confidence region around the
Bartlett spectrum.


Files to open: A: ASCII, B: BIN, D: ASCII or BIN depends on target, fmt in
                namelist &param
D     21 - signal 'X' input
D     22 ~ signal 'Y' to subtract from 'X'
B     25 ~ PEF filter input for confidence limits
B     41 - PEF filter output                           1)
A     42 - PSP output
A iun_wf ~ Wiener filter input
B iun_wh - PEF whitening filter output / re-input      1)
          ~ open optional
          - open compulsory

1) If unit 25 is not opened, 41 must be opened since a new filter bank
    will be produced.

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

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)

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]

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

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.

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 ?