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:
ier=int_prompt_s18
&('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
qsavesp=.true.
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)
Example:
...
42 B ${ODIR:[o]}/grav-1h.psp
43 B ${ODIR:[o]}/grav-1h.memsp
...
¶m
...
iuspout=42, iumemspout=43
...
&end
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
GRDSPS>
to enter W filename, example
W owf/grav.ra.psp
3.2 Automatic, same example:
In the namelist, set
qsp_autowf=.true., wfname='owf/grav.ra.psp'
You can specify environment parameters with
defaults, like
wfname='${ODIR:[o]}/grav.ra.psp'
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
owf/grav.ra.psp.cfd
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 ¶m
----------------------------------------------------------------------------
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 qsavesp=.true.
B 52 - Power spectrum, optional (code change necessary: don't if not open)
B $w ~ weights to restore effects of data gaps in
covariance 2)
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
lvl_debug - Debug level, relevant also for output to unit 99
0 - only little is
printed,
1 - medium, with /cufios/ qpr print option .true.
2 - much, with /cufios/ qpr,qdbg options
.true.
[0]
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)
nmems - Length of spectrum
domain for MEM; note ...
qonly_mem - If .true., stop after the MEM spectrum has been
computed.
nmems will be the length of the spectrum then.
If .false., the length of the MEM spectrum is taken
from 2 * l_section, so that the MEM spectrum can be used
to plot the uncertainty corridor around the Bartlett
spectrum.
[.false.]
rec_mrs &
rec_mrs_y - Symbolic value for missing records in
input file [2*999.9]
dat_mrs -
" "
" " internal data
qremdc - Remove DC-level from
input
[.true.]
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
qspout_frsp -
if .true. under
qsavesp=.true. and unit 51 opened,
output the MEM spectra after diff-filter restoration,
if .false. ..., before.
[.false.]
iu_bartlett_cov_reweight
-
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.
Use `]` as the target string delimiter.
Specify `BIN´ for
binary files
[2*'>> ]']
fmt &
fmt_y - (char*64) Data
format. Use `BIN´ for binary access.
fmt only: environment variables will be replaced.
Specify `L<@´
if a label is to be read from the
comment part of unit 21 in the open-file block.
['(3i2,1x,i2,4x,f15.0)']
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
start_hour - (real*8) Start at this hour
q_cut_at_epoch - (logical) Must be set .true. that the start
parameters
take
effect
[start at first sample]
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
_____________________________________________________________________________
Environment parameters
Presently just a few which, when non-empty strings
are assigned, will supersede
the namelist settings.
SASM03_GRA - don't use graphics if the
value starts with `N´ `n´ `F´ or `f´
also `.f´
and `.F´ will work. I.e. "no" or "false".
SASM03_DFF1 - The value v of the difference
filter (z - v) in z-transform lingo.
SASM03_PEFORDERS -
Alt.1: #pef_order_show(1),#pef_show,#pef_order_show(#pef_show),AUTO
The PEF-orders shown will step down a harmonic ladder.
Alt.2: #pef_order_show(1),#pef_order_show(2),...
Max.
20 orders can be shown.
_____________________________________________________________________________
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.
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.
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:
dff1=0.d0
scale=1.d0
q_bartlett_keepdc=.true.
acv_file='onoff.acv'
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
iu_bartlett_cov_reweight=59
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/g100201-130515-1h-lapww-M6.ra
to override the default.
Also fmt is scanned for environment parameters.
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
Q
¶m
lvllen=-1, qshlvl=.false.
rec_mrs=-99999.9, dat_mrs=-99999.9
l_section=1024, logxax=.true.
target='BIN', fmt='${noiseexp:[U]}'
units='[mm]$'
mx_miss=0, dff1=0.d0, ibeg_burg=24, qrfsp=.false.
qrepair=.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
qgra=.true.
cutlvl=1.d-8
scale=1.0d0
time_scale=1.d0
q_tsf_edit=.false., tsf_edit_name='LOPASS'
iuspout=42
qonly_mem=.false., nmems=16384, iumemspout=43
&end
TSF EDIT LOPASS
FILTER D:WD KAIS 100 2.1 0.0 0.02 0.5 0.0
END
Example: generated by pspd +T -Emodulate.tse,R +
o/g090701-140131-1h-boa-nk.pt.ts -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
o/g090701-140131-1h-boa-nk.pt.pwrsp.cfd
tsf-edit is switched off here; it was invoked in the pre-tslist
step.
21 ^ /data1/hgs/scratch/ts/tsftmp.ts
41 < ${ODIR:[o]}/g090701-140131-1h-boa-nk.pt.pef
42 B ${ODIR:[o]}/g090701-140131-1h-boa-nk.pt.psp
43 B ${ODIR:[o]}/g090701-140131-1h-boa-nk.pt.memsp
99 * nothing
C setenv OTHERTERM 'B /dev/pts#'
Q
¶m
Display=' '
lvllen=-1, qshlvl=.false.
rec_mrs=-99999.9, dat_mrs=-99999.9
l_section=2048, logxax=.true.
target='BIN', fmt='U'
units='[mm]$'
mx_miss=0
dff1=0.0d0, qrfsp=.false.
qrepair=.true.
ibeg_burg=24
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/g090701-140131-1h-boa-nk.pt.pwrsp'
cutlvl=1.d-8
scale=1.0d0
time_scale=1.d0
q_tsf_edit=.false., tsf_edit_name='LOPASS'
iuspout=42, iumemspout=43
qonly_mem=.false., nmems=16384
&end