sasm03  - Power Spectrum Estimation with MEM and Bartlett-method

USAGE:
    sasm031 @ <ins-file> [:target]

PURPOSE:
Batch processing, Bartlett power spectrum and/or MEM-method incl. Prediction Error Filter (PEF) assembly

INSTRUCTION FILE:
target>
Open-file block
Namelist  &param

FILES TO OPEN:

Data description
 PSP / MEM
 File type
 I / O
 File unit
Time series, ts-file or mc-file
 PSP + MEM
 BIN
 I
 21
Autocovariance
 PSP
 ASCII
 O
 iuo_acv
Windowed Autocovariance
 PSP
 ASCII
 O
 iuo_wacv
Power spectrum
 PSP
 ASCII
 O
 iuo_psp
PEF-bank for prefilter
 PSP+MEM
 BIN
 I
 iui_pef
PEF-bank for MEM
 MEM
 BIN
 I
 iui_pef
New PEF-bank
 MEM
 BIN
 O
 iuo_pef
MEM spectrum
 MEM
 ASCII
 O
 iuo_memsp
Diagnostic print 1)
 PSP
 ASCII
 O
 99
________________________________________________________
1) If not opened, a file fort.99 will be produced.
    In the open-file block you can add
    99 ${OTHERTERM:[B tmp/sasm031.dgn]}
    or 
    99 ${OTHERTERM:[* ]}
    and issue  
    setenv OTHERTERM "B /dev/pts/n"

NAMELIST &param

Parameter
 ____type___
 description
 default
PSP




l_section
 integer
 Bartlett section length
 1024

l_window
 integer
 Window length, could be l_section/2
 l_section

t_window
 char*4
 Window type ('HANN' 'KAIS' ..., cf ~/sas/p/window.f
 'KAIS'

p_window
 real*8
 Window parameter, depends on type, default for 'KAIS':
 1.8,

qbartlett_psp
 logical
 Do PSP
 .true.

qpsp_db
 logical
 Output is in decibels
 .true.

frq_scale
 char*16
 'mHz'  'Hz'  'cyc/h' 'cyc/d'  'cyc/y'
 'cyc/h'

iuo_psp
 integer
 Output file unit for PSP, -1 means no output
 -1

iuo_acv
 integer
 ... for Autocovariance
 -1

iuo_acv integer ... for windowed Autocovariance -1





MEM



PEF
 ibeg_burg
 integer
Where in the time series to start
 0

l_burg
 integer
 Length, 0 means all
 0

qnew_pef
 logical
 Produce new PEF
 .true.

iui_pef
 integer
 if not qnew_pef, input PEF from this file unit
 -1





SP
 maxpef_memsp
 integer
 Maximum length of new PEF
 21

npef_memsp
 integer
 Length of PEF to use in MEM spectrum
 10

l_memsp
 integer
 Length of MEM spectrum domain
 1024

qmemsp
 logical
 Do MEM spectrum
 .false.

qpsp_db
 logical
 identical as in PSP-mode

frq_scale
 char*16
 identical as in PSP-mode

iuo_memsp
 ingeger
 Output file unit for MEMSP, -1 means no output -1





INPUT




fmt
 char*64
 Input ascii: the format expression to read_fuf
 'U'



'U' - binary data, unlabeled



'L:label' - binary data (mc-file), the label

target
 char*16
 The target string to read_fuf, 'BIN' for binary files
 'BIN'

hms_pos
 integer
 0 - 4, ascii input, depth of time of day, h ... ms
 0

ibeg, iend
 integer
 range of input data, sample numbers, 0,0 means all
 0, 0

qrepair
 logical
 repair gaps by replacement with DC-level, occurs after diff-filter
 .true.

dff1
 real*8
 Difference filter (1, dff1)
 (1., -1.)

qrfsp
 logical
 Compensate diff-filter attenuation in output
 .true.

pef_order 1)
 integer
 Pre-filter type PEF, if iui_pef > 0
 3
________________________________________________________
1) This filter cannot be compensated in output - yet. T.b.done as need arises....
What to do with PSP values that are negative?
It happens with too pure sinusoids. It leads to NaN when taking the logarithm.
First, run with qpsp_db=.false. and check the values.
There will be an option to choose among: Suppress, use Absolute value, substitute a Default, substitute the Mean of the neighbours, ignore and print NaN
  psp_repair_option='c'  , c  € { S A D M N }

EXAMPLE:
cd ~/TD
sasm031 @ sasm031.ins :M

M>
21 ^ d/G1_garb_160102-1s.mc
31 > o/G1_garb_160102-1s.pef
41 B o/G1_garb_160102-1s.memsp
   q
 &param
 qbartlett_psp=.false.
 fmt='L:G|B'
 iuo_pef=31
 qnew_pef=.true.
 qmemsp=.true., l_memsp=65536, npef_memsp=20, maxpef_memsp=40
 iuo_memsp=41
 qrfsp=.true.
 frq_scale='mHz'
 &end


sasm031 @ sasm031.ins :B

B>
21 ^ d/G1_garb_160102-1s.mc
31 B o/G1_garb_160102-1s.acv
41 B o/G1_garb_160102-1s.psp
   q
 &param
 fmt='L:G|B'
 l_section=2048
 l_window=1024
 p_window=2.2
 iuo_acv=31
 iuo_psp=41
 frq_scale='mHz'
 &end

More examples needed (and t.b. tested): Pre-filter with a PEF, construct a new PEF

ddssdd

.bye