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 ¶m
----------------------------------------------------------------------------
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
operations)
'X' or any other char: compute PSP of time series X
['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.
Specify:
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 ?
____________________________________________________________________________